Fixed segfault error with the information/statistics display
[mldemos:mldemos.git] / _AlgorithmsPlugins / KPCA / Eigen / src / Sparse / SparseSelfAdjointView.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.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_SPARSE_SELFADJOINTVIEW_H
26 #define EIGEN_SPARSE_SELFADJOINTVIEW_H
27
28 /** \class SparseSelfAdjointView
29   *
30   *
31   * \brief Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix.
32   *
33   * \param MatrixType the type of the dense matrix storing the coefficients
34   * \param UpLo can be either \c Lower or \c Upper
35   *
36   * This class is an expression of a sefladjoint matrix from a triangular part of a matrix
37   * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView()
38   * and most of the time this is the only way that it is used.
39   *
40   * \sa SparseMatrixBase::selfAdjointView()
41   */
42 template<typename Lhs, typename Rhs, int UpLo>
43 class SparseSelfAdjointTimeDenseProduct;
44
45 template<typename Lhs, typename Rhs, int UpLo>
46 class DenseTimeSparseSelfAdjointProduct;
47
48 template<typename MatrixType,int UpLo>
49 class SparseSymmetricPermutationProduct;
50
51 namespace internal {
52   
53 template<typename MatrixType, unsigned int UpLo>
54 struct traits<SparseSelfAdjointView<MatrixType,UpLo> > : traits<MatrixType> {
55 };
56
57 template<int SrcUpLo,int DstUpLo,typename MatrixType,int DestOrder>
58 void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm = 0);
59
60 template<int UpLo,typename MatrixType,int DestOrder>
61 void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm = 0);
62
63 }
64
65 template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView
66   : public EigenBase<SparseSelfAdjointView<MatrixType,UpLo> >
67 {
68   public:
69
70     typedef typename MatrixType::Scalar Scalar;
71     typedef typename MatrixType::Index Index;
72     typedef Matrix<Index,Dynamic,1> VectorI;
73     typedef typename MatrixType::Nested MatrixTypeNested;
74     typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
75
76     inline SparseSelfAdjointView(const MatrixType& matrix) : m_matrix(matrix)
77     {
78       eigen_assert(rows()==cols() && "SelfAdjointView is only for squared matrices");
79     }
80
81     inline Index rows() const { return m_matrix.rows(); }
82     inline Index cols() const { return m_matrix.cols(); }
83
84     /** \internal \returns a reference to the nested matrix */
85     const _MatrixTypeNested& matrix() const { return m_matrix; }
86     _MatrixTypeNested& matrix() { return m_matrix.const_cast_derived(); }
87
88     /** Efficient sparse self-adjoint matrix times dense vector/matrix product */
89     template<typename OtherDerived>
90     SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>
91     operator*(const MatrixBase<OtherDerived>& rhs) const
92     {
93       return SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>(m_matrix, rhs.derived());
94     }
95
96     /** Efficient dense vector/matrix times sparse self-adjoint matrix product */
97     template<typename OtherDerived> friend
98     DenseTimeSparseSelfAdjointProduct<OtherDerived,MatrixType,UpLo>
99     operator*(const MatrixBase<OtherDerived>& lhs, const SparseSelfAdjointView& rhs)
100     {
101       return DenseTimeSparseSelfAdjointProduct<OtherDerived,_MatrixTypeNested,UpLo>(lhs.derived(), rhs.m_matrix);
102     }
103
104     /** Perform a symmetric rank K update of the selfadjoint matrix \c *this:
105       * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix.
106       *
107       * \returns a reference to \c *this
108       *
109       * Note that it is faster to set alpha=0 than initializing the matrix to zero
110       * and then keep the default value alpha=1.
111       *
112       * To perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply
113       * call this function with u.adjoint().
114       */
115     template<typename DerivedU>
116     SparseSelfAdjointView& rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
117     
118     /** \internal triggered by sparse_matrix = SparseSelfadjointView; */
119     template<typename DestScalar> void evalTo(SparseMatrix<DestScalar>& _dest) const
120     {
121       internal::permute_symm_to_fullsymm<UpLo>(m_matrix, _dest);
122     }
123     
124     template<typename DestScalar> void evalTo(DynamicSparseMatrix<DestScalar>& _dest) const
125     {
126       // TODO directly evaluate into _dest;
127       SparseMatrix<DestScalar> tmp(_dest.rows(),_dest.cols());
128       internal::permute_symm_to_fullsymm<UpLo>(m_matrix, tmp);
129       _dest = tmp;
130     }
131     
132     /** \returns an expression of P^-1 H P */
133     SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo> twistedBy(const PermutationMatrix<Dynamic>& perm) const
134     {
135       return SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo>(m_matrix, perm);
136     }
137     
138     template<typename SrcMatrixType,int SrcUpLo>
139     SparseSelfAdjointView& operator=(const SparseSymmetricPermutationProduct<SrcMatrixType,SrcUpLo>& permutedMatrix)
140     {
141       permutedMatrix.evalTo(*this);
142       return *this;
143     }
144     
145
146     // const SparseLLT<PlainObject, UpLo> llt() const;
147     // const SparseLDLT<PlainObject, UpLo> ldlt() const;
148
149   protected:
150
151     const typename MatrixType::Nested m_matrix;
152     mutable VectorI m_countPerRow;
153     mutable VectorI m_countPerCol;
154 };
155
156 /***************************************************************************
157 * Implementation of SparseMatrixBase methods
158 ***************************************************************************/
159
160 template<typename Derived>
161 template<unsigned int UpLo>
162 const SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView() const
163 {
164   return derived();
165 }
166
167 template<typename Derived>
168 template<unsigned int UpLo>
169 SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView()
170 {
171   return derived();
172 }
173
174 /***************************************************************************
175 * Implementation of SparseSelfAdjointView methods
176 ***************************************************************************/
177
178 template<typename MatrixType, unsigned int UpLo>
179 template<typename DerivedU>
180 SparseSelfAdjointView<MatrixType,UpLo>&
181 SparseSelfAdjointView<MatrixType,UpLo>::rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha)
182 {
183   SparseMatrix<Scalar,MatrixType::Flags&RowMajorBit?RowMajor:ColMajor> tmp = u * u.adjoint();
184   if(alpha==Scalar(0))
185     m_matrix.const_cast_derived() = tmp.template triangularView<UpLo>();
186   else
187     m_matrix.const_cast_derived() += alpha * tmp.template triangularView<UpLo>();
188
189   return *this;
190 }
191
192 /***************************************************************************
193 * Implementation of sparse self-adjoint time dense matrix
194 ***************************************************************************/
195
196 namespace internal {
197 template<typename Lhs, typename Rhs, int UpLo>
198 struct traits<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo> >
199  : traits<ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
200 {
201   typedef Dense StorageKind;
202 };
203 }
204
205 template<typename Lhs, typename Rhs, int UpLo>
206 class SparseSelfAdjointTimeDenseProduct
207   : public ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
208 {
209   public:
210     EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseSelfAdjointTimeDenseProduct)
211
212     SparseSelfAdjointTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
213     {}
214
215     template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const
216     {
217       // TODO use alpha
218       eigen_assert(alpha==Scalar(1) && "alpha != 1 is not implemented yet, sorry");
219       typedef typename internal::remove_all<Lhs>::type _Lhs;
220       typedef typename internal::remove_all<Rhs>::type _Rhs;
221       typedef typename _Lhs::InnerIterator LhsInnerIterator;
222       enum {
223         LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit,
224         ProcessFirstHalf =
225                  ((UpLo&(Upper|Lower))==(Upper|Lower))
226               || ( (UpLo&Upper) && !LhsIsRowMajor)
227               || ( (UpLo&Lower) && LhsIsRowMajor),
228         ProcessSecondHalf = !ProcessFirstHalf
229       };
230       for (Index j=0; j<m_lhs.outerSize(); ++j)
231       {
232         LhsInnerIterator i(m_lhs,j);
233         if (ProcessSecondHalf && i && (i.index()==j))
234         {
235           dest.row(j) += i.value() * m_rhs.row(j);
236           ++i;
237         }
238         Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest.row(LhsIsRowMajor ? j : 0));
239         for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
240         {
241           Index a = LhsIsRowMajor ? j : i.index();
242           Index b = LhsIsRowMajor ? i.index() : j;
243           typename Lhs::Scalar v = i.value();
244           dest.row(a) += (v) * m_rhs.row(b);
245           dest.row(b) += internal::conj(v) * m_rhs.row(a);
246         }
247         if (ProcessFirstHalf && i && (i.index()==j))
248           dest.row(j) += i.value() * m_rhs.row(j);
249       }
250     }
251
252   private:
253     SparseSelfAdjointTimeDenseProduct& operator=(const SparseSelfAdjointTimeDenseProduct&);
254 };
255
256 namespace internal {
257 template<typename Lhs, typename Rhs, int UpLo>
258 struct traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> >
259  : traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
260 {};
261 }
262
263 template<typename Lhs, typename Rhs, int UpLo>
264 class DenseTimeSparseSelfAdjointProduct
265   : public ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
266 {
267   public:
268     EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseSelfAdjointProduct)
269
270     DenseTimeSparseSelfAdjointProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
271     {}
272
273     template<typename Dest> void scaleAndAddTo(Dest& /*dest*/, Scalar /*alpha*/) const
274     {
275       // TODO
276     }
277
278   private:
279     DenseTimeSparseSelfAdjointProduct& operator=(const DenseTimeSparseSelfAdjointProduct&);
280 };
281
282 /***************************************************************************
283 * Implementation of symmetric copies and permutations
284 ***************************************************************************/
285 namespace internal {
286   
287 template<typename MatrixType, int UpLo>
288 struct traits<SparseSymmetricPermutationProduct<MatrixType,UpLo> > : traits<MatrixType> {
289 };
290
291 template<int UpLo,typename MatrixType,int DestOrder>
292 void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm)
293 {
294   typedef typename MatrixType::Index Index;
295   typedef typename MatrixType::Scalar Scalar;
296   typedef SparseMatrix<Scalar,DestOrder,Index> Dest;
297   typedef Matrix<Index,Dynamic,1> VectorI;
298   
299   Dest& dest(_dest.derived());
300   enum {
301     StorageOrderMatch = int(Dest::IsRowMajor) == int(MatrixType::IsRowMajor)
302   };
303   eigen_assert(perm==0);
304   Index size = mat.rows();
305   VectorI count;
306   count.resize(size);
307   count.setZero();
308   dest.resize(size,size);
309   for(Index j = 0; j<size; ++j)
310   {
311     Index jp = perm ? perm[j] : j;
312     for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
313     {
314       Index i = it.index();
315       Index ip = perm ? perm[i] : i;
316       if(i==j)
317         count[ip]++;
318       else if((UpLo==Lower && i>j) || (UpLo==Upper && i<j))
319       {
320         count[ip]++;
321         count[jp]++;
322       }
323     }
324   }
325   Index nnz = count.sum();
326   
327   // reserve space
328   dest.reserve(nnz);
329   dest._outerIndexPtr()[0] = 0;
330   for(Index j=0; j<size; ++j)
331     dest._outerIndexPtr()[j+1] = dest._outerIndexPtr()[j] + count[j];
332   for(Index j=0; j<size; ++j)
333     count[j] = dest._outerIndexPtr()[j];
334   
335   // copy data
336   for(Index j = 0; j<size; ++j)
337   {
338     Index jp = perm ? perm[j] : j;
339     for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
340     {
341       Index i = it.index();
342       Index ip = perm ? perm[i] : i;
343       if(i==j)
344       {
345         int k = count[ip]++;
346         dest._innerIndexPtr()[k] = ip;
347         dest._valuePtr()[k] = it.value();
348       }
349       else if((UpLo==Lower && i>j) || (UpLo==Upper && i<j))
350       {
351         int k = count[jp]++;
352         dest._innerIndexPtr()[k] = ip;
353         dest._valuePtr()[k] = it.value();
354         k = count[ip]++;
355         dest._innerIndexPtr()[k] = jp;
356         dest._valuePtr()[k] = internal::conj(it.value());
357       }
358     }
359   }
360 }
361
362 template<int SrcUpLo,int DstUpLo,typename MatrixType,int DestOrder>
363 void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm)
364 {
365   typedef typename MatrixType::Index Index;
366   typedef typename MatrixType::Scalar Scalar;
367   typedef SparseMatrix<Scalar,DestOrder,Index> Dest;
368   Dest& dest(_dest.derived());
369   typedef Matrix<Index,Dynamic,1> VectorI;
370   //internal::conj_if<SrcUpLo!=DstUpLo> cj;
371   
372   Index size = mat.rows();
373   VectorI count(size);
374   count.setZero();
375   dest.resize(size,size);
376   for(Index j = 0; j<size; ++j)
377   {
378     Index jp = perm ? perm[j] : j;
379     for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
380     {
381       Index i = it.index();
382       if((SrcUpLo==Lower && i<j) || (SrcUpLo==Upper && i>j))
383         continue;
384                   
385       Index ip = perm ? perm[i] : i;
386       count[DstUpLo==Lower ? std::min(ip,jp) : std::max(ip,jp)]++;
387     }
388   }
389   dest._outerIndexPtr()[0] = 0;
390   for(Index j=0; j<size; ++j)
391     dest._outerIndexPtr()[j+1] = dest._outerIndexPtr()[j] + count[j];
392   dest.resizeNonZeros(dest._outerIndexPtr()[size]);
393   for(Index j=0; j<size; ++j)
394     count[j] = dest._outerIndexPtr()[j];
395   
396   for(Index j = 0; j<size; ++j)
397   {
398     Index jp = perm ? perm[j] : j;
399     for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
400     {
401       Index i = it.index();
402       if((SrcUpLo==Lower && i<j) || (SrcUpLo==Upper && i>j))
403         continue;
404                   
405       Index ip = perm? perm[i] : i;
406       Index k = count[DstUpLo==Lower ? std::min(ip,jp) : std::max(ip,jp)]++;
407       dest._innerIndexPtr()[k] = DstUpLo==Lower ? std::max(ip,jp) : std::min(ip,jp);
408       
409       if((DstUpLo==Lower && ip<jp) || (DstUpLo==Upper && ip>jp))
410         dest._valuePtr()[k] = conj(it.value());
411       else
412         dest._valuePtr()[k] = it.value();
413     }
414   }
415 }
416
417 }
418
419 template<typename MatrixType,int UpLo>
420 class SparseSymmetricPermutationProduct
421   : public EigenBase<SparseSymmetricPermutationProduct<MatrixType,UpLo> >
422 {
423     typedef PermutationMatrix<Dynamic> Perm;
424   public:
425     typedef typename MatrixType::Scalar Scalar;
426     typedef typename MatrixType::Index Index;
427     typedef Matrix<Index,Dynamic,1> VectorI;
428     typedef typename MatrixType::Nested MatrixTypeNested;
429     typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
430     
431     SparseSymmetricPermutationProduct(const MatrixType& mat, const Perm& perm)
432       : m_matrix(mat), m_perm(perm)
433     {}
434     
435     inline Index rows() const { return m_matrix.rows(); }
436     inline Index cols() const { return m_matrix.cols(); }
437     
438     template<typename DestScalar> void evalTo(SparseMatrix<DestScalar>& _dest) const
439     {
440       internal::permute_symm_to_fullsymm<UpLo>(m_matrix,_dest,m_perm.indices().data());
441     }
442     
443     template<typename DestType,unsigned int DestUpLo> void evalTo(SparseSelfAdjointView<DestType,DestUpLo>& dest) const
444     {
445       internal::permute_symm_to_symm<UpLo,DestUpLo>(m_matrix,dest.matrix(),m_perm.indices().data());
446     }
447     
448   protected:
449     const MatrixTypeNested m_matrix;
450     const Perm& m_perm;
451
452 };
453
454 #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H