/[ff3d]/ff3d/algebra/TinyMatrix.hpp
ViewVC logotype

Diff of /ff3d/algebra/TinyMatrix.hpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.4 by delpinux, Sun Mar 9 19:59:45 2003 UTC revision 1.5 by delpinux, Tue Mar 11 17:43:27 2003 UTC
# Line 265  public: Line 265  public:
265    {    {
266      TinyMatrix<N,N,T> M;      TinyMatrix<N,N,T> M;
267      TinyMatrix<N,N,T> lu;      TinyMatrix<N,N,T> lu;
268      LU(*this, lu);      TinyMatrix<N,N,T> P;
269        LU(*this, lu, P);
270    
271      TinyVector<N,T> v;      TinyVector<N,T> v;
272      TinyVector<N,T> u;      TinyVector<N,T> u;
# Line 274  public: Line 275  public:
275        for (size_t j=0; j<N; ++j)        for (size_t j=0; j<N; ++j)
276          v[j] = (i==j)?1:0;          v[j] = (i==j)?1:0;
277    
278        u = lu.applyLU(v);        u = lu.applyLU(P*v);
279        for (size_t j=0; j<N; ++j)        for (size_t j=0; j<N; ++j)
280          M(j,i) = u[j];          M(j,i) = u[j];
281      }      }
# Line 520  public: Line 521  public:
521      }      }
522    }    }
523    
524    /*! Returns the LU factorisation of a matrix    /*! Computes P, L and U such that PA=LU.
525        L and U are stored in a single matrix. L non zero coefficients are i<j.
526      \remark As only sens for square matrix      \remark As only sens for square matrix
527    */    */
528    friend void LU(const TinyMatrix<N,N,T>& A,    friend void LU(const TinyMatrix<N,N,T>& A,
529                   TinyMatrix<N,N,T>& lu)                   TinyMatrix<N,N,T>& lu,
530                     TinyMatrix<N,N,T>& P)
531    {    {
532        // Reinitialize P to Identity
533        for (size_t i=0; i<N; i++)
534          for (size_t j=0; j<N; j++) {
535            P(i,j) = (i==j);
536          }
537    
538      lu = A;      lu = A;
539            
540      for (int k=0; k<N-1; k++) {      for (int k=0; k<N-1; k++) {
541  #warning Do not proceed to swaping anymore when needed.        // Checkes for Pivot.
542  //       int maxLine = k;        int maxLine = k;
543  //       for (int i=k+1; i<N; i++) {        for (int i=k+1; i<N; i++) {
544  //      maxLine = (std::abs(lu(i,k))<std::abs(lu(maxLine,k)))?maxLine:i;          maxLine = (std::abs(lu(i,k))<std::abs(lu(maxLine,k)))?maxLine:i;
545  //       }        }
546  //       lu.swapLines(k,maxLine);        lu.swapLines(k,maxLine);
547          P.swapLines(k,maxLine);
548    
549        for (int i=k+1; i<N; i++) {        for (int i=k+1; i<N; i++) {
550          lu(i,k) /= lu(k,k);          lu(i,k) /= lu(k,k);
551          for (int j=k+1; j<N; j++)          for (int j=k+1; j<N; j++)
# Line 575  public: Line 586  public:
586                                     const TinyMatrix<N,N,T>& A)                                     const TinyMatrix<N,N,T>& A)
587    {    {
588      TinyMatrix<N,N,T> lu;      TinyMatrix<N,N,T> lu;
589      LU(A,lu);      TinyMatrix<N,N,T> P;
590        LU(A,lu,P);
591    
592      return lu.applyLU(b);      return lu.applyLU(P*b);
593    }    }
594    
595    //! Constructor does nothing.    //! Constructor does nothing.

Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26