IT++ Logo

gf2mat.cpp

Go to the documentation of this file.
00001 
00030 #include <itpp/base/gf2mat.h>
00031 #include <itpp/base/specmat.h>
00032 #include <itpp/base/matfunc.h>
00033 #include <itpp/base/converters.h>
00034 #include <iostream>
00035 
00036 namespace itpp
00037 {
00038 
00039 // ====== IMPLEMENTATION OF THE ALIST CLASS ==========
00040 
00041 GF2mat_sparse_alist::GF2mat_sparse_alist(const std::string &fname)
00042     : data_ok(false)
00043 {
00044   read(fname);
00045 }
00046 
00047 void GF2mat_sparse_alist::read(const std::string &fname)
00048 {
00049   std::ifstream file;
00050   std::string line;
00051   std::stringstream ss;
00052 
00053   file.open(fname.c_str());
00054   it_assert(file.is_open(),
00055             "GF2mat_sparse_alist::read(): Could not open file \""
00056             << fname << "\" for reading");
00057 
00058   // parse N and M:
00059   getline(file, line);
00060   ss << line;
00061   ss >> N >> M;
00062   it_assert(!ss.fail(),
00063             "GF2mat_sparse_alist::read(): Wrong alist data (N or M)");
00064   it_assert((N > 0) && (M > 0),
00065             "GF2mat_sparse_alist::read(): Wrong alist data");
00066   ss.seekg(0, std::ios::end);
00067   ss.clear();
00068 
00069   // parse max_num_n and max_num_m
00070   getline(file, line);
00071   ss << line;
00072   ss >> max_num_n >> max_num_m;
00073   it_assert(!ss.fail(),
00074             "GF2mat_sparse_alist::read(): Wrong alist data (max_num_{n,m})");
00075   it_assert((max_num_n >= 0) && (max_num_n <= N) &&
00076             (max_num_m >= 0) && (max_num_m <= M),
00077             "GF2mat_sparse_alist::read(): Wrong alist data");
00078   ss.seekg(0, std::ios::end);
00079   ss.clear();
00080 
00081   // parse weight of each column n
00082   num_nlist.set_size(N);
00083   num_nlist.clear();
00084   getline(file, line);
00085   ss << line;
00086   for (int i = 0; i < N; i++) {
00087     ss >> num_nlist(i);
00088     it_assert(!ss.fail(),
00089               "GF2mat_sparse_alist::read(): Wrong alist data (num_nlist("
00090               << i << "))");
00091     it_assert((num_nlist(i) >= 0) && (num_nlist(i) <= M),
00092               "GF2mat_sparse_alist::read(): Wrong alist data (num_nlist("
00093               << i << "))");
00094   }
00095   ss.seekg(0, std::ios::end);
00096   ss.clear();
00097 
00098   // parse weight of each row m
00099   num_mlist.set_size(M);
00100   num_mlist.clear();
00101   getline(file, line);
00102   ss << line;
00103   for (int i = 0; i < M; i++) {
00104     ss >> num_mlist(i);
00105     it_assert(!ss.fail(),
00106               "GF2mat_sparse_alist::read(): Wrong alist data (num_mlist("
00107               << i << "))");
00108     it_assert((num_mlist(i) >= 0) && (num_mlist(i) <= N),
00109               "GF2mat_sparse_alist::read(): Wrong alist data (num_mlist("
00110               << i << "))");
00111   }
00112   ss.seekg(0, std::ios::end);
00113   ss.clear();
00114 
00115   // parse coordinates in the n direction with non-zero entries
00116   nlist.set_size(N, max_num_n);
00117   nlist.clear();
00118   for (int i = 0; i < N; i++) {
00119     getline(file, line);
00120     ss << line;
00121     for (int j = 0; j < num_nlist(i); j++) {
00122       ss >> nlist(i, j);
00123       it_assert(!ss.fail(),
00124                 "GF2mat_sparse_alist::read(): Wrong alist data (nlist("
00125                 << i << "," << j << "))");
00126       it_assert((nlist(i, j) >= 0) && (nlist(i, j) <= M),
00127                 "GF2mat_sparse_alist::read(): Wrong alist data (nlist("
00128                 << i << "," << j << "))");
00129     }
00130     ss.seekg(0, std::ios::end);
00131     ss.clear();
00132   }
00133 
00134   // parse coordinates in the m direction with non-zero entries
00135   mlist.set_size(M, max_num_m);
00136   mlist.clear();
00137   for (int i = 0; i < M; i++) {
00138     getline(file, line);
00139     ss << line;
00140     for (int j = 0; j < num_mlist(i); j++) {
00141       ss >> mlist(i, j);
00142       it_assert(!ss.fail(),
00143                 "GF2mat_sparse_alist::read(): Wrong alist data (mlist("
00144                 << i << "," << j << "))");
00145       it_assert((mlist(i, j) >= 0) && (mlist(i, j) <= N),
00146                 "GF2mat_sparse_alist::read(): Wrong alist data (mlist("
00147                 << i << "," << j << "))");
00148     }
00149     ss.seekg(0, std::ios::end);
00150     ss.clear();
00151   }
00152 
00153   file.close();
00154   data_ok = true;
00155 }
00156 
00157 void GF2mat_sparse_alist::write(const std::string &fname) const
00158 {
00159   it_assert(data_ok,
00160             "GF2mat_sparse_alist::write(): alist data not ready for writing");
00161 
00162   std::ofstream file(fname.c_str(), std::ofstream::out);
00163   it_assert(file.is_open(),
00164             "GF2mat_sparse_alist::write(): Could not open file \""
00165             << fname << "\" for writing");
00166 
00167   file << N << " " << M << std::endl;
00168   file << max_num_n << " " << max_num_m << std::endl;
00169 
00170   for (int i = 0; i < num_nlist.length() - 1; i++)
00171     file << num_nlist(i) << " ";
00172   file << num_nlist(num_nlist.length() - 1) << std::endl;
00173 
00174   for (int i = 0; i < num_mlist.length() - 1; i++)
00175     file << num_mlist(i) << " ";
00176   file << num_mlist(num_mlist.length() - 1) << std::endl;
00177 
00178   for (int i = 0; i < N; i++) {
00179     for (int j = 0; j < num_nlist(i) - 1; j++)
00180       file << nlist(i, j) << " ";
00181     file << nlist(i, num_nlist(i) - 1) << std::endl;
00182   }
00183 
00184   for (int i = 0; i < M; i++) {
00185     for (int j = 0; j < num_mlist(i) - 1; j++)
00186       file << mlist(i, j) << " ";
00187     file << mlist(i, num_mlist(i) - 1) << std::endl;
00188   }
00189 
00190   file.close();
00191 }
00192 
00193 
00194 GF2mat_sparse GF2mat_sparse_alist::to_sparse(bool transpose) const
00195 {
00196   GF2mat_sparse sbmat(M, N, max_num_m);
00197 
00198   for (int i = 0; i < M; i++) {
00199     for (int j = 0; j < num_mlist(i); j++) {
00200       sbmat.set_new(i, mlist(i, j) - 1, bin(1));
00201     }
00202   }
00203   sbmat.compact();
00204 
00205   if (transpose) {
00206     return sbmat.transpose();
00207   }
00208   else {
00209     return sbmat;
00210   }
00211 }
00212 
00213 
00214 // ----------------------------------------------------------------------
00215 // WARNING: This method is very slow. Sparse_Mat has to be extended with
00216 // some extra functions to improve the performance of this.
00217 // ----------------------------------------------------------------------
00218 void GF2mat_sparse_alist::from_sparse(const GF2mat_sparse &sbmat,
00219                                       bool transpose)
00220 {
00221   if (transpose) {
00222     from_sparse(sbmat.transpose(), false);
00223   }
00224   else {
00225     // check matrix dimension
00226     M = sbmat.rows();
00227     N = sbmat.cols();
00228 
00229     num_mlist.set_size(M);
00230     num_nlist.set_size(N);
00231 
00232     // fill mlist matrix, num_mlist vector and max_num_m
00233     mlist.set_size(0, 0);
00234     int tmp_cols = 0; // initial number of allocated columns
00235     for (int i = 0; i < M; i++) {
00236       ivec temp_row(0);
00237       for (int j = 0; j < N; j++) {
00238         if (sbmat(i, j) == bin(1)) {
00239           temp_row = concat(temp_row, j + 1);
00240         }
00241       }
00242       int trs = temp_row.size();
00243       if (trs > tmp_cols) {
00244         tmp_cols = trs;
00245         mlist.set_size(M, tmp_cols, true);
00246       }
00247       else if (trs < tmp_cols) {
00248         temp_row.set_size(tmp_cols, true);
00249       }
00250       mlist.set_row(i, temp_row);
00251       num_mlist(i) = trs;
00252     }
00253     max_num_m = max(num_mlist);
00254 
00255     // fill nlist matrix, num_nlist vector and max_num_n
00256     nlist.set_size(0, 0);
00257     tmp_cols = 0; // initial number of allocated columns
00258     for (int j = 0; j < N; j++) {
00259       ivec temp_row = sbmat.get_col(j).get_nz_indices() + 1;
00260       int trs = temp_row.size();
00261       if (trs > tmp_cols) {
00262         tmp_cols = trs;
00263         nlist.set_size(N, tmp_cols, true);
00264       }
00265       else if (trs < tmp_cols) {
00266         temp_row.set_size(tmp_cols, true);
00267       }
00268       nlist.set_row(j, temp_row);
00269       num_nlist(j) = trs;
00270     }
00271     max_num_n = max(num_nlist);
00272 
00273     data_ok = true;
00274   }
00275 }
00276 
00277 
00278 // ----------------------------------------------------------------------
00279 // Implementation of a dense GF2 matrix class
00280 // ----------------------------------------------------------------------
00281 
00282 GF2mat::GF2mat(int i, int j): nrows(i), ncols(j),
00283     nwords((j >> shift_divisor) + 1)
00284 {
00285   data.set_size(nrows, nwords);
00286   data.clear();
00287 }
00288 
00289 GF2mat::GF2mat(): nrows(1), ncols(1), nwords(1)
00290 {
00291   data.set_size(nrows, nwords);
00292   data.clear();
00293 }
00294 
00295 GF2mat::GF2mat(const bvec &x, bool is_column)
00296 {
00297   if (is_column) {     // create column vector
00298     nrows = length(x);
00299     ncols = 1;
00300     nwords = 1;
00301     data.set_size(nrows, nwords);
00302     data.clear();
00303     for (int i = 0; i < nrows; i++) {
00304       set(i, 0, x(i));
00305     }
00306   }
00307   else {    // create row vector
00308     nrows = 1;
00309     ncols = length(x);
00310     nwords = (ncols >> shift_divisor) + 1;
00311     data.set_size(nrows, nwords);
00312     data.clear();
00313     for (int i = 0; i < ncols; i++) {
00314       set(0, i, x(i));
00315     }
00316   }
00317 }
00318 
00319 
00320 GF2mat::GF2mat(const bmat &X): nrows(X.rows()), ncols(X.cols())
00321 {
00322   nwords = (ncols >> shift_divisor) + 1;
00323   data.set_size(nrows, nwords);
00324   data.clear();
00325   for (int i = 0; i < nrows; i++) {
00326     for (int j = 0; j < ncols; j++) {
00327       set(i, j, X(i, j));
00328     }
00329   }
00330 }
00331 
00332 
00333 GF2mat::GF2mat(const GF2mat_sparse &X)
00334 {
00335   nrows = X.rows();
00336   ncols = X.cols();
00337   nwords = (ncols >> shift_divisor) + 1;
00338   data.set_size(nrows, nwords);
00339   for (int i = 0; i < nrows; i++) {
00340     for (int j = 0; j < nwords; j++) {
00341       data(i, j) = 0;
00342     }
00343   }
00344 
00345   for (int j = 0; j < ncols; j++) {
00346     for (int i = 0; i < X.get_col(j).nnz(); i++)   {
00347       bin b = X.get_col(j).get_nz_data(i);
00348       set(X.get_col(j).get_nz_index(i), j, b);
00349     }
00350   }
00351 }
00352 
00353 GF2mat::GF2mat(const GF2mat_sparse &X, int m1, int n1, int m2, int n2)
00354 {
00355   it_assert(X.rows() > m2, "GF2mat(): indexes out of range");
00356   it_assert(X.cols() > n2, "GF2mat(): indexes out of range");
00357   it_assert(m1 >= 0 && n1 >= 0 && m2 >= m1 && n2 >= n1,
00358             "GF2mat::GF2mat(): indexes out of range");
00359 
00360   nrows = m2 - m1 + 1;
00361   ncols = n2 - n1 + 1;
00362   nwords = (ncols >> shift_divisor) + 1;
00363   data.set_size(nrows, nwords);
00364 
00365   for (int i = 0; i < nrows; i++) {
00366     for (int j = 0; j < nwords; j++) {
00367       data(i, j) = 0;
00368     }
00369   }
00370 
00371   for (int i = 0; i < nrows; i++) {
00372     for (int j = 0; j < ncols; j++)   {
00373       bin b = X(i + m1, j + n1);
00374       set(i, j, b);
00375     }
00376   }
00377 }
00378 
00379 
00380 GF2mat::GF2mat(const GF2mat_sparse &X, const ivec &columns)
00381 {
00382   it_assert(X.cols() > max(columns),
00383             "GF2mat::GF2mat(): index out of range");
00384   it_assert(min(columns) >= 0,
00385             "GF2mat::GF2mat(): column index must be positive");
00386 
00387   nrows = X.rows();
00388   ncols = length(columns);
00389   nwords = (ncols >> shift_divisor) + 1;
00390   data.set_size(nrows, nwords);
00391 
00392   for (int i = 0; i < nrows; i++) {
00393     for (int j = 0; j < nwords; j++) {
00394       data(i, j) = 0;
00395     }
00396   }
00397 
00398   for (int j = 0; j < ncols; j++) {
00399     for (int i = 0; i < X.get_col(columns(j)).nnz(); i++)   {
00400       bin b = X.get_col(columns(j)).get_nz_data(i);
00401       set(X.get_col(columns(j)).get_nz_index(i), j, b);
00402     }
00403   }
00404 }
00405 
00406 
00407 void GF2mat::set_size(int m, int n, bool copy)
00408 {
00409   nrows = m;
00410   ncols = n;
00411   nwords = (ncols >> shift_divisor) + 1;
00412   data.set_size(nrows, nwords, copy);
00413   if (!copy)
00414     data.clear();
00415 }
00416 
00417 
00418 GF2mat_sparse GF2mat::sparsify() const
00419 {
00420   GF2mat_sparse Z(nrows, ncols);
00421   for (int i = 0; i < nrows; i++) {
00422     for (int j = 0; j < ncols; j++) {
00423       if (get(i, j) == 1) {
00424         Z.set(i, j, 1);
00425       }
00426     }
00427   }
00428 
00429   return Z;
00430 }
00431 
00432 bvec GF2mat::bvecify() const
00433 {
00434   it_assert(nrows == 1 || ncols == 1,
00435             "GF2mat::bvecify() matrix must be a vector");
00436   int n = (nrows == 1 ? ncols : nrows);
00437   bvec result(n);
00438   if (nrows == 1) {
00439     for (int i = 0; i < n; i++) {
00440       result(i) = get(0, i);
00441     }
00442   }
00443   else {
00444     for (int i = 0; i < n; i++) {
00445       result(i) = get(i, 0);
00446     }
00447   }
00448   return result;
00449 }
00450 
00451 
00452 void GF2mat::set_row(int i, bvec x)
00453 {
00454   it_assert(length(x) == ncols,
00455             "GF2mat::set_row(): dimension mismatch");
00456   for (int j = 0; j < ncols; j++) {
00457     set(i, j, x(j));
00458   }
00459 }
00460 
00461 void GF2mat::set_col(int j, bvec x)
00462 {
00463   it_assert(length(x) == nrows,
00464             "GF2mat::set_col(): dimension mismatch");
00465   for (int i = 0; i < nrows; i++) {
00466     set(i, j, x(i));
00467   }
00468 }
00469 
00470 
00471 GF2mat gf2dense_eye(int m)
00472 {
00473   GF2mat Z(m, m);
00474   for (int i = 0; i < m; i++) {
00475     Z.set(i, i, 1);
00476   }
00477   return Z;
00478 }
00479 
00480 GF2mat GF2mat::get_submatrix(int m1, int n1, int m2, int n2) const
00481 {
00482   it_assert(m1 >= 0 && n1 >= 0 && m2 >= m1 && n2 >= n1
00483             && m2 < nrows && n2 < ncols,
00484             "GF2mat::get_submatrix() index out of range");
00485   GF2mat result(m2 - m1 + 1, n2 - n1 + 1);
00486 
00487   for (int i = m1; i <= m2; i++) {
00488     for (int j = n1; j <= n2; j++) {
00489       result.set(i - m1, j - n1, get(i, j));
00490     }
00491   }
00492 
00493   return result;
00494 }
00495 
00496 
00497 GF2mat GF2mat::concatenate_vertical(const GF2mat &X) const
00498 {
00499   it_assert(X.ncols == ncols,
00500             "GF2mat::concatenate_vertical(): dimension mismatch");
00501   it_assert(X.nwords == nwords,
00502             "GF2mat::concatenate_vertical(): dimension mismatch");
00503 
00504   GF2mat result(nrows + X.nrows, ncols);
00505   for (int i = 0; i < nrows; i++) {
00506     for (int j = 0; j < nwords; j++) {
00507       result.data(i, j) = data(i, j);
00508     }
00509   }
00510 
00511   for (int i = 0; i < X.nrows; i++) {
00512     for (int j = 0; j < nwords; j++) {
00513       result.data(i + nrows, j) = X.data(i, j);
00514     }
00515   }
00516 
00517   return result;
00518 }
00519 
00520 GF2mat GF2mat::concatenate_horizontal(const GF2mat &X) const
00521 {
00522   it_assert(X.nrows == nrows,
00523             "GF2mat::concatenate_horizontal(): dimension mismatch");
00524 
00525   GF2mat result(nrows, X.ncols + ncols);
00526   for (int i = 0; i < nrows; i++) {
00527     for (int j = 0; j < ncols; j++) {
00528       result.set(i, j, get(i, j));
00529     }
00530   }
00531 
00532   for (int i = 0; i < nrows; i++) {
00533     for (int j = 0; j < X.ncols; j++) {
00534       result.set(i, j + ncols, X.get(i, j));
00535     }
00536   }
00537 
00538   return result;
00539 }
00540 
00541 bvec GF2mat::get_row(int i) const
00542 {
00543   bvec result(ncols);
00544   for (int j = 0; j < ncols; j++) {
00545     result(j) = get(i, j);
00546   }
00547 
00548   return result;
00549 }
00550 
00551 bvec GF2mat::get_col(int j) const
00552 {
00553   bvec result(nrows);
00554   for (int i = 0; i < nrows; i++) {
00555     result(i) = get(i, j);
00556   }
00557 
00558   return result;
00559 }
00560 
00561 
00562 int GF2mat::T_fact(GF2mat &T, GF2mat &U, ivec &perm) const
00563 {
00564   T = gf2dense_eye(nrows);
00565   U = *this;
00566 
00567   perm = zeros_i(ncols);
00568   for (int i = 0; i < ncols; i++) {
00569     perm(i) = i;
00570   }
00571 
00572   if (nrows > 250) {   // avoid cluttering output ...
00573     it_info_debug("Performing T-factorization of GF(2) matrix...  rows: "
00574                   << nrows << " cols: "  << ncols << " .... " << std::endl);
00575   }
00576   int pdone = 0;
00577   for (int j = 0; j < nrows; j++) {
00578     // Now working on diagonal element j,j
00579     // First try find a row with a 1 in column i
00580     int i1, j1;
00581     for (i1 = j; i1 < nrows; i1++) {
00582       for (j1 = j; j1 < ncols; j1++) {
00583         if (U.get(i1, j1) == 1) { goto found; }
00584       }
00585     }
00586 
00587     return j;
00588 
00589   found:
00590     U.swap_rows(i1, j);
00591     T.swap_rows(i1, j);
00592     U.swap_cols(j1, j);
00593 
00594     int temp = perm(j);
00595     perm(j) = perm(j1);
00596     perm(j1) = temp;
00597 
00598     // now subtract row i from remaining rows
00599     for (int i1 = j + 1; i1 < nrows; i1++)  {
00600       if (U.get(i1, j) == 1) {
00601         int ptemp = floor_i(100.0 * (i1 + j * nrows) / (nrows * nrows));
00602         if (nrows > 250 && ptemp > pdone + 10) {
00603           it_info_debug(ptemp << "% done.");
00604           pdone = ptemp;
00605         }
00606         U.add_rows(i1, j);
00607         T.add_rows(i1, j);
00608       }
00609     }
00610   }
00611   return nrows;
00612 }
00613 
00614 
00615 int GF2mat::T_fact_update_bitflip(GF2mat &T, GF2mat &U,
00616                                   ivec &perm, int rank,
00617                                   int r, int c) const
00618 {
00619   // First, update U (before re-triangulization)
00620   int c0 = c;
00621   for (c = 0; c < ncols; c++) {
00622     if (perm(c) == c0) {
00623       goto foundidx;
00624     }
00625   }
00626   it_error("GF2mat::T_fact_update_bitflip() - internal error");
00627 
00628 foundidx:
00629   for (int i = 0; i < nrows; i++) {
00630     if (T.get(i, r) == 1) {
00631       U.addto_element(i, c, 1);
00632     }
00633   }
00634 
00635   // first move column c to the end
00636   bvec lastcol = U.get_col(c);
00637   int temp_perm = perm(c);
00638   for (int j = c; j < ncols - 1; j++) {
00639     U.set_col(j, U.get_col(j + 1));
00640     perm(j) = perm(j + 1);
00641   }
00642   U.set_col(ncols - 1, lastcol);
00643   perm(ncols - 1) = temp_perm;
00644 
00645   // then, if the matrix is tall, also move row c to the end
00646   if (nrows >= ncols) {
00647     bvec lastrow_U = U.get_row(c);
00648     bvec lastrow_T = T.get_row(c);
00649     for (int i = c; i < nrows - 1; i++) {
00650       U.set_row(i, U.get_row(i + 1));
00651       T.set_row(i, T.get_row(i + 1));
00652     }
00653     U.set_row(nrows - 1, lastrow_U);
00654     T.set_row(nrows - 1, lastrow_T);
00655 
00656     // Do Gaussian elimination on the last row
00657     for (int j = c; j < ncols; j++)  {
00658       if (U.get(nrows - 1, j) == 1) {
00659         U.add_rows(nrows - 1, j);
00660         T.add_rows(nrows - 1, j);
00661       }
00662     }
00663   }
00664 
00665   // Now, continue T-factorization from the point (rank-1,rank-1)
00666   for (int j = rank - 1; j < nrows; j++) {
00667     int i1, j1;
00668     for (i1 = j; i1 < nrows; i1++) {
00669       for (j1 = j; j1 < ncols; j1++) {
00670         if (U.get(i1, j1) == 1) {
00671           goto found;
00672         }
00673       }
00674     }
00675 
00676     return j;
00677 
00678   found:
00679     U.swap_rows(i1, j);
00680     T.swap_rows(i1, j);
00681     U.swap_cols(j1, j);
00682 
00683     int temp = perm(j);
00684     perm(j) = perm(j1);
00685     perm(j1) = temp;
00686 
00687     for (int i1 = j + 1; i1 < nrows; i1++)  {
00688       if (U.get(i1, j) == 1) {
00689         U.add_rows(i1, j);
00690         T.add_rows(i1, j);
00691       }
00692     }
00693   }
00694 
00695   return nrows;
00696 }
00697 
00698 bool GF2mat::T_fact_update_addcol(GF2mat &T, GF2mat &U,
00699                                   ivec &perm, bvec newcol) const
00700 {
00701   int i0 = T.rows();
00702   int j0 = U.cols();
00703   it_assert(j0 > 0, "GF2mat::T_fact_update_addcol(): dimension mismatch");
00704   it_assert(i0 == T.cols(),
00705             "GF2mat::T_fact_update_addcol(): dimension mismatch");
00706   it_assert(i0 == U.rows(),
00707             "GF2mat::T_fact_update_addcol(): dimension mismatch");
00708   it_assert(length(perm) == j0,
00709             "GF2mat::T_fact_update_addcol(): dimension mismatch");
00710   it_assert(U.get(j0 - 1, j0 - 1) == 1,
00711             "GF2mat::T_fact_update_addcol(): dimension mismatch");
00712   // The following test is VERY TIME-CONSUMING
00713   it_assert_debug(U.row_rank() == j0,
00714                   "GF2mat::T_fact_update_addcol(): factorization has incorrect rank");
00715 
00716   bvec z = T * newcol;
00717   GF2mat Utemp = U.concatenate_horizontal(GF2mat(z, 1));
00718 
00719   // start working on position (j0,j0)
00720   int i;
00721   for (i = j0; i < i0; i++) {
00722     if (Utemp.get(i, j0) == 1) {
00723       goto found;
00724     }
00725   }
00726   return (false); // adding the new column would not improve the rank
00727 
00728 found:
00729   perm.set_length(j0 + 1, true);
00730   perm(j0) = j0;
00731 
00732   Utemp.swap_rows(i, j0);
00733   T.swap_rows(i, j0);
00734 
00735   for (int i1 = j0 + 1; i1 < i0; i1++)  {
00736     if (Utemp.get(i1, j0) == 1) {
00737       Utemp.add_rows(i1, j0);
00738       T.add_rows(i1, j0);
00739     }
00740   }
00741 
00742   U = Utemp;
00743   return (true); // the new column was successfully added
00744 }
00745 
00746 
00747 
00748 
00749 GF2mat GF2mat::inverse() const
00750 {
00751   it_assert(nrows == ncols, "GF2mat::inverse(): Matrix must be square");
00752 
00753   // first compute the T-factorization
00754   GF2mat T, U;
00755   ivec perm;
00756   int rank = T_fact(T, U, perm);
00757   it_assert(rank == ncols, "GF2mat::inverse(): Matrix is not full rank");
00758 
00759   // backward substitution
00760   for (int i = ncols - 2; i >= 0; i--) {
00761     for (int j = ncols - 1; j > i; j--) {
00762       if (U.get(i, j) == 1) {
00763         U.add_rows(i, j);
00764         T.add_rows(i, j);
00765       }
00766     }
00767   }
00768   T.permute_rows(perm, 1);
00769   return T;
00770 }
00771 
00772 int GF2mat::row_rank() const
00773 {
00774   GF2mat T, U;
00775   ivec perm;
00776   int rank = T_fact(T, U, perm);
00777   return rank;
00778 }
00779 
00780 bool GF2mat::is_zero() const
00781 {
00782   for (int i = 0; i < nrows; i++) {
00783     for (int j = 0; j < nwords; j++) {
00784       if (data(i, j) != 0) {
00785         return false;
00786       }
00787     }
00788   }
00789   return true;
00790 }
00791 
00792 bool GF2mat::operator==(const GF2mat &X) const
00793 {
00794   if (X.nrows != nrows) { return false; }
00795   if (X.ncols != ncols) { return false; }
00796   it_assert(X.nwords == nwords, "GF2mat::operator==() dimension mismatch");
00797 
00798   for (int i = 0; i < nrows; i++) {
00799     for (int j = 0; j < nwords; j++) {
00800       //      if (X.get(i,j)!=get(i,j)) {
00801       if (X.data(i, j) != data(i, j)) {
00802         return false;
00803       }
00804     }
00805   }
00806   return true;
00807 }
00808 
00809 
00810 void GF2mat::add_rows(int i, int j)
00811 {
00812   it_assert(i >= 0 && i < nrows, "GF2mat::add_rows(): out of range");
00813   it_assert(j >= 0 && j < nrows, "GF2mat::add_rows(): out of range");
00814   for (int k = 0; k < nwords; k++)   {
00815     data(i, k) ^= data(j, k);
00816   }
00817 }
00818 
00819 void GF2mat::swap_rows(int i, int j)
00820 {
00821   it_assert(i >= 0 && i < nrows, "GF2mat::swap_rows(): index out of range");
00822   it_assert(j >= 0 && j < nrows, "GF2mat::swap_rows(): index out of range");
00823   for (int k = 0; k < nwords; k++) {
00824     unsigned char temp = data(i, k);
00825     data(i, k) = data(j, k);
00826     data(j, k) = temp;
00827   }
00828 }
00829 
00830 void GF2mat::swap_cols(int i, int j)
00831 {
00832   it_assert(i >= 0 && i < ncols, "GF2mat::swap_cols(): index out of range");
00833   it_assert(j >= 0 && j < ncols, "GF2mat::swap_cols(): index out of range");
00834   bvec temp = get_col(i);
00835   set_col(i, get_col(j));
00836   set_col(j, temp);
00837 }
00838 
00839 
00840 void GF2mat::operator=(const GF2mat &X)
00841 {
00842   nrows = X.nrows;
00843   ncols = X.ncols;
00844   nwords = X.nwords;
00845   data = X.data;
00846 }
00847 
00848 GF2mat operator*(const GF2mat &X, const GF2mat &Y)
00849 {
00850   it_assert(X.ncols == Y.nrows, "GF2mat::operator*(): dimension mismatch");
00851   it_assert(X.nwords > 0, "Gfmat::operator*(): dimension mismatch");
00852   it_assert(Y.nwords > 0, "Gfmat::operator*(): dimension mismatch");
00853 
00854   /*
00855   // this can be done more efficiently?
00856   GF2mat result(X.nrows,Y.ncols);
00857   for (int i=0; i<X.nrows; i++) {
00858   for (int j=0; j<Y.ncols; j++) {
00859   bin b=0;
00860   for (int k=0; k<X.ncols; k++) {
00861   bin x = X.get(i,k);
00862   bin y = Y.get(k,j);
00863   b ^= (x&y);
00864   }
00865   result.set(i,j,b);
00866   }
00867   }
00868   return result; */
00869 
00870   // is this better?
00871   return mult_trans(X, Y.transpose());
00872 }
00873 
00874 bvec operator*(const GF2mat &X, const bvec &y)
00875 {
00876   it_assert(length(y) == X.ncols, "GF2mat::operator*(): dimension mismatch");
00877   it_assert(X.nwords > 0, "Gfmat::operator*(): dimension mismatch");
00878 
00879   /*
00880   // this can be done more efficiently?
00881   bvec result = zeros_b(X.nrows);
00882   for (int i=0; i<X.nrows; i++) {
00883   // do the nwords-1 data columns first
00884   for (int j=0; j<X.nwords-1; j++) {
00885   int ind = j<<shift_divisor;
00886   unsigned char r=X.data(i,j);
00887   while (r) {
00888   result(i) ^= (r & y(ind));
00889   r >>= 1;
00890   ind++;
00891   }
00892   }
00893   // do the last column separately
00894   for (int j=(X.nwords-1)<<shift_divisor; j<X.ncols; j++) {
00895   result(i) ^= (X.get(i,j) & y(j));
00896   }
00897   }
00898   return result; */
00899 
00900   // is this better?
00901   return (mult_trans(X, GF2mat(y, 0))).bvecify();
00902 }
00903 
00904 GF2mat mult_trans(const GF2mat &X, const GF2mat &Y)
00905 {
00906   it_assert(X.ncols == Y.ncols, "GF2mat::mult_trans(): dimension mismatch");
00907   it_assert(X.nwords > 0, "GF2mat::mult_trans(): dimension mismatch");
00908   it_assert(Y.nwords > 0, "GF2mat::mult_trans(): dimension mismatch");
00909   it_assert(X.nwords == Y.nwords, "GF2mat::mult_trans(): dimension mismatch");
00910 
00911   GF2mat result(X.nrows, Y.nrows);
00912 
00913   for (int i = 0; i < X.nrows; i++) {
00914     for (int j = 0; j < Y.nrows; j++) {
00915       bin b = 0;
00916       int kindx = i;
00917       int kindy = j;
00918       for (int k = 0; k < X.nwords; k++) {
00919         unsigned char r = X.data(kindx) & Y.data(kindy);
00920         /* The following can be speeded up by using a small lookup
00921            table for the number of ones and shift only a few times (or
00922            not at all if table is large) */
00923         while (r) {
00924           b ^= r & 1;
00925           r >>= 1;
00926         };
00927         kindx += X.nrows;
00928         kindy += Y.nrows;
00929       }
00930       result.set(i, j, b);
00931     }
00932   }
00933   return result;
00934 }
00935 
00936 GF2mat GF2mat::transpose() const
00937 {
00938   // CAN BE SPEEDED UP
00939   GF2mat result(ncols, nrows);
00940 
00941   for (int i = 0; i < nrows; i++) {
00942     for (int j = 0; j < ncols; j++) {
00943       result.set(j, i, get(i, j));
00944     }
00945   }
00946   return result;
00947 }
00948 
00949 GF2mat operator+(const GF2mat &X, const GF2mat &Y)
00950 {
00951   it_assert(X.nrows == Y.nrows, "GF2mat::operator+(): dimension mismatch");
00952   it_assert(X.ncols == Y.ncols, "GF2mat::operator+(): dimension mismatch");
00953   it_assert(X.nwords == Y.nwords, "GF2mat::operator+(): dimension mismatch");
00954   GF2mat result(X.nrows, X.ncols);
00955 
00956   for (int i = 0; i < X.nrows; i++) {
00957     for (int j = 0; j < X.nwords; j++) {
00958       result.data(i, j) = X.data(i, j) ^ Y.data(i, j);
00959     }
00960   }
00961 
00962   return result;
00963 }
00964 
00965 void GF2mat::permute_cols(ivec &perm, bool I)
00966 {
00967   it_assert(length(perm) == ncols,
00968             "GF2mat::permute_cols(): dimensions do not match");
00969 
00970   GF2mat temp = (*this);
00971   for (int j = 0; j < ncols; j++) {
00972     if (I == 0) {
00973       set_col(j, temp.get_col(perm(j)));
00974     }
00975     else {
00976       set_col(perm(j), temp.get_col(j));
00977     }
00978   }
00979 }
00980 
00981 void GF2mat::permute_rows(ivec &perm, bool I)
00982 {
00983   it_assert(length(perm) == nrows,
00984             "GF2mat::permute_rows(): dimensions do not match");
00985 
00986   GF2mat temp = (*this);
00987   for (int i = 0; i < nrows; i++) {
00988     if (I == 0) {
00989       for (int j = 0; j < nwords; j++) {
00990         data(i, j) = temp.data(perm(i), j);
00991       }
00992     }
00993     else {
00994       for (int j = 0; j < nwords; j++) {
00995         data(perm(i), j) = temp.data(i, j);
00996       }
00997     }
00998   }
00999 }
01000 
01001 
01002 std::ostream &operator<<(std::ostream &os, const GF2mat &X)
01003 {
01004   int i, j;
01005   os << "---- GF(2) matrix of dimension " << X.nrows << "*" << X.ncols
01006   << " -- Density: " << X.density() << " ----" << std::endl;
01007 
01008   for (i = 0; i < X.nrows; i++) {
01009     os << "      ";
01010     for (j = 0; j < X.ncols; j++) {
01011       os << X.get(i, j) << " ";
01012     }
01013     os << std::endl;
01014   }
01015 
01016   return os;
01017 }
01018 
01019 double GF2mat::density() const
01020 {
01021   int no_of_ones = 0;
01022 
01023   for (int i = 0; i < nrows; i++) {
01024     for (int j = 0; j < ncols; j++) {
01025       no_of_ones += (get(i, j) == 1 ? 1 : 0);
01026     }
01027   }
01028 
01029   return ((double) no_of_ones) / (nrows*ncols);
01030 }
01031 
01032 
01033 it_file &operator<<(it_file &f, const GF2mat &X)
01034 {
01035   // 3 64-bit unsigned words for: nrows, ncols and nwords + rest for char data
01036   uint64_t bytecount = 3 * sizeof(uint64_t)
01037                        + X.nrows * X.nwords * sizeof(char);
01038   f.write_data_header("GF2mat", bytecount);
01039 
01040   f.low_level_write(static_cast<uint64_t>(X.nrows));
01041   f.low_level_write(static_cast<uint64_t>(X.ncols));
01042   f.low_level_write(static_cast<uint64_t>(X.nwords));
01043   for (int i = 0; i < X.nrows; i++) {
01044     for (int j = 0; j < X.nwords; j++) {
01045       f.low_level_write(static_cast<char>(X.data(i, j)));
01046     }
01047   }
01048   return f;
01049 }
01050 
01051 it_ifile &operator>>(it_ifile &f, GF2mat &X)
01052 {
01053   it_file::data_header h;
01054 
01055   f.read_data_header(h);
01056   if (h.type == "GF2mat") {
01057     uint64_t tmp;
01058     f.low_level_read(tmp);
01059     X.nrows = static_cast<int>(tmp);
01060     f.low_level_read(tmp);
01061     X.ncols = static_cast<int>(tmp);
01062     f.low_level_read(tmp);
01063     X.nwords = static_cast<int>(tmp);
01064     X.data.set_size(X.nrows, X.nwords);
01065     for (int i = 0; i < X.nrows; i++) {
01066       for (int j = 0; j < X.nwords; j++) {
01067         char r;
01068         f.low_level_read(r);
01069         X.data(i, j) = static_cast<unsigned char>(r);
01070       }
01071     }
01072   }
01073   else {
01074     it_error("it_ifile &operator>>() - internal error");
01075   }
01076 
01077   return f;
01078 }
01079 
01080 } // namespace itpp
01081 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Wed Jan 20 23:03:03 2010 for IT++ by Doxygen 1.6.2