00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_PASTIXSUPPORT_H
00026 #define EIGEN_PASTIXSUPPORT_H
00027
00028 namespace Eigen {
00029
00039 template<typename _MatrixType, bool IsStrSym = false> class PastixLU;
00040 template<typename _MatrixType, int Options> class PastixLLT;
00041 template<typename _MatrixType, int Options> class PastixLDLT;
00042
00043 namespace internal
00044 {
00045
00046 template<class Pastix> struct pastix_traits;
00047
00048 template<typename _MatrixType>
00049 struct pastix_traits< PastixLU<_MatrixType> >
00050 {
00051 typedef _MatrixType MatrixType;
00052 typedef typename _MatrixType::Scalar Scalar;
00053 typedef typename _MatrixType::RealScalar RealScalar;
00054 typedef typename _MatrixType::Index Index;
00055 };
00056
00057 template<typename _MatrixType, int Options>
00058 struct pastix_traits< PastixLLT<_MatrixType,Options> >
00059 {
00060 typedef _MatrixType MatrixType;
00061 typedef typename _MatrixType::Scalar Scalar;
00062 typedef typename _MatrixType::RealScalar RealScalar;
00063 typedef typename _MatrixType::Index Index;
00064 };
00065
00066 template<typename _MatrixType, int Options>
00067 struct pastix_traits< PastixLDLT<_MatrixType,Options> >
00068 {
00069 typedef _MatrixType MatrixType;
00070 typedef typename _MatrixType::Scalar Scalar;
00071 typedef typename _MatrixType::RealScalar RealScalar;
00072 typedef typename _MatrixType::Index Index;
00073 };
00074
00075 void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, float *vals, int *perm, int * invp, float *x, int nbrhs, int *iparm, double *dparm)
00076 {
00077 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
00078 if (nbrhs == 0) x = NULL;
00079 s_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm);
00080 }
00081
00082 void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, double *vals, int *perm, int * invp, double *x, int nbrhs, int *iparm, double *dparm)
00083 {
00084 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
00085 if (nbrhs == 0) x = NULL;
00086 d_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm);
00087 }
00088
00089 void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex<float> *vals, int *perm, int * invp, std::complex<float> *x, int nbrhs, int *iparm, double *dparm)
00090 {
00091 c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<COMPLEX*>(vals), perm, invp, reinterpret_cast<COMPLEX*>(x), nbrhs, iparm, dparm);
00092 }
00093
00094 void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex<double> *vals, int *perm, int * invp, std::complex<double> *x, int nbrhs, int *iparm, double *dparm)
00095 {
00096 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
00097 if (nbrhs == 0) x = NULL;
00098 z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<DCOMPLEX*>(vals), perm, invp, reinterpret_cast<DCOMPLEX*>(x), nbrhs, iparm, dparm);
00099 }
00100
00101
00102 template <typename MatrixType>
00103 void EigenToFortranNumbering (MatrixType& mat)
00104 {
00105 if ( !(mat.outerIndexPtr()[0]) )
00106 {
00107 int i;
00108 for(i = 0; i <= mat.rows(); ++i)
00109 ++mat.outerIndexPtr()[i];
00110 for(i = 0; i < mat.nonZeros(); ++i)
00111 ++mat.innerIndexPtr()[i];
00112 }
00113 }
00114
00115
00116 template <typename MatrixType>
00117 void EigenToCNumbering (MatrixType& mat)
00118 {
00119
00120 if ( mat.outerIndexPtr()[0] == 1 )
00121 {
00122 int i;
00123 for(i = 0; i <= mat.rows(); ++i)
00124 --mat.outerIndexPtr()[i];
00125 for(i = 0; i < mat.nonZeros(); ++i)
00126 --mat.innerIndexPtr()[i];
00127 }
00128 }
00129
00130
00131
00132
00133
00134
00135
00136
00137 template <typename MatrixType>
00138 void EigenSymmetrizeMatrixGraph (const MatrixType& In, MatrixType& Out, MatrixType& StrMatTrans, bool& hasTranspose)
00139 {
00140 eigen_assert(In.cols()==In.rows() && " Can only symmetrize the graph of a square matrix");
00141 if (!hasTranspose)
00142 {
00143 StrMatTrans = In.transpose();
00144
00145 for (int i = 0; i < StrMatTrans.rows(); i++)
00146 {
00147 for (typename MatrixType::InnerIterator it(StrMatTrans, i); it; ++it)
00148 it.valueRef() = 0.0;
00149 }
00150 hasTranspose = true;
00151 }
00152 Out = (StrMatTrans + In).eval();
00153 }
00154
00155 }
00156
00157
00158
00159 template <class Derived>
00160 class PastixBase
00161 {
00162 public:
00163 typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType;
00164 typedef _MatrixType MatrixType;
00165 typedef typename MatrixType::Scalar Scalar;
00166 typedef typename MatrixType::RealScalar RealScalar;
00167 typedef typename MatrixType::Index Index;
00168 typedef Matrix<Scalar,Dynamic,1> Vector;
00169
00170 public:
00171
00172 PastixBase():m_initisOk(false),m_analysisIsOk(false),m_factorizationIsOk(false),m_isInitialized(false)
00173 {
00174 m_pastixdata = 0;
00175 m_hasTranspose = false;
00176 PastixInit();
00177 }
00178
00179 ~PastixBase()
00180 {
00181 PastixDestroy();
00182 }
00183
00184
00185 void PastixInit();
00186
00187
00188 Derived& analyzePattern (MatrixType& mat);
00189
00190
00191 Derived& factorize (MatrixType& mat);
00192
00197 template<typename Rhs>
00198 inline const internal::solve_retval<PastixBase, Rhs>
00199 solve(const MatrixBase<Rhs>& b) const
00200 {
00201 eigen_assert(m_isInitialized && "Pastix solver is not initialized.");
00202 eigen_assert(rows()==b.rows()
00203 && "PastixBase::solve(): invalid number of rows of the right hand side matrix b");
00204 return internal::solve_retval<PastixBase, Rhs>(*this, b.derived());
00205 }
00206
00207 template<typename Rhs,typename Dest>
00208 bool _solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const;
00209
00211 template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
00212 void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
00213 {
00214 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00215 eigen_assert(rows()==b.rows());
00216
00217
00218 static const int NbColsAtOnce = 1;
00219 int rhsCols = b.cols();
00220 int size = b.rows();
00221 Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols);
00222 for(int k=0; k<rhsCols; k+=NbColsAtOnce)
00223 {
00224 int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
00225 tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
00226 tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols));
00227 dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
00228 }
00229 }
00230
00231 Derived& derived()
00232 {
00233 return *static_cast<Derived*>(this);
00234 }
00235 const Derived& derived() const
00236 {
00237 return *static_cast<const Derived*>(this);
00238 }
00239
00245 Array<Index,IPARM_SIZE,1>& iparm()
00246 {
00247 return m_iparm;
00248 }
00249
00254 int& iparm(int idxparam)
00255 {
00256 return m_iparm(idxparam);
00257 }
00258
00263 Array<RealScalar,IPARM_SIZE,1>& dparm()
00264 {
00265 return m_dparm;
00266 }
00267
00268
00273 double& dparm(int idxparam)
00274 {
00275 return m_dparm(idxparam);
00276 }
00277
00278 inline Index cols() const { return m_size; }
00279 inline Index rows() const { return m_size; }
00280
00289 ComputationInfo info() const
00290 {
00291 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00292 return m_info;
00293 }
00294
00299 template<typename Rhs>
00300 inline const internal::sparse_solve_retval<PastixBase, Rhs>
00301 solve(const SparseMatrixBase<Rhs>& b) const
00302 {
00303 eigen_assert(m_isInitialized && "Pastix LU, LLT or LDLT is not initialized.");
00304 eigen_assert(rows()==b.rows()
00305 && "PastixBase::solve(): invalid number of rows of the right hand side matrix b");
00306 return internal::sparse_solve_retval<PastixBase, Rhs>(*this, b.derived());
00307 }
00308
00309 protected:
00310
00311 void PastixDestroy()
00312 {
00313 eigen_assert(m_initisOk && "The Pastix structure should be allocated first");
00314 m_iparm(IPARM_START_TASK) = API_TASK_CLEAN;
00315 m_iparm(IPARM_END_TASK) = API_TASK_CLEAN;
00316 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(),
00317 m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 1, m_iparm.data(), m_dparm.data());
00318 }
00319
00320 Derived& compute (MatrixType& mat);
00321
00322 int m_initisOk;
00323 int m_analysisIsOk;
00324 int m_factorizationIsOk;
00325 bool m_isInitialized;
00326 mutable ComputationInfo m_info;
00327 mutable pastix_data_t *m_pastixdata;
00328 mutable SparseMatrix<Scalar, ColMajor> m_mat_null;
00329 mutable Matrix<Scalar, Dynamic,1> m_vec_null;
00330 mutable SparseMatrix<Scalar, ColMajor> m_StrMatTrans;
00331 mutable bool m_hasTranspose;
00332 mutable int m_comm;
00333 mutable Matrix<Index,IPARM_SIZE,1> m_iparm;
00334 mutable Matrix<double,DPARM_SIZE,1> m_dparm;
00335 mutable Matrix<Index,Dynamic,1> m_perm;
00336 mutable Matrix<Index,Dynamic,1> m_invp;
00337 mutable int m_ordering;
00338 mutable int m_amalgamation;
00339 mutable int m_size;
00340
00341 private:
00342 PastixBase(PastixBase& ) {}
00343
00344 };
00345
00350 template <class Derived>
00351 void PastixBase<Derived>::PastixInit()
00352 {
00353 m_size = 0;
00354 m_iparm.resize(IPARM_SIZE);
00355 m_dparm.resize(DPARM_SIZE);
00356
00357 m_iparm(IPARM_MODIFY_PARAMETER) = API_NO;
00358 if(m_pastixdata)
00359 {
00360
00361 PastixDestroy();
00362 m_pastixdata = 0;
00363 m_iparm(IPARM_MODIFY_PARAMETER) = API_YES;
00364 m_hasTranspose = false;
00365 }
00366
00367 m_iparm(IPARM_START_TASK) = API_TASK_INIT;
00368 m_iparm(IPARM_END_TASK) = API_TASK_INIT;
00369 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
00370 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(),
00371 m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 1, m_iparm.data(), m_dparm.data());
00372
00373 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
00374
00375
00376 if(m_iparm(IPARM_ERROR_NUMBER)) {
00377 m_info = InvalidInput;
00378 m_initisOk = false;
00379 }
00380 else {
00381 m_info = Success;
00382 m_initisOk = true;
00383 }
00384 }
00385
00386 template <class Derived>
00387 Derived& PastixBase<Derived>::compute(MatrixType& mat)
00388 {
00389 eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared");
00390 typedef typename MatrixType::Scalar Scalar;
00391
00392
00393 m_size = mat.rows();
00394
00395 internal::EigenToFortranNumbering(mat);
00396 analyzePattern(mat);
00397 factorize(mat);
00398 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
00399 if (m_factorizationIsOk) m_isInitialized = true;
00400
00401
00402 internal::EigenToCNumbering(mat);
00403
00404 return derived();
00405 }
00406
00407
00408 template <class Derived>
00409 Derived& PastixBase<Derived>::analyzePattern(MatrixType& mat)
00410 {
00411 eigen_assert(m_initisOk && "PastixInit should be called first to set the default parameters");
00412 m_size = mat.rows();
00413 m_perm.resize(m_size);
00414 m_invp.resize(m_size);
00415
00416
00417 internal::EigenToFortranNumbering(mat);
00418
00419 m_iparm(IPARM_START_TASK) = API_TASK_ORDERING;
00420 m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE;
00421
00422 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
00423 mat.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 0, m_iparm.data(), m_dparm.data());
00424
00425
00426 if(m_iparm(IPARM_ERROR_NUMBER)) {
00427 m_info = NumericalIssue;
00428 m_analysisIsOk = false;
00429 }
00430 else {
00431 m_info = Success;
00432 m_analysisIsOk = true;
00433 }
00434 return derived();
00435 }
00436
00437 template <class Derived>
00438 Derived& PastixBase<Derived>::factorize(MatrixType& mat)
00439 {
00440 eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase");
00441 m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT;
00442 m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT;
00443 m_size = mat.rows();
00444
00445
00446 internal::EigenToFortranNumbering(mat);
00447
00448 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
00449 mat.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 0, m_iparm.data(), m_dparm.data());
00450
00451
00452 if(m_iparm(IPARM_ERROR_NUMBER)) {
00453 m_info = NumericalIssue;
00454 m_factorizationIsOk = false;
00455 m_isInitialized = false;
00456 }
00457 else {
00458 m_info = Success;
00459 m_factorizationIsOk = true;
00460 m_isInitialized = true;
00461 }
00462 return derived();
00463 }
00464
00465
00466 template<typename Base>
00467 template<typename Rhs,typename Dest>
00468 bool PastixBase<Base>::_solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const
00469 {
00470 eigen_assert(m_isInitialized && "The matrix should be factorized first");
00471 EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
00472 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00473 int rhs = 1;
00474
00475 x = b;
00476
00477 for (int i = 0; i < b.cols(); i++){
00478 m_iparm(IPARM_START_TASK) = API_TASK_SOLVE;
00479 m_iparm(IPARM_END_TASK) = API_TASK_REFINE;
00480 m_iparm(IPARM_RHS_MAKING) = API_RHS_B;
00481 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, x.rows(), m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(),
00482 m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data());
00483 }
00484
00485 if(m_iparm(IPARM_ERROR_NUMBER)) {
00486 m_info = NumericalIssue;
00487 return false;
00488 }
00489 else {
00490 return true;
00491 }
00492 }
00493
00513 template<typename _MatrixType, bool IsStrSym>
00514 class PastixLU : public PastixBase< PastixLU<_MatrixType> >
00515 {
00516 public:
00517 typedef _MatrixType MatrixType;
00518 typedef PastixBase<PastixLU<MatrixType> > Base;
00519 typedef typename MatrixType::Scalar Scalar;
00520 typedef SparseMatrix<Scalar, ColMajor> PaStiXType;
00521
00522 public:
00523 PastixLU():Base() {}
00524
00525 PastixLU(const MatrixType& matrix):Base()
00526 {
00527 compute(matrix);
00528 }
00534 void compute (const MatrixType& matrix)
00535 {
00536
00537 Base::PastixInit();
00538 PaStiXType temp(matrix.rows(), matrix.cols());
00539
00540 if (IsStrSym)
00541 temp = matrix;
00542 else
00543 {
00544 internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans, m_hasTranspose);
00545 }
00546 m_iparm[IPARM_SYM] = API_SYM_NO;
00547 m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
00548 Base::compute(temp);
00549 }
00555 void analyzePattern(const MatrixType& matrix)
00556 {
00557
00558 Base::PastixInit();
00559
00560 SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00561
00562 if (IsStrSym)
00563 temp = matrix;
00564 else
00565 {
00566 internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans,m_hasTranspose);
00567 }
00568
00569 m_iparm(IPARM_SYM) = API_SYM_NO;
00570 m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
00571 Base::analyzePattern(temp);
00572 }
00573
00579 void factorize(const MatrixType& matrix)
00580 {
00581
00582 SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00583
00584 if (IsStrSym)
00585 temp = matrix;
00586 else
00587 {
00588 internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans,m_hasTranspose);
00589 }
00590 m_iparm(IPARM_SYM) = API_SYM_NO;
00591 m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
00592 Base::factorize(temp);
00593 }
00594 protected:
00595 using Base::m_iparm;
00596 using Base::m_dparm;
00597 using Base::m_StrMatTrans;
00598 using Base::m_hasTranspose;
00599
00600 private:
00601 PastixLU(PastixLU& ) {}
00602 };
00603
00618 template<typename _MatrixType, int _UpLo>
00619 class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> >
00620 {
00621 public:
00622 typedef _MatrixType MatrixType;
00623 typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base;
00624 typedef typename MatrixType::Scalar Scalar;
00625 typedef typename MatrixType::Index Index;
00626
00627 public:
00628 enum { UpLo = _UpLo };
00629 PastixLLT():Base() {}
00630
00631 PastixLLT(const MatrixType& matrix):Base()
00632 {
00633 compute(matrix);
00634 }
00635
00639 void compute (const MatrixType& matrix)
00640 {
00641
00642 Base::PastixInit();
00643 SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00644 PermutationMatrix<Dynamic,Dynamic,Index> pnull;
00645 temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
00646 m_iparm(IPARM_SYM) = API_SYM_YES;
00647 m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
00648 Base::compute(temp);
00649 }
00650
00655 void analyzePattern(const MatrixType& matrix)
00656 {
00657 Base::PastixInit();
00658
00659 SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00660 PermutationMatrix<Dynamic,Dynamic,Index> pnull;
00661 temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
00662 m_iparm(IPARM_SYM) = API_SYM_YES;
00663 m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
00664 Base::analyzePattern(temp);
00665 }
00669 void factorize(const MatrixType& matrix)
00670 {
00671
00672 SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00673 PermutationMatrix<Dynamic,Dynamic,Index> pnull;
00674 temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
00675 m_iparm(IPARM_SYM) = API_SYM_YES;
00676 m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
00677 Base::factorize(temp);
00678 }
00679 protected:
00680 using Base::m_iparm;
00681
00682 private:
00683 PastixLLT(PastixLLT& ) {}
00684 };
00685
00700 template<typename _MatrixType, int _UpLo>
00701 class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> >
00702 {
00703 public:
00704 typedef _MatrixType MatrixType;
00705 typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base;
00706 typedef typename MatrixType::Scalar Scalar;
00707 typedef typename MatrixType::Index Index;
00708
00709 public:
00710 enum { UpLo = _UpLo };
00711 PastixLDLT():Base() {}
00712
00713 PastixLDLT(const MatrixType& matrix):Base()
00714 {
00715 compute(matrix);
00716 }
00717
00721 void compute (const MatrixType& matrix)
00722 {
00723 Base::PastixInit();
00724
00725 SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00726 PermutationMatrix<Dynamic,Dynamic,Index> pnull;
00727 temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
00728 m_iparm(IPARM_SYM) = API_SYM_YES;
00729 m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
00730 Base::compute(temp);
00731 }
00732
00737 void analyzePattern(const MatrixType& matrix)
00738 {
00739 Base::PastixInit();
00740
00741 SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00742 PermutationMatrix<Dynamic,Dynamic,Index> pnull;
00743 temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
00744
00745 m_iparm(IPARM_SYM) = API_SYM_YES;
00746 m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
00747 Base::analyzePattern(temp);
00748 }
00752 void factorize(const MatrixType& matrix)
00753 {
00754
00755 SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00756 PermutationMatrix<Dynamic,Dynamic,Index> pnull;
00757 temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
00758
00759 m_iparm(IPARM_SYM) = API_SYM_YES;
00760 m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
00761 Base::factorize(temp);
00762 }
00763
00764 protected:
00765 using Base::m_iparm;
00766
00767 private:
00768 PastixLDLT(PastixLDLT& ) {}
00769 };
00770
00771 namespace internal {
00772
00773 template<typename _MatrixType, typename Rhs>
00774 struct solve_retval<PastixBase<_MatrixType>, Rhs>
00775 : solve_retval_base<PastixBase<_MatrixType>, Rhs>
00776 {
00777 typedef PastixBase<_MatrixType> Dec;
00778 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00779
00780 template<typename Dest> void evalTo(Dest& dst) const
00781 {
00782 dec()._solve(rhs(),dst);
00783 }
00784 };
00785
00786 template<typename _MatrixType, typename Rhs>
00787 struct sparse_solve_retval<PastixBase<_MatrixType>, Rhs>
00788 : sparse_solve_retval_base<PastixBase<_MatrixType>, Rhs>
00789 {
00790 typedef PastixBase<_MatrixType> Dec;
00791 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00792
00793 template<typename Dest> void evalTo(Dest& dst) const
00794 {
00795 dec()._solve_sparse(rhs(),dst);
00796 }
00797 };
00798
00799 }
00800
00801 }
00802
00803 #endif