← Back to branch summary

~ubuntu-branches/ubuntu/trusty/blender/trusty

« back to all changes in this revision

Viewing changes to extern/Eigen3/Eigen/src/SuperLUSupport/SuperLUSupport.h

  • Committer: Package Import Robot
  • Author(s): Jeremy Bicha
  • Date: 2013-03-06 12:08:47 UTC
  • mfrom: (1.5.1) (14.1.8 experimental)
  • Revision ID: package-import@ubuntu.com-20130306120847-frjfaryb2zrotwcg
Tags: 2.66a-1ubuntu1
https://launchpad.net/bugs/1076930
https://launchpad.net/bugs/1089256
https://launchpad.net/bugs/1052743
https://launchpad.net/bugs/999024
* Resynchronize with Debian (LP: #1076930, #1089256, #1052743, #999024,
  #1122888, #1147084)
* debian/control:
  - Lower build-depends on libavcodec-dev since we're not
    doing the libav9 transition in Ubuntu yet

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/Eigen3/Eigen/src/SuperLUSupport/SuperLUSupport.h
 
1
// This file is part of Eigen, a lightweight C++ template library
 
2
// for linear algebra.
 
3
//
 
4
// Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
 
5
//
 
6
// This Source Code Form is subject to the terms of the Mozilla
 
7
// Public License v. 2.0. If a copy of the MPL was not distributed
 
8
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
9
 
 
10
#ifndef EIGEN_SUPERLUSUPPORT_H
 
11
#define EIGEN_SUPERLUSUPPORT_H
 
12
 
 
13
namespace Eigen { 
 
14
 
 
15
#define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE)            \
 
16
    extern "C" {                                                                                          \
 
17
      typedef struct { FLOATTYPE for_lu; FLOATTYPE total_needed; int expansions; } PREFIX##mem_usage_t;   \
 
18
      extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *,                  \
 
19
                                char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *,           \
 
20
                                void *, int, SuperMatrix *, SuperMatrix *,                                \
 
21
                                FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *,                       \
 
22
                                PREFIX##mem_usage_t *, SuperLUStat_t *, int *);                           \
 
23
    }                                                                                                     \
 
24
    inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A,                                \
 
25
         int *perm_c, int *perm_r, int *etree, char *equed,                                               \
 
26
         FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L,                                                      \
 
27
         SuperMatrix *U, void *work, int lwork,                                                           \
 
28
         SuperMatrix *B, SuperMatrix *X,                                                                  \
 
29
         FLOATTYPE *recip_pivot_growth,                                                                   \
 
30
         FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr,                                              \
 
31
         SuperLUStat_t *stats, int *info, KEYTYPE) {                                                      \
 
32
    PREFIX##mem_usage_t mem_usage;                                                                        \
 
33
    PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L,                                      \
 
34
         U, work, lwork, B, X, recip_pivot_growth, rcond,                                                 \
 
35
         ferr, berr, &mem_usage, stats, info);                                                            \
 
36
    return mem_usage.for_lu; /* bytes used by the factor storage */                                       \
 
37
  }
 
38
 
 
39
DECL_GSSVX(s,float,float)
 
40
DECL_GSSVX(c,float,std::complex<float>)
 
41
DECL_GSSVX(d,double,double)
 
42
DECL_GSSVX(z,double,std::complex<double>)
 
43
 
 
44
#ifdef MILU_ALPHA
 
45
#define EIGEN_SUPERLU_HAS_ILU
 
46
#endif
 
47
 
 
48
#ifdef EIGEN_SUPERLU_HAS_ILU
 
49
 
 
50
// similarly for the incomplete factorization using gsisx
 
51
#define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE)                                                    \
 
52
    extern "C" {                                                                                \
 
53
      extern void PREFIX##gsisx(superlu_options_t *, SuperMatrix *, int *, int *, int *,        \
 
54
                         char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *,        \
 
55
                         void *, int, SuperMatrix *, SuperMatrix *, FLOATTYPE *, FLOATTYPE *,   \
 
56
                         PREFIX##mem_usage_t *, SuperLUStat_t *, int *);                        \
 
57
    }                                                                                           \
 
58
    inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A,                      \
 
59
         int *perm_c, int *perm_r, int *etree, char *equed,                                     \
 
60
         FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L,                                            \
 
61
         SuperMatrix *U, void *work, int lwork,                                                 \
 
62
         SuperMatrix *B, SuperMatrix *X,                                                        \
 
63
         FLOATTYPE *recip_pivot_growth,                                                         \
 
64
         FLOATTYPE *rcond,                                                                      \
 
65
         SuperLUStat_t *stats, int *info, KEYTYPE) {                                            \
 
66
    PREFIX##mem_usage_t mem_usage;                                                              \
 
67
    PREFIX##gsisx(options, A, perm_c, perm_r, etree, equed, R, C, L,                            \
 
68
         U, work, lwork, B, X, recip_pivot_growth, rcond,                                       \
 
69
         &mem_usage, stats, info);                                                              \
 
70
    return mem_usage.for_lu; /* bytes used by the factor storage */                             \
 
71
  }
 
72
 
 
73
DECL_GSISX(s,float,float)
 
74
DECL_GSISX(c,float,std::complex<float>)
 
75
DECL_GSISX(d,double,double)
 
76
DECL_GSISX(z,double,std::complex<double>)
 
77
 
 
78
#endif
 
79
 
 
80
template<typename MatrixType>
 
81
struct SluMatrixMapHelper;
 
82
 
 
83
/** \internal
 
84
  *
 
85
  * A wrapper class for SuperLU matrices. It supports only compressed sparse matrices
 
86
  * and dense matrices. Supernodal and other fancy format are not supported by this wrapper.
 
87
  *
 
88
  * This wrapper class mainly aims to avoids the need of dynamic allocation of the storage structure.
 
89
  */
 
90
struct SluMatrix : SuperMatrix
 
91
{
 
92
  SluMatrix()
 
93
  {
 
94
    Store = &storage;
 
95
  }
 
96
 
 
97
  SluMatrix(const SluMatrix& other)
 
98
    : SuperMatrix(other)
 
99
  {
 
100
    Store = &storage;
 
101
    storage = other.storage;
 
102
  }
 
103
 
 
104
  SluMatrix& operator=(const SluMatrix& other)
 
105
  {
 
106
    SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
 
107
    Store = &storage;
 
108
    storage = other.storage;
 
109
    return *this;
 
110
  }
 
111
 
 
112
  struct
 
113
  {
 
114
    union {int nnz;int lda;};
 
115
    void *values;
 
116
    int *innerInd;
 
117
    int *outerInd;
 
118
  } storage;
 
119
 
 
120
  void setStorageType(Stype_t t)
 
121
  {
 
122
    Stype = t;
 
123
    if (t==SLU_NC || t==SLU_NR || t==SLU_DN)
 
124
      Store = &storage;
 
125
    else
 
126
    {
 
127
      eigen_assert(false && "storage type not supported");
 
128
      Store = 0;
 
129
    }
 
130
  }
 
131
 
 
132
  template<typename Scalar>
 
133
  void setScalarType()
 
134
  {
 
135
    if (internal::is_same<Scalar,float>::value)
 
136
      Dtype = SLU_S;
 
137
    else if (internal::is_same<Scalar,double>::value)
 
138
      Dtype = SLU_D;
 
139
    else if (internal::is_same<Scalar,std::complex<float> >::value)
 
140
      Dtype = SLU_C;
 
141
    else if (internal::is_same<Scalar,std::complex<double> >::value)
 
142
      Dtype = SLU_Z;
 
143
    else
 
144
    {
 
145
      eigen_assert(false && "Scalar type not supported by SuperLU");
 
146
    }
 
147
  }
 
148
 
 
149
  template<typename MatrixType>
 
150
  static SluMatrix Map(MatrixBase<MatrixType>& _mat)
 
151
  {
 
152
    MatrixType& mat(_mat.derived());
 
153
    eigen_assert( ((MatrixType::Flags&RowMajorBit)!=RowMajorBit) && "row-major dense matrices are not supported by SuperLU");
 
154
    SluMatrix res;
 
155
    res.setStorageType(SLU_DN);
 
156
    res.setScalarType<typename MatrixType::Scalar>();
 
157
    res.Mtype     = SLU_GE;
 
158
 
 
159
    res.nrow      = mat.rows();
 
160
    res.ncol      = mat.cols();
 
161
 
 
162
    res.storage.lda       = MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride();
 
163
    res.storage.values    = mat.data();
 
164
    return res;
 
165
  }
 
166
 
 
167
  template<typename MatrixType>
 
168
  static SluMatrix Map(SparseMatrixBase<MatrixType>& mat)
 
169
  {
 
170
    SluMatrix res;
 
171
    if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
 
172
    {
 
173
      res.setStorageType(SLU_NR);
 
174
      res.nrow      = mat.cols();
 
175
      res.ncol      = mat.rows();
 
176
    }
 
177
    else
 
178
    {
 
179
      res.setStorageType(SLU_NC);
 
180
      res.nrow      = mat.rows();
 
181
      res.ncol      = mat.cols();
 
182
    }
 
183
 
 
184
    res.Mtype       = SLU_GE;
 
185
 
 
186
    res.storage.nnz       = mat.nonZeros();
 
187
    res.storage.values    = mat.derived().valuePtr();
 
188
    res.storage.innerInd  = mat.derived().innerIndexPtr();
 
189
    res.storage.outerInd  = mat.derived().outerIndexPtr();
 
190
 
 
191
    res.setScalarType<typename MatrixType::Scalar>();
 
192
 
 
193
    // FIXME the following is not very accurate
 
194
    if (MatrixType::Flags & Upper)
 
195
      res.Mtype = SLU_TRU;
 
196
    if (MatrixType::Flags & Lower)
 
197
      res.Mtype = SLU_TRL;
 
198
 
 
199
    eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
 
200
 
 
201
    return res;
 
202
  }
 
203
};
 
204
 
 
205
template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
 
206
struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,Options,MRows,MCols> >
 
207
{
 
208
  typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType;
 
209
  static void run(MatrixType& mat, SluMatrix& res)
 
210
  {
 
211
    eigen_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU");
 
212
    res.setStorageType(SLU_DN);
 
213
    res.setScalarType<Scalar>();
 
214
    res.Mtype     = SLU_GE;
 
215
 
 
216
    res.nrow      = mat.rows();
 
217
    res.ncol      = mat.cols();
 
218
 
 
219
    res.storage.lda       = mat.outerStride();
 
220
    res.storage.values    = mat.data();
 
221
  }
 
222
};
 
223
 
 
224
template<typename Derived>
 
225
struct SluMatrixMapHelper<SparseMatrixBase<Derived> >
 
226
{
 
227
  typedef Derived MatrixType;
 
228
  static void run(MatrixType& mat, SluMatrix& res)
 
229
  {
 
230
    if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
 
231
    {
 
232
      res.setStorageType(SLU_NR);
 
233
      res.nrow      = mat.cols();
 
234
      res.ncol      = mat.rows();
 
235
    }
 
236
    else
 
237
    {
 
238
      res.setStorageType(SLU_NC);
 
239
      res.nrow      = mat.rows();
 
240
      res.ncol      = mat.cols();
 
241
    }
 
242
 
 
243
    res.Mtype       = SLU_GE;
 
244
 
 
245
    res.storage.nnz       = mat.nonZeros();
 
246
    res.storage.values    = mat.valuePtr();
 
247
    res.storage.innerInd  = mat.innerIndexPtr();
 
248
    res.storage.outerInd  = mat.outerIndexPtr();
 
249
 
 
250
    res.setScalarType<typename MatrixType::Scalar>();
 
251
 
 
252
    // FIXME the following is not very accurate
 
253
    if (MatrixType::Flags & Upper)
 
254
      res.Mtype = SLU_TRU;
 
255
    if (MatrixType::Flags & Lower)
 
256
      res.Mtype = SLU_TRL;
 
257
 
 
258
    eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
 
259
  }
 
260
};
 
261
 
 
262
namespace internal {
 
263
 
 
264
template<typename MatrixType>
 
265
SluMatrix asSluMatrix(MatrixType& mat)
 
266
{
 
267
  return SluMatrix::Map(mat);
 
268
}
 
269
 
 
270
/** View a Super LU matrix as an Eigen expression */
 
271
template<typename Scalar, int Flags, typename Index>
 
272
MappedSparseMatrix<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat)
 
273
{
 
274
  eigen_assert((Flags&RowMajor)==RowMajor && sluMat.Stype == SLU_NR
 
275
         || (Flags&ColMajor)==ColMajor && sluMat.Stype == SLU_NC);
 
276
 
 
277
  Index outerSize = (Flags&RowMajor)==RowMajor ? sluMat.ncol : sluMat.nrow;
 
278
 
 
279
  return MappedSparseMatrix<Scalar,Flags,Index>(
 
280
    sluMat.nrow, sluMat.ncol, sluMat.storage.outerInd[outerSize],
 
281
    sluMat.storage.outerInd, sluMat.storage.innerInd, reinterpret_cast<Scalar*>(sluMat.storage.values) );
 
282
}
 
283
 
 
284
} // end namespace internal
 
285
 
 
286
/** \ingroup SuperLUSupport_Module
 
287
  * \class SuperLUBase
 
288
  * \brief The base class for the direct and incomplete LU factorization of SuperLU
 
289
  */
 
290
template<typename _MatrixType, typename Derived>
 
291
class SuperLUBase : internal::noncopyable
 
292
{
 
293
  public:
 
294
    typedef _MatrixType MatrixType;
 
295
    typedef typename MatrixType::Scalar Scalar;
 
296
    typedef typename MatrixType::RealScalar RealScalar;
 
297
    typedef typename MatrixType::Index Index;
 
298
    typedef Matrix<Scalar,Dynamic,1> Vector;
 
299
    typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
 
300
    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;    
 
301
    typedef SparseMatrix<Scalar> LUMatrixType;
 
302
 
 
303
  public:
 
304
 
 
305
    SuperLUBase() {}
 
306
 
 
307
    ~SuperLUBase()
 
308
    {
 
309
      clearFactors();
 
310
    }
 
311
    
 
312
    Derived& derived() { return *static_cast<Derived*>(this); }
 
313
    const Derived& derived() const { return *static_cast<const Derived*>(this); }
 
314
    
 
315
    inline Index rows() const { return m_matrix.rows(); }
 
316
    inline Index cols() const { return m_matrix.cols(); }
 
317
    
 
318
    /** \returns a reference to the Super LU option object to configure the  Super LU algorithms. */
 
319
    inline superlu_options_t& options() { return m_sluOptions; }
 
320
    
 
321
    /** \brief Reports whether previous computation was successful.
 
322
      *
 
323
      * \returns \c Success if computation was succesful,
 
324
      *          \c NumericalIssue if the matrix.appears to be negative.
 
325
      */
 
326
    ComputationInfo info() const
 
327
    {
 
328
      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
 
329
      return m_info;
 
330
    }
 
331
 
 
332
    /** Computes the sparse Cholesky decomposition of \a matrix */
 
333
    void compute(const MatrixType& matrix)
 
334
    {
 
335
      derived().analyzePattern(matrix);
 
336
      derived().factorize(matrix);
 
337
    }
 
338
    
 
339
    /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
 
340
      *
 
341
      * \sa compute()
 
342
      */
 
343
    template<typename Rhs>
 
344
    inline const internal::solve_retval<SuperLUBase, Rhs> solve(const MatrixBase<Rhs>& b) const
 
345
    {
 
346
      eigen_assert(m_isInitialized && "SuperLU is not initialized.");
 
347
      eigen_assert(rows()==b.rows()
 
348
                && "SuperLU::solve(): invalid number of rows of the right hand side matrix b");
 
349
      return internal::solve_retval<SuperLUBase, Rhs>(*this, b.derived());
 
350
    }
 
351
    
 
352
    /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
 
353
      *
 
354
      * \sa compute()
 
355
      */
 
356
//     template<typename Rhs>
 
357
//     inline const internal::sparse_solve_retval<SuperLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const
 
358
//     {
 
359
//       eigen_assert(m_isInitialized && "SuperLU is not initialized.");
 
360
//       eigen_assert(rows()==b.rows()
 
361
//                 && "SuperLU::solve(): invalid number of rows of the right hand side matrix b");
 
362
//       return internal::sparse_solve_retval<SuperLU, Rhs>(*this, b.derived());
 
363
//     }
 
364
    
 
365
    /** Performs a symbolic decomposition on the sparcity of \a matrix.
 
366
      *
 
367
      * This function is particularly useful when solving for several problems having the same structure.
 
368
      * 
 
369
      * \sa factorize()
 
370
      */
 
371
    void analyzePattern(const MatrixType& /*matrix*/)
 
372
    {
 
373
      m_isInitialized = true;
 
374
      m_info = Success;
 
375
      m_analysisIsOk = true;
 
376
      m_factorizationIsOk = false;
 
377
    }
 
378
    
 
379
    template<typename Stream>
 
380
    void dumpMemory(Stream& s)
 
381
    {}
 
382
    
 
383
  protected:
 
384
    
 
385
    void initFactorization(const MatrixType& a)
 
386
    {
 
387
      set_default_options(&this->m_sluOptions);
 
388
      
 
389
      const int size = a.rows();
 
390
      m_matrix = a;
 
391
 
 
392
      m_sluA = internal::asSluMatrix(m_matrix);
 
393
      clearFactors();
 
394
 
 
395
      m_p.resize(size);
 
396
      m_q.resize(size);
 
397
      m_sluRscale.resize(size);
 
398
      m_sluCscale.resize(size);
 
399
      m_sluEtree.resize(size);
 
400
 
 
401
      // set empty B and X
 
402
      m_sluB.setStorageType(SLU_DN);
 
403
      m_sluB.setScalarType<Scalar>();
 
404
      m_sluB.Mtype          = SLU_GE;
 
405
      m_sluB.storage.values = 0;
 
406
      m_sluB.nrow           = 0;
 
407
      m_sluB.ncol           = 0;
 
408
      m_sluB.storage.lda    = size;
 
409
      m_sluX                = m_sluB;
 
410
      
 
411
      m_extractedDataAreDirty = true;
 
412
    }
 
413
    
 
414
    void init()
 
415
    {
 
416
      m_info = InvalidInput;
 
417
      m_isInitialized = false;
 
418
      m_sluL.Store = 0;
 
419
      m_sluU.Store = 0;
 
420
    }
 
421
    
 
422
    void extractData() const;
 
423
 
 
424
    void clearFactors()
 
425
    {
 
426
      if(m_sluL.Store)
 
427
        Destroy_SuperNode_Matrix(&m_sluL);
 
428
      if(m_sluU.Store)
 
429
        Destroy_CompCol_Matrix(&m_sluU);
 
430
 
 
431
      m_sluL.Store = 0;
 
432
      m_sluU.Store = 0;
 
433
 
 
434
      memset(&m_sluL,0,sizeof m_sluL);
 
435
      memset(&m_sluU,0,sizeof m_sluU);
 
436
    }
 
437
 
 
438
    // cached data to reduce reallocation, etc.
 
439
    mutable LUMatrixType m_l;
 
440
    mutable LUMatrixType m_u;
 
441
    mutable IntColVectorType m_p;
 
442
    mutable IntRowVectorType m_q;
 
443
 
 
444
    mutable LUMatrixType m_matrix;  // copy of the factorized matrix
 
445
    mutable SluMatrix m_sluA;
 
446
    mutable SuperMatrix m_sluL, m_sluU;
 
447
    mutable SluMatrix m_sluB, m_sluX;
 
448
    mutable SuperLUStat_t m_sluStat;
 
449
    mutable superlu_options_t m_sluOptions;
 
450
    mutable std::vector<int> m_sluEtree;
 
451
    mutable Matrix<RealScalar,Dynamic,1> m_sluRscale, m_sluCscale;
 
452
    mutable Matrix<RealScalar,Dynamic,1> m_sluFerr, m_sluBerr;
 
453
    mutable char m_sluEqued;
 
454
 
 
455
    mutable ComputationInfo m_info;
 
456
    bool m_isInitialized;
 
457
    int m_factorizationIsOk;
 
458
    int m_analysisIsOk;
 
459
    mutable bool m_extractedDataAreDirty;
 
460
    
 
461
  private:
 
462
    SuperLUBase(SuperLUBase& ) { }
 
463
};
 
464
 
 
465
 
 
466
/** \ingroup SuperLUSupport_Module
 
467
  * \class SuperLU
 
468
  * \brief A sparse direct LU factorization and solver based on the SuperLU library
 
469
  *
 
470
  * This class allows to solve for A.X = B sparse linear problems via a direct LU factorization
 
471
  * using the SuperLU library. The sparse matrix A must be squared and invertible. The vectors or matrices
 
472
  * X and B can be either dense or sparse.
 
473
  *
 
474
  * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
 
475
  *
 
476
  * \sa \ref TutorialSparseDirectSolvers
 
477
  */
 
478
template<typename _MatrixType>
 
479
class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> >
 
480
{
 
481
  public:
 
482
    typedef SuperLUBase<_MatrixType,SuperLU> Base;
 
483
    typedef _MatrixType MatrixType;
 
484
    typedef typename Base::Scalar Scalar;
 
485
    typedef typename Base::RealScalar RealScalar;
 
486
    typedef typename Base::Index Index;
 
487
    typedef typename Base::IntRowVectorType IntRowVectorType;
 
488
    typedef typename Base::IntColVectorType IntColVectorType;    
 
489
    typedef typename Base::LUMatrixType LUMatrixType;
 
490
    typedef TriangularView<LUMatrixType, Lower|UnitDiag>  LMatrixType;
 
491
    typedef TriangularView<LUMatrixType,  Upper>           UMatrixType;
 
492
 
 
493
  public:
 
494
 
 
495
    SuperLU() : Base() { init(); }
 
496
 
 
497
    SuperLU(const MatrixType& matrix) : Base()
 
498
    {
 
499
      Base::init();
 
500
      compute(matrix);
 
501
    }
 
502
 
 
503
    ~SuperLU()
 
504
    {
 
505
    }
 
506
    
 
507
    /** Performs a symbolic decomposition on the sparcity of \a matrix.
 
508
      *
 
509
      * This function is particularly useful when solving for several problems having the same structure.
 
510
      * 
 
511
      * \sa factorize()
 
512
      */
 
513
    void analyzePattern(const MatrixType& matrix)
 
514
    {
 
515
      m_info = InvalidInput;
 
516
      m_isInitialized = false;
 
517
      Base::analyzePattern(matrix);
 
518
    }
 
519
    
 
520
    /** Performs a numeric decomposition of \a matrix
 
521
      *
 
522
      * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
 
523
      *
 
524
      * \sa analyzePattern()
 
525
      */
 
526
    void factorize(const MatrixType& matrix);
 
527
    
 
528
    #ifndef EIGEN_PARSED_BY_DOXYGEN
 
529
    /** \internal */
 
530
    template<typename Rhs,typename Dest>
 
531
    void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
 
532
    #endif // EIGEN_PARSED_BY_DOXYGEN
 
533
    
 
534
    inline const LMatrixType& matrixL() const
 
535
    {
 
536
      if (m_extractedDataAreDirty) this->extractData();
 
537
      return m_l;
 
538
    }
 
539
 
 
540
    inline const UMatrixType& matrixU() const
 
541
    {
 
542
      if (m_extractedDataAreDirty) this->extractData();
 
543
      return m_u;
 
544
    }
 
545
 
 
546
    inline const IntColVectorType& permutationP() const
 
547
    {
 
548
      if (m_extractedDataAreDirty) this->extractData();
 
549
      return m_p;
 
550
    }
 
551
 
 
552
    inline const IntRowVectorType& permutationQ() const
 
553
    {
 
554
      if (m_extractedDataAreDirty) this->extractData();
 
555
      return m_q;
 
556
    }
 
557
    
 
558
    Scalar determinant() const;
 
559
    
 
560
  protected:
 
561
    
 
562
    using Base::m_matrix;
 
563
    using Base::m_sluOptions;
 
564
    using Base::m_sluA;
 
565
    using Base::m_sluB;
 
566
    using Base::m_sluX;
 
567
    using Base::m_p;
 
568
    using Base::m_q;
 
569
    using Base::m_sluEtree;
 
570
    using Base::m_sluEqued;
 
571
    using Base::m_sluRscale;
 
572
    using Base::m_sluCscale;
 
573
    using Base::m_sluL;
 
574
    using Base::m_sluU;
 
575
    using Base::m_sluStat;
 
576
    using Base::m_sluFerr;
 
577
    using Base::m_sluBerr;
 
578
    using Base::m_l;
 
579
    using Base::m_u;
 
580
    
 
581
    using Base::m_analysisIsOk;
 
582
    using Base::m_factorizationIsOk;
 
583
    using Base::m_extractedDataAreDirty;
 
584
    using Base::m_isInitialized;
 
585
    using Base::m_info;
 
586
    
 
587
    void init()
 
588
    {
 
589
      Base::init();
 
590
      
 
591
      set_default_options(&this->m_sluOptions);
 
592
      m_sluOptions.PrintStat        = NO;
 
593
      m_sluOptions.ConditionNumber  = NO;
 
594
      m_sluOptions.Trans            = NOTRANS;
 
595
      m_sluOptions.ColPerm          = COLAMD;
 
596
    }
 
597
    
 
598
    
 
599
  private:
 
600
    SuperLU(SuperLU& ) { }
 
601
};
 
602
 
 
603
template<typename MatrixType>
 
604
void SuperLU<MatrixType>::factorize(const MatrixType& a)
 
605
{
 
606
  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
 
607
  if(!m_analysisIsOk)
 
608
  {
 
609
    m_info = InvalidInput;
 
610
    return;
 
611
  }
 
612
  
 
613
  this->initFactorization(a);
 
614
  
 
615
  int info = 0;
 
616
  RealScalar recip_pivot_growth, rcond;
 
617
  RealScalar ferr, berr;
 
618
 
 
619
  StatInit(&m_sluStat);
 
620
  SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
 
621
                &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
 
622
                &m_sluL, &m_sluU,
 
623
                NULL, 0,
 
624
                &m_sluB, &m_sluX,
 
625
                &recip_pivot_growth, &rcond,
 
626
                &ferr, &berr,
 
627
                &m_sluStat, &info, Scalar());
 
628
  StatFree(&m_sluStat);
 
629
 
 
630
  m_extractedDataAreDirty = true;
 
631
 
 
632
  // FIXME how to better check for errors ???
 
633
  m_info = info == 0 ? Success : NumericalIssue;
 
634
  m_factorizationIsOk = true;
 
635
}
 
636
 
 
637
template<typename MatrixType>
 
638
template<typename Rhs,typename Dest>
 
639
void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
 
640
{
 
641
  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
 
642
 
 
643
  const int size = m_matrix.rows();
 
644
  const int rhsCols = b.cols();
 
645
  eigen_assert(size==b.rows());
 
646
 
 
647
  m_sluOptions.Trans = NOTRANS;
 
648
  m_sluOptions.Fact = FACTORED;
 
649
  m_sluOptions.IterRefine = NOREFINE;
 
650
  
 
651
 
 
652
  m_sluFerr.resize(rhsCols);
 
653
  m_sluBerr.resize(rhsCols);
 
654
  m_sluB = SluMatrix::Map(b.const_cast_derived());
 
655
  m_sluX = SluMatrix::Map(x.derived());
 
656
  
 
657
  typename Rhs::PlainObject b_cpy;
 
658
  if(m_sluEqued!='N')
 
659
  {
 
660
    b_cpy = b;
 
661
    m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());  
 
662
  }
 
663
 
 
664
  StatInit(&m_sluStat);
 
665
  int info = 0;
 
666
  RealScalar recip_pivot_growth, rcond;
 
667
  SuperLU_gssvx(&m_sluOptions, &m_sluA,
 
668
                m_q.data(), m_p.data(),
 
669
                &m_sluEtree[0], &m_sluEqued,
 
670
                &m_sluRscale[0], &m_sluCscale[0],
 
671
                &m_sluL, &m_sluU,
 
672
                NULL, 0,
 
673
                &m_sluB, &m_sluX,
 
674
                &recip_pivot_growth, &rcond,
 
675
                &m_sluFerr[0], &m_sluBerr[0],
 
676
                &m_sluStat, &info, Scalar());
 
677
  StatFree(&m_sluStat);
 
678
  m_info = info==0 ? Success : NumericalIssue;
 
679
}
 
680
 
 
681
// the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
 
682
//
 
683
//  Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
 
684
//
 
685
//  THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
 
686
//  EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
 
687
//
 
688
template<typename MatrixType, typename Derived>
 
689
void SuperLUBase<MatrixType,Derived>::extractData() const
 
690
{
 
691
  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
 
692
  if (m_extractedDataAreDirty)
 
693
  {
 
694
    int         upper;
 
695
    int         fsupc, istart, nsupr;
 
696
    int         lastl = 0, lastu = 0;
 
697
    SCformat    *Lstore = static_cast<SCformat*>(m_sluL.Store);
 
698
    NCformat    *Ustore = static_cast<NCformat*>(m_sluU.Store);
 
699
    Scalar      *SNptr;
 
700
 
 
701
    const int size = m_matrix.rows();
 
702
    m_l.resize(size,size);
 
703
    m_l.resizeNonZeros(Lstore->nnz);
 
704
    m_u.resize(size,size);
 
705
    m_u.resizeNonZeros(Ustore->nnz);
 
706
 
 
707
    int* Lcol = m_l.outerIndexPtr();
 
708
    int* Lrow = m_l.innerIndexPtr();
 
709
    Scalar* Lval = m_l.valuePtr();
 
710
 
 
711
    int* Ucol = m_u.outerIndexPtr();
 
712
    int* Urow = m_u.innerIndexPtr();
 
713
    Scalar* Uval = m_u.valuePtr();
 
714
 
 
715
    Ucol[0] = 0;
 
716
    Ucol[0] = 0;
 
717
 
 
718
    /* for each supernode */
 
719
    for (int k = 0; k <= Lstore->nsuper; ++k)
 
720
    {
 
721
      fsupc   = L_FST_SUPC(k);
 
722
      istart  = L_SUB_START(fsupc);
 
723
      nsupr   = L_SUB_START(fsupc+1) - istart;
 
724
      upper   = 1;
 
725
 
 
726
      /* for each column in the supernode */
 
727
      for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
 
728
      {
 
729
        SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
 
730
 
 
731
        /* Extract U */
 
732
        for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
 
733
        {
 
734
          Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
 
735
          /* Matlab doesn't like explicit zero. */
 
736
          if (Uval[lastu] != 0.0)
 
737
            Urow[lastu++] = U_SUB(i);
 
738
        }
 
739
        for (int i = 0; i < upper; ++i)
 
740
        {
 
741
          /* upper triangle in the supernode */
 
742
          Uval[lastu] = SNptr[i];
 
743
          /* Matlab doesn't like explicit zero. */
 
744
          if (Uval[lastu] != 0.0)
 
745
            Urow[lastu++] = L_SUB(istart+i);
 
746
        }
 
747
        Ucol[j+1] = lastu;
 
748
 
 
749
        /* Extract L */
 
750
        Lval[lastl] = 1.0; /* unit diagonal */
 
751
        Lrow[lastl++] = L_SUB(istart + upper - 1);
 
752
        for (int i = upper; i < nsupr; ++i)
 
753
        {
 
754
          Lval[lastl] = SNptr[i];
 
755
          /* Matlab doesn't like explicit zero. */
 
756
          if (Lval[lastl] != 0.0)
 
757
            Lrow[lastl++] = L_SUB(istart+i);
 
758
        }
 
759
        Lcol[j+1] = lastl;
 
760
 
 
761
        ++upper;
 
762
      } /* for j ... */
 
763
 
 
764
    } /* for k ... */
 
765
 
 
766
    // squeeze the matrices :
 
767
    m_l.resizeNonZeros(lastl);
 
768
    m_u.resizeNonZeros(lastu);
 
769
 
 
770
    m_extractedDataAreDirty = false;
 
771
  }
 
772
}
 
773
 
 
774
template<typename MatrixType>
 
775
typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const
 
776
{
 
777
  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for computing the determinant, you must first call either compute() or analyzePattern()/factorize()");
 
778
  
 
779
  if (m_extractedDataAreDirty)
 
780
    this->extractData();
 
781
 
 
782
  Scalar det = Scalar(1);
 
783
  for (int j=0; j<m_u.cols(); ++j)
 
784
  {
 
785
    if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0)
 
786
    {
 
787
      int lastId = m_u.outerIndexPtr()[j+1]-1;
 
788
      eigen_assert(m_u.innerIndexPtr()[lastId]<=j);
 
789
      if (m_u.innerIndexPtr()[lastId]==j)
 
790
        det *= m_u.valuePtr()[lastId];
 
791
    }
 
792
  }
 
793
  if(m_sluEqued!='N')
 
794
    return det/m_sluRscale.prod()/m_sluCscale.prod();
 
795
  else
 
796
    return det;
 
797
}
 
798
 
 
799
#ifdef EIGEN_PARSED_BY_DOXYGEN
 
800
#define EIGEN_SUPERLU_HAS_ILU
 
801
#endif
 
802
 
 
803
#ifdef EIGEN_SUPERLU_HAS_ILU
 
804
 
 
805
/** \ingroup SuperLUSupport_Module
 
806
  * \class SuperILU
 
807
  * \brief A sparse direct \b incomplete LU factorization and solver based on the SuperLU library
 
808
  *
 
809
  * This class allows to solve for an approximate solution of A.X = B sparse linear problems via an incomplete LU factorization
 
810
  * using the SuperLU library. This class is aimed to be used as a preconditioner of the iterative linear solvers.
 
811
  *
 
812
  * \warning This class requires SuperLU 4 or later.
 
813
  *
 
814
  * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
 
815
  *
 
816
  * \sa \ref TutorialSparseDirectSolvers, class ConjugateGradient, class BiCGSTAB
 
817
  */
 
818
 
 
819
template<typename _MatrixType>
 
820
class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> >
 
821
{
 
822
  public:
 
823
    typedef SuperLUBase<_MatrixType,SuperILU> Base;
 
824
    typedef _MatrixType MatrixType;
 
825
    typedef typename Base::Scalar Scalar;
 
826
    typedef typename Base::RealScalar RealScalar;
 
827
    typedef typename Base::Index Index;
 
828
 
 
829
  public:
 
830
 
 
831
    SuperILU() : Base() { init(); }
 
832
 
 
833
    SuperILU(const MatrixType& matrix) : Base()
 
834
    {
 
835
      init();
 
836
      compute(matrix);
 
837
    }
 
838
 
 
839
    ~SuperILU()
 
840
    {
 
841
    }
 
842
    
 
843
    /** Performs a symbolic decomposition on the sparcity of \a matrix.
 
844
      *
 
845
      * This function is particularly useful when solving for several problems having the same structure.
 
846
      * 
 
847
      * \sa factorize()
 
848
      */
 
849
    void analyzePattern(const MatrixType& matrix)
 
850
    {
 
851
      Base::analyzePattern(matrix);
 
852
    }
 
853
    
 
854
    /** Performs a numeric decomposition of \a matrix
 
855
      *
 
856
      * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
 
857
      *
 
858
      * \sa analyzePattern()
 
859
      */
 
860
    void factorize(const MatrixType& matrix);
 
861
    
 
862
    #ifndef EIGEN_PARSED_BY_DOXYGEN
 
863
    /** \internal */
 
864
    template<typename Rhs,typename Dest>
 
865
    void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
 
866
    #endif // EIGEN_PARSED_BY_DOXYGEN
 
867
    
 
868
  protected:
 
869
    
 
870
    using Base::m_matrix;
 
871
    using Base::m_sluOptions;
 
872
    using Base::m_sluA;
 
873
    using Base::m_sluB;
 
874
    using Base::m_sluX;
 
875
    using Base::m_p;
 
876
    using Base::m_q;
 
877
    using Base::m_sluEtree;
 
878
    using Base::m_sluEqued;
 
879
    using Base::m_sluRscale;
 
880
    using Base::m_sluCscale;
 
881
    using Base::m_sluL;
 
882
    using Base::m_sluU;
 
883
    using Base::m_sluStat;
 
884
    using Base::m_sluFerr;
 
885
    using Base::m_sluBerr;
 
886
    using Base::m_l;
 
887
    using Base::m_u;
 
888
    
 
889
    using Base::m_analysisIsOk;
 
890
    using Base::m_factorizationIsOk;
 
891
    using Base::m_extractedDataAreDirty;
 
892
    using Base::m_isInitialized;
 
893
    using Base::m_info;
 
894
 
 
895
    void init()
 
896
    {
 
897
      Base::init();
 
898
      
 
899
      ilu_set_default_options(&m_sluOptions);
 
900
      m_sluOptions.PrintStat        = NO;
 
901
      m_sluOptions.ConditionNumber  = NO;
 
902
      m_sluOptions.Trans            = NOTRANS;
 
903
      m_sluOptions.ColPerm          = MMD_AT_PLUS_A;
 
904
      
 
905
      // no attempt to preserve column sum
 
906
      m_sluOptions.ILU_MILU = SILU;
 
907
      // only basic ILU(k) support -- no direct control over memory consumption
 
908
      // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
 
909
      // and set ILU_FillFactor to max memory growth
 
910
      m_sluOptions.ILU_DropRule = DROP_BASIC;
 
911
      m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision()*10;
 
912
    }
 
913
    
 
914
  private:
 
915
    SuperILU(SuperILU& ) { }
 
916
};
 
917
 
 
918
template<typename MatrixType>
 
919
void SuperILU<MatrixType>::factorize(const MatrixType& a)
 
920
{
 
921
  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
 
922
  if(!m_analysisIsOk)
 
923
  {
 
924
    m_info = InvalidInput;
 
925
    return;
 
926
  }
 
927
  
 
928
  this->initFactorization(a);
 
929
 
 
930
  int info = 0;
 
931
  RealScalar recip_pivot_growth, rcond;
 
932
 
 
933
  StatInit(&m_sluStat);
 
934
  SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
 
935
                &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
 
936
                &m_sluL, &m_sluU,
 
937
                NULL, 0,
 
938
                &m_sluB, &m_sluX,
 
939
                &recip_pivot_growth, &rcond,
 
940
                &m_sluStat, &info, Scalar());
 
941
  StatFree(&m_sluStat);
 
942
 
 
943
  // FIXME how to better check for errors ???
 
944
  m_info = info == 0 ? Success : NumericalIssue;
 
945
  m_factorizationIsOk = true;
 
946
}
 
947
 
 
948
template<typename MatrixType>
 
949
template<typename Rhs,typename Dest>
 
950
void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
 
951
{
 
952
  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
 
953
 
 
954
  const int size = m_matrix.rows();
 
955
  const int rhsCols = b.cols();
 
956
  eigen_assert(size==b.rows());
 
957
 
 
958
  m_sluOptions.Trans = NOTRANS;
 
959
  m_sluOptions.Fact = FACTORED;
 
960
  m_sluOptions.IterRefine = NOREFINE;
 
961
 
 
962
  m_sluFerr.resize(rhsCols);
 
963
  m_sluBerr.resize(rhsCols);
 
964
  m_sluB = SluMatrix::Map(b.const_cast_derived());
 
965
  m_sluX = SluMatrix::Map(x.derived());
 
966
 
 
967
  typename Rhs::PlainObject b_cpy;
 
968
  if(m_sluEqued!='N')
 
969
  {
 
970
    b_cpy = b;
 
971
    m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());  
 
972
  }
 
973
  
 
974
  int info = 0;
 
975
  RealScalar recip_pivot_growth, rcond;
 
976
 
 
977
  StatInit(&m_sluStat);
 
978
  SuperLU_gsisx(&m_sluOptions, &m_sluA,
 
979
                m_q.data(), m_p.data(),
 
980
                &m_sluEtree[0], &m_sluEqued,
 
981
                &m_sluRscale[0], &m_sluCscale[0],
 
982
                &m_sluL, &m_sluU,
 
983
                NULL, 0,
 
984
                &m_sluB, &m_sluX,
 
985
                &recip_pivot_growth, &rcond,
 
986
                &m_sluStat, &info, Scalar());
 
987
  StatFree(&m_sluStat);
 
988
 
 
989
  m_info = info==0 ? Success : NumericalIssue;
 
990
}
 
991
#endif
 
992
 
 
993
namespace internal {
 
994
  
 
995
template<typename _MatrixType, typename Derived, typename Rhs>
 
996
struct solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
 
997
  : solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
 
998
{
 
999
  typedef SuperLUBase<_MatrixType,Derived> Dec;
 
1000
  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
 
1001
 
 
1002
  template<typename Dest> void evalTo(Dest& dst) const
 
1003
  {
 
1004
    dec().derived()._solve(rhs(),dst);
 
1005
  }
 
1006
};
 
1007
 
 
1008
template<typename _MatrixType, typename Derived, typename Rhs>
 
1009
struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
 
1010
  : sparse_solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
 
1011
{
 
1012
  typedef SuperLUBase<_MatrixType,Derived> Dec;
 
1013
  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
 
1014
 
 
1015
  template<typename Dest> void evalTo(Dest& dst) const
 
1016
  {
 
1017
    dec().derived()._solve(rhs(),dst);
 
1018
  }
 
1019
};
 
1020
 
 
1021
} // end namespace internal
 
1022
 
 
1023
} // end namespace Eigen
 
1024
 
 
1025
#endif // EIGEN_SUPERLUSUPPORT_H