Bugzilla – Attachment 215 Details for
Bug 334
Implement LU, ILU0 and solve in SparseExtra
Home
|
New
|
Browse
|
Search
|
[?]
|
Reports
|
Requests
|
Help
|
Log In
[x]
|
Forgot Password
Login:
[x]
This bugzilla service is closed. All entries have been migrated to
https://gitlab.com/libeigen/eigen
[patch]
SparseLU.h
SparseLU_path.diff (text/plain), 4.41 KB, created by
Donald McGillavry
on 2011-09-30 12:42:55 UTC
(
hide
)
Description:
SparseLU.h
Filename:
MIME Type:
Creator:
Donald McGillavry
Created:
2011-09-30 12:42:55 UTC
Size:
4.41 KB
patch
obsolete
>--- SparseLU_original.h 2011-05-30 19:15:37.000000000 +0600 >+++ SparseLU.h 2011-09-30 17:43:34.199611500 +0600 >@@ -123,7 +123,7 @@ > > template<typename BDerived, typename XDerived> > bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, >- const int transposed = SvNoTrans) const; >+ const int transposed=SvNoTrans, const int CompleteSolve=0) const; > > /** \returns true if the factorization succeeded */ > inline bool succeeded(void) const { return m_succeeded; } >@@ -133,15 +133,64 @@ > int m_flags; > mutable int m_status; > bool m_succeeded; >+ mutable SparseMatrix<Scalar,ColMajor> m_matrix; > }; > > /** Computes / recomputes the LU decomposition of matrix \a a > * using the default algorithm. > */ > template<typename _MatrixType, typename Backend> >-void SparseLU<_MatrixType,Backend>::compute(const _MatrixType& ) >+void SparseLU<_MatrixType,Backend>::compute(const _MatrixType& A) > { >- eigen_assert(false && "not implemented yet"); >+ if (A.cols()!=A.rows()) {m_succeeded=false; return;} >+ //Zero Fill-in ILU (ILU(0)) >+ //algorithm jki-version >+ //Saad Y. Iterative Methods for Sparse Linear Systems. (SIAM, Philadelphia, 2003) >+ //optimized for Compressed Sparse Column-major format >+ if (m_flags&IncompleteFactorization){ >+ m_matrix=A; >+ //Size of matrix >+ const int n=m_matrix.cols(); >+ int *indexDiagonal=new int[n]; >+ int *outerIndex=m_matrix._outerIndexPtr(),*innerIndex=m_matrix._innerIndexPtr(); >+ Scalar *matrixValue=m_matrix._valuePtr(); >+ //Search outer indexes of diagonal elements >+ for (int j=0;j<n;j++){ >+ int i=outerIndex[j],endColJ=outerIndex[j+1]; >+ if(i==endColJ) {m_succeeded=false; return;} >+ while (i<endColJ && innerIndex[i]<j)i++; >+ if (innerIndex[i]==j)indexDiagonal[j]=i; >+ else {m_succeeded=false; return;} >+ } >+ //A(i,0)/=A(0,0) >+ Scalar diagonalElement=matrixValue[indexDiagonal[0]]; >+ for (int i=indexDiagonal[0]+1; i<outerIndex[1]; i++) >+ matrixValue[i]/=diagonalElement; >+ //A(i,j)-=A(i,k)*A(k,j) >+ for (int j=1;j<n;j++){ >+ int endColJ=outerIndex[j+1]; >+ for (int k=outerIndex[j];k<indexDiagonal[j];k++){ >+ int colJ=k+1, colK=outerIndex[innerIndex[k]],\ >+ endColK=outerIndex[innerIndex[k]+1]; >+ while (colK<endColK && colJ<endColJ){ >+ if(innerIndex[colK]==innerIndex[colJ]){ >+ matrixValue[colJ]-=matrixValue[k]*matrixValue[colK]; >+ colJ++; >+ colK++; >+ } >+ else if(innerIndex[colK]>innerIndex[colJ])colJ++; >+ else colK++; >+ } >+ } >+ //A(i,j)/=A(j,j) >+ diagonalElement=matrixValue[indexDiagonal[j]]; >+ for (int i=indexDiagonal[j]+1; i<outerIndex[j+1]; i++) >+ matrixValue[i]/=diagonalElement; >+ } >+ delete [] indexDiagonal; >+ m_succeeded=true; >+ } >+ else eigen_assert(false && "not implemented yet"); > } > > /** Computes *x = U^-1 L^-1 b >@@ -154,10 +203,36 @@ > */ > template<typename _MatrixType, typename Backend> > template<typename BDerived, typename XDerived> >-bool SparseLU<_MatrixType,Backend>::solve(const MatrixBase<BDerived> &, MatrixBase<XDerived>* , const int ) const >+bool SparseLU<_MatrixType,Backend>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> *x, const int transposed, const int CompleteSolve) const > { > eigen_assert(false && "not implemented yet"); > return false; >+ if (m_succeeded==false) return false; >+ //If the ILU(0) factorization >+ if (m_flags&IncompleteFactorization) >+ { >+ *x=b; >+ if (transposed==SvNoTrans){ >+ // Solve L*y=b >+ if (CompleteSolve==0 || CompleteSolve&Lower) *x=m_matrix.triangularView<UnitLower>().solve(*x); >+ // Solve U*x=y >+ if (CompleteSolve==0 || CompleteSolve&Upper) *x=m_matrix.triangularView<Upper>().solve(*x); >+ } >+ if (transposed==SvTranspose){ >+ // Solve L^T*y=b >+ if (CompleteSolve==0 || CompleteSolve&Lower) *x=m_matrix.transpose().triangularView<UnitLower>().solve(*x); >+ // Solve U^T*x=y >+ if (CompleteSolve==0 || CompleteSolve&Upper) *x=m_matrix.transpose().triangularView<Upper>().solve(*x); >+ } >+ if (transposed==SvAdjoint) { >+ // Solve L^T*y=b >+ if (CompleteSolve==0 || CompleteSolve&Lower) *x=m_matrix.adjoint().triangularView<UnitLower>().solve(*x); >+ // Solve U^T*x=y >+ if (CompleteSolve==0 || CompleteSolve&Upper) *x=m_matrix.adjoint().triangularView<Upper>().solve(*x); >+ } >+ return true; >+ } >+ else {eigen_assert(false && "not implemented yet"); return false;} > } > > #endif // EIGEN_SPARSELU_H
You cannot view the attachment while viewing its details because your browser does not support IFRAMEs.
View the attachment on a separate page
.
View Attachment As Diff
View Attachment As Raw
Actions:
View
|
Diff
Attachments on
bug 334
: 215