Changeset 4073


Ignore:
Timestamp:
11/19/11 10:20:35 (3 years ago)
Author:
saunders
Message:

revision of 2local smith form due to memory issues in blas-matrix

Location:
trunk/linbox
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/linbox/linbox/algorithms/smith-form-local2.inl

    r4013 r4073  
    3434        { 
    3535        public: 
    36                 typedef Local2_32 LocalPID; 
    37                 typedef LocalPID::Element Elt; 
     36                typedef Local2_32 LocalPIR; 
     37                typedef LocalPIR::Element Elt; 
    3838 
    3939                template<class Matrix> 
    40                 std::list<Elt>& operator()(std::list<Elt>& L, Matrix& A, const LocalPID& R) 
     40                std::list<Elt>& operator()(std::list<Elt>& L, Matrix& A, const LocalPIR& R) 
    4141                {   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); 
    4344                } 
    4445 
    45                 template<class Matrix> 
    4646                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) 
    4848                { 
    49                         if ( A.rowdim() == 0 || A.coldim() == 0 ) 
     49                        if ( m == 0 || n == 0 ) 
    5050                                return L; 
    5151 
    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) 
    5657                        { 
    57                                 for (q = p->begin(); q != p->end(); ++q) 
     58                                for (j = 0; j != n; ++j) 
    5859                                { 
    59                                         R.gcdin(g, *q); 
     60                                        R.gcdin(g, Ap[i*stride + j]); 
    6061                                        if ( R.isUnit(g) ) break; 
    6162                                } 
    6263                                if ( R.isUnit(g) ) break; 
    6364                        } 
    64                         //std::cout << "g = " << (int)g <<"\n"; 
    6565                        if ( R.isZero(g) ) 
    6666                        { 
    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); 
    7368                                return L; 
    7469                        } 
    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) 
    7890                                { 
    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] ); 
    10794                                } 
     95                                L.push_back(d); 
     96                                return smithStep(L, d, Ap + stride+1,m-1, n-1, stride, R); 
     97                        } 
    10898                        else 
    10999                        { 
    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); 
    114106                        } 
    115107                } 
  • trunk/linbox/tests/test-dense.C

    r4051 r4073  
    5454// using long on purpose 
    5555template <class Field> 
    56 static bool testIdentity (Field &F, long n, int iterations) 
     56static bool testIdentity (Field &F, size_t n, int iterations) 
    5757{ 
    5858        typedef typename Vector<Field>::Dense Vector; 
     
    7474        F.init (one, 1); 
    7575 
    76         for (long i = 0; i < n; i++) 
     76        for (size_t i = 0; i < n; i++) 
    7777                I.setEntry (i, i, one); 
    7878 
     
    8787                iter_passed = true; 
    8888 
    89                 for (long j = 0; j < n; j++) 
     89                for (size_t j = 0; j < n; j++) 
    9090                        r.random (v[j]); 
    9191 
Note: See TracChangeset for help on using the changeset viewer.