aboutsummaryrefslogtreecommitdiffhomepage
path: root/Eigen
diff options
context:
space:
mode:
Diffstat (limited to 'Eigen')
-rw-r--r--Eigen/src/LU/LU.h5
-rw-r--r--Eigen/src/Sparse/SparseLLT.h2
-rw-r--r--Eigen/src/Sparse/SparseMatrix.h4
-rw-r--r--Eigen/src/Sparse/SuperLUSupport.h191
-rw-r--r--Eigen/src/Sparse/UmfPackSupport.h146
5 files changed, 301 insertions, 47 deletions
diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h
index 554d8bd3c..fda369499 100644
--- a/Eigen/src/LU/LU.h
+++ b/Eigen/src/LU/LU.h
@@ -190,10 +190,7 @@ template<typename MatrixType> class LU
* \sa MatrixBase::solveTriangular(), kernel(), computeKernel(), inverse(), computeInverse()
*/
template<typename OtherDerived, typename ResultType>
- bool solve(
- const MatrixBase<OtherDerived>& b,
- ResultType *result
- ) const;
+ bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
/** \returns the determinant of the matrix of which
* *this is the LU decomposition. It has only linear complexity
diff --git a/Eigen/src/Sparse/SparseLLT.h b/Eigen/src/Sparse/SparseLLT.h
index b7d4f5bbd..f9eb459dc 100644
--- a/Eigen/src/Sparse/SparseLLT.h
+++ b/Eigen/src/Sparse/SparseLLT.h
@@ -160,7 +160,7 @@ void SparseLLT<MatrixType,Backend>::compute(const MatrixType& a)
{
Scalar y = it.value();
x -= ei_abs2(y);
- ++it; // skip j-th element, and process remaing column coefficients
+ ++it; // skip j-th element, and process remaining column coefficients
tempVector.restart();
for (; it; ++it)
{
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index 74e81f4bb..1346b310a 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -189,6 +189,10 @@ class SparseMatrix
m_outerSize = outerSize;
}
}
+ void resizeNonZeros(int size)
+ {
+ m_data.resize(size);
+ }
inline SparseMatrix()
: m_outerSize(0), m_innerSize(0), m_outerIndex(0)
diff --git a/Eigen/src/Sparse/SuperLUSupport.h b/Eigen/src/Sparse/SuperLUSupport.h
index 2ccfe2267..aeab835aa 100644
--- a/Eigen/src/Sparse/SuperLUSupport.h
+++ b/Eigen/src/Sparse/SuperLUSupport.h
@@ -211,6 +211,10 @@ class SparseLU<MatrixType,SuperLU> : public SparseLU<MatrixType>
typedef typename Base::Scalar Scalar;
typedef typename Base::RealScalar RealScalar;
typedef Matrix<Scalar,Dynamic,1> Vector;
+ typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
+ typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
+ typedef SparseMatrix<Scalar,Lower|UnitDiagBit> LMatrixType;
+ typedef SparseMatrix<Scalar,Upper> UMatrixType;
using Base::m_flags;
using Base::m_status;
@@ -231,23 +235,59 @@ class SparseLU<MatrixType,SuperLU> : public SparseLU<MatrixType>
{
}
+ inline const LMatrixType& matrixL() const
+ {
+ if (m_extractedDataAreDirty) extractData();
+ return m_l;
+ }
+
+ inline const UMatrixType& matrixU() const
+ {
+ if (m_extractedDataAreDirty) extractData();
+ return m_u;
+ }
+
+ inline const IntColVectorType& permutationP() const
+ {
+ if (m_extractedDataAreDirty) extractData();
+ return m_p;
+ }
+
+ inline const IntRowVectorType& permutationQ() const
+ {
+ if (m_extractedDataAreDirty) extractData();
+ return m_q;
+ }
+
+ Scalar determinant() const;
+
template<typename BDerived, typename XDerived>
bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const;
void compute(const MatrixType& matrix);
protected:
- // cached data to reduce reallocation:
+
+ void extractData() const;
+
+ protected:
+ // cached data to reduce reallocation, etc.
+ mutable LMatrixType m_l;
+ mutable UMatrixType m_u;
+ mutable IntColVectorType m_p;
+ mutable IntRowVectorType m_q;
+
mutable SparseMatrix<Scalar> m_matrix;
mutable SluMatrix m_sluA;
- mutable SuperMatrix m_sluL, m_sluU,;
+ mutable SuperMatrix m_sluL, m_sluU;
mutable SluMatrix m_sluB, m_sluX;
mutable SuperLUStat_t m_sluStat;
mutable superlu_options_t m_sluOptions;
- mutable std::vector<int> m_sluEtree, m_sluPermR, m_sluPermC;
+ mutable std::vector<int> m_sluEtree;
mutable std::vector<RealScalar> m_sluRscale, m_sluCscale;
mutable std::vector<RealScalar> m_sluFerr, m_sluBerr;
mutable char m_sluEqued;
+ mutable bool m_extractedDataAreDirty;
};
template<typename MatrixType>
@@ -261,6 +301,7 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
m_sluOptions.PrintStat = NO;
m_sluOptions.ConditionNumber = NO;
m_sluOptions.Trans = NOTRANS;
+ // m_sluOptions.Equil = NO;
switch (Base::orderingMethod())
{
@@ -279,8 +320,8 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
m_sluEqued = 'B';
int info = 0;
- m_sluPermR.resize(size);
- m_sluPermC.resize(size);
+ m_p.resize(size);
+ m_q.resize(size);
m_sluRscale.resize(size);
m_sluCscale.resize(size);
m_sluEtree.resize(size);
@@ -298,7 +339,7 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
m_sluX = m_sluB;
StatInit(&m_sluStat);
- SuperLU_gssvx(&m_sluOptions, &m_sluA, &m_sluPermC[0], &m_sluPermR[0], &m_sluEtree[0],
+ SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
&m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
&m_sluL, &m_sluU,
NULL, 0,
@@ -308,26 +349,12 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
&m_sluStat, &info, Scalar());
StatFree(&m_sluStat);
+ m_extractedDataAreDirty = true;
+
// FIXME how to better check for errors ???
Base::m_succeeded = (info == 0);
}
-// template<typename MatrixType>
-// inline const MatrixType&
-// SparseLU<MatrixType,SuperLU>::matrixL() const
-// {
-// ei_assert(false && "matrixL() is Not supported by the SuperLU backend");
-// return m_matrix;
-// }
-//
-// template<typename MatrixType>
-// inline const MatrixType&
-// SparseLU<MatrixType,SuperLU>::matrixU() const
-// {
-// ei_assert(false && "matrixU() is Not supported by the SuperLU backend");
-// return m_matrix;
-// }
-
template<typename MatrixType>
template<typename BDerived,typename XDerived>
bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> *x) const
@@ -349,7 +376,7 @@ bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, MatrixBa
RealScalar recip_pivot_gross, rcond;
SuperLU_gssvx(
&m_sluOptions, &m_sluA,
- &m_sluPermC[0], &m_sluPermR[0],
+ m_q.data(), m_p.data(),
&m_sluEtree[0], &m_sluEqued,
&m_sluRscale[0], &m_sluCscale[0],
&m_sluL, &m_sluU,
@@ -363,4 +390,122 @@ bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, MatrixBa
return info==0;
}
+//
+// the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
+//
+// Copyright (c) 1994 by Xerox Corporation. All rights reserved.
+//
+// THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
+// EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
+//
+template<typename MatrixType>
+void SparseLU<MatrixType,SuperLU>::extractData() const
+{
+ if (m_extractedDataAreDirty)
+ {
+ int upper;
+ int fsupc, istart, nsupr;
+ int lastl = 0, lastu = 0;
+ SCformat *Lstore = static_cast<SCformat*>(m_sluL.Store);
+ NCformat *Ustore = static_cast<NCformat*>(m_sluU.Store);
+ Scalar *SNptr;
+
+ const int size = m_matrix.rows();
+ m_l.resize(size,size);
+ m_l.resizeNonZeros(Lstore->nnz);
+ m_u.resize(size,size);
+ m_u.resizeNonZeros(Ustore->nnz);
+
+ int* Lcol = m_l._outerIndexPtr();
+ int* Lrow = m_l._innerIndexPtr();
+ Scalar* Lval = m_l._valuePtr();
+
+ int* Ucol = m_u._outerIndexPtr();
+ int* Urow = m_u._innerIndexPtr();
+ Scalar* Uval = m_u._valuePtr();
+
+ Ucol[0] = 0;
+ Ucol[0] = 0;
+
+ /* for each supernode */
+ for (int k = 0; k <= Lstore->nsuper; ++k)
+ {
+ fsupc = L_FST_SUPC(k);
+ istart = L_SUB_START(fsupc);
+ nsupr = L_SUB_START(fsupc+1) - istart;
+ upper = 1;
+
+ /* for each column in the supernode */
+ for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
+ {
+ SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
+
+ /* Extract U */
+ for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
+ {
+ Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
+ /* Matlab doesn't like explicit zero. */
+ if (Uval[lastu] != 0.0)
+ Urow[lastu++] = U_SUB(i);
+ }
+ for (int i = 0; i < upper; ++i)
+ {
+ /* upper triangle in the supernode */
+ Uval[lastu] = SNptr[i];
+ /* Matlab doesn't like explicit zero. */
+ if (Uval[lastu] != 0.0)
+ Urow[lastu++] = L_SUB(istart+i);
+ }
+ Ucol[j+1] = lastu;
+
+ /* Extract L */
+ Lval[lastl] = 1.0; /* unit diagonal */
+ Lrow[lastl++] = L_SUB(istart + upper - 1);
+ for (int i = upper; i < nsupr; ++i)
+ {
+ Lval[lastl] = SNptr[i];
+ /* Matlab doesn't like explicit zero. */
+ if (Lval[lastl] != 0.0)
+ Lrow[lastl++] = L_SUB(istart+i);
+ }
+ Lcol[j+1] = lastl;
+
+ ++upper;
+ } /* for j ... */
+
+ } /* for k ... */
+
+ // squeeze the matrices :
+ m_l.resizeNonZeros(lastl);
+ m_u.resizeNonZeros(lastu);
+
+ m_extractedDataAreDirty = false;
+ }
+}
+
+template<typename MatrixType>
+typename SparseLU<MatrixType,SuperLU>::Scalar SparseLU<MatrixType,SuperLU>::determinant() const
+{
+ if (m_extractedDataAreDirty)
+ extractData();
+
+ // TODO this code coule be moved to the default/base backend
+ // FIXME perhaps we have to take into account the scale factors m_sluRscale and m_sluCscale ???
+ Scalar det = Scalar(1);
+ for (int j=0; j<m_u.cols(); ++j)
+ {
+ if (m_u._outerIndexPtr()[j+1]-m_u._outerIndexPtr()[j] > 0)
+ {
+ int lastId = m_u._outerIndexPtr()[j+1]-1;
+ ei_assert(m_u._innerIndexPtr()[lastId]<=j);
+ if (m_u._innerIndexPtr()[lastId]==j)
+ {
+ det *= m_u._valuePtr()[lastId];
+ }
+ }
+ // std::cout << m_sluRscale[j] << " " << m_sluCscale[j] << " ";
+ }
+ return det;
+}
+
#endif // EIGEN_SUPERLUSUPPORT_H
diff --git a/Eigen/src/Sparse/UmfPackSupport.h b/Eigen/src/Sparse/UmfPackSupport.h
index 751e390f6..00fa7e57d 100644
--- a/Eigen/src/Sparse/UmfPackSupport.h
+++ b/Eigen/src/Sparse/UmfPackSupport.h
@@ -25,7 +25,21 @@
#ifndef EIGEN_UMFPACKSUPPORT_H
#define EIGEN_UMFPACKSUPPORT_H
-/* TODO extract L, extrac U, compute det, etc... */
+/* TODO extract L, extract U, compute det, etc... */
+
+// generic double/complex<double> wrapper functions:
+
+inline void umfpack_free_numeric(void **Numeric, double)
+{ umfpack_di_free_numeric(Numeric); }
+
+inline void umfpack_free_numeric(void **Numeric, std::complex<double>)
+{ umfpack_zi_free_numeric(Numeric); }
+
+inline void umfpack_free_symbolic(void **Symbolic, double)
+{ umfpack_di_free_symbolic(Symbolic); }
+
+inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>)
+{ umfpack_zi_free_symbolic(Symbolic); }
inline int umfpack_symbolic(int n_row,int n_col,
const int Ap[], const int Ai[], const double Ax[], void **Symbolic,
@@ -69,6 +83,39 @@ inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::co
return umfpack_zi_solve(sys,Ap,Ai,&Ax[0].real(),0,&X[0].real(),0,&B[0].real(),0,Numeric,Control,Info);
}
+inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double)
+{
+ return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
+}
+
+inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, std::complex<double>)
+{
+ return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
+}
+
+inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[],
+ int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric)
+{
+ return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
+}
+
+inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], int Up[], int Ui[], std::complex<double> Ux[],
+ int P[], int Q[], std::complex<double> Dx[], int *do_recip, double Rs[], void *Numeric)
+{
+ return umfpack_zi_get_numeric(Lp,Lj,Lx?&Lx[0].real():0,0,Up,Ui,Ux?&Ux[0].real():0,0,P,Q,
+ Dx?&Dx[0].real():0,0,do_recip,Rs,Numeric);
+}
+
+inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
+{
+ return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info);
+}
+
+inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
+{
+ return umfpack_zi_get_determinant(&Mx->real(),0,Ex,NumericHandle,User_Info);
+}
+
template<typename MatrixType>
class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType>
@@ -78,6 +125,10 @@ class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType>
typedef typename Base::Scalar Scalar;
typedef typename Base::RealScalar RealScalar;
typedef Matrix<Scalar,Dynamic,1> Vector;
+ typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
+ typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
+ typedef SparseMatrix<Scalar,Lower|UnitDiagBit> LMatrixType;
+ typedef SparseMatrix<Scalar,Upper> UMatrixType;
using Base::m_flags;
using Base::m_status;
@@ -97,18 +148,53 @@ class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType>
~SparseLU()
{
if (m_numeric)
- umfpack_di_free_numeric(&m_numeric);
+ umfpack_free_numeric(&m_numeric,Scalar());
+ }
+
+ inline const LMatrixType& matrixL() const
+ {
+ if (m_extractedDataAreDirty) extractData();
+ return m_l;
}
+ inline const UMatrixType& matrixU() const
+ {
+ if (m_extractedDataAreDirty) extractData();
+ return m_u;
+ }
+
+ inline const IntColVectorType& permutationP() const
+ {
+ if (m_extractedDataAreDirty) extractData();
+ return m_p;
+ }
+
+ inline const IntRowVectorType& permutationQ() const
+ {
+ if (m_extractedDataAreDirty) extractData();
+ return m_q;
+ }
+
+ Scalar determinant() const;
+
template<typename BDerived, typename XDerived>
bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const;
void compute(const MatrixType& matrix);
protected:
+
+ void extractData() const;
+
+ protected:
// cached data:
void* m_numeric;
const MatrixType* m_matrixRef;
+ mutable LMatrixType m_l;
+ mutable UMatrixType m_u;
+ mutable IntColVectorType m_p;
+ mutable IntRowVectorType m_q;
+ mutable bool m_extractedDataAreDirty;
};
template<typename MatrixType>
@@ -121,7 +207,7 @@ void SparseLU<MatrixType,UmfPack>::compute(const MatrixType& a)
m_matrixRef = &a;
if (m_numeric)
- umfpack_di_free_numeric(&m_numeric);
+ umfpack_free_numeric(&m_numeric,Scalar());
void* symbolic;
int errorCode = 0;
@@ -131,26 +217,48 @@ void SparseLU<MatrixType,UmfPack>::compute(const MatrixType& a)
errorCode = umfpack_numeric(a._outerIndexPtr(), a._innerIndexPtr(), a._valuePtr(),
symbolic, &m_numeric, 0, 0);
- umfpack_di_free_symbolic(&symbolic);
+ umfpack_free_symbolic(&symbolic,Scalar());
+
+ m_extractedDataAreDirty = true;
Base::m_succeeded = (errorCode==0);
}
-// template<typename MatrixType>
-// inline const MatrixType&
-// SparseLU<MatrixType,SuperLU>::matrixL() const
-// {
-// ei_assert(false && "matrixL() is Not supported by the SuperLU backend");
-// return m_matrix;
-// }
-//
-// template<typename MatrixType>
-// inline const MatrixType&
-// SparseLU<MatrixType,SuperLU>::matrixU() const
-// {
-// ei_assert(false && "matrixU() is Not supported by the SuperLU backend");
-// return m_matrix;
-// }
+template<typename MatrixType>
+void SparseLU<MatrixType,UmfPack>::extractData() const
+{
+ if (m_extractedDataAreDirty)
+ {
+ // get size of the data
+ int lnz, unz, rows, cols, nz_udiag;
+ umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar());
+
+ // allocate data
+ m_l.resize(rows,std::min(rows,cols));
+ m_l.resizeNonZeros(lnz);
+
+ m_u.resize(std::min(rows,cols),cols);
+ m_u.resizeNonZeros(unz);
+
+ m_p.resize(rows);
+ m_q.resize(cols);
+
+ // extract
+ umfpack_get_numeric(m_l._outerIndexPtr(), m_l._innerIndexPtr(), m_l._valuePtr(),
+ m_u._outerIndexPtr(), m_u._innerIndexPtr(), m_u._valuePtr(),
+ m_p.data(), m_q.data(), 0, 0, 0, m_numeric);
+
+ m_extractedDataAreDirty = false;
+ }
+}
+
+template<typename MatrixType>
+typename SparseLU<MatrixType,UmfPack>::Scalar SparseLU<MatrixType,UmfPack>::determinant() const
+{
+ Scalar det;
+ umfpack_get_determinant(&det, 0, m_numeric, 0);
+ return det;
+}
template<typename MatrixType>
template<typename BDerived,typename XDerived>
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-29 15:34:31 +0000