mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-01-18 14:34:17 +08:00
clean a bit AMD and SimplicialCholesky and add support for partly stored selfadjoint matrices
This commit is contained in:
Gael Guennebaud
parent
1618df55df
commit
3e99356b59
@ -51,7 +51,7 @@ enum {
|
||||
OrderingMask = 0x0f00
|
||||
};
|
||||
|
||||
#include "src/misc/Solve.h"
|
||||
#include "../../Eigen/src/misc/Solve.h"
|
||||
|
||||
#include "src/SparseExtra/RandomSetter.h"
|
||||
#include "src/SparseExtra/Solve.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
|
||||
|
||||
@ -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)
|
||||
@ -352,6 +353,10 @@ class CholmodDecomposition
|
||||
}
|
||||
#endif // EIGEN_PARSED_BY_DOXYGEN
|
||||
|
||||
template<typename Stream>
|
||||
void dumpMemory(Stream& s)
|
||||
{}
|
||||
|
||||
protected:
|
||||
mutable cholmod_common m_cholmod;
|
||||
cholmod_factor* m_cholmodFactor;
|
||||
|
||||
@ -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);
|
||||
void analyzePattern(const MatrixType& a);
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
/** 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);
|
||||
void factorize(const MatrixType& a);
|
||||
|
||||
const Index* Lp = m_matrix._outerIndexPtr();
|
||||
Index* Li = m_matrix._innerIndexPtr();
|
||||
Scalar* Lx = m_matrix._valuePtr();
|
||||
/** \returns the permutation P
|
||||
* \sa permutationPinv() */
|
||||
const PermutationMatrix<Dynamic>& permutationP() const
|
||||
{ return m_P; }
|
||||
|
||||
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;
|
||||
}
|
||||
/** \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 */
|
||||
@ -409,7 +250,28 @@ 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>
|
||||
|
||||
@ -96,30 +96,30 @@ template<typename Scalar> void sparse_ldlt(int rows, int cols)
|
||||
x = SimplicialCholesky<SparseMatrix<Scalar>, Upper>().setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt).compute(m3).solve(b);
|
||||
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, upper, single dense rhs");
|
||||
|
||||
// x = SimplicialCholesky<SparseMatrix<Scalar>, Lower>(m3_lo).solve(b);
|
||||
// VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, lower only, single dense rhs");
|
||||
x = SimplicialCholesky<SparseMatrix<Scalar>, Lower>(m3_lo).solve(b);
|
||||
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, lower only, single dense rhs");
|
||||
|
||||
// x = SimplicialCholesky<SparseMatrix<Scalar>, Upper>(m3_up).solve(b);
|
||||
// VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, upper only, single dense rhs");
|
||||
x = SimplicialCholesky<SparseMatrix<Scalar>, Upper>(m3_up).solve(b);
|
||||
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, upper only, single dense rhs");
|
||||
|
||||
|
||||
// with multiple rhs
|
||||
ref_X = refMat3.template selfadjointView<Lower>().llt().solve(B);
|
||||
|
||||
X = SimplicialCholesky<SparseMatrix<Scalar>, Lower>()/*.setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt)*/.compute(m3).solve(B);
|
||||
X = SimplicialCholesky<SparseMatrix<Scalar>, Lower>().setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt).compute(m3).solve(B);
|
||||
VERIFY(ref_X.isApprox(X,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, lower, multiple dense rhs");
|
||||
|
||||
// X = SimplicialCholesky<SparseMatrix<Scalar>, Upper>().setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt).compute(m3).solve(B);
|
||||
// VERIFY(ref_X.isApprox(X,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, upper, multiple dense rhs");
|
||||
X = SimplicialCholesky<SparseMatrix<Scalar>, Upper>().setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt).compute(m3).solve(B);
|
||||
VERIFY(ref_X.isApprox(X,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, upper, multiple dense rhs");
|
||||
|
||||
|
||||
// // with a sparse rhs
|
||||
// with a sparse rhs
|
||||
// SparseMatrix<Scalar> spB(rows,cols), spX(rows,cols);
|
||||
// B.diagonal().array() += 1;
|
||||
// spB = B.sparseView(0.5,1);
|
||||
//
|
||||
// ref_X = refMat3.template selfadjointView<Lower>().llt().solve(DenseMatrix(spB));
|
||||
|
||||
//
|
||||
// spX = SimplicialCholesky<SparseMatrix<Scalar>, Lower>(m3).solve(spB);
|
||||
// VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: cholmod solve, multiple sparse rhs");
|
||||
//
|
||||
|
||||
Loading…
Reference in New Issue
Block a user