Changeset 4073
 Timestamp:
 11/19/11 10:20:35 (3 years ago)
 Location:
 trunk/linbox
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/linbox/linbox/algorithms/smithformlocal2.inl
r4013 r4073 34 34 { 35 35 public: 36 typedef Local2_32 LocalPI D;37 typedef LocalPI D::Element Elt;36 typedef Local2_32 LocalPIR; 37 typedef LocalPIR::Element Elt; 38 38 39 39 template<class Matrix> 40 std::list<Elt>& operator()(std::list<Elt>& L, Matrix& A, const LocalPI D& R)40 std::list<Elt>& operator()(std::list<Elt>& L, Matrix& A, const LocalPIR& R) 41 41 { Elt d; R.init(d, 1); 42 return smithStep(L, d, A, R); 42 Elt *p = &(A[0][0]); 43 return smithStep(L, d, p, A.rowdim(), A.coldim(), A.getStride(), R); 43 44 } 44 45 45 template<class Matrix>46 46 std::list<Elt>& 47 smithStep(std::list<Elt>& L, Elt& d, Matrix& A, const LocalPID& R)47 smithStep(std::list<Elt>& L, Elt& d, Elt* Ap, size_t m, size_t n, size_t stride, const LocalPIR& R) 48 48 { 49 if ( A.rowdim() == 0  A.coldim()== 0 )49 if ( m == 0  n == 0 ) 50 50 return L; 51 51 52 LocalPID::Exponent g = LocalPID::Exponent(32); //R.init(g, 0); // must change to 2^31 maybe. 53 typename Matrix::RowIterator p; 54 typename Matrix::Row::iterator q, r; 55 for ( p = A.rowBegin(); p != A.rowEnd(); ++p) 52 LocalPIR::Exponent g = LocalPIR::Exponent(32); //R.init(g, 0); // must change to 2^31 maybe. 53 size_t i, j, k; 54 /* Arguably this search order should be reversed to increase the likelyhood of no col swap, 55 assuming row swaps cheaper. Not so, however on my example. bds 11Nov */ 56 for ( i = 0; i != m; ++i) 56 57 { 57 for ( q = p>begin(); q != p>end(); ++q)58 for (j = 0; j != n; ++j) 58 59 { 59 R.gcdin(g, *q);60 R.gcdin(g, Ap[i*stride + j]); 60 61 if ( R.isUnit(g) ) break; 61 62 } 62 63 if ( R.isUnit(g) ) break; 63 64 } 64 //std::cout << "g = " << (int)g <<"\n";65 65 if ( R.isZero(g) ) 66 66 { 67 // std::cout << " R.isZero(g) is used\n"; 68 // std::cout << A.rowdim() << " " << A.coldim() << "\n"; 69 L.insert(L.end(), 70 (A.rowdim() < A.coldim()) ? A.rowdim() : A.coldim(), 71 0 72 ); 67 L.insert(L.end(), (m < n) ? m : n, 0); 73 68 return L; 74 69 } 75 if ( p != A.rowEnd() ) // g is a unit and, 76 // because this is a local ring, value at which this first happened 77 // also is a unit. 70 if ( i != m ) // g is a unit and, because this is a local ring, 71 // value at which this first happened also is a unit. 72 { // put pivot in 0,0 position 73 if ( i != 0 ) // swap rows 74 std::swap_ranges(Ap, Ap+n, Ap + i*stride); 75 if ( j != 0 ) // swap cols 76 for(k = 0; k != m; ++k) 77 std::swap(Ap[k*stride + 0], Ap[k*stride + j]); 78 79 // elimination step  crude and for dense only  fix later 80 // Want to use a block method or "left looking" elimination. 81 Elt f; R.inv(f, Ap[0*stride + 0] ); 82 R.negin(f); 83 84 // normalize first row to 1, ... 85 for ( j = 0; j != n; ++j) 86 R.mulin(Ap[0*stride + j], f); 87 88 // eliminate in subsequent rows 89 for ( i = 1; i != m; ++i) 78 90 { 79 if ( p != A.rowBegin() ) 80 swap_ranges(A.rowBegin()>begin(), A.rowBegin()>end(), 81 p>begin()); 82 if ( q != p>begin() ) 83 swap_ranges(A.colBegin()>begin(), A.colBegin()>end(), 84 (A.colBegin() + (int)(q  p>begin()))>begin()); 85 86 // eliminate step  crude and for dense only  fix later 87 // Want to use a block method or "left looking" elimination. 88 //std::cout << " Value of A[0][0]: " << *(A.rowBegin() > begin()) <<"\n"; 89 Elt f; R.inv(f, *(A.rowBegin()>begin() ) ); 90 R.negin(f); 91 // normalize first row to 1, ... 92 //std::cout << "f = " << f << "\n"; 93 //A.write(std::cout); 94 95 for ( q = A.rowBegin()>begin() /*+ 1*/; q != A.rowBegin()>end(); ++q) 96 R.mulin(*q, f); 97 // 98 // eliminate in subsequent rows 99 for ( p = A.rowBegin() + 1; p != A.rowEnd(); ++p) 100 for ( q = p>begin() + 1, r = A.rowBegin()>begin() + 1, f = *(p > begin()); 101 q != p>end(); ++q, ++r ) 102 R.axpyin( *q, f, *r ); 103 104 BlasMatrix<LocalPID> Ap(A, 1, 1, A.rowdim()  1, A.coldim()  1); 105 L.push_back(d); 106 return smithStep(L, d, Ap, R); 91 f = Ap[i*stride + 0]; 92 for ( j = 0; j != n; ++j) 93 R.axpyin( Ap[i*stride +j], f, Ap[0*stride +j] ); 107 94 } 95 L.push_back(d); 96 return smithStep(L, d, Ap + stride+1,m1, n1, stride, R); 97 } 108 98 else 109 99 { 110 typename Matrix::Iterator p_it; 111 for (p_it = A.Begin(); p_it != A.End(); ++p_it) 112 R.divin(*p_it, g); 113 return smithStep(L, R.mulin(d, g), A, R); 100 for ( i = 0; i != m; ++i) 101 for ( j = 0; j != n; ++j) 102 { 103 R.divin(Ap[i*stride + j], g); 104 } 105 return smithStep(L, R.mulin(d, g), Ap, m, n, stride, R); 114 106 } 115 107 } 
trunk/linbox/tests/testdense.C
r4051 r4073 54 54 // using long on purpose 55 55 template <class Field> 56 static bool testIdentity (Field &F, longn, int iterations)56 static bool testIdentity (Field &F, size_t n, int iterations) 57 57 { 58 58 typedef typename Vector<Field>::Dense Vector; … … 74 74 F.init (one, 1); 75 75 76 for ( longi = 0; i < n; i++)76 for (size_t i = 0; i < n; i++) 77 77 I.setEntry (i, i, one); 78 78 … … 87 87 iter_passed = true; 88 88 89 for ( longj = 0; j < n; j++)89 for (size_t j = 0; j < n; j++) 90 90 r.random (v[j]); 91 91
Note: See TracChangeset
for help on using the changeset viewer.