← Back to branch summary

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

« back to all changes in this revision

Viewing changes to extern/Eigen3/Eigen/src/PardisoSupport/PardisoSupport.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/PardisoSupport/PardisoSupport.h
 
1
/*
 
2
 Copyright (c) 2011, Intel Corporation. All rights reserved.
 
3
 
 
4
 Redistribution and use in source and binary forms, with or without modification,
 
5
 are permitted provided that the following conditions are met:
 
6
 
 
7
 * Redistributions of source code must retain the above copyright notice, this
 
8
   list of conditions and the following disclaimer.
 
9
 * Redistributions in binary form must reproduce the above copyright notice,
 
10
   this list of conditions and the following disclaimer in the documentation
 
11
   and/or other materials provided with the distribution.
 
12
 * Neither the name of Intel Corporation nor the names of its contributors may
 
13
   be used to endorse or promote products derived from this software without
 
14
   specific prior written permission.
 
15
 
 
16
 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 
17
 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 
18
 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 
19
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 
20
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 
21
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 
22
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 
23
 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 
24
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 
25
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
26
 
 
27
 ********************************************************************************
 
28
 *   Content : Eigen bindings to Intel(R) MKL PARDISO
 
29
 ********************************************************************************
 
30
*/
 
31
 
 
32
#ifndef EIGEN_PARDISOSUPPORT_H
 
33
#define EIGEN_PARDISOSUPPORT_H
 
34
 
 
35
namespace Eigen { 
 
36
 
 
37
template<typename _MatrixType> class PardisoLU;
 
38
template<typename _MatrixType, int Options=Upper> class PardisoLLT;
 
39
template<typename _MatrixType, int Options=Upper> class PardisoLDLT;
 
40
 
 
41
namespace internal
 
42
{
 
43
  template<typename Index>
 
44
  struct pardiso_run_selector
 
45
  {
 
46
    static Index run( _MKL_DSS_HANDLE_t pt, Index maxfct, Index mnum, Index type, Index phase, Index n, void *a,
 
47
                      Index *ia, Index *ja, Index *perm, Index nrhs, Index *iparm, Index msglvl, void *b, void *x)
 
48
    {
 
49
      Index error = 0;
 
50
      ::pardiso(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
 
51
      return error;
 
52
    }
 
53
  };
 
54
  template<>
 
55
  struct pardiso_run_selector<long long int>
 
56
  {
 
57
    typedef long long int Index;
 
58
    static Index run( _MKL_DSS_HANDLE_t pt, Index maxfct, Index mnum, Index type, Index phase, Index n, void *a,
 
59
                      Index *ia, Index *ja, Index *perm, Index nrhs, Index *iparm, Index msglvl, void *b, void *x)
 
60
    {
 
61
      Index error = 0;
 
62
      ::pardiso_64(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
 
63
      return error;
 
64
    }
 
65
  };
 
66
 
 
67
  template<class Pardiso> struct pardiso_traits;
 
68
 
 
69
  template<typename _MatrixType>
 
70
  struct pardiso_traits< PardisoLU<_MatrixType> >
 
71
  {
 
72
    typedef _MatrixType MatrixType;
 
73
    typedef typename _MatrixType::Scalar Scalar;
 
74
    typedef typename _MatrixType::RealScalar RealScalar;
 
75
    typedef typename _MatrixType::Index Index;
 
76
  };
 
77
 
 
78
  template<typename _MatrixType, int Options>
 
79
  struct pardiso_traits< PardisoLLT<_MatrixType, Options> >
 
80
  {
 
81
    typedef _MatrixType MatrixType;
 
82
    typedef typename _MatrixType::Scalar Scalar;
 
83
    typedef typename _MatrixType::RealScalar RealScalar;
 
84
    typedef typename _MatrixType::Index Index;
 
85
  };
 
86
 
 
87
  template<typename _MatrixType, int Options>
 
88
  struct pardiso_traits< PardisoLDLT<_MatrixType, Options> >
 
89
  {
 
90
    typedef _MatrixType MatrixType;
 
91
    typedef typename _MatrixType::Scalar Scalar;
 
92
    typedef typename _MatrixType::RealScalar RealScalar;
 
93
    typedef typename _MatrixType::Index Index;    
 
94
  };
 
95
 
 
96
}
 
97
 
 
98
template<class Derived>
 
99
class PardisoImpl
 
100
{
 
101
    typedef internal::pardiso_traits<Derived> Traits;
 
102
  public:
 
103
    typedef typename Traits::MatrixType MatrixType;
 
104
    typedef typename Traits::Scalar Scalar;
 
105
    typedef typename Traits::RealScalar RealScalar;
 
106
    typedef typename Traits::Index Index;
 
107
    typedef SparseMatrix<Scalar,RowMajor,Index> SparseMatrixType;
 
108
    typedef Matrix<Scalar,Dynamic,1> VectorType;
 
109
    typedef Matrix<Index, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
 
110
    typedef Matrix<Index, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
 
111
    enum {
 
112
      ScalarIsComplex = NumTraits<Scalar>::IsComplex
 
113
    };
 
114
 
 
115
    PardisoImpl()
 
116
    {
 
117
      eigen_assert((sizeof(Index) >= sizeof(_INTEGER_t) && sizeof(Index) <= 8) && "Non-supported index type");
 
118
      m_iparm.setZero();
 
119
      m_msglvl = 0; // No output
 
120
      m_initialized = false;
 
121
    }
 
122
 
 
123
    ~PardisoImpl()
 
124
    {
 
125
      pardisoRelease();
 
126
    }
 
127
 
 
128
    inline Index cols() const { return m_size; }
 
129
    inline Index rows() const { return m_size; }
 
130
  
 
131
    /** \brief Reports whether previous computation was successful.
 
132
      *
 
133
      * \returns \c Success if computation was succesful,
 
134
      *          \c NumericalIssue if the matrix appears to be negative.
 
135
      */
 
136
    ComputationInfo info() const
 
137
    {
 
138
      eigen_assert(m_initialized && "Decomposition is not initialized.");
 
139
      return m_info;
 
140
    }
 
141
 
 
142
    /** \warning for advanced usage only.
 
143
      * \returns a reference to the parameter array controlling PARDISO.
 
144
      * See the PARDISO manual to know how to use it. */
 
145
    Array<Index,64,1>& pardisoParameterArray()
 
146
    {
 
147
      return m_iparm;
 
148
    }
 
149
    
 
150
    /** Performs a symbolic decomposition on the sparcity of \a matrix.
 
151
      *
 
152
      * This function is particularly useful when solving for several problems having the same structure.
 
153
      * 
 
154
      * \sa factorize()
 
155
      */
 
156
    Derived& analyzePattern(const MatrixType& matrix);
 
157
    
 
158
    /** Performs a numeric decomposition of \a matrix
 
159
      *
 
160
      * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
 
161
      *
 
162
      * \sa analyzePattern()
 
163
      */
 
164
    Derived& factorize(const MatrixType& matrix);
 
165
 
 
166
    Derived& compute(const MatrixType& matrix);
 
167
    
 
168
    /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
 
169
      *
 
170
      * \sa compute()
 
171
      */
 
172
    template<typename Rhs>
 
173
    inline const internal::solve_retval<PardisoImpl, Rhs>
 
174
    solve(const MatrixBase<Rhs>& b) const
 
175
    {
 
176
      eigen_assert(m_initialized && "Pardiso solver is not initialized.");
 
177
      eigen_assert(rows()==b.rows()
 
178
                && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
 
179
      return internal::solve_retval<PardisoImpl, Rhs>(*this, b.derived());
 
180
    }
 
181
 
 
182
    /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
 
183
      *
 
184
      * \sa compute()
 
185
      */
 
186
    template<typename Rhs>
 
187
    inline const internal::sparse_solve_retval<PardisoImpl, Rhs>
 
188
    solve(const SparseMatrixBase<Rhs>& b) const
 
189
    {
 
190
      eigen_assert(m_initialized && "Pardiso solver is not initialized.");
 
191
      eigen_assert(rows()==b.rows()
 
192
                && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
 
193
      return internal::sparse_solve_retval<PardisoImpl, Rhs>(*this, b.derived());
 
194
    }
 
195
 
 
196
    Derived& derived()
 
197
    {
 
198
      return *static_cast<Derived*>(this);
 
199
    }
 
200
    const Derived& derived() const
 
201
    {
 
202
      return *static_cast<const Derived*>(this);
 
203
    }
 
204
 
 
205
    template<typename BDerived, typename XDerived>
 
206
    bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const;
 
207
 
 
208
    /** \internal */
 
209
    template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
 
210
    void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
 
211
    {
 
212
      eigen_assert(m_size==b.rows());
 
213
 
 
214
      // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
 
215
      static const int NbColsAtOnce = 4;
 
216
      int rhsCols = b.cols();
 
217
      int size = b.rows();
 
218
      // Pardiso cannot solve in-place,
 
219
      // so we need two temporaries
 
220
      Eigen::Matrix<DestScalar,Dynamic,Dynamic,ColMajor> tmp_rhs(size,rhsCols);
 
221
      Eigen::Matrix<DestScalar,Dynamic,Dynamic,ColMajor> tmp_res(size,rhsCols);
 
222
      for(int k=0; k<rhsCols; k+=NbColsAtOnce)
 
223
      {
 
224
        int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
 
225
        tmp_rhs.leftCols(actualCols) = b.middleCols(k,actualCols);
 
226
        tmp_res.leftCols(actualCols) = derived().solve(tmp_rhs.leftCols(actualCols));
 
227
        dest.middleCols(k,actualCols) = tmp_res.leftCols(actualCols).sparseView();
 
228
      }
 
229
    }
 
230
 
 
231
  protected:
 
232
    void pardisoRelease()
 
233
    {
 
234
      if(m_initialized) // Factorization ran at least once
 
235
      {
 
236
        internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, -1, m_size, 0, 0, 0, m_perm.data(), 0,
 
237
                                                   m_iparm.data(), m_msglvl, 0, 0);
 
238
      }
 
239
    }
 
240
 
 
241
    void pardisoInit(int type)
 
242
    {
 
243
      m_type = type;
 
244
      bool symmetric = abs(m_type) < 10;
 
245
      m_iparm[0] = 1;   // No solver default
 
246
      m_iparm[1] = 3;   // use Metis for the ordering
 
247
      m_iparm[2] = 1;   // Numbers of processors, value of OMP_NUM_THREADS
 
248
      m_iparm[3] = 0;   // No iterative-direct algorithm
 
249
      m_iparm[4] = 0;   // No user fill-in reducing permutation
 
250
      m_iparm[5] = 0;   // Write solution into x
 
251
      m_iparm[6] = 0;   // Not in use
 
252
      m_iparm[7] = 2;   // Max numbers of iterative refinement steps
 
253
      m_iparm[8] = 0;   // Not in use
 
254
      m_iparm[9] = 13;  // Perturb the pivot elements with 1E-13
 
255
      m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS
 
256
      m_iparm[11] = 0;  // Not in use
 
257
      m_iparm[12] = symmetric ? 0 : 1;  // Maximum weighted matching algorithm is switched-off (default for symmetric).
 
258
                                        // Try m_iparm[12] = 1 in case of inappropriate accuracy
 
259
      m_iparm[13] = 0;  // Output: Number of perturbed pivots
 
260
      m_iparm[14] = 0;  // Not in use
 
261
      m_iparm[15] = 0;  // Not in use
 
262
      m_iparm[16] = 0;  // Not in use
 
263
      m_iparm[17] = -1; // Output: Number of nonzeros in the factor LU
 
264
      m_iparm[18] = -1; // Output: Mflops for LU factorization
 
265
      m_iparm[19] = 0;  // Output: Numbers of CG Iterations
 
266
      
 
267
      m_iparm[20] = 0;  // 1x1 pivoting
 
268
      m_iparm[26] = 0;  // No matrix checker
 
269
      m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0;
 
270
      m_iparm[34] = 1;  // C indexing
 
271
      m_iparm[59] = 1;  // Automatic switch between In-Core and Out-of-Core modes
 
272
    }
 
273
 
 
274
  protected:
 
275
    // cached data to reduce reallocation, etc.
 
276
    
 
277
    void manageErrorCode(Index error)
 
278
    {
 
279
      switch(error)
 
280
      {
 
281
        case 0:
 
282
          m_info = Success;
 
283
          break;
 
284
        case -4:
 
285
        case -7:
 
286
          m_info = NumericalIssue;
 
287
          break;
 
288
        default:
 
289
          m_info = InvalidInput;
 
290
      }
 
291
    }
 
292
 
 
293
    mutable SparseMatrixType m_matrix;
 
294
    ComputationInfo m_info;
 
295
    bool m_initialized, m_analysisIsOk, m_factorizationIsOk;
 
296
    Index m_type, m_msglvl;
 
297
    mutable void *m_pt[64];
 
298
    mutable Array<Index,64,1> m_iparm;
 
299
    mutable IntColVectorType m_perm;
 
300
    Index m_size;
 
301
    
 
302
  private:
 
303
    PardisoImpl(PardisoImpl &) {}
 
304
};
 
305
 
 
306
template<class Derived>
 
307
Derived& PardisoImpl<Derived>::compute(const MatrixType& a)
 
308
{
 
309
  m_size = a.rows();
 
310
  eigen_assert(a.rows() == a.cols());
 
311
 
 
312
  pardisoRelease();
 
313
  memset(m_pt, 0, sizeof(m_pt));
 
314
  m_perm.setZero(m_size);
 
315
  derived().getMatrix(a);
 
316
  
 
317
  Index error;
 
318
  error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 12, m_size,
 
319
                                                     m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
 
320
                                                     m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
 
321
 
 
322
  manageErrorCode(error);
 
323
  m_analysisIsOk = true;
 
324
  m_factorizationIsOk = true;
 
325
  m_initialized = true;
 
326
  return derived();
 
327
}
 
328
 
 
329
template<class Derived>
 
330
Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a)
 
331
{
 
332
  m_size = a.rows();
 
333
  eigen_assert(m_size == a.cols());
 
334
 
 
335
  pardisoRelease();
 
336
  memset(m_pt, 0, sizeof(m_pt));
 
337
  m_perm.setZero(m_size);
 
338
  derived().getMatrix(a);
 
339
  
 
340
  Index error;
 
341
  error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 11, m_size,
 
342
                                                     m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
 
343
                                                     m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
 
344
  
 
345
  manageErrorCode(error);
 
346
  m_analysisIsOk = true;
 
347
  m_factorizationIsOk = false;
 
348
  m_initialized = true;
 
349
  return derived();
 
350
}
 
351
 
 
352
template<class Derived>
 
353
Derived& PardisoImpl<Derived>::factorize(const MatrixType& a)
 
354
{
 
355
  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
 
356
  eigen_assert(m_size == a.rows() && m_size == a.cols());
 
357
  
 
358
  derived().getMatrix(a);
 
359
 
 
360
  Index error;  
 
361
  error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 22, m_size,
 
362
                                                     m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
 
363
                                                     m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
 
364
  
 
365
  manageErrorCode(error);
 
366
  m_factorizationIsOk = true;
 
367
  return derived();
 
368
}
 
369
 
 
370
template<class Base>
 
371
template<typename BDerived,typename XDerived>
 
372
bool PardisoImpl<Base>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const
 
373
{
 
374
  if(m_iparm[0] == 0) // Factorization was not computed
 
375
    return false;
 
376
 
 
377
  //Index n = m_matrix.rows();
 
378
  Index nrhs = Index(b.cols());
 
379
  eigen_assert(m_size==b.rows());
 
380
  eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported");
 
381
  eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported");
 
382
  eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
 
383
 
 
384
 
 
385
//  switch (transposed) {
 
386
//    case SvNoTrans    : m_iparm[11] = 0 ; break;
 
387
//    case SvTranspose  : m_iparm[11] = 2 ; break;
 
388
//    case SvAdjoint    : m_iparm[11] = 1 ; break;
 
389
//    default:
 
390
//      //std::cerr << "Eigen: transposition  option \"" << transposed << "\" not supported by the PARDISO backend\n";
 
391
//      m_iparm[11] = 0;
 
392
//  }
 
393
 
 
394
  Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data());
 
395
  Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp;
 
396
  
 
397
  // Pardiso cannot solve in-place
 
398
  if(rhs_ptr == x.derived().data())
 
399
  {
 
400
    tmp = b;
 
401
    rhs_ptr = tmp.data();
 
402
  }
 
403
  
 
404
  Index error;
 
405
  error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 33, m_size,
 
406
                                                     m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
 
407
                                                     m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
 
408
                                                     rhs_ptr, x.derived().data());
 
409
 
 
410
  return error==0;
 
411
}
 
412
 
 
413
 
 
414
/** \ingroup PardisoSupport_Module
 
415
  * \class PardisoLU
 
416
  * \brief A sparse direct LU factorization and solver based on the PARDISO library
 
417
  *
 
418
  * This class allows to solve for A.X = B sparse linear problems via a direct LU factorization
 
419
  * using the Intel MKL PARDISO library. The sparse matrix A must be squared and invertible.
 
420
  * The vectors or matrices X and B can be either dense or sparse.
 
421
  *
 
422
  * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
 
423
  *
 
424
  * \sa \ref TutorialSparseDirectSolvers
 
425
  */
 
426
template<typename MatrixType>
 
427
class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> >
 
428
{
 
429
  protected:
 
430
    typedef PardisoImpl< PardisoLU<MatrixType> > Base;
 
431
    typedef typename Base::Scalar Scalar;
 
432
    typedef typename Base::RealScalar RealScalar;
 
433
    using Base::pardisoInit;
 
434
    using Base::m_matrix;
 
435
    friend class PardisoImpl< PardisoLU<MatrixType> >;
 
436
 
 
437
  public:
 
438
 
 
439
    using Base::compute;
 
440
    using Base::solve;
 
441
 
 
442
    PardisoLU()
 
443
      : Base()
 
444
    {
 
445
      pardisoInit(Base::ScalarIsComplex ? 13 : 11);
 
446
    }
 
447
 
 
448
    PardisoLU(const MatrixType& matrix)
 
449
      : Base()
 
450
    {
 
451
      pardisoInit(Base::ScalarIsComplex ? 13 : 11);
 
452
      compute(matrix);
 
453
    }
 
454
  protected:
 
455
    void getMatrix(const MatrixType& matrix)
 
456
    {
 
457
      m_matrix = matrix;
 
458
    }
 
459
    
 
460
  private:
 
461
    PardisoLU(PardisoLU& ) {}
 
462
};
 
463
 
 
464
/** \ingroup PardisoSupport_Module
 
465
  * \class PardisoLLT
 
466
  * \brief A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library
 
467
  *
 
468
  * This class allows to solve for A.X = B sparse linear problems via a LL^T Cholesky factorization
 
469
  * using the Intel MKL PARDISO library. The sparse matrix A must be selfajoint and positive definite.
 
470
  * The vectors or matrices X and B can be either dense or sparse.
 
471
  *
 
472
  * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
 
473
  * \tparam UpLo can be any bitwise combination of Upper, Lower. The default is Upper, meaning only the upper triangular part has to be used.
 
474
  *         Upper|Lower can be used to tell both triangular parts can be used as input.
 
475
  *
 
476
  * \sa \ref TutorialSparseDirectSolvers
 
477
  */
 
478
template<typename MatrixType, int _UpLo>
 
479
class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
 
480
{
 
481
  protected:
 
482
    typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
 
483
    typedef typename Base::Scalar Scalar;
 
484
    typedef typename Base::Index Index;
 
485
    typedef typename Base::RealScalar RealScalar;
 
486
    using Base::pardisoInit;
 
487
    using Base::m_matrix;
 
488
    friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
 
489
 
 
490
  public:
 
491
 
 
492
    enum { UpLo = _UpLo };
 
493
    using Base::compute;
 
494
    using Base::solve;
 
495
 
 
496
    PardisoLLT()
 
497
      : Base()
 
498
    {
 
499
      pardisoInit(Base::ScalarIsComplex ? 4 : 2);
 
500
    }
 
501
 
 
502
    PardisoLLT(const MatrixType& matrix)
 
503
      : Base()
 
504
    {
 
505
      pardisoInit(Base::ScalarIsComplex ? 4 : 2);
 
506
      compute(matrix);
 
507
    }
 
508
    
 
509
  protected:
 
510
    
 
511
    void getMatrix(const MatrixType& matrix)
 
512
    {
 
513
      // PARDISO supports only upper, row-major matrices
 
514
      PermutationMatrix<Dynamic,Dynamic,Index> p_null;
 
515
      m_matrix.resize(matrix.rows(), matrix.cols());
 
516
      m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
 
517
    }
 
518
    
 
519
  private:
 
520
    PardisoLLT(PardisoLLT& ) {}
 
521
};
 
522
 
 
523
/** \ingroup PardisoSupport_Module
 
524
  * \class PardisoLDLT
 
525
  * \brief A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library
 
526
  *
 
527
  * This class allows to solve for A.X = B sparse linear problems via a LDL^T Cholesky factorization
 
528
  * using the Intel MKL PARDISO library. The sparse matrix A is assumed to be selfajoint and positive definite.
 
529
  * For complex matrices, A can also be symmetric only, see the \a Options template parameter.
 
530
  * The vectors or matrices X and B can be either dense or sparse.
 
531
  *
 
532
  * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
 
533
  * \tparam Options can be any bitwise combination of Upper, Lower, and Symmetric. The default is Upper, meaning only the upper triangular part has to be used.
 
534
  *         Symmetric can be used for symmetric, non-selfadjoint complex matrices, the default being to assume a selfadjoint matrix.
 
535
  *         Upper|Lower can be used to tell both triangular parts can be used as input.
 
536
  *
 
537
  * \sa \ref TutorialSparseDirectSolvers
 
538
  */
 
539
template<typename MatrixType, int Options>
 
540
class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> >
 
541
{
 
542
  protected:
 
543
    typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
 
544
    typedef typename Base::Scalar Scalar;
 
545
    typedef typename Base::Index Index;
 
546
    typedef typename Base::RealScalar RealScalar;
 
547
    using Base::pardisoInit;
 
548
    using Base::m_matrix;
 
549
    friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
 
550
 
 
551
  public:
 
552
 
 
553
    using Base::compute;
 
554
    using Base::solve;
 
555
    enum { UpLo = Options&(Upper|Lower) };
 
556
 
 
557
    PardisoLDLT()
 
558
      : Base()
 
559
    {
 
560
      pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
 
561
    }
 
562
 
 
563
    PardisoLDLT(const MatrixType& matrix)
 
564
      : Base()
 
565
    {
 
566
      pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
 
567
      compute(matrix);
 
568
    }
 
569
    
 
570
    void getMatrix(const MatrixType& matrix)
 
571
    {
 
572
      // PARDISO supports only upper, row-major matrices
 
573
      PermutationMatrix<Dynamic,Dynamic,Index> p_null;
 
574
      m_matrix.resize(matrix.rows(), matrix.cols());
 
575
      m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
 
576
    }
 
577
    
 
578
  private:
 
579
    PardisoLDLT(PardisoLDLT& ) {}
 
580
};
 
581
 
 
582
namespace internal {
 
583
  
 
584
template<typename _Derived, typename Rhs>
 
585
struct solve_retval<PardisoImpl<_Derived>, Rhs>
 
586
  : solve_retval_base<PardisoImpl<_Derived>, Rhs>
 
587
{
 
588
  typedef PardisoImpl<_Derived> Dec;
 
589
  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
 
590
 
 
591
  template<typename Dest> void evalTo(Dest& dst) const
 
592
  {
 
593
    dec()._solve(rhs(),dst);
 
594
  }
 
595
};
 
596
 
 
597
template<typename Derived, typename Rhs>
 
598
struct sparse_solve_retval<PardisoImpl<Derived>, Rhs>
 
599
  : sparse_solve_retval_base<PardisoImpl<Derived>, Rhs>
 
600
{
 
601
  typedef PardisoImpl<Derived> Dec;
 
602
  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
 
603
 
 
604
  template<typename Dest> void evalTo(Dest& dst) const
 
605
  {
 
606
    dec().derived()._solve_sparse(rhs(),dst);
 
607
  }
 
608
};
 
609
 
 
610
} // end namespace internal
 
611
 
 
612
} // end namespace Eigen
 
613
 
 
614
#endif // EIGEN_PARDISOSUPPORT_H