diff options
author | Gael Guennebaud <g.gael@free.fr> | 2010-11-18 10:30:52 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2010-11-18 10:30:52 +0100 |
commit | 3e99356b59ea194599cb80adb7b10989d6f4a54a (patch) | |
tree | 6b41dbaf5b3cc2cf1734a2a40c200d3d39babfd9 /unsupported/Eigen/src/SparseExtra | |
parent | 1618df55dfc971113a974841b52e86391f445ff2 (diff) |
clean a bit AMD and SimplicialCholesky and add support for partly stored selfadjoint matrices
Diffstat (limited to 'unsupported/Eigen/src/SparseExtra')
-rw-r--r-- | unsupported/Eigen/src/SparseExtra/Amd.h | 67 | ||||
-rw-r--r-- | unsupported/Eigen/src/SparseExtra/CholmodSupport.h | 5 | ||||
-rw-r--r-- | unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h | 374 |
3 files changed, 220 insertions, 226 deletions
diff --git a/unsupported/Eigen/src/SparseExtra/Amd.h b/unsupported/Eigen/src/SparseExtra/Amd.h index a716305fd..f52183232 100644 --- a/unsupported/Eigen/src/SparseExtra/Amd.h +++ b/unsupported/Eigen/src/SparseExtra/Amd.h @@ -96,18 +96,14 @@ Index cs_tdfs(Index j, Index k, Index *head, const Index *next, Index *post, Ind return k; } -/** keeps off-diagonal entries; drops diagonal entries */ -template<typename Index, typename Scalar> -struct keep_diag { - inline bool operator() (const Index& row, const Index& col, const Scalar&) const - { - return row!=col; - } -}; -/** p = amd(A+A') if symmetric is true, or amd(A'A) otherwise */ -template<typename Scalar, int Options, typename Index> -int *minimum_degree_ordering(int order, const SparseMatrix<Scalar,Options,Index>& A) /* order 0:natural, 1:Chol, 2:LU, 3:QR */ +/** \internal + * Approximate minimum degree ordering algorithm. + * \returns the permutation P reducing the fill-in of the input matrix \a C + * The input matrix \a C must be a selfadjoint compressed column major SparseMatrix object. Both the upper and lower parts have to be stored, but the diagonal entries are optional. + * On exit the values of C are destroyed */ +template<typename Scalar, typename Index> +void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, PermutationMatrix<Dynamic>& perm) { typedef SparseMatrix<Scalar,ColMajor,Index> CCS; @@ -115,49 +111,13 @@ int *minimum_degree_ordering(int order, const SparseMatrix<Scalar,Options,Index> k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi, ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, t; unsigned int h; - /* --- Construct matrix C ----------------------------------------------- */ - if(order <= 0 || order > 3) return (NULL); /* check */ - Index m = A.rows(); - Index n = A.cols(); + Index n = C.cols(); dense = std::max<Index> (16, 10 * sqrt ((double) n)); /* find dense threshold */ dense = std::min<Index> (n-2, dense); - CCS C; - if(order == 1 && n == m) // Cholesky - { - C = A + SparseMatrix<Scalar,Options,Index>(A.adjoint()); - } - else if(order == 2) // LU - { - CCS AT = A.adjoint(); - // drop dense columns from AT - Index* ATp = AT._outerIndexPtr(); - Index* ATi = AT._innerIndexPtr(); - Index p2; - Index j; - for(p2 = 0, j = 0; j < m; j++) - { - Index p = ATp[j]; // column j of AT starts here - ATp[j] = p2; // new column j starts here - if(ATp[j+1] - p > dense) continue; // skip dense col j - for(; p < ATp[j+1]; p++) - ATi[p2++] = ATi[p]; - } - ATp[m] = p2; // finalize AT - // TODO this could be implemented using a sparse filter expression - // TODO do a cheap selfadjoint rank update - C = AT * AT.adjoint(); // C=A'*A with no dense rows - } - else // QR - { - C = A.adjoint() * A; - } - - C.prune(keep_diag<Index,Scalar>()); - - Index cnz = A.nonZeros(); - Index* P = new Index[n+1]; /* allocate result */ + Index cnz = C.nonZeros(); + perm.resize(n+1); t = cnz + cnz/5 + 2*n; /* add elbow room to C */ C.resizeNonZeros(t); @@ -170,7 +130,7 @@ int *minimum_degree_ordering(int order, const SparseMatrix<Scalar,Options,Index> Index* degree = W + 5*(n+1); Index* w = W + 6*(n+1); Index* hhead = W + 7*(n+1); - Index* last = P; /* use P as workspace for last */ + Index* last = perm.indices().data(); /* use P as workspace for last */ /* --- Initialize quotient graph ---------------------------------------- */ Index* Cp = C._outerIndexPtr(); @@ -475,11 +435,12 @@ int *minimum_degree_ordering(int order, const SparseMatrix<Scalar,Options,Index> } for(k = 0, i = 0; i <= n; i++) /* postorder the assembly tree */ { - if(Cp[i] == -1) k = cs_tdfs (i, k, head, next, P, w); + if(Cp[i] == -1) k = cs_tdfs (i, k, head, next, perm.indices().data(), w); } + + perm.indices().conservativeResize(n); delete[] W; - return P; } } // namespace internal diff --git a/unsupported/Eigen/src/SparseExtra/CholmodSupport.h b/unsupported/Eigen/src/SparseExtra/CholmodSupport.h index 14f18d30f..3e502c0aa 100644 --- a/unsupported/Eigen/src/SparseExtra/CholmodSupport.h +++ b/unsupported/Eigen/src/SparseExtra/CholmodSupport.h @@ -177,6 +177,7 @@ class CholmodDecomposition : m_cholmodFactor(0), m_info(Success), m_isInitialized(false) { cholmod_start(&m_cholmod); + setMode(CholmodLDLt); } CholmodDecomposition(const MatrixType& matrix) @@ -351,6 +352,10 @@ class CholmodDecomposition cholmod_free_sparse(&x_cs, &m_cholmod); } #endif // EIGEN_PARSED_BY_DOXYGEN + + template<typename Stream> + void dumpMemory(Stream& s) + {} protected: mutable cholmod_common m_cholmod; diff --git a/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h b/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h index 6b22d1502..302be1cce 100644 --- a/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h +++ b/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h @@ -71,7 +71,7 @@ enum SimplicialCholeskyMode { /** \brief A direct sparse Cholesky factorization * * This class allows to solve for A.X = B sparse linear problems via a LL^T or LDL^T Cholesky factorization. - * The sparse matrix A must be selfajoint and positive definite. The vectors or matrices + * The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices * X and B can be either dense or sparse. * * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> @@ -180,73 +180,8 @@ class SimplicialCholesky * * \sa factorize() */ - void analyzePattern(const MatrixType& a) - { - eigen_assert(a.rows()==a.cols()); - const Index size = a.rows(); - m_matrix.resize(size, size); - m_parent.resize(size); - m_nonZerosPerCol.resize(size); - - Index* tags = ei_aligned_stack_new(Index, size); - - // TODO allows to configure the permutation - const Index* P = internal::minimum_degree_ordering(1, a); - const Index* Pinv = 0; - if(P) - { - m_P.indices() = VectorXi::Map(P,size); - m_Pinv = m_P.inverse(); - Pinv = m_Pinv.indices().data(); - } - else - { - m_P.resize(0); - m_Pinv.resize(0); - } - - for(Index k = 0; k < size; ++k) - { - /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */ - m_parent[k] = -1; /* parent of k is not yet known */ - tags[k] = k; /* mark node k as visited */ - m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */ - Index kk = P ? P[k] : k; /* kth original, or permuted, column */ - for(typename MatrixType::InnerIterator it(a,kk); it; ++it) - { - /* A (i,k) is nonzero (original or permuted A) */ - Index i = Pinv ? Pinv[it.index()] : it.index(); - if(i < k) - { - /* follow path from i to root of etree, stop at flagged node */ - for(; tags[i] != k; i = m_parent[i]) - { - /* find parent of i if not yet determined */ - if (m_parent[i] == -1) - m_parent[i] = k; - ++m_nonZerosPerCol[i]; /* L (k,i) is nonzero */ - tags[i] = k; /* mark i as visited */ - } - } - } - } - - // release worspace - ei_aligned_stack_delete(Index, tags, size); - - /* construct Lp index array from m_nonZerosPerCol column counts */ - Index* Lp = m_matrix._outerIndexPtr(); - Lp[0] = 0; - for(Index k = 0; k < size; ++k) - Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (m_LDLt ? 0 : 1); - - m_matrix.resizeNonZeros(Lp[size]); - - m_isInitialized = true; - m_info = Success; - m_analysisIsOk = true; - m_factorizationIsOk = false; - } + void analyzePattern(const MatrixType& a); + /** Performs a numeric decomposition of \a matrix * @@ -254,111 +189,17 @@ class SimplicialCholesky * * \sa analyzePattern() */ - void factorize(const MatrixType& a) - { - eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); - eigen_assert(a.rows()==a.cols()); - const Index size = a.rows(); - eigen_assert(m_parent.size()==size); - eigen_assert(m_nonZerosPerCol.size()==size); - - const Index* Lp = m_matrix._outerIndexPtr(); - Index* Li = m_matrix._innerIndexPtr(); - Scalar* Lx = m_matrix._valuePtr(); - - Scalar* y = ei_aligned_stack_new(Scalar, size); - Index* pattern = ei_aligned_stack_new(Index, size); - Index* tags = ei_aligned_stack_new(Index, size); - - Index* P = 0; - Index* Pinv = 0; - - if(m_P.size()==size) - { - P = m_P.indices().data(); - Pinv = m_Pinv.indices().data(); - } - - bool ok = true; - - m_diag.resize(m_LDLt ? size : 0); - - for(Index k = 0; k < size; ++k) - { - /* compute nonzero pattern of kth row of L, in topological order */ - y[k] = 0.0; /* Y(0:k) is now all zero */ - Index top = size; /* stack for pattern is empty */ - tags[k] = k; /* mark node k as visited */ - m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */ - Index kk = (P) ? (P[k]) : (k); /* kth original, or permuted, column */ - for(typename MatrixType::InnerIterator it(a,kk); it; ++it) - { - Index i = Pinv ? Pinv[it.index()] : it.index(); - if(i <= k) - { - y[i] += internal::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */ - Index len; - for(len = 0; tags[i] != k; i = m_parent[i]) - { - pattern[len++] = i; /* L(k,i) is nonzero */ - tags[i] = k; /* mark i as visited */ - } - while(len > 0) - pattern[--top] = pattern[--len]; - } - } - - /* compute numerical values kth row of L (a sparse triangular solve) */ - Scalar d = y[k]; // get D(k,k) and clear Y(k) - y[k] = 0.0; - for(; top < size; ++top) - { - if(1) - { - Index i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */ - Scalar yi = y[i]; /* get and clear Y(i) */ - y[i] = 0.0; - - /* the nonzero entry L(k,i) */ - Scalar l_ki; - if(m_LDLt) - l_ki = yi / m_diag[i]; - else - yi = l_ki = yi / Lx[Lp[i]]; - - Index p2 = Lp[i] + m_nonZerosPerCol[i]; - Index p; - for(p = Lp[i] + (m_LDLt ? 0 : 1); p < p2; ++p) - y[Li[p]] -= internal::conj(Lx[p]) * yi; - d -= l_ki * internal::conj(yi); - Li[p] = k; /* store L(k,i) in column form of L */ - Lx[p] = l_ki; - ++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */ - } - } - if(m_LDLt) - m_diag[k] = d; - else - { - Index p = Lp[k]+m_nonZerosPerCol[k]++; - Li[p] = k ; /* store L(k,k) = sqrt (d) in column k */ - Lx[p] = internal::sqrt(d) ; - } - if(d == Scalar(0)) - { - ok = false; /* failure, D(k,k) is zero */ - break; - } - } - - // release workspace - ei_aligned_stack_delete(Scalar, y, size); - ei_aligned_stack_delete(Index, pattern, size); - ei_aligned_stack_delete(Index, tags, size); - - m_info = ok ? Success : NumericalIssue; - m_factorizationIsOk = true; - } + void factorize(const MatrixType& a); + + /** \returns the permutation P + * \sa permutationPinv() */ + const PermutationMatrix<Dynamic>& permutationP() const + { return m_P; } + + /** \returns the inverse P^-1 of the permutation P + * \sa permutationP() */ + const PermutationMatrix<Dynamic>& permutationPinv() const + { return m_Pinv; } #ifndef EIGEN_PARSED_BY_DOXYGEN /** \internal */ @@ -408,8 +249,29 @@ class SimplicialCholesky } */ #endif // EIGEN_PARSED_BY_DOXYGEN + + template<typename Stream> + void dumpMemory(Stream& s) + { + int total = 0; + s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n"; + s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n"; + s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n"; + s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n"; + s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n"; + s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n"; + s << " TOTAL: " << (total>> 20) << "Mb" << "\n"; + } protected: + /** keeps off-diagonal entries; drops diagonal entries */ + struct keep_diag { + inline bool operator() (const Index& row, const Index& col, const Scalar&) const + { + return row!=col; + } + }; + mutable ComputationInfo m_info; bool m_isInitialized; bool m_factorizationIsOk; @@ -424,6 +286,172 @@ class SimplicialCholesky PermutationMatrix<Dynamic> m_Pinv; // the inverse permutation }; +template<typename _MatrixType, int _UpLo> +void SimplicialCholesky<_MatrixType,_UpLo>::analyzePattern(const MatrixType& a) +{ + eigen_assert(a.rows()==a.cols()); + const Index size = a.rows(); + m_matrix.resize(size, size); + m_parent.resize(size); + m_nonZerosPerCol.resize(size); + + Index* tags = ei_aligned_stack_new(Index, size); + + // TODO allows to configure the permutation + { + CholMatrixType C; + C = a.template selfadjointView<UpLo>(); + // remove diagonal entries: + C.prune(keep_diag()); + internal::minimum_degree_ordering(C, m_P); + } + + if(m_P.size()>0) + m_Pinv = m_P.inverse(); + else + m_Pinv.resize(0); + + SparseMatrix<Scalar,ColMajor,Index> ap(size,size); + ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv); + + for(Index k = 0; k < size; ++k) + { + /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */ + m_parent[k] = -1; /* parent of k is not yet known */ + tags[k] = k; /* mark node k as visited */ + m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */ + for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it) + { + Index i = it.index(); + if(i < k) + { + /* follow path from i to root of etree, stop at flagged node */ + for(; tags[i] != k; i = m_parent[i]) + { + /* find parent of i if not yet determined */ + if (m_parent[i] == -1) + m_parent[i] = k; + m_nonZerosPerCol[i]++; /* L (k,i) is nonzero */ + tags[i] = k; /* mark i as visited */ + } + } + } + } + + // release workspace + ei_aligned_stack_delete(Index, tags, size); + + /* construct Lp index array from m_nonZerosPerCol column counts */ + Index* Lp = m_matrix._outerIndexPtr(); + Lp[0] = 0; + for(Index k = 0; k < size; ++k) + Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (m_LDLt ? 0 : 1); + + m_matrix.resizeNonZeros(Lp[size]); + + m_isInitialized = true; + m_info = Success; + m_analysisIsOk = true; + m_factorizationIsOk = false; +} + + +template<typename _MatrixType, int _UpLo> +void SimplicialCholesky<_MatrixType,_UpLo>::factorize(const MatrixType& a) +{ + eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); + eigen_assert(a.rows()==a.cols()); + const Index size = a.rows(); + eigen_assert(m_parent.size()==size); + eigen_assert(m_nonZerosPerCol.size()==size); + + const Index* Lp = m_matrix._outerIndexPtr(); + Index* Li = m_matrix._innerIndexPtr(); + Scalar* Lx = m_matrix._valuePtr(); + + Scalar* y = ei_aligned_stack_new(Scalar, size); + Index* pattern = ei_aligned_stack_new(Index, size); + Index* tags = ei_aligned_stack_new(Index, size); + + SparseMatrix<Scalar,ColMajor,Index> ap(size,size); + ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv); + + bool ok = true; + m_diag.resize(m_LDLt ? size : 0); + + for(Index k = 0; k < size; ++k) + { + // compute nonzero pattern of kth row of L, in topological order + y[k] = 0.0; // Y(0:k) is now all zero + Index top = size; // stack for pattern is empty + tags[k] = k; // mark node k as visited + m_nonZerosPerCol[k] = 0; // count of nonzeros in column k of L + for(typename MatrixType::InnerIterator it(ap,k); it; ++it) + { + Index i = it.index(); + if(i <= k) + { + y[i] += internal::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */ + Index len; + for(len = 0; tags[i] != k; i = m_parent[i]) + { + pattern[len++] = i; /* L(k,i) is nonzero */ + tags[i] = k; /* mark i as visited */ + } + while(len > 0) + pattern[--top] = pattern[--len]; + } + } + + /* compute numerical values kth row of L (a sparse triangular solve) */ + Scalar d = y[k]; // get D(k,k) and clear Y(k) + y[k] = 0.0; + for(; top < size; ++top) + { + Index i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */ + Scalar yi = y[i]; /* get and clear Y(i) */ + y[i] = 0.0; + + /* the nonzero entry L(k,i) */ + Scalar l_ki; + if(m_LDLt) + l_ki = yi / m_diag[i]; + else + yi = l_ki = yi / Lx[Lp[i]]; + + Index p2 = Lp[i] + m_nonZerosPerCol[i]; + Index p; + for(p = Lp[i] + (m_LDLt ? 0 : 1); p < p2; ++p) + y[Li[p]] -= internal::conj(Lx[p]) * yi; + d -= l_ki * internal::conj(yi); + Li[p] = k; /* store L(k,i) in column form of L */ + Lx[p] = l_ki; + ++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */ + } + if(m_LDLt) + m_diag[k] = d; + else + { + Index p = Lp[k]+m_nonZerosPerCol[k]++; + Li[p] = k ; /* store L(k,k) = sqrt (d) in column k */ + Lx[p] = internal::sqrt(d) ; + } + if(d == Scalar(0)) + { + ok = false; /* failure, D(k,k) is zero */ + break; + } + } + + // release workspace + ei_aligned_stack_delete(Scalar, y, size); + ei_aligned_stack_delete(Index, pattern, size); + ei_aligned_stack_delete(Index, tags, size); + + m_info = ok ? Success : NumericalIssue; + m_factorizationIsOk = true; +} + namespace internal { template<typename _MatrixType, int _UpLo, typename Rhs> |