/usr/include/eigen3/Eigen/src/SparseLU/SparseLU.h is in libeigen3-dev 3.2.2-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 | // This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_SPARSE_LU_H
#define EIGEN_SPARSE_LU_H
namespace Eigen {
template <typename _MatrixType, typename _OrderingType = COLAMDOrdering<typename _MatrixType::Index> > class SparseLU;
template <typename MappedSparseMatrixType> struct SparseLUMatrixLReturnType;
template <typename MatrixLType, typename MatrixUType> struct SparseLUMatrixUReturnType;
/** \ingroup SparseLU_Module
* \class SparseLU
*
* \brief Sparse supernodal LU factorization for general matrices
*
* This class implements the supernodal LU factorization for general matrices.
* It uses the main techniques from the sequential SuperLU package
* (http://crd-legacy.lbl.gov/~xiaoye/SuperLU/). It handles transparently real
* and complex arithmetics with single and double precision, depending on the
* scalar type of your input matrix.
* The code has been optimized to provide BLAS-3 operations during supernode-panel updates.
* It benefits directly from the built-in high-performant Eigen BLAS routines.
* Moreover, when the size of a supernode is very small, the BLAS calls are avoided to
* enable a better optimization from the compiler. For best performance,
* you should compile it with NDEBUG flag to avoid the numerous bounds checking on vectors.
*
* An important parameter of this class is the ordering method. It is used to reorder the columns
* (and eventually the rows) of the matrix to reduce the number of new elements that are created during
* numerical factorization. The cheapest method available is COLAMD.
* See \link OrderingMethods_Module the OrderingMethods module \endlink for the list of
* built-in and external ordering methods.
*
* Simple example with key steps
* \code
* VectorXd x(n), b(n);
* SparseMatrix<double, ColMajor> A;
* SparseLU<SparseMatrix<scalar, ColMajor>, COLAMDOrdering<Index> > solver;
* // fill A and b;
* // Compute the ordering permutation vector from the structural pattern of A
* solver.analyzePattern(A);
* // Compute the numerical factorization
* solver.factorize(A);
* //Use the factors to solve the linear system
* x = solver.solve(b);
* \endcode
*
* \warning The input matrix A should be in a \b compressed and \b column-major form.
* Otherwise an expensive copy will be made. You can call the inexpensive makeCompressed() to get a compressed matrix.
*
* \note Unlike the initial SuperLU implementation, there is no step to equilibrate the matrix.
* For badly scaled matrices, this step can be useful to reduce the pivoting during factorization.
* If this is the case for your matrices, you can try the basic scaling method at
* "unsupported/Eigen/src/IterativeSolvers/Scaling.h"
*
* \tparam _MatrixType The type of the sparse matrix. It must be a column-major SparseMatrix<>
* \tparam _OrderingType The ordering method to use, either AMD, COLAMD or METIS. Default is COLMAD
*
*
* \sa \ref TutorialSparseDirectSolvers
* \sa \ref OrderingMethods_Module
*/
template <typename _MatrixType, typename _OrderingType>
class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typename _MatrixType::Index>
{
public:
typedef _MatrixType MatrixType;
typedef _OrderingType OrderingType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index;
typedef SparseMatrix<Scalar,ColMajor,Index> NCMatrix;
typedef internal::MappedSuperNodalMatrix<Scalar, Index> SCMatrix;
typedef Matrix<Scalar,Dynamic,1> ScalarVector;
typedef Matrix<Index,Dynamic,1> IndexVector;
typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
typedef internal::SparseLUImpl<Scalar, Index> Base;
public:
SparseLU():m_isInitialized(true),m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1)
{
initperfvalues();
}
SparseLU(const MatrixType& matrix):m_isInitialized(true),m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1)
{
initperfvalues();
compute(matrix);
}
~SparseLU()
{
// Free all explicit dynamic pointers
}
void analyzePattern (const MatrixType& matrix);
void factorize (const MatrixType& matrix);
void simplicialfactorize(const MatrixType& matrix);
/**
* Compute the symbolic and numeric factorization of the input sparse matrix.
* The input matrix should be in column-major storage.
*/
void compute (const MatrixType& matrix)
{
// Analyze
analyzePattern(matrix);
//Factorize
factorize(matrix);
}
inline Index rows() const { return m_mat.rows(); }
inline Index cols() const { return m_mat.cols(); }
/** Indicate that the pattern of the input matrix is symmetric */
void isSymmetric(bool sym)
{
m_symmetricmode = sym;
}
/** \returns an expression of the matrix L, internally stored as supernodes
* The only operation available with this expression is the triangular solve
* \code
* y = b; matrixL().solveInPlace(y);
* \endcode
*/
SparseLUMatrixLReturnType<SCMatrix> matrixL() const
{
return SparseLUMatrixLReturnType<SCMatrix>(m_Lstore);
}
/** \returns an expression of the matrix U,
* The only operation available with this expression is the triangular solve
* \code
* y = b; matrixU().solveInPlace(y);
* \endcode
*/
SparseLUMatrixUReturnType<SCMatrix,MappedSparseMatrix<Scalar,ColMajor,Index> > matrixU() const
{
return SparseLUMatrixUReturnType<SCMatrix, MappedSparseMatrix<Scalar,ColMajor,Index> >(m_Lstore, m_Ustore);
}
/**
* \returns a reference to the row matrix permutation \f$ P_r \f$ such that \f$P_r A P_c^T = L U\f$
* \sa colsPermutation()
*/
inline const PermutationType& rowsPermutation() const
{
return m_perm_r;
}
/**
* \returns a reference to the column matrix permutation\f$ P_c^T \f$ such that \f$P_r A P_c^T = L U\f$
* \sa rowsPermutation()
*/
inline const PermutationType& colsPermutation() const
{
return m_perm_c;
}
/** Set the threshold used for a diagonal entry to be an acceptable pivot. */
void setPivotThreshold(const RealScalar& thresh)
{
m_diagpivotthresh = thresh;
}
/** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A.
*
* \warning the destination matrix X in X = this->solve(B) must be colmun-major.
*
* \sa compute()
*/
template<typename Rhs>
inline const internal::solve_retval<SparseLU, Rhs> solve(const MatrixBase<Rhs>& B) const
{
eigen_assert(m_factorizationIsOk && "SparseLU is not initialized.");
eigen_assert(rows()==B.rows()
&& "SparseLU::solve(): invalid number of rows of the right hand side matrix B");
return internal::solve_retval<SparseLU, Rhs>(*this, B.derived());
}
/** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A.
*
* \sa compute()
*/
template<typename Rhs>
inline const internal::sparse_solve_retval<SparseLU, Rhs> solve(const SparseMatrixBase<Rhs>& B) const
{
eigen_assert(m_factorizationIsOk && "SparseLU is not initialized.");
eigen_assert(rows()==B.rows()
&& "SparseLU::solve(): invalid number of rows of the right hand side matrix B");
return internal::sparse_solve_retval<SparseLU, Rhs>(*this, B.derived());
}
/** \brief Reports whether previous computation was successful.
*
* \returns \c Success if computation was succesful,
* \c NumericalIssue if the LU factorization reports a problem, zero diagonal for instance
* \c InvalidInput if the input matrix is invalid
*
* \sa iparm()
*/
ComputationInfo info() const
{
eigen_assert(m_isInitialized && "Decomposition is not initialized.");
return m_info;
}
/**
* \returns A string describing the type of error
*/
std::string lastErrorMessage() const
{
return m_lastError;
}
template<typename Rhs, typename Dest>
bool _solve(const MatrixBase<Rhs> &B, MatrixBase<Dest> &X_base) const
{
Dest& X(X_base.derived());
eigen_assert(m_factorizationIsOk && "The matrix should be factorized first");
EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
// Permute the right hand side to form X = Pr*B
// on return, X is overwritten by the computed solution
X.resize(B.rows(),B.cols());
// this ugly const_cast_derived() helps to detect aliasing when applying the permutations
for(Index j = 0; j < B.cols(); ++j)
X.col(j) = rowsPermutation() * B.const_cast_derived().col(j);
//Forward substitution with L
this->matrixL().solveInPlace(X);
this->matrixU().solveInPlace(X);
// Permute back the solution
for (Index j = 0; j < B.cols(); ++j)
X.col(j) = colsPermutation().inverse() * X.col(j);
return true;
}
/**
* \returns the absolute value of the determinant of the matrix of which
* *this is the QR decomposition.
*
* \warning a determinant can be very big or small, so for matrices
* of large enough dimension, there is a risk of overflow/underflow.
* One way to work around that is to use logAbsDeterminant() instead.
*
* \sa logAbsDeterminant(), signDeterminant()
*/
Scalar absDeterminant()
{
eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
// Initialize with the determinant of the row matrix
Scalar det = Scalar(1.);
//Note that the diagonal blocks of U are stored in supernodes,
// which are available in the L part :)
for (Index j = 0; j < this->cols(); ++j)
{
for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
{
if(it.row() < j) continue;
if(it.row() == j)
{
det *= (std::abs)(it.value());
break;
}
}
}
return det;
}
/** \returns the natural log of the absolute value of the determinant of the matrix
* of which **this is the QR decomposition
*
* \note This method is useful to work around the risk of overflow/underflow that's
* inherent to the determinant computation.
*
* \sa absDeterminant(), signDeterminant()
*/
Scalar logAbsDeterminant() const
{
eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
Scalar det = Scalar(0.);
for (Index j = 0; j < this->cols(); ++j)
{
for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
{
if(it.row() < j) continue;
if(it.row() == j)
{
det += (std::log)((std::abs)(it.value()));
break;
}
}
}
return det;
}
/** \returns A number representing the sign of the determinant
*
* \sa absDeterminant(), logAbsDeterminant()
*/
Scalar signDeterminant()
{
eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
return Scalar(m_detPermR);
}
protected:
// Functions
void initperfvalues()
{
m_perfv.panel_size = 1;
m_perfv.relax = 1;
m_perfv.maxsuper = 128;
m_perfv.rowblk = 16;
m_perfv.colblk = 8;
m_perfv.fillfactor = 20;
}
// Variables
mutable ComputationInfo m_info;
bool m_isInitialized;
bool m_factorizationIsOk;
bool m_analysisIsOk;
std::string m_lastError;
NCMatrix m_mat; // The input (permuted ) matrix
SCMatrix m_Lstore; // The lower triangular matrix (supernodal)
MappedSparseMatrix<Scalar,ColMajor,Index> m_Ustore; // The upper triangular matrix
PermutationType m_perm_c; // Column permutation
PermutationType m_perm_r ; // Row permutation
IndexVector m_etree; // Column elimination tree
typename Base::GlobalLU_t m_glu;
// SparseLU options
bool m_symmetricmode;
// values for performance
internal::perfvalues<Index> m_perfv;
RealScalar m_diagpivotthresh; // Specifies the threshold used for a diagonal entry to be an acceptable pivot
Index m_nnzL, m_nnzU; // Nonzeros in L and U factors
Index m_detPermR; // Determinant of the coefficient matrix
private:
// Disable copy constructor
SparseLU (const SparseLU& );
}; // End class SparseLU
// Functions needed by the anaysis phase
/**
* Compute the column permutation to minimize the fill-in
*
* - Apply this permutation to the input matrix -
*
* - Compute the column elimination tree on the permuted matrix
*
* - Postorder the elimination tree and the column permutation
*
*/
template <typename MatrixType, typename OrderingType>
void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
{
//TODO It is possible as in SuperLU to compute row and columns scaling vectors to equilibrate the matrix mat.
OrderingType ord;
ord(mat,m_perm_c);
// Apply the permutation to the column of the input matrix
//First copy the whole input matrix.
m_mat = mat;
if (m_perm_c.size()) {
m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. FIXME : This vector is filled but not subsequently used.
//Then, permute only the column pointers
const Index * outerIndexPtr;
if (mat.isCompressed()) outerIndexPtr = mat.outerIndexPtr();
else
{
Index *outerIndexPtr_t = new Index[mat.cols()+1];
for(Index i = 0; i <= mat.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i];
outerIndexPtr = outerIndexPtr_t;
}
for (Index i = 0; i < mat.cols(); i++)
{
m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i];
}
if(!mat.isCompressed()) delete[] outerIndexPtr;
}
// Compute the column elimination tree of the permuted matrix
IndexVector firstRowElt;
internal::coletree(m_mat, m_etree,firstRowElt);
// In symmetric mode, do not do postorder here
if (!m_symmetricmode) {
IndexVector post, iwork;
// Post order etree
internal::treePostorder(m_mat.cols(), m_etree, post);
// Renumber etree in postorder
Index m = m_mat.cols();
iwork.resize(m+1);
for (Index i = 0; i < m; ++i) iwork(post(i)) = post(m_etree(i));
m_etree = iwork;
// Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree
PermutationType post_perm(m);
for (Index i = 0; i < m; i++)
post_perm.indices()(i) = post(i);
// Combine the two permutations : postorder the permutation for future use
if(m_perm_c.size()) {
m_perm_c = post_perm * m_perm_c;
}
} // end postordering
m_analysisIsOk = true;
}
// Functions needed by the numerical factorization phase
/**
* - Numerical factorization
* - Interleaved with the symbolic factorization
* On exit, info is
*
* = 0: successful factorization
*
* > 0: if info = i, and i is
*
* <= A->ncol: U(i,i) is exactly zero. The factorization has
* been completed, but the factor U is exactly singular,
* and division by zero will occur if it is used to solve a
* system of equations.
*
* > A->ncol: number of bytes allocated when memory allocation
* failure occurred, plus A->ncol. If lwork = -1, it is
* the estimated amount of space needed, plus A->ncol.
*/
template <typename MatrixType, typename OrderingType>
void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
{
using internal::emptyIdxLU;
eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
eigen_assert((matrix.rows() == matrix.cols()) && "Only for squared matrices");
typedef typename IndexVector::Scalar Index;
// Apply the column permutation computed in analyzepattern()
// m_mat = matrix * m_perm_c.inverse();
m_mat = matrix;
if (m_perm_c.size())
{
m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers.
//Then, permute only the column pointers
const Index * outerIndexPtr;
if (matrix.isCompressed()) outerIndexPtr = matrix.outerIndexPtr();
else
{
Index* outerIndexPtr_t = new Index[matrix.cols()+1];
for(Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i];
outerIndexPtr = outerIndexPtr_t;
}
for (Index i = 0; i < matrix.cols(); i++)
{
m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i];
}
if(!matrix.isCompressed()) delete[] outerIndexPtr;
}
else
{ //FIXME This should not be needed if the empty permutation is handled transparently
m_perm_c.resize(matrix.cols());
for(Index i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i;
}
Index m = m_mat.rows();
Index n = m_mat.cols();
Index nnz = m_mat.nonZeros();
Index maxpanel = m_perfv.panel_size * m;
// Allocate working storage common to the factor routines
Index lwork = 0;
Index info = Base::memInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu);
if (info)
{
m_lastError = "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ;
m_factorizationIsOk = false;
return ;
}
// Set up pointers for integer working arrays
IndexVector segrep(m); segrep.setZero();
IndexVector parent(m); parent.setZero();
IndexVector xplore(m); xplore.setZero();
IndexVector repfnz(maxpanel);
IndexVector panel_lsub(maxpanel);
IndexVector xprune(n); xprune.setZero();
IndexVector marker(m*internal::LUNoMarker); marker.setZero();
repfnz.setConstant(-1);
panel_lsub.setConstant(-1);
// Set up pointers for scalar working arrays
ScalarVector dense;
dense.setZero(maxpanel);
ScalarVector tempv;
tempv.setZero(internal::LUnumTempV(m, m_perfv.panel_size, m_perfv.maxsuper, /*m_perfv.rowblk*/m) );
// Compute the inverse of perm_c
PermutationType iperm_c(m_perm_c.inverse());
// Identify initial relaxed snodes
IndexVector relax_end(n);
if ( m_symmetricmode == true )
Base::heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
else
Base::relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
m_perm_r.resize(m);
m_perm_r.indices().setConstant(-1);
marker.setConstant(-1);
m_detPermR = 1; // Record the determinant of the row permutation
m_glu.supno(0) = emptyIdxLU; m_glu.xsup.setConstant(0);
m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) = Index(0);
// Work on one 'panel' at a time. A panel is one of the following :
// (a) a relaxed supernode at the bottom of the etree, or
// (b) panel_size contiguous columns, <panel_size> defined by the user
Index jcol;
IndexVector panel_histo(n);
Index pivrow; // Pivotal row number in the original row matrix
Index nseg1; // Number of segments in U-column above panel row jcol
Index nseg; // Number of segments in each U-column
Index irep;
Index i, k, jj;
for (jcol = 0; jcol < n; )
{
// Adjust panel size so that a panel won't overlap with the next relaxed snode.
Index panel_size = m_perfv.panel_size; // upper bound on panel width
for (k = jcol + 1; k < (std::min)(jcol+panel_size, n); k++)
{
if (relax_end(k) != emptyIdxLU)
{
panel_size = k - jcol;
break;
}
}
if (k == n)
panel_size = n - jcol;
// Symbolic outer factorization on a panel of columns
Base::panel_dfs(m, panel_size, jcol, m_mat, m_perm_r.indices(), nseg1, dense, panel_lsub, segrep, repfnz, xprune, marker, parent, xplore, m_glu);
// Numeric sup-panel updates in topological order
Base::panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_glu);
// Sparse LU within the panel, and below the panel diagonal
for ( jj = jcol; jj< jcol + panel_size; jj++)
{
k = (jj - jcol) * m; // Column index for w-wide arrays
nseg = nseg1; // begin after all the panel segments
//Depth-first-search for the current column
VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m);
VectorBlock<IndexVector> repfnz_k(repfnz, k, m);
info = Base::column_dfs(m, jj, m_perm_r.indices(), m_perfv.maxsuper, nseg, panel_lsubk, segrep, repfnz_k, xprune, marker, parent, xplore, m_glu);
if ( info )
{
m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() ";
m_info = NumericalIssue;
m_factorizationIsOk = false;
return;
}
// Numeric updates to this column
VectorBlock<ScalarVector> dense_k(dense, k, m);
VectorBlock<IndexVector> segrep_k(segrep, nseg1, m-nseg1);
info = Base::column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu);
if ( info )
{
m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() ";
m_info = NumericalIssue;
m_factorizationIsOk = false;
return;
}
// Copy the U-segments to ucol(*)
info = Base::copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu);
if ( info )
{
m_lastError = "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() ";
m_info = NumericalIssue;
m_factorizationIsOk = false;
return;
}
// Form the L-segment
info = Base::pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu);
if ( info )
{
m_lastError = "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT ";
std::ostringstream returnInfo;
returnInfo << info;
m_lastError += returnInfo.str();
m_info = NumericalIssue;
m_factorizationIsOk = false;
return;
}
// Update the determinant of the row permutation matrix
if (pivrow != jj) m_detPermR *= -1;
// Prune columns (0:jj-1) using column jj
Base::pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu);
// Reset repfnz for this column
for (i = 0; i < nseg; i++)
{
irep = segrep(i);
repfnz_k(irep) = emptyIdxLU;
}
} // end SparseLU within the panel
jcol += panel_size; // Move to the next panel
} // end for -- end elimination
// Count the number of nonzeros in factors
Base::countnz(n, m_nnzL, m_nnzU, m_glu);
// Apply permutation to the L subscripts
Base::fixupL(n, m_perm_r.indices(), m_glu);
// Create supernode matrix L
m_Lstore.setInfos(m, n, m_glu.lusup, m_glu.xlusup, m_glu.lsub, m_glu.xlsub, m_glu.supno, m_glu.xsup);
// Create the column major upper sparse matrix U;
new (&m_Ustore) MappedSparseMatrix<Scalar, ColMajor, Index> ( m, n, m_nnzU, m_glu.xusub.data(), m_glu.usub.data(), m_glu.ucol.data() );
m_info = Success;
m_factorizationIsOk = true;
}
template<typename MappedSupernodalType>
struct SparseLUMatrixLReturnType : internal::no_assignment_operator
{
typedef typename MappedSupernodalType::Index Index;
typedef typename MappedSupernodalType::Scalar Scalar;
SparseLUMatrixLReturnType(const MappedSupernodalType& mapL) : m_mapL(mapL)
{ }
Index rows() { return m_mapL.rows(); }
Index cols() { return m_mapL.cols(); }
template<typename Dest>
void solveInPlace( MatrixBase<Dest> &X) const
{
m_mapL.solveInPlace(X);
}
const MappedSupernodalType& m_mapL;
};
template<typename MatrixLType, typename MatrixUType>
struct SparseLUMatrixUReturnType : internal::no_assignment_operator
{
typedef typename MatrixLType::Index Index;
typedef typename MatrixLType::Scalar Scalar;
SparseLUMatrixUReturnType(const MatrixLType& mapL, const MatrixUType& mapU)
: m_mapL(mapL),m_mapU(mapU)
{ }
Index rows() { return m_mapL.rows(); }
Index cols() { return m_mapL.cols(); }
template<typename Dest> void solveInPlace(MatrixBase<Dest> &X) const
{
Index nrhs = X.cols();
Index n = X.rows();
// Backward solve with U
for (Index k = m_mapL.nsuper(); k >= 0; k--)
{
Index fsupc = m_mapL.supToCol()[k];
Index lda = m_mapL.colIndexPtr()[fsupc+1] - m_mapL.colIndexPtr()[fsupc]; // leading dimension
Index nsupc = m_mapL.supToCol()[k+1] - fsupc;
Index luptr = m_mapL.colIndexPtr()[fsupc];
if (nsupc == 1)
{
for (Index j = 0; j < nrhs; j++)
{
X(fsupc, j) /= m_mapL.valuePtr()[luptr];
}
}
else
{
Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) );
Map< Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
U = A.template triangularView<Upper>().solve(U);
}
for (Index j = 0; j < nrhs; ++j)
{
for (Index jcol = fsupc; jcol < fsupc + nsupc; jcol++)
{
typename MatrixUType::InnerIterator it(m_mapU, jcol);
for ( ; it; ++it)
{
Index irow = it.index();
X(irow, j) -= X(jcol, j) * it.value();
}
}
}
} // End For U-solve
}
const MatrixLType& m_mapL;
const MatrixUType& m_mapU;
};
namespace internal {
template<typename _MatrixType, typename Derived, typename Rhs>
struct solve_retval<SparseLU<_MatrixType,Derived>, Rhs>
: solve_retval_base<SparseLU<_MatrixType,Derived>, Rhs>
{
typedef SparseLU<_MatrixType,Derived> Dec;
EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
template<typename Dest> void evalTo(Dest& dst) const
{
dec()._solve(rhs(),dst);
}
};
template<typename _MatrixType, typename Derived, typename Rhs>
struct sparse_solve_retval<SparseLU<_MatrixType,Derived>, Rhs>
: sparse_solve_retval_base<SparseLU<_MatrixType,Derived>, Rhs>
{
typedef SparseLU<_MatrixType,Derived> Dec;
EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
template<typename Dest> void evalTo(Dest& dst) const
{
this->defaultEvalTo(dst);
}
};
} // end namespace internal
} // End namespace Eigen
#endif
|