Repository navigation
Expand file tree
/
Copy pathSparseLU.h
More file actions
775 lines (688 loc) · 27.3 KB
/
Copy pathSparseLU.h
File metadata and controls
775 lines (688 loc) · 27.3 KB
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
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
// 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-2014 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::StorageIndex> > 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
*
* \implsparsesolverconcept
*
* \sa \ref TutorialSparseSolverConcept
* \sa \ref OrderingMethods_Module
*/
template <typename _MatrixType, typename _OrderingType>
class SparseLU : public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >, public internal::SparseLUImpl<typename _MatrixType::Scalar, typename _MatrixType::StorageIndex>
{
protected:
typedef SparseSolverBase<SparseLU<_MatrixType,_OrderingType> > APIBase;
using APIBase::m_isInitialized;
public:
using APIBase::_solve_impl;
typedef _MatrixType MatrixType;
typedef _OrderingType OrderingType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::StorageIndex StorageIndex;
typedef SparseMatrix<Scalar,ColMajor,StorageIndex> NCMatrix;
typedef internal::MappedSuperNodalMatrix<Scalar, StorageIndex> SCMatrix;
typedef Matrix<Scalar,Dynamic,1> ScalarVector;
typedef Matrix<StorageIndex,Dynamic,1> IndexVector;
typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
typedef internal::SparseLUImpl<Scalar, StorageIndex> Base;
enum {
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
public:
SparseLU():m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1)
{
initperfvalues();
}
explicit SparseLU(const MatrixType& matrix)
: 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,StorageIndex> > matrixU() const
{
return SparseLUMatrixUReturnType<SCMatrix, MappedSparseMatrix<Scalar,ColMajor,StorageIndex> >(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;
}
#ifdef EIGEN_PARSED_BY_DOXYGEN
/** \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 Solve<SparseLU, Rhs> solve(const MatrixBase<Rhs>& B) const;
#endif // EIGEN_PARSED_BY_DOXYGEN
/** \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_impl(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()
{
using std::abs;
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.index() == j)
{
det *= 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
{
using std::log;
using std::abs;
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 += log(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.");
// Initialize with the determinant of the row matrix
Index det = 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.index() == j)
{
if(it.value()<0)
det = -det;
else if(it.value()==0)
return 0;
break;
}
}
}
return det * m_detPermR * m_detPermC;
}
/** \returns The determinant of the matrix.
*
* \sa absDeterminant(), logAbsDeterminant()
*/
Scalar determinant()
{
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.index() == j)
{
det *= it.value();
break;
}
}
}
return (m_detPermR * m_detPermC) > 0 ? det : -det;
}
protected:
// Functions
void initperfvalues()
{
m_perfv.panel_size = 16;
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_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,StorageIndex> 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 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, m_detPermC; // Determinants of the permutation matrices
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.
// Firstly, copy the whole input matrix.
m_mat = mat;
// Compute fill-in ordering
OrderingType ord;
ord(m_mat,m_perm_c);
// Apply the permutation to the column of the input matrix
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
ei_declare_aligned_stack_constructed_variable(StorageIndex,outerIndexPtr,mat.cols()+1,mat.isCompressed()?const_cast<StorageIndex*>(mat.outerIndexPtr()):0);
// If the input matrix 'mat' is uncompressed, then the outer-indices do not match the ones of m_mat, and a copy is thus needed.
if(!mat.isCompressed())
IndexVector::Map(outerIndexPtr, mat.cols()+1) = IndexVector::Map(m_mat.outerIndexPtr(),mat.cols()+1);
// Apply the permutation and compute the nnz per column.
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];
}
}
// 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(StorageIndex(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 StorageIndex;
m_isInitialized = true;
// 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 StorageIndex * outerIndexPtr;
if (matrix.isCompressed()) outerIndexPtr = matrix.outerIndexPtr();
else
{
StorageIndex* outerIndexPtr_t = new StorageIndex[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(StorageIndex 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
// FIXME: the following test is not correct, we should probably take iperm_c into account and pivrow is not directly the row pivot.
if (pivrow != jj) m_detPermR = -m_detPermR;
// 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
m_detPermR = m_perm_r.determinant();
m_detPermC = m_perm_c.determinant();
// 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, StorageIndex> ( 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::Scalar Scalar;
explicit 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::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, ColMajor>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) );
Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 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;
};
} // End namespace Eigen
#endif