Cycles: svn merge -r41225:41232 ^/trunk/blender
[blender.git] / extern / Eigen2 / Eigen / src / Sparse / SuperLUSupport.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra. Eigen itself is part of the KDE project.
3 //
4 // Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24
25 #ifndef EIGEN_SUPERLUSUPPORT_H
26 #define EIGEN_SUPERLUSUPPORT_H
27
28 // declaration of gssvx taken from GMM++
29 #define DECL_GSSVX(NAMESPACE,FNAME,FLOATTYPE,KEYTYPE) \
30     inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A,  \
31          int *perm_c, int *perm_r, int *etree, char *equed,  \
32          FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L,         \
33          SuperMatrix *U, void *work, int lwork,              \
34          SuperMatrix *B, SuperMatrix *X,                     \
35          FLOATTYPE *recip_pivot_growth,                      \
36          FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \
37          SuperLUStat_t *stats, int *info, KEYTYPE) {         \
38     using namespace NAMESPACE; \
39     mem_usage_t mem_usage;                                    \
40     NAMESPACE::FNAME(options, A, perm_c, perm_r, etree, equed, R, C, L,  \
41          U, work, lwork, B, X, recip_pivot_growth, rcond,    \
42          ferr, berr, &mem_usage, stats, info);               \
43     return mem_usage.for_lu; /* bytes used by the factor storage */     \
44   }
45
46 DECL_GSSVX(SuperLU_S,sgssvx,float,float)
47 DECL_GSSVX(SuperLU_C,cgssvx,float,std::complex<float>)
48 DECL_GSSVX(SuperLU_D,dgssvx,double,double)
49 DECL_GSSVX(SuperLU_Z,zgssvx,double,std::complex<double>)
50
51 template<typename MatrixType>
52 struct SluMatrixMapHelper;
53
54 /** \internal
55   *
56   * A wrapper class for SuperLU matrices. It supports only compressed sparse matrices
57   * and dense matrices. Supernodal and other fancy format are not supported by this wrapper.
58   *
59   * This wrapper class mainly aims to avoids the need of dynamic allocation of the storage structure.
60   */
61 struct SluMatrix : SuperMatrix
62 {
63   SluMatrix()
64   {
65     Store = &storage;
66   }
67
68   SluMatrix(const SluMatrix& other)
69     : SuperMatrix(other)
70   {
71     Store = &storage;
72     storage = other.storage;
73   }
74   
75   SluMatrix& operator=(const SluMatrix& other)
76   {
77     SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
78     Store = &storage;
79     storage = other.storage;
80     return *this;
81   }
82
83   struct
84   {
85     union {int nnz;int lda;};
86     void *values;
87     int *innerInd;
88     int *outerInd;
89   } storage;
90
91   void setStorageType(Stype_t t)
92   {
93     Stype = t;
94     if (t==SLU_NC || t==SLU_NR || t==SLU_DN)
95       Store = &storage;
96     else
97     {
98       ei_assert(false && "storage type not supported");
99       Store = 0;
100     }
101   }
102
103   template<typename Scalar>
104   void setScalarType()
105   {
106     if (ei_is_same_type<Scalar,float>::ret)
107       Dtype = SLU_S;
108     else if (ei_is_same_type<Scalar,double>::ret)
109       Dtype = SLU_D;
110     else if (ei_is_same_type<Scalar,std::complex<float> >::ret)
111       Dtype = SLU_C;
112     else if (ei_is_same_type<Scalar,std::complex<double> >::ret)
113       Dtype = SLU_Z;
114     else
115     {
116       ei_assert(false && "Scalar type not supported by SuperLU");
117     }
118   }
119
120   template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
121   static SluMatrix Map(Matrix<Scalar,Rows,Cols,Options,MRows,MCols>& mat)
122   {
123     typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType;
124     ei_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU");
125     SluMatrix res;
126     res.setStorageType(SLU_DN);
127     res.setScalarType<Scalar>();
128     res.Mtype     = SLU_GE;
129
130     res.nrow      = mat.rows();
131     res.ncol      = mat.cols();
132
133     res.storage.lda       = mat.stride();
134     res.storage.values    = mat.data();
135     return res;
136   }
137
138   template<typename MatrixType>
139   static SluMatrix Map(SparseMatrixBase<MatrixType>& mat)
140   {
141     SluMatrix res;
142     if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
143     {
144       res.setStorageType(SLU_NR);
145       res.nrow      = mat.cols();
146       res.ncol      = mat.rows();
147     }
148     else
149     {
150       res.setStorageType(SLU_NC);
151       res.nrow      = mat.rows();
152       res.ncol      = mat.cols();
153     }
154
155     res.Mtype     = SLU_GE;
156
157     res.storage.nnz       = mat.nonZeros();
158     res.storage.values    = mat.derived()._valuePtr();
159     res.storage.innerInd  = mat.derived()._innerIndexPtr();
160     res.storage.outerInd  = mat.derived()._outerIndexPtr();
161
162     res.setScalarType<typename MatrixType::Scalar>();
163
164     // FIXME the following is not very accurate
165     if (MatrixType::Flags & UpperTriangular)
166       res.Mtype = SLU_TRU;
167     if (MatrixType::Flags & LowerTriangular)
168       res.Mtype = SLU_TRL;
169     if (MatrixType::Flags & SelfAdjoint)
170       ei_assert(false && "SelfAdjoint matrix shape not supported by SuperLU");
171     return res;
172   }
173 };
174
175 template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
176 struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,Options,MRows,MCols> >
177 {
178   typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType;
179   static void run(MatrixType& mat, SluMatrix& res)
180   {
181     ei_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU");
182     res.setStorageType(SLU_DN);
183     res.setScalarType<Scalar>();
184     res.Mtype     = SLU_GE;
185
186     res.nrow      = mat.rows();
187     res.ncol      = mat.cols();
188
189     res.storage.lda       = mat.stride();
190     res.storage.values    = mat.data();
191   }
192 };
193
194 template<typename Derived>
195 struct SluMatrixMapHelper<SparseMatrixBase<Derived> >
196 {
197   typedef Derived MatrixType;
198   static void run(MatrixType& mat, SluMatrix& res)
199   {
200     if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
201     {
202       res.setStorageType(SLU_NR);
203       res.nrow      = mat.cols();
204       res.ncol      = mat.rows();
205     }
206     else
207     {
208       res.setStorageType(SLU_NC);
209       res.nrow      = mat.rows();
210       res.ncol      = mat.cols();
211     }
212
213     res.Mtype     = SLU_GE;
214
215     res.storage.nnz       = mat.nonZeros();
216     res.storage.values    = mat._valuePtr();
217     res.storage.innerInd  = mat._innerIndexPtr();
218     res.storage.outerInd  = mat._outerIndexPtr();
219
220     res.setScalarType<typename MatrixType::Scalar>();
221
222     // FIXME the following is not very accurate
223     if (MatrixType::Flags & UpperTriangular)
224       res.Mtype = SLU_TRU;
225     if (MatrixType::Flags & LowerTriangular)
226       res.Mtype = SLU_TRL;
227     if (MatrixType::Flags & SelfAdjoint)
228       ei_assert(false && "SelfAdjoint matrix shape not supported by SuperLU");
229   }
230 };
231
232 template<typename Derived>
233 SluMatrix SparseMatrixBase<Derived>::asSluMatrix()
234 {
235   return SluMatrix::Map(derived());
236 }
237
238 /** View a Super LU matrix as an Eigen expression */
239 template<typename Scalar, int Flags>
240 MappedSparseMatrix<Scalar,Flags>::MappedSparseMatrix(SluMatrix& sluMat)
241 {
242   if ((Flags&RowMajorBit)==RowMajorBit)
243   {
244     assert(sluMat.Stype == SLU_NR);
245     m_innerSize   = sluMat.ncol;
246     m_outerSize   = sluMat.nrow;
247   }
248   else
249   {
250     assert(sluMat.Stype == SLU_NC);
251     m_innerSize   = sluMat.nrow;
252     m_outerSize   = sluMat.ncol;
253   }
254   m_outerIndex = sluMat.storage.outerInd;
255   m_innerIndices = sluMat.storage.innerInd;
256   m_values = reinterpret_cast<Scalar*>(sluMat.storage.values);
257   m_nnz = sluMat.storage.outerInd[m_outerSize];
258 }
259
260 template<typename MatrixType>
261 class SparseLU<MatrixType,SuperLU> : public SparseLU<MatrixType>
262 {
263   protected:
264     typedef SparseLU<MatrixType> Base;
265     typedef typename Base::Scalar Scalar;
266     typedef typename Base::RealScalar RealScalar;
267     typedef Matrix<Scalar,Dynamic,1> Vector;
268     typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
269     typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
270     typedef SparseMatrix<Scalar,LowerTriangular|UnitDiagBit> LMatrixType;
271     typedef SparseMatrix<Scalar,UpperTriangular> UMatrixType;
272     using Base::m_flags;
273     using Base::m_status;
274
275   public:
276
277     SparseLU(int flags = NaturalOrdering)
278       : Base(flags)
279     {
280     }
281
282     SparseLU(const MatrixType& matrix, int flags = NaturalOrdering)
283       : Base(flags)
284     {
285       compute(matrix);
286     }
287
288     ~SparseLU()
289     {
290     }
291
292     inline const LMatrixType& matrixL() const
293     {
294       if (m_extractedDataAreDirty) extractData();
295       return m_l;
296     }
297
298     inline const UMatrixType& matrixU() const
299     {
300       if (m_extractedDataAreDirty) extractData();
301       return m_u;
302     }
303
304     inline const IntColVectorType& permutationP() const
305     {
306       if (m_extractedDataAreDirty) extractData();
307       return m_p;
308     }
309
310     inline const IntRowVectorType& permutationQ() const
311     {
312       if (m_extractedDataAreDirty) extractData();
313       return m_q;
314     }
315
316     Scalar determinant() const;
317
318     template<typename BDerived, typename XDerived>
319     bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const;
320
321     void compute(const MatrixType& matrix);
322
323   protected:
324
325     void extractData() const;
326
327   protected:
328     // cached data to reduce reallocation, etc.
329     mutable LMatrixType m_l;
330     mutable UMatrixType m_u;
331     mutable IntColVectorType m_p;
332     mutable IntRowVectorType m_q;
333
334     mutable SparseMatrix<Scalar> m_matrix;
335     mutable SluMatrix m_sluA;
336     mutable SuperMatrix m_sluL, m_sluU;
337     mutable SluMatrix m_sluB, m_sluX;
338     mutable SuperLUStat_t m_sluStat;
339     mutable superlu_options_t m_sluOptions;
340     mutable std::vector<int> m_sluEtree;
341     mutable std::vector<RealScalar> m_sluRscale, m_sluCscale;
342     mutable std::vector<RealScalar> m_sluFerr, m_sluBerr;
343     mutable char m_sluEqued;
344     mutable bool m_extractedDataAreDirty;
345 };
346
347 template<typename MatrixType>
348 void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
349 {
350   const int size = a.rows();
351   m_matrix = a;
352
353   set_default_options(&m_sluOptions);
354   m_sluOptions.ColPerm = NATURAL;
355   m_sluOptions.PrintStat = NO;
356   m_sluOptions.ConditionNumber = NO;
357   m_sluOptions.Trans = NOTRANS;
358   // m_sluOptions.Equil = NO;
359
360   switch (Base::orderingMethod())
361   {
362       case NaturalOrdering          : m_sluOptions.ColPerm = NATURAL; break;
363       case MinimumDegree_AT_PLUS_A  : m_sluOptions.ColPerm = MMD_AT_PLUS_A; break;
364       case MinimumDegree_ATA        : m_sluOptions.ColPerm = MMD_ATA; break;
365       case ColApproxMinimumDegree   : m_sluOptions.ColPerm = COLAMD; break;
366       default:
367         std::cerr << "Eigen: ordering method \"" << Base::orderingMethod() << "\" not supported by the SuperLU backend\n";
368         m_sluOptions.ColPerm = NATURAL;
369   };
370
371   m_sluA = m_matrix.asSluMatrix();
372   memset(&m_sluL,0,sizeof m_sluL);
373   memset(&m_sluU,0,sizeof m_sluU);
374   m_sluEqued = 'B';
375   int info = 0;
376
377   m_p.resize(size);
378   m_q.resize(size);
379   m_sluRscale.resize(size);
380   m_sluCscale.resize(size);
381   m_sluEtree.resize(size);
382
383   RealScalar recip_pivot_gross, rcond;
384   RealScalar ferr, berr;
385
386   // set empty B and X
387   m_sluB.setStorageType(SLU_DN);
388   m_sluB.setScalarType<Scalar>();
389   m_sluB.Mtype = SLU_GE;
390   m_sluB.storage.values = 0;
391   m_sluB.nrow = m_sluB.ncol = 0;
392   m_sluB.storage.lda = size;
393   m_sluX = m_sluB;
394
395   StatInit(&m_sluStat);
396   SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
397           &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
398           &m_sluL, &m_sluU,
399           NULL, 0,
400           &m_sluB, &m_sluX,
401           &recip_pivot_gross, &rcond,
402           &ferr, &berr,
403           &m_sluStat, &info, Scalar());
404   StatFree(&m_sluStat);
405
406   m_extractedDataAreDirty = true;
407
408   // FIXME how to better check for errors ???
409   Base::m_succeeded = (info == 0);
410 }
411
412 template<typename MatrixType>
413 template<typename BDerived,typename XDerived>
414 bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> *x) const
415 {
416   const int size = m_matrix.rows();
417   const int rhsCols = b.cols();
418   ei_assert(size==b.rows());
419
420   m_sluOptions.Fact = FACTORED;
421   m_sluOptions.IterRefine = NOREFINE;
422
423   m_sluFerr.resize(rhsCols);
424   m_sluBerr.resize(rhsCols);
425   m_sluB = SluMatrix::Map(b.const_cast_derived());
426   m_sluX = SluMatrix::Map(x->derived());
427
428   StatInit(&m_sluStat);
429   int info = 0;
430   RealScalar recip_pivot_gross, rcond;
431   SuperLU_gssvx(
432     &m_sluOptions, &m_sluA,
433     m_q.data(), m_p.data(),
434     &m_sluEtree[0], &m_sluEqued,
435     &m_sluRscale[0], &m_sluCscale[0],
436     &m_sluL, &m_sluU,
437     NULL, 0,
438     &m_sluB, &m_sluX,
439     &recip_pivot_gross, &rcond,
440     &m_sluFerr[0], &m_sluBerr[0],
441     &m_sluStat, &info, Scalar());
442   StatFree(&m_sluStat);
443
444   return info==0;
445 }
446
447 //
448 // the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
449 //
450 //  Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
451 //
452 //  THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
453 //  EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
454 //
455 template<typename MatrixType>
456 void SparseLU<MatrixType,SuperLU>::extractData() const
457 {
458   if (m_extractedDataAreDirty)
459   {
460     int         upper;
461     int         fsupc, istart, nsupr;
462     int         lastl = 0, lastu = 0;
463     SCformat    *Lstore = static_cast<SCformat*>(m_sluL.Store);
464     NCformat    *Ustore = static_cast<NCformat*>(m_sluU.Store);
465     Scalar      *SNptr;
466
467     const int size = m_matrix.rows();
468     m_l.resize(size,size);
469     m_l.resizeNonZeros(Lstore->nnz);
470     m_u.resize(size,size);
471     m_u.resizeNonZeros(Ustore->nnz);
472
473     int* Lcol = m_l._outerIndexPtr();
474     int* Lrow = m_l._innerIndexPtr();
475     Scalar* Lval = m_l._valuePtr();
476
477     int* Ucol = m_u._outerIndexPtr();
478     int* Urow = m_u._innerIndexPtr();
479     Scalar* Uval = m_u._valuePtr();
480
481     Ucol[0] = 0;
482     Ucol[0] = 0;
483
484     /* for each supernode */
485     for (int k = 0; k <= Lstore->nsuper; ++k)
486     {
487       fsupc   = L_FST_SUPC(k);
488       istart  = L_SUB_START(fsupc);
489       nsupr   = L_SUB_START(fsupc+1) - istart;
490       upper   = 1;
491
492       /* for each column in the supernode */
493       for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
494       {
495         SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
496
497         /* Extract U */
498         for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
499         {
500           Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
501           /* Matlab doesn't like explicit zero. */
502           if (Uval[lastu] != 0.0)
503             Urow[lastu++] = U_SUB(i);
504         }
505         for (int i = 0; i < upper; ++i)
506         {
507           /* upper triangle in the supernode */
508           Uval[lastu] = SNptr[i];
509           /* Matlab doesn't like explicit zero. */
510           if (Uval[lastu] != 0.0)
511             Urow[lastu++] = L_SUB(istart+i);
512         }
513         Ucol[j+1] = lastu;
514
515         /* Extract L */
516         Lval[lastl] = 1.0; /* unit diagonal */
517         Lrow[lastl++] = L_SUB(istart + upper - 1);
518         for (int i = upper; i < nsupr; ++i)
519         {
520           Lval[lastl] = SNptr[i];
521           /* Matlab doesn't like explicit zero. */
522           if (Lval[lastl] != 0.0)
523             Lrow[lastl++] = L_SUB(istart+i);
524         }
525         Lcol[j+1] = lastl;
526
527         ++upper;
528       } /* for j ... */
529
530     } /* for k ... */
531
532     // squeeze the matrices :
533     m_l.resizeNonZeros(lastl);
534     m_u.resizeNonZeros(lastu);
535
536     m_extractedDataAreDirty = false;
537   }
538 }
539
540 template<typename MatrixType>
541 typename SparseLU<MatrixType,SuperLU>::Scalar SparseLU<MatrixType,SuperLU>::determinant() const
542 {
543   if (m_extractedDataAreDirty)
544     extractData();
545
546   // TODO this code coule be moved to the default/base backend
547   // FIXME perhaps we have to take into account the scale factors m_sluRscale and m_sluCscale ???
548   Scalar det = Scalar(1);
549   for (int j=0; j<m_u.cols(); ++j)
550   {
551     if (m_u._outerIndexPtr()[j+1]-m_u._outerIndexPtr()[j] > 0)
552     {
553       int lastId = m_u._outerIndexPtr()[j+1]-1;
554       ei_assert(m_u._innerIndexPtr()[lastId]<=j);
555       if (m_u._innerIndexPtr()[lastId]==j)
556       {
557         det *= m_u._valuePtr()[lastId];
558       }
559     }
560     // std::cout << m_sluRscale[j] << " " << m_sluCscale[j] << "   ";
561   }
562   return det;
563 }
564
565 #endif // EIGEN_SUPERLUSUPPORT_H