Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members   Examples  

gauss.h

Go to the documentation of this file.
00001 //N.Yu.Zolotykh 1999, 2000
00002 //University of Nizhni Novgorod, Russia
00003 
00004 #include"matrices.h"
00005 #include"gcd.h"
00006 #include"powerest.h"
00007 
00008 #ifndef GAUSS_H_
00009 #define GAUSS_H_
00010 
00017 
00028 
00029 #ifdef no_default_arguments
00030 template <class field_item>
00031   void rref(const matrix<field_item>& A,
00032     matrix<field_item>& B, matrix<field_item>& Q,
00033     vector<index>& basis, field_item& det,
00034     int movie, ostr& movie_stream);
00035 #else
00036 template <class field_item>
00037   void rref(const matrix<field_item>& A,
00038     matrix<field_item>& B, matrix<field_item>& Q,
00039     vector<index>& basis, field_item& det,
00040     int movie = 0, ostr& movie_stream = cout);
00041 #endif
00042 
00047 
00048 template <class field_item>
00049   matrix<field_item> inv(const matrix<field_item>& A);
00050 
00055 
00056 template <class field_item>
00057   field_item det(const matrix<field_item>& A);
00058 
00059 
00069 
00070 #ifdef no_default_arguments
00071 template <class int_item>
00072   void rref_int(const matrix<int_item>& A,
00073     matrix<int_item>& B, matrix<int_item>& Q,
00074     vector<index>& basis, int_item& det,
00075     int movie, ostr& movie_stream);
00076 #else
00077 template <class int_item>
00078   void rref_int(const matrix<int_item>& A,
00079     matrix<int_item>& B, matrix<int_item>& Q,
00080     vector<index>& basis, int_item& det,
00081     int movie = 0, ostr& movie_stream = cout);
00082 #endif
00083 
00084 
00085 
00089 
00090 template <class int_item>
00091   int_item det_int(const matrix<int_item>& A);
00092 
00093 
00094 
00104 
00105 #ifdef no_default_arguments
00106 template <class int_item>
00107   void smith(const matrix<int_item>& A,
00108     matrix<int_item>& B, matrix<int_item>& P,
00109     matrix<int_item>& Q, index& rank, int_item& det,
00110     int movie, ostr& movie_stream);
00111 #else
00112 template <class int_item>
00113   void smith(const matrix<int_item>& A,
00114     matrix<int_item>& B, matrix<int_item>& P,
00115     matrix<int_item>& Q, index& rank, int_item& det,
00116     int movie = 0, ostr& movie_stream = cout);
00117 #endif
00118 
00119 // I M P L E M E N T A T I O N
00120 
00121 
00122 #ifdef no_default_arguments
00123 template <class field_item>
00124   void rref(const matrix<field_item>& A,
00125     matrix<field_item>& B, matrix<field_item>& Q,
00126     vector<index>& basis, field_item& det,
00127     int movie, ostr& movie_stream)
00128 #else
00129 template <class field_item>
00130   void rref(const matrix<field_item>& A,
00131     matrix<field_item>& B, matrix<field_item>& Q,
00132     vector<index>& basis, field_item& det,
00133     int movie, ostr& movie_stream)
00134 #endif
00135 
00136 
00137 {
00138   det = 1;
00139   index m = A.get_m();
00140   index n = A.get_n();
00141   Q = eye(field_item(1), m);
00142   B = A;
00143   basis.resize(0);
00144 
00145   if (movie)
00146   {
00147     movie_stream << "Reduced row echelon form\n";
00148     print2(movie_stream, B, Q);
00149   }
00150 
00151   for (index i = 0, j = 0; i < m && j < n; j++)
00152   {
00153     // main loop:
00154     // for any cols j of B
00155     field_item max = 0;
00156     index i_max = i;
00157 
00158     // find the biggest (in absolute value) entry in the column j
00159     if (movie)
00160       movie_stream
00161         << "Find the biggest (in absolute value) entry in the column "
00162         << j << endl;
00163     for (index ii = i; ii < m; ii++)
00164     {
00165       if (absolute_value(B(ii, j)) > max)
00166       {
00167             max = absolute_value(B(ii, j));
00168             i_max = ii;
00169       }
00170     }
00171     if (max == field_item(0))
00172     {
00173       // of main loop: column j is negligible
00174       if (movie)
00175         movie_stream << "Column " << j << " is negligible\n";
00176       det = 0;
00177       continue;
00178     }
00179 
00180     if (movie)
00181       movie_stream << "Pivot item: (" << i_max << ", " << j << ")\n";
00182 
00183     if (i != i_max)
00184     {
00185       B.swap_rows(i, i_max);
00186       Q.swap_rows(i, i_max);
00187       det = -det;
00188       if (movie)
00189       {
00190         movie_stream << "Swap rows: " << i << ", " << i_max << endl;
00191         print2(movie_stream, B, Q);
00192       }
00193     }
00194 
00195     field_item pivot = B(i, j);
00196     B.div_row(i, pivot);
00197     Q.div_row(i, pivot);
00198     det = det * pivot;
00199 
00200     for (index ii = 0; ii < m; ii++)  // eliminate in column j
00201       if (ii != i)
00202       {
00203         field_item alpha = B(ii, j);
00204         B.add_row(ii, i, -alpha);
00205         Q.add_row(ii, i, -alpha);
00206       }
00207     if (movie)
00208     {
00209       movie_stream << "Elimination in column " << j << endl;
00210       print2(movie_stream, B, Q);
00211     }
00212 
00213     basis.join(j);
00214     i++;
00215   } // end of main loop
00216 }
00217 
00218 template <class field_item>
00219   field_item det(const matrix<field_item>& A)
00220 {
00221   // precondition: A must be square
00222   if (A.get_m() != A.get_n())
00223   {
00224     matrix_error (st_matrix_must_be_square, "matrix must be square");
00225     return 0;
00226   }
00227 
00228   matrix<field_item> B;
00229   matrix<field_item> Q;
00230   vector<index> basis;
00231 
00232   field_item det;
00233 
00234   rref(A, B, Q, basis, det, 0, cout);
00235 
00236   return det;
00237 }
00238 
00239 
00240 
00241 
00242 template <class field_item>
00243   matrix<field_item> inv(const matrix<field_item>& A)
00244 {
00245   // precondition: A must be square
00246   if (A.get_m() != A.get_n())
00247   {
00248     matrix_error (st_matrix_must_be_square, "matrix must be square");
00249     return 0;
00250   }
00251 
00252   matrix<field_item> B;
00253   matrix<field_item> Q;
00254   vector<index> basis;
00255 
00256   field_item det;
00257 
00258   rref(A, B, Q, basis, det, 0, cout);
00259 
00260   if (det == 0)
00261     matrix_error (st_matrix_is_singular, "matrix is singular");
00262 
00263   return Q;
00264 }
00265 
00266 
00267 #ifdef no_default_arguments
00268 template <class int_item>
00269   void rref_int(const matrix<int_item>& A,
00270     matrix<int_item>& B, matrix<int_item>& Q,
00271     vector<index>& basis, int_item& det,
00272     int movie, ostr& movie_stream)
00273 #else
00274 template <class int_item>
00275   void rref_int(const matrix<int_item>& A,
00276     matrix<int_item>& B, matrix<int_item>& Q,
00277     vector<index>& basis, int_item& det,
00278     int movie, ostr& movie_stream)
00279 #endif
00280 {
00281   det = 1;
00282   int_item det_denom = 1; // differs from field version
00283   index m = A.get_m();
00284   index n = A.get_n();
00285   Q = eye(int_item(1), m);
00286   B = A;
00287   basis.resize(0);
00288 
00289   if (movie)
00290   {
00291     movie_stream << "Reduced row echelon form for integral matrix\n";
00292     print2(movie_stream, B, Q);
00293   }
00294 
00295   for (index i = 0, j = 0; i < m && j < n; j++)
00296   {
00297     // main loop:
00298     // for any cols j of B
00299     int_item max = 0;
00300     index i_max = i;
00301 
00302     // find non-zero entry in the column j
00303     if (movie)
00304       movie_stream
00305         << "Find non-zero entry in the column "
00306         << j << endl;
00307     for (index ii = i; ii < m; ii++)
00308     {
00309       if (absolute_value(B(ii, j)) > max)
00310       {
00311         max = absolute_value(B(ii, j));
00312         i_max = ii;
00313       }
00314     }
00315     if (max == 0)
00316     {
00317       // of main loop: column j is negligible
00318       if (movie)
00319         movie_stream << "Column " << j << " is negligible\n";
00320       continue;
00321     }
00322 
00323     if (movie)
00324       movie_stream << "Pivot item: (" << i_max << ", " << j << ")\n";
00325 
00326     if (i != i_max)
00327     {
00328       B.swap_rows(i, i_max);
00329       Q.swap_rows(i, i_max);
00330       det = -det;
00331       if (movie)
00332       {
00333         movie_stream << "Swap rows: " << i << ", " << i_max << endl;
00334         print2(movie_stream, B, Q);
00335       }
00336     }
00337 
00338     int_item pivot = B(i, j);
00339 
00340     if (pivot < 0)      // differs from field version
00341     {                   // differs from field version
00342       B.mult_row(i, -1);// differs from field version
00343       Q.mult_row(i, -1);// differs from field version
00344       pivot = -pivot;   // differs from field version
00345       det = -det;       // differs from field version
00346     }                   // differs from field version
00347 
00348     for (index ii = 0; ii < m; ii++)    // eliminate in column j
00349       if (ii != i)
00350       {
00351         int_item alpha = B(ii, j);
00352         int_item delta = gcd(pivot, alpha);// differs from field version
00353         int_item b_ii = -alpha / delta;  // differs from field version
00354         int_item b_i = pivot / delta;  // differs from field version
00355         B.mult_row(ii, b_i);  // differs from field version
00356         B.add_row(ii, i, b_ii); // differs from field version
00357         Q.mult_row(ii, b_i);  // differs from field version
00358         Q.add_row(ii, i, b_ii); // differs from field version
00359         det_denom = det_denom * b_i;  // differs from field version
00360         alpha = gcd(gcd(B.row(ii)), gcd(Q.row(ii)));  // differs from field version
00361         B.div_row(ii, alpha);   // differs from field version
00362         Q.div_row(ii, alpha); // differs from field version
00363         det = det * alpha;  // differs from field version
00364       }
00365     if (movie)
00366     {
00367       movie_stream << "Elimination in column " << j << endl;
00368       print2(movie_stream, B, Q);
00369     }
00370 
00371 
00372     basis.join(j);
00373     i++;
00374   } // end of main loop
00375 
00376   for (index i = 0; i < basis.get_n(); i++) // differs from field version
00377     det = det * B(i, basis[i]); // differs from field version
00378   det = det/det_denom;  // differs from field version
00379 }
00380 
00381 
00382 template <class int_item>
00383   int_item det_int(const matrix<int_item>& A)
00384 {
00385   // precondition: A must be square
00386   if (A.get_m() != A.get_n())
00387   {
00388     matrix_error (st_matrix_must_be_square,
00389       "matrix must be square");
00390     return 0;
00391   }
00392 
00393   matrix<int_item> B;
00394   matrix<int_item> Q;
00395   vector<index> basis;
00396   int_item d;
00397 
00398   rref_int(A, B, Q, basis, d, 0, cout);
00399 
00400   return d;
00401 }
00402 
00403 #ifdef no_default_arguments
00404 template <class int_item>
00405   void smith(const matrix<int_item>& A,
00406     matrix<int_item>& B, matrix<int_item>& P,
00407     matrix<int_item>& Q, index& rank, int_item& det,
00408     int movie, ostr& movie_stream)
00409 #else
00410 template <class int_item>
00411   void smith(const matrix<int_item>& A,
00412     matrix<int_item>& B, matrix<int_item>& P,
00413     matrix<int_item>& Q, index& rank, int_item& det,
00414     int movie, ostr& movie_stream)
00415 #endif
00416 {
00417   det = 1;
00418   index m = A.get_m();
00419   index n = A.get_n();
00420   Q = eye(int_item(1), n);
00421   P = eye(int_item(1), m);
00422   B = A;
00423 
00424   if (movie)
00425   {
00426     movie_stream << "Smith's normal diagonal form\n";
00427     print3(movie_stream, Q, P, B);
00428   }
00429 
00430   for (index corner = 0; corner < m && corner < n; corner++)
00431   {
00432     index new_i, new_j, i, j;
00433     int_item pivot_item;
00434 
00435     if (!find_pivot_item(B, corner, new_i, new_j))
00436       // the smallest (in absolute value) non-zero entry in the matrix
00437       break; // all others entries are 0
00438 
00439     if (movie)
00440       movie_stream << "The smallest (in absolute value) non-zero entry: ("
00441         << new_i << ", " << new_j << ")\n";
00442 
00443     while(1)
00444     {
00445       do{
00446         i = new_i;
00447         j = new_j;
00448         pivot_item = B(i, j);
00449 
00450         if (movie)
00451           movie_stream << "Pivoting item: (" << i << ", " << j << ")\n";
00452 
00453         for (index jj = 0; jj < n; jj++)
00454           if (jj != j)
00455           {
00456             int_item p = quotient(B(i, jj), pivot_item);
00457             B.add_col(jj, j, -p);
00458             Q.add_col(jj, j, -p);
00459           }
00460         for (index ii = 0; ii < m; ii++)
00461           if (ii != i)
00462           {
00463             int_item p = quotient(B(ii, j), pivot_item);
00464             B.add_row(ii, i, -p);
00465             P.add_row(ii, i, -p);
00466           }
00467 
00468         if (movie)
00469         {
00470           movie_stream << "After pivoting step:\n";
00471           print3(movie_stream, Q, P, B);
00472         }
00473 
00474         new_j = find_pivot_item_in_row(B, corner, i);
00475         if(new_j == j)
00476           new_i = find_pivot_item_in_col(B, corner, j);
00477         else
00478           new_i = i;
00479       }while (i != new_i || j != new_j);
00480         // while pivot entry is changing,
00481         // i.e. until others entries in row i and column j are 0
00482 
00483       index k, l;
00484       if (!find_non_divisor_of_item(B, corner, pivot_item, k, l))
00485         break;
00486 
00487       B.add_row(i, k, 1);
00488       P.add_row(i, k, 1);
00489 
00490       if (movie)
00491       {
00492         movie_stream << "Non-divisor entry: (" << k << ", " << l << ")\n";
00493         movie_stream << "Add the " << k << "-th row to the " << i
00494           << "-th row\n";
00495         print3(movie_stream, Q, P, B);
00496       }
00497 
00498     }
00499 
00500     if (i != corner)
00501     {
00502       B.swap_rows(i, corner);
00503       P.swap_rows(i, corner);
00504       det = -det;
00505     }
00506     if (j != corner)
00507     {
00508       B.swap_cols(j, corner);
00509       Q.swap_cols(j, corner);
00510       det = -det;
00511     }
00512     if (B(corner, corner) < 0)
00513     {
00514       B.mult_row(corner, -1);
00515       P.mult_row(corner, -1);
00516       det = -det;
00517     }
00518     if (movie)
00519       print3(movie_stream, Q, P, B);
00520 
00521   }
00522 
00523   rank = 0;
00524 
00525   for (index i = 0; i < m && i < n; i++)
00526     if(B(i, i) != 0)
00527     {
00528       det = det * B(i, i);
00529       rank++;
00530     }
00531     else
00532       break;
00533 }
00534 
00535 template <class int_item>
00536   int find_pivot_item(matrix<int_item>& B, index corner,
00537     index& min_i, index& min_j)
00538 {
00539   int_item min = 0;
00540   for (index i = corner; i < B.get_m(); i++)
00541     for (index j = corner; j < B.get_n(); j++)
00542       if (B(i, j) != 0 && (min == 0 || absolute_value(B(i, j)) < min))
00543       {
00544         min = absolute_value (B(i, j));
00545         min_i = i;
00546         min_j = j;
00547       }
00548   if (min == 0)
00549     return 0;
00550   return 1;
00551 }
00552 
00553 
00554 template <class int_item>
00555   index find_pivot_item_in_row(matrix<int_item>& B, index corner,
00556     index i)
00557 {
00558   int_item min = 0;
00559   index min_j = 0;
00560 
00561   for (index j = corner; j < B.get_n(); j++)
00562     if (B(i, j) != 0 && (min == 0 || absolute_value(B(i, j)) < min))
00563     {
00564       min = absolute_value (B(i, j));
00565       min_j = j;
00566     }
00567 
00568   return min_j;
00569 }
00570 
00571 template <class int_item>
00572   index find_pivot_item_in_col(matrix<int_item>& B, index corner,
00573     index j)
00574 {
00575   int_item min = 0;
00576   index min_i = 0;
00577 
00578   for (index i = corner; i < B.get_m(); i++)
00579     if (B(i, j) != 0 && (min == 0 || absolute_value(B(i, j)) < min))
00580     {
00581       min = absolute_value (B(i, j));
00582       min_i = i;
00583     }
00584 
00585   return min_i;
00586 }
00587 
00588 template <class int_item>
00589   int find_non_divisor_of_item(matrix<int_item>& B, index corner,
00590     int_item item, index& k, index& l)
00591 {
00592   for (index i = corner; i < B.get_m(); i++)
00593     for (index j = corner; j < B.get_n(); j++)
00594       if (B(i, j) % item != 0)
00595       {
00596         k = i;
00597         l = j;
00598         return 1;
00599       }
00600   return 0;
00601 }
00602 
00603 
00604 
00605 #endif GAUSS_H_

Generated at Tue Jan 22 20:37:04 2002 for Arageli by doxygen1.2.9.1 written by Dimitri van Heesch, © 1997-2001