daa7564777de14b09bae98dd5fd090475d4bdd15
[scilab.git] / scilab / modules / ast / src / cpp / types / sparse.cpp
1 /*
2 *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 *  Copyright (C) 2010 - DIGITEO - Bernard HUGUENEY
4 *
5 *  This file must be used under the terms of the CeCILL.
6 *  This source file is licensed as described in the file COPYING, which
7 *  you should have received as part of this distribution.  The terms
8 *  are also available at
9 *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 *
11 */
12
13 #include <sstream>
14 #include <math.h>
15 #include <Eigen/Sparse>
16 #include <complex>
17 #include <iterator>
18 #include <algorithm>
19
20 #include <Eigen/Core>
21 #include <Eigen/IterativeLinearSolvers>
22 #include <Eigen/SparseCholesky>
23
24 #include "sparse.hxx"
25 #include "types.hxx"
26 #include "tostring_common.hxx"
27 #include "double.hxx"
28 #include "matrixiterator.hxx"
29 #include "types_subtraction.hxx"
30 #include "types_addition.hxx"
31 #include "types_multiplication.hxx"
32 #include "configvariable.hxx"
33 #include "scilabWrite.hxx"
34
35 #include "sparseOp.hxx"
36
37 extern "C"
38 {
39 #include "elem_common.h"
40 }
41 namespace
42 {
43
44 /* used for debuging output
45 */
46 template<typename Os, typename In, typename Sz> Os& writeData(wchar_t const* title, In beg, Sz n, Os& os)
47 {
48     os << title;
49     /* TODO: use tostring_common (with a kind of std::boolalpha for boolean output)
50     */
51     mycopy_n(beg, n, std::ostream_iterator<typename std::iterator_traits<In>::value_type, char>(os, L" "));
52     os << std::endl;
53     return os;
54 }
55
56 struct Printer
57 {
58     Printer (int precision) : p(precision)
59     {
60     }
61     template<typename T>
62     std::wstring emptyName( /* */) const
63     {
64         return L" zero";
65     }
66
67     template<typename T>
68     std::wstring operator()(T const& t) const
69     {
70         //never call ?
71         std::wostringstream ostr;
72         ostr.precision(p);
73         ostr << t;
74         return ostr.str();
75     }
76     int p;
77 };
78
79 template<>
80 std::wstring Printer::operator()(bool const& b) const
81 {
82     if (b)
83     {
84         return L"T";
85     }
86     else
87     {
88         return L"F";
89     }
90 }
91
92 template<>
93 std::wstring Printer::operator()(double const& d) const
94 {
95     std::wostringstream ostr;
96     DoubleFormat df;
97     getDoubleFormat(d, &df);
98     addDoubleValue(&ostr, d, &df);
99     return ostr.str();
100 }
101
102 template<>
103 std::wstring Printer::operator()(std::complex<double > const& c) const
104 {
105     std::wostringstream ostr;
106     int iLen = 0;
107     DoubleFormat dfR, dfI;
108     getComplexFormat(c.real(), c.imag(), &iLen, &dfR, &dfI);
109     addDoubleComplexValue(&ostr, c.real(), c.imag(), iLen, &dfR, &dfI);
110     return ostr.str();
111 }
112
113 template<>
114 std::wstring Printer::emptyName<bool>() const
115 {
116     return L"False";
117 }
118
119
120 template<typename T> std::wstring toString(T const& m, int precision)
121 {
122     std::wostringstream ostr;
123
124     int iWidthRows  = 0;
125     int iWidthCols  = 0;
126     getSignedIntFormat(m.rows(), &iWidthRows);
127     getSignedIntFormat(m.cols(), &iWidthCols);
128
129     ostr << L"(" ;
130     addUnsignedIntValue<unsigned long long>(&ostr, m.rows(), iWidthRows);
131     ostr << ",";
132     addUnsignedIntValue<unsigned long long>(&ostr, m.cols(), iWidthCols);
133     ostr << L")";
134
135     Printer p(precision);
136     if (!m.nonZeros())
137     {
138         ostr << ( p.emptyName<typename Eigen::internal::traits<T>::Scalar>());
139     }
140     ostr << L" sparse matrix\n\n";
141
142     const typename Eigen::internal::traits<T>::Index* pIColPos      = m.innerIndexPtr();
143     const typename Eigen::internal::traits<T>::Index* pINbItemByRow = m.outerIndexPtr();
144
145     int iPos = 0;
146
147     for (size_t j = 1 ; j < m.rows() + 1 ; j++)
148     {
149         for (size_t i = pINbItemByRow[j - 1] ; i < pINbItemByRow[j] ; i++)
150         {
151             ostr << L"(";
152             addUnsignedIntValue<unsigned long long>(&ostr, (int)j, iWidthRows);
153             ostr << L",";
154             addUnsignedIntValue<unsigned long long>(&ostr, pIColPos[iPos] + 1, iWidthCols);
155             ostr << L")\t" << p(m.valuePtr()[iPos]) << std::endl;
156
157             iPos++;
158         }
159     }
160
161     return ostr.str();
162 }
163
164 /** utility function to compare two Eigen::Sparse matrices to equality
165 */
166 template<typename T> bool equal(T const& s1, T const& s2)
167 {
168     bool res(true);
169     // only compares elts when both inner iterators are "defined", so we assert that we compared all the non zero values
170     // i.e. the inner iterators where defined for the same values
171     std::size_t nbElts(0);
172
173     for (int k = 0; res && k != s1.outerSize(); ++k)
174     {
175         for (typename T::InnerIterator it1(s1, k), it2(s2, k); res && it1 && it2 ; ++it1, ++it2, ++nbElts)
176         {
177             res = (it1.value() == it2.value()
178                    && it1.row() == it2.row()
179                    && it1.col() == it2.col());
180         }
181     }
182     return res && (nbElts == s1.nonZeros()) && (nbElts == s2.nonZeros());
183 }
184 /**
185 utility function to set non zero values of an Eigen::Sparse matrix to a fixed values
186 @param s : sparse matrix to modify
187 @param v : value to set (default to 1.)
188 */
189 template<typename T> bool setNonZero(T& s, typename Eigen::internal::traits<T>::Scalar v = 1.)
190 {
191     for (typename Eigen::internal::traits<T>::Index j = 0; j < s.outerSize(); ++j)
192     {
193         for (typename T::InnerIterator it(s, j); it; ++it)
194         {
195             it.valueRef() = v;
196         }
197     }
198     return true;
199 }
200
201
202
203 template<typename Src, typename Sp>
204 void doAppend(Src SPARSE_CONST& src, int r, int c, Sp& dest)
205 {
206     typedef typename Eigen::internal::traits<Sp>::Scalar data_t;
207     mycopy_n(makeMatrixIterator<data_t>(src, makeNonZerosIterator(src)), nonZeros(src)
208              , makeMatrixIterator<data_t>(dest, makeTranslatedIterator(makeNonZerosIterator(src), Coords2D(r, c))));
209 }
210
211 template<typename Scalar1, typename Scalar2>
212 void doAppend(Eigen::SparseMatrix<Scalar1, Eigen::RowMajor> SPARSE_CONST& src, int r, int c, Eigen::SparseMatrix<Scalar2, Eigen::RowMajor>& dest)
213 {
214     typedef typename Eigen::SparseMatrix<Scalar1, Eigen::RowMajor>::InnerIterator srcIt_t;
215     for (std::size_t k = 0; k != src.outerSize(); ++k)
216     {
217         for (srcIt_t it(src, (int)k); it; ++it)
218         {
219             dest.insert( it.row() + r, it.col() + c) =  it.value();
220         }
221     }
222 }
223 /*
224 Sp is an Eigen::SparseMatrix
225 */
226 template<typename Sp, typename M>
227 void cwiseInPlaceProduct(Sp& sp, M SPARSE_CONST& m)
228 {
229     // should be a transform_n() over makeNonZerosIterator(src)
230     for (std::size_t k = 0; k != sp.outerSize(); ++k)
231     {
232         for (typename Sp::InnerIterator it(sp, k); it; ++it)
233         {
234             it.valueRef() *= get<typename Eigen::internal::traits<Sp>::Scalar >(m, it.row(), it.col());
235         }
236     }
237
238 }
239 }
240 namespace types
241 {
242
243 template<typename T, typename Arg>
244 T* create_new(Arg const& a)
245 {
246     return 0;
247 }
248
249 template<>
250 Double* create_new(double const& d)
251 {
252     Double* res(new Double(1, 1, false));
253     res->set(0, 0, d);
254     return res;
255 }
256
257 template<>
258 Double* create_new(std::complex<double>const& c)
259 {
260     Double* res(new Double(1, 1, true));
261     res->set(0, 0, c.real());
262     res->setImg(0, 0, c.imag());
263     return res;
264 }
265
266 template<>
267 Double* create_new(Sparse const& s)
268 {
269     Sparse& cs(const_cast<Sparse&>(s)); // inherited member functions are not const-correct
270     Double* res(new Double(cs.getRows(), cs.getCols(), cs.isComplex()));
271     const_cast<Sparse&>(s).fill(*res);
272     return res;
273 }
274
275
276 Sparse::~Sparse()
277 {
278     delete matrixReal;
279     delete matrixCplx;
280 #ifndef NDEBUG
281     Inspector::removeItem(this);
282 #endif
283 }
284
285 Sparse::Sparse(Sparse const& src)
286     : matrixReal(src.matrixReal ? new RealSparse_t(*src.matrixReal) : 0)
287     , matrixCplx(src.matrixCplx ? new CplxSparse_t(*src.matrixCplx) : 0)
288
289 {
290     m_iRows = const_cast<Sparse*>(&src)->getRows();
291     m_iCols = const_cast<Sparse*>(&src)->getCols();
292     m_iSize = m_iRows * m_iCols;
293     m_iDims = 2;
294     m_piDims[0] = m_iRows;
295     m_piDims[1] = m_iCols;
296 #ifndef NDEBUG
297     Inspector::addItem(this);
298 #endif
299 }
300
301 Sparse::Sparse(int _iRows, int _iCols, bool cplx)
302     : matrixReal(cplx ? 0 : new RealSparse_t(_iRows, _iCols))
303     , matrixCplx(cplx ? new CplxSparse_t(_iRows, _iCols) : 0)
304 {
305     m_iRows = _iRows;
306     m_iCols = _iCols;
307     m_iSize = _iRows * _iCols;
308     m_iDims = 2;
309     m_piDims[0] = _iRows;
310     m_piDims[1] = _iCols;
311 #ifndef NDEBUG
312     Inspector::addItem(this);
313 #endif
314 }
315
316 Sparse::Sparse(Double SPARSE_CONST& src)
317 {
318     //compute idx
319     int size = src.getSize();
320     int row = src.getRows();
321     Double* idx = new Double(src.getSize(), 2);
322     double* p = idx->get();
323     for (int i = 0; i < size; ++i)
324     {
325         p[i]        = (double)(i % row) + 1;
326         p[i + size] = (double)(i / row) + 1;
327     }
328     create2(src.getRows(), src.getCols(), src, *idx);
329     idx->killMe();
330 #ifndef NDEBUG
331     Inspector::addItem(this);
332 #endif
333 }
334
335 Sparse::Sparse(Double SPARSE_CONST& src, Double SPARSE_CONST& idx)
336 {
337     int idxrow = idx.getRows();
338     int rows = static_cast<int>(*std::max_element(idx.get(), idx.get() + idxrow));
339     int cols = static_cast<int>(*std::max_element(idx.get() + idxrow, idx.get() + idxrow * 2));
340
341     create2(rows, cols, src, idx);
342 #ifndef NDEBUG
343     Inspector::removeItem(this);
344 #endif
345 }
346
347 Sparse::Sparse(Double SPARSE_CONST& src, Double SPARSE_CONST& idx, Double SPARSE_CONST& dims)
348 {
349     create2(static_cast<int>(dims.get(0)), static_cast<int>(dims.get(1)), src, idx);
350 #ifndef NDEBUG
351     Inspector::addItem(this);
352 #endif
353 }
354
355 Sparse::Sparse(RealSparse_t* realSp, CplxSparse_t* cplxSp):  matrixReal(realSp), matrixCplx(cplxSp)
356 {
357     if (realSp)
358     {
359         m_iCols = realSp->cols();
360         m_iRows = realSp->rows();
361     }
362     else
363     {
364         m_iCols = cplxSp->cols();
365         m_iRows = cplxSp->rows();
366     }
367     m_iSize = m_iCols * m_iRows;
368     m_iDims = 2;
369     m_piDims[0] = m_iRows;
370     m_piDims[1] = m_iCols;
371 #ifndef NDEBUG
372     Inspector::addItem(this);
373 #endif
374 }
375
376 Sparse::Sparse(Double SPARSE_CONST& xadj, Double SPARSE_CONST& adjncy, Double SPARSE_CONST& src, std::size_t r, std::size_t c)
377 {
378     Adjacency a(xadj.get(), adjncy.get());
379     create(static_cast<int>(r), static_cast<int>(c), src, makeIteratorFromVar(a), src.getSize());
380 #ifndef NDEBUG
381     Inspector::addItem(this);
382 #endif
383 }
384
385 Sparse::Sparse(int rows, int cols, int nonzeros, int* inner, int* outer, double* real, double* img)
386 {
387     int* out = nullptr;
388     int* in = nullptr;
389
390     if (img)
391     {
392         matrixCplx = new CplxSparse_t(rows, cols);
393         matrixCplx->reserve((int)nonzeros);
394         out = matrixCplx->outerIndexPtr();
395         in = matrixCplx->innerIndexPtr();
396         matrixReal = nullptr;
397     }
398     else
399     {
400         matrixReal = new RealSparse_t(rows, cols);
401         matrixReal->reserve((int)nonzeros);
402         out = matrixReal->outerIndexPtr();
403         in = matrixReal->innerIndexPtr();
404         matrixCplx = nullptr;
405     }
406
407     //update outerIndexPtr
408     memcpy(out, outer, sizeof(int) * (rows + 1));
409     //update innerIndexPtr
410     memcpy(in, inner, sizeof(int) * nonzeros);
411
412     if (img)
413     {
414         std::complex<double>* data = matrixCplx->valuePtr();
415         for (int i = 0; i < nonzeros; ++i)
416         {
417             data[i] = std::complex<double>(real[i], img[i]);
418         }
419     }
420     else
421     {
422         double* data = matrixReal->valuePtr();
423         for (int i = 0; i < nonzeros; ++i)
424         {
425             data[i] = real[i];
426         }
427
428     }
429
430     m_iCols = cols;
431     m_iRows = rows;
432     m_iSize = cols * rows;
433     m_iDims = 2;
434     m_piDims[0] = m_iRows;
435     m_piDims[1] = m_iCols;
436
437     matrixCplx ? matrixCplx->resizeNonZeros(nonzeros) : matrixReal->resizeNonZeros(nonzeros);
438     //finalize();
439 }
440
441 template<typename DestIter>
442 void Sparse::create(int rows, int cols, Double SPARSE_CONST& src, DestIter o, std::size_t n)
443 {
444     m_iCols = cols;
445     m_iRows = rows;
446     m_iSize = cols * rows;
447     m_iDims = 2;
448     m_piDims[0] = m_iRows;
449     m_piDims[1] = m_iCols;
450
451     if (src.isComplex())
452     {
453         matrixReal = 0;
454         matrixCplx = new CplxSparse_t(rows, cols);
455         matrixCplx->reserve((int)n);
456         mycopy_n(makeMatrixIterator<std::complex<double> >(src, RowWiseFullIterator(src.getRows(), src.getCols())), n, makeMatrixIterator<std::complex<double> >(*matrixCplx, o));
457     }
458     else
459     {
460         matrixReal = new RealSparse_t(rows, cols);
461         matrixReal->reserve((int)n);
462         matrixCplx = 0;
463         mycopy_n(makeMatrixIterator<double >(src, RowWiseFullIterator(src.getRows(), src.getCols())), n
464                  , makeMatrixIterator<double>(*matrixReal, o));
465     }
466     finalize();
467 }
468
469 void Sparse::create2(int rows, int cols, Double SPARSE_CONST& src, Double SPARSE_CONST& idx)
470 {
471     int nnz = src.getSize();
472     double* i = idx.get();
473     double* j = i + idx.getRows();
474     double* valR = src.get();
475
476     if (src.isComplex())
477     {
478         matrixReal = 0;
479
480         typedef Eigen::Triplet<std::complex<double> > T;
481         std::vector<T> tripletList;
482         tripletList.reserve((int)nnz);
483
484         double* valI = src.getImg();
485
486         for (int k = 0; k < nnz; ++k)
487         {
488             tripletList.push_back(T(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, std::complex<double>(valR[k], valI[k])));
489         }
490
491         matrixCplx = new CplxSparse_t(rows, cols);
492         matrixCplx->setFromTriplets(tripletList.begin(), tripletList.end());
493         m_iRows = matrixCplx->rows();
494         m_iCols = matrixCplx->cols();
495     }
496     else
497     {
498         matrixCplx = 0;
499
500         typedef Eigen::Triplet<double> T;
501         std::vector<T> tripletList;
502         tripletList.reserve((int)nnz);
503
504         for (int k = 0; k < nnz; ++k)
505         {
506             tripletList.push_back(T(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, valR[k]));
507         }
508
509         matrixReal = new RealSparse_t(rows, cols);
510         matrixReal->setFromTriplets(tripletList.begin(), tripletList.end());
511
512         m_iRows = matrixReal->rows();
513         m_iCols = matrixReal->cols();
514     }
515
516     m_iSize = m_iCols * m_iRows;
517     m_iDims = 2;
518     m_piDims[0] = m_iRows;
519     m_piDims[1] = m_iCols;
520     finalize();
521 }
522
523 void Sparse::fill(Double& dest, int r, int c) SPARSE_CONST
524 {
525     Sparse & cthis(const_cast<Sparse&>(*this));
526     if (isComplex())
527     {
528         mycopy_n(makeMatrixIterator<std::complex<double> >(*matrixCplx, RowWiseFullIterator(cthis.getRows(), cthis.getCols())), cthis.getSize()
529         , makeMatrixIterator<std::complex<double> >(dest, RowWiseFullIterator(dest.getRows(), dest.getCols(), r, c)));
530     }
531     else
532     {
533         mycopy_n( makeMatrixIterator<double>(*matrixReal,  RowWiseFullIterator(cthis.getRows(), cthis.getCols())), cthis.getSize()
534         , makeMatrixIterator<double >(dest, RowWiseFullIterator(dest.getRows(), dest.getCols(), r, c)));
535     }
536 }
537
538 bool Sparse::set(int _iRows, int _iCols, std::complex<double> v, bool _bFinalize)
539 {
540     if (_iRows >= getRows() || _iCols >= getCols())
541     {
542         return false;
543     }
544
545     if (matrixReal)
546     {
547         matrixReal->coeffRef(_iRows, _iCols) = v.real();
548     }
549     else
550     {
551         matrixCplx->coeffRef(_iRows, _iCols) = v;
552     }
553
554     if (_bFinalize)
555     {
556         finalize();
557     }
558     return true;
559 }
560
561 bool Sparse::set(int _iRows, int _iCols, double _dblReal, bool _bFinalize)
562 {
563     if (_iRows >= getRows() || _iCols >= getCols())
564     {
565         return false;
566     }
567
568     if (matrixReal)
569     {
570         matrixReal->coeffRef(_iRows, _iCols) = _dblReal;
571     }
572     else
573     {
574         matrixCplx->coeffRef(_iRows, _iCols) = std::complex<double>(_dblReal, 0);
575     }
576
577
578     if (_bFinalize)
579     {
580         finalize();
581     }
582     return true;
583 }
584
585 void Sparse::finalize()
586 {
587     if (isComplex())
588     {
589         matrixCplx->prune(&keepForSparse<std::complex<double> >);
590         matrixCplx->finalize();
591     }
592     else
593     {
594         matrixReal->prune(&keepForSparse<double>);
595         matrixReal->finalize();
596     }
597
598 }
599
600 bool Sparse::neg(InternalType *& out)
601 {
602     SparseBool * _out = new SparseBool(getRows(), getCols());
603     types::neg(getRows(), getCols(), matrixReal, _out->matrixBool);
604     out = _out;
605
606     return true;
607 }
608
609
610 bool Sparse::isComplex() const
611 {
612     return static_cast<bool>(matrixCplx != NULL);
613 }
614
615 // TODO: should have both a bounds checking and a non-checking interface to elt access
616 double* Sparse::get()
617 {
618     if (isComplex() == false)
619     {
620         return matrixReal->valuePtr();
621     }
622
623     return nullptr;
624 }
625
626 double  Sparse::get(int _iRows, int _iCols) const
627 {
628     return getReal(_iRows, _iCols);
629 }
630
631 double Sparse::getReal(int _iRows, int _iCols) const
632 {
633     double res = 0;
634     if (matrixReal)
635     {
636         res = matrixReal->coeff(_iRows, _iCols);
637     }
638     else
639     {
640         res = matrixCplx->coeff(_iRows, _iCols).real();
641     }
642     return res;
643 }
644
645 std::complex<double>* Sparse::getImg()
646 {
647     if (isComplex())
648     {
649         return matrixCplx->valuePtr();
650     }
651
652     return nullptr;
653 }
654
655 std::complex<double> Sparse::getImg(int _iRows, int _iCols) const
656 {
657     std::complex<double> res;
658     if (matrixCplx)
659     {
660         res = matrixCplx->coeff(_iRows, _iCols);
661     }
662     else
663     {
664         res = std::complex<double>(matrixReal->coeff(_iRows, _iCols), 0.);
665     }
666
667     return res;
668 }
669
670 void Sparse::whoAmI() SPARSE_CONST
671 {
672     std::cout << "types::Sparse";
673 }
674
675 Sparse* Sparse::clone(void) const
676 {
677     return new Sparse(*this);
678 }
679
680 bool Sparse::zero_set()
681 {
682     if (matrixReal)
683     {
684         matrixReal->setZero();
685     }
686     else
687     {
688         matrixCplx->setZero();
689     }
690
691     return true;
692 }
693
694 // TODO: handle precision and line length
695 bool Sparse::toString(std::wostringstream& ostr) const
696 {
697     int iPrecision = ConfigVariable::getFormatSize();
698     std::wstring res;
699     if (matrixReal)
700     {
701         res = ::toString(*matrixReal, iPrecision);
702     }
703     else
704     {
705         res = ::toString(*matrixCplx, iPrecision);
706     }
707
708     ostr << res;
709     return true;
710 }
711
712 bool Sparse::resize(int _iNewRows, int _iNewCols)
713 {
714     if (_iNewRows <= getRows() && _iNewCols <= getCols())
715     {
716         //nothing to do: hence we do NOT fail
717         return true;
718     }
719
720     bool res = false;
721     try
722     {
723         if (matrixReal)
724         {
725             //item count
726             size_t iNonZeros = nonZeros();
727             RealSparse_t *newReal = new RealSparse_t(_iNewRows, _iNewCols);
728             newReal->reserve((int)iNonZeros);
729
730
731             //coords
732             int* pRows = new int[iNonZeros * 2];
733             outputRowCol(pRows);
734             int* pCols = pRows + iNonZeros;
735
736             //values
737             double* pNonZeroR = new double[iNonZeros];
738             double* pNonZeroI = new double[iNonZeros];
739             outputValues(pNonZeroR, pNonZeroI);
740
741             typedef Eigen::Triplet<double> triplet;
742             std::vector<triplet> tripletList;
743
744             for (size_t i = 0 ; i < iNonZeros ; i++)
745             {
746                 tripletList.push_back(triplet((int)pRows[i] - 1, (int)pCols[i] - 1, pNonZeroR[i]));
747             }
748
749             newReal->setFromTriplets(tripletList.begin(), tripletList.end());
750
751             delete matrixReal;
752             matrixReal = newReal;
753             delete[] pRows;
754             delete[] pNonZeroR;
755             delete[] pNonZeroI;
756         }
757         else
758         {
759             //item count
760             size_t iNonZeros = nonZeros();
761             CplxSparse_t *newCplx = new CplxSparse_t(_iNewRows, _iNewCols);
762             newCplx->reserve((int)iNonZeros);
763
764             //coords
765             int* pRows = new int[iNonZeros * 2];
766             outputRowCol(pRows);
767             int* pCols = pRows + iNonZeros;
768
769             //values
770             double* pNonZeroR = new double[iNonZeros];
771             double* pNonZeroI = new double[iNonZeros];
772             outputValues(pNonZeroR, pNonZeroI);
773
774             typedef Eigen::Triplet<std::complex<double> > triplet;
775             std::vector<triplet> tripletList;
776
777             for (size_t i = 0 ; i < iNonZeros ; i++)
778             {
779                 tripletList.push_back(triplet((int)pRows[i] - 1, (int)pCols[i] - 1, std::complex<double>(pNonZeroR[i], pNonZeroI[i])));
780             }
781
782             newCplx->setFromTriplets(tripletList.begin(), tripletList.end());
783
784
785             delete matrixCplx;
786             matrixCplx = newCplx;
787             delete[] pRows;
788             delete[] pNonZeroR;
789             delete[] pNonZeroI;
790         }
791
792         m_iRows = _iNewRows;
793         m_iCols = _iNewCols;
794         m_iSize = _iNewRows * _iNewCols;
795         m_piDims[0] = m_iRows;
796         m_piDims[1] = m_iCols;
797
798         res = true;
799     }
800     catch (...)
801     {
802         res = false;
803     }
804     return res;
805 }
806 // TODO decide if a complex matrix with 0 imag can be == to a real matrix
807 // not true for dense (cf double.cpp)
808 bool Sparse::operator==(const InternalType& it) SPARSE_CONST
809 {
810     Sparse* otherSparse = const_cast<Sparse*>(dynamic_cast<Sparse const*>(&it));/* types::GenericType is not const-correct :( */
811     Sparse & cthis (const_cast<Sparse&>(*this));
812
813     if (otherSparse == NULL)
814     {
815         return false;
816     }
817
818     if (otherSparse->getRows() != cthis.getRows())
819     {
820         return false;
821     }
822
823     if (otherSparse->getCols() != cthis.getCols())
824     {
825         return false;
826     }
827
828     if (otherSparse->isComplex() != isComplex())
829     {
830         return false;
831     }
832
833     if (isComplex())
834     {
835         return equal(*matrixCplx, *otherSparse->matrixCplx);
836     }
837     else
838     {
839         return equal(*matrixReal, *otherSparse->matrixReal);
840     }
841 }
842
843 bool Sparse::one_set()
844 {
845     if (isComplex())
846     {
847         return setNonZero(*matrixCplx);
848     }
849     else
850     {
851         return setNonZero(*matrixReal);
852     }
853 }
854
855 void Sparse::toComplex()
856 {
857     if (!isComplex())
858     {
859         try
860         {
861             matrixCplx = new CplxSparse_t(matrixReal->cast<std::complex<double> >());
862             delete matrixReal;
863             matrixReal = NULL;
864         }
865         catch (...)
866         {
867             delete matrixCplx;
868             matrixCplx = NULL;
869             throw;
870         }
871     }
872 }
873
874 InternalType* Sparse::insertNew(typed_list* _pArgs, InternalType* _pSource)
875 {
876     typed_list pArg;
877     InternalType *pOut  = NULL;
878     Sparse* pSource = _pSource->getAs<Sparse>();
879
880     int iDims           = (int)_pArgs->size();
881     int* piMaxDim       = new int[iDims];
882     int* piCountDim     = new int[iDims];
883     bool bComplex       = pSource->isComplex();
884     bool bUndefine      = false;
885
886     //evaluate each argument and replace by appropriate value and compute the count of combinations
887     int iSeqCount = checkIndexesArguments(NULL, _pArgs, &pArg, piMaxDim, piCountDim);
888     if (iSeqCount == 0)
889     {
890         //free pArg content
891         cleanIndexesArguments(_pArgs, &pArg);
892         return createEmptyDouble();
893     }
894
895     if (iSeqCount < 0)
896     {
897         iSeqCount = -iSeqCount;
898         bUndefine = true;
899     }
900
901     if (bUndefine)
902     {
903         //manage : and $ in creation by insertion
904         int iSource         = 0;
905         int *piSourceDims   = pSource->getDimsArray();
906
907         for (int i = 0 ; i < iDims ; i++)
908         {
909             if (pArg[i] == NULL)
910             {
911                 //undefine value
912                 if (pSource->isScalar())
913                 {
914                     piMaxDim[i]     = 1;
915                     piCountDim[i]   = 1;
916                 }
917                 else
918                 {
919                     piMaxDim[i]     = piSourceDims[iSource];
920                     piCountDim[i]   = piSourceDims[iSource];
921                 }
922                 iSource++;
923                 //replace pArg value by the new one
924                 pArg[i] = createDoubleVector(piMaxDim[i]);
925             }
926             //else
927             //{
928             //    piMaxDim[i] = piCountDim[i];
929             //}
930         }
931     }
932
933     //remove last dimension at size 1
934     //remove last dimension if are == 1
935     for (int i = (iDims - 1) ; i >= 2 ; i--)
936     {
937         if (piMaxDim[i] == 1)
938         {
939             iDims--;
940             pArg.pop_back();
941         }
942         else
943         {
944             break;
945         }
946     }
947
948     if (checkArgValidity(pArg) == false)
949     {
950         //free pArg content
951         cleanIndexesArguments(_pArgs, &pArg);
952         //contain bad index, like <= 0, ...
953         return NULL;
954     }
955
956     if (iDims == 1)
957     {
958         if (pSource->getCols() == 1)
959         {
960             pOut = new Sparse(piCountDim[0], 1, bComplex);
961         }
962         else
963         {
964             //rows == 1
965             pOut = new Sparse(1, piCountDim[0], bComplex);
966         }
967     }
968     else
969     {
970         pOut = new Sparse(piMaxDim[0], piMaxDim[1], bComplex);
971         //pOut = pSource->createEmpty(iDims, piMaxDim, bComplex);
972     }
973
974     //fill with null item
975     Sparse* pSpOut = pOut->getAs<Sparse>();
976
977     //insert values in new matrix
978     InternalType* pOut2 = pSpOut->insert(&pArg, pSource);
979     if (pOut != pOut2)
980     {
981         delete pOut;
982     }
983
984     //free pArg content
985     cleanIndexesArguments(_pArgs, &pArg);
986
987     return pOut2;
988 }
989
990 Sparse* Sparse::insert(typed_list* _pArgs, InternalType* _pSource)
991 {
992     bool bNeedToResize  = false;
993     int iDims           = (int)_pArgs->size();
994     if (iDims > 2)
995     {
996         //sparse are only in 2 dims
997         return NULL;
998     }
999
1000     typed_list pArg;
1001
1002     int piMaxDim[2];
1003     int piCountDim[2];
1004
1005     //on case of resize
1006     int iNewRows    = 0;
1007     int iNewCols    = 0;
1008     Double* pSource = _pSource->getAs<Double>();
1009
1010     //evaluate each argument and replace by appropriate value and compute the count of combinations
1011     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1012     if (iSeqCount == 0)
1013     {
1014         //free pArg content
1015         cleanIndexesArguments(_pArgs, &pArg);
1016         return this;
1017     }
1018
1019     if (iDims < 2)
1020     {
1021         //see as vector
1022         if (getRows() == 1 || getCols() == 1)
1023         {
1024             //vector or scalar
1025             if (getSize() < piMaxDim[0])
1026             {
1027                 bNeedToResize = true;
1028
1029                 //need to enlarge sparse dimensions
1030                 if (getCols() == 1 || getSize() == 0)
1031                 {
1032                     //column vector
1033                     iNewRows    = piMaxDim[0];
1034                     iNewCols    = 1;
1035                 }
1036                 else if (getRows() == 1)
1037                 {
1038                     //row vector
1039                     iNewRows    = 1;
1040                     iNewCols    = piMaxDim[0];
1041                 }
1042             }
1043         }
1044         else if (getSize() < piMaxDim[0])
1045         {
1046             //free pArg content
1047             cleanIndexesArguments(_pArgs, &pArg);
1048             //out of range
1049             return NULL;
1050         }
1051     }
1052     else
1053     {
1054         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
1055         {
1056             bNeedToResize = true;
1057             iNewRows = std::max(getRows(), piMaxDim[0]);
1058             iNewCols = std::max(getCols(), piMaxDim[1]);
1059         }
1060     }
1061
1062     //check number of insertion
1063     if (pSource->isScalar() == false && pSource->getSize() != iSeqCount)
1064     {
1065         //free pArg content
1066         cleanIndexesArguments(_pArgs, &pArg);
1067         return NULL;
1068     }
1069
1070     //now you are sure to be able to insert values
1071     if (bNeedToResize)
1072     {
1073         if (resize(iNewRows, iNewCols) == false)
1074         {
1075             //free pArg content
1076             cleanIndexesArguments(_pArgs, &pArg);
1077             return NULL;
1078         }
1079     }
1080
1081     //update complexity
1082     if (pSource->isComplex() && isComplex() == false)
1083     {
1084         toComplex();
1085     }
1086
1087
1088     if (iDims == 1)
1089     {
1090         double* pIdx = pArg[0]->getAs<Double>()->get();
1091         for (int i = 0 ; i < iSeqCount ; i++)
1092         {
1093             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
1094             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
1095             if (pSource->isScalar())
1096             {
1097                 if (pSource->isComplex())
1098                 {
1099                     set(iRow, iCol, std::complex<double>(pSource->get(0), pSource->getImg(0)), false);
1100                 }
1101                 else
1102                 {
1103                     set(iRow, iCol, pSource->get(0), false);
1104                 }
1105             }
1106             else
1107             {
1108                 if (pSource->isComplex())
1109                 {
1110                     set(iRow, iCol, std::complex<double>(pSource->get(i), pSource->getImg(i)), false);
1111                 }
1112                 else
1113                 {
1114                     set(iRow, iCol, pSource->get(i), false);
1115                 }
1116             }
1117         }
1118     }
1119     else
1120     {
1121         double* pIdxRow = pArg[0]->getAs<Double>()->get();
1122         int iRowSize    = pArg[0]->getAs<Double>()->getSize();
1123         double* pIdxCol = pArg[1]->getAs<Double>()->get();
1124
1125         for (int i = 0 ; i < iSeqCount ; i++)
1126         {
1127             if (pSource->isScalar())
1128             {
1129                 if (pSource->isComplex())
1130                 {
1131                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, std::complex<double>(pSource->get(0), pSource->getImg(0)), false);
1132                 }
1133                 else
1134                 {
1135                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(0), false);
1136                 }
1137             }
1138             else
1139             {
1140                 int iRowOrig = i % pSource->getRows();
1141                 int iColOrig = i / pSource->getRows();
1142
1143                 if (pSource->isComplex())
1144                 {
1145                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, std::complex<double>(pSource->get(iRowOrig, iColOrig), pSource->getImg(iRowOrig, iColOrig)), false);
1146                 }
1147                 else
1148                 {
1149                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(iRowOrig, iColOrig), false);
1150                 }
1151             }
1152         }
1153     }
1154
1155     finalize();
1156
1157     //free pArg content
1158     cleanIndexesArguments(_pArgs, &pArg);
1159
1160     return this;
1161 }
1162
1163 Sparse* Sparse::insert(typed_list* _pArgs, Sparse* _pSource)
1164 {
1165     bool bNeedToResize  = false;
1166     int iDims           = (int)_pArgs->size();
1167     if (iDims > 2)
1168     {
1169         //sparse are only in 2 dims
1170         return NULL;
1171     }
1172
1173     typed_list pArg;
1174
1175     int piMaxDim[2];
1176     int piCountDim[2];
1177
1178     //on case of resize
1179     int iNewRows    = 0;
1180     int iNewCols    = 0;
1181
1182     //evaluate each argument and replace by appropriate value and compute the count of combinations
1183     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1184     if (iSeqCount == 0)
1185     {
1186         //free pArg content
1187         cleanIndexesArguments(_pArgs, &pArg);
1188         return this;
1189     }
1190
1191     if (iDims < 2)
1192     {
1193         //see as vector
1194         if (getRows() == 1 || getCols() == 1)
1195         {
1196             //vector or scalar
1197             bNeedToResize = true;
1198             if (getSize() < piMaxDim[0])
1199             {
1200                 //need to enlarge sparse dimensions
1201                 if (getCols() == 1 || getSize() == 0)
1202                 {
1203                     //column vector
1204                     iNewRows    = piMaxDim[0];
1205                     iNewCols    = 1;
1206                 }
1207                 else if (getRows() == 1)
1208                 {
1209                     //row vector
1210                     iNewRows    = 1;
1211                     iNewCols    = piMaxDim[0];
1212                 }
1213             }
1214         }
1215         else if (getSize() < piMaxDim[0])
1216         {
1217             //free pArg content
1218             cleanIndexesArguments(_pArgs, &pArg);
1219             //out of range
1220             return NULL;
1221         }
1222     }
1223     else
1224     {
1225         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
1226         {
1227             bNeedToResize = true;
1228             iNewRows = std::max(getRows(), piMaxDim[0]);
1229             iNewCols = std::max(getCols(), piMaxDim[1]);
1230         }
1231     }
1232
1233     //check number of insertion
1234     if (_pSource->isScalar() == false && _pSource->getSize() != iSeqCount)
1235     {
1236         //free pArg content
1237         cleanIndexesArguments(_pArgs, &pArg);
1238         return NULL;
1239     }
1240
1241     //now you are sure to be able to insert values
1242     if (bNeedToResize)
1243     {
1244         if (resize(iNewRows, iNewCols) == false)
1245         {
1246             //free pArg content
1247             cleanIndexesArguments(_pArgs, &pArg);
1248             return NULL;
1249         }
1250     }
1251
1252     //update complexity
1253     if (_pSource->isComplex() && isComplex() == false)
1254     {
1255         toComplex();
1256     }
1257
1258     if (iDims == 1)
1259     {
1260         double* pIdx = pArg[0]->getAs<Double>()->get();
1261         for (int i = 0 ; i < iSeqCount ; i++)
1262         {
1263             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
1264             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
1265
1266             if (_pSource->isScalar())
1267             {
1268                 if (_pSource->isComplex())
1269                 {
1270                     set(iRow, iCol, _pSource->getImg(0, 0), false);
1271                 }
1272                 else
1273                 {
1274                     set(iRow, iCol, _pSource->get(0, 0), false);
1275                 }
1276             }
1277             else
1278             {
1279                 int iRowOrig = i % _pSource->getRows();
1280                 int iColOrig = i / _pSource->getRows();
1281                 if (_pSource->isComplex())
1282                 {
1283                     set(iRow, iCol, _pSource->getImg(iRowOrig, iColOrig), false);
1284                 }
1285                 else
1286                 {
1287                     set(iRow, iCol, _pSource->get(iRowOrig, iColOrig), false);
1288                 }
1289             }
1290         }
1291     }
1292     else
1293     {
1294         double* pIdxRow = pArg[0]->getAs<Double>()->get();
1295         int iRowSize    = pArg[0]->getAs<Double>()->getSize();
1296         double* pIdxCol = pArg[1]->getAs<Double>()->get();
1297
1298         for (int i = 0 ; i < iSeqCount ; i++)
1299         {
1300             if (_pSource->isScalar())
1301             {
1302                 if (_pSource->isComplex())
1303                 {
1304                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->getImg(0, 0), false);
1305                 }
1306                 else
1307                 {
1308                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(0, 0), false);
1309                 }
1310             }
1311             else
1312             {
1313                 int iRowOrig = i % _pSource->getRows();
1314                 int iColOrig = i / _pSource->getRows();
1315                 if (_pSource->isComplex())
1316                 {
1317                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->getImg(iRowOrig, iColOrig), false);
1318                 }
1319                 else
1320                 {
1321                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(iRowOrig, iColOrig), false);
1322                 }
1323             }
1324         }
1325     }
1326
1327     finalize();
1328
1329     //free pArg content
1330     cleanIndexesArguments(_pArgs, &pArg);
1331
1332     return this;
1333 }
1334
1335 Sparse* Sparse::remove(typed_list* _pArgs)
1336 {
1337     Sparse* pOut = NULL;
1338     int iDims = (int)_pArgs->size();
1339     if (iDims > 2)
1340     {
1341         //sparse are only in 2 dims
1342         return NULL;
1343     }
1344
1345     typed_list pArg;
1346
1347     int piMaxDim[2];
1348     int piCountDim[2];
1349
1350     //evaluate each argument and replace by appropriate value and compute the count of combinations
1351     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1352     if (iSeqCount == 0)
1353     {
1354         //free pArg content
1355         cleanIndexesArguments(_pArgs, &pArg);
1356         return this;
1357     }
1358
1359     bool* pbFull = new bool[iDims];
1360     //coord must represent all values on a dimension
1361     for (int i = 0 ; i < iDims ; i++)
1362     {
1363         pbFull[i]       = false;
1364         int iDimToCheck = getVarMaxDim(i, iDims);
1365         int iIndexSize  = pArg[i]->getAs<GenericType>()->getSize();
1366
1367         //we can have index more than once
1368         if (iIndexSize >= iDimToCheck)
1369         {
1370             //size is good, now check datas
1371             double* pIndexes = getDoubleArrayFromDouble(pArg[i]);
1372             for (int j = 0 ; j < iDimToCheck ; j++)
1373             {
1374                 bool bFind = false;
1375                 for (int k = 0 ; k < iIndexSize ; k++)
1376                 {
1377                     if ((int)pIndexes[k] == j + 1)
1378                     {
1379                         bFind = true;
1380                         break;
1381                     }
1382                 }
1383                 pbFull[i]  = bFind;
1384             }
1385         }
1386     }
1387
1388     //only one dims can be not full/entire
1389     bool bNotEntire = false;
1390     int iNotEntire  = 0;
1391     bool bTooMuchNotEntire = false;
1392     for (int i = 0 ; i < iDims ; i++)
1393     {
1394         if (pbFull[i] == false)
1395         {
1396             if (bNotEntire == false)
1397             {
1398                 bNotEntire = true;
1399                 iNotEntire = i;
1400             }
1401             else
1402             {
1403                 bTooMuchNotEntire = true;
1404                 break;
1405             }
1406         }
1407     }
1408
1409     if (bTooMuchNotEntire == true)
1410     {
1411         //free pArg content
1412         cleanIndexesArguments(_pArgs, &pArg);
1413         return NULL;
1414     }
1415
1416     delete[] pbFull;
1417
1418     //find index to keep
1419     int iNotEntireSize          = pArg[iNotEntire]->getAs<GenericType>()->getSize();
1420     double* piNotEntireIndex    = getDoubleArrayFromDouble(pArg[iNotEntire]);
1421     int iKeepSize               = getVarMaxDim(iNotEntire, iDims);
1422     bool* pbKeep                = new bool[iKeepSize];
1423
1424     //fill pbKeep with true value
1425     for (int i = 0 ; i < iKeepSize ; i++)
1426     {
1427         pbKeep[i] = true;
1428     }
1429
1430     for (int i = 0 ; i < iNotEntireSize ; i++)
1431     {
1432         int idx = (int)piNotEntireIndex[i] - 1;
1433
1434         //don't care of value out of bounds
1435         if (idx < iKeepSize)
1436         {
1437             pbKeep[idx] = false;
1438         }
1439     }
1440
1441     int iNewDimSize = 0;
1442     for (int i = 0 ; i < iKeepSize ; i++)
1443     {
1444         if (pbKeep[i] == true)
1445         {
1446             iNewDimSize++;
1447         }
1448     }
1449     delete[] pbKeep;
1450
1451     int* piNewDims = new int[iDims];
1452     for (int i = 0 ; i < iDims ; i++)
1453     {
1454         if (i == iNotEntire)
1455         {
1456             piNewDims[i] = iNewDimSize;
1457         }
1458         else
1459         {
1460             piNewDims[i] = getVarMaxDim(i, iDims);
1461         }
1462     }
1463
1464     //remove last dimension if are == 1
1465     int iOrigDims = iDims;
1466     for (int i = (iDims - 1) ; i >= 2 ; i--)
1467     {
1468         if (piNewDims[i] == 1)
1469         {
1470             iDims--;
1471         }
1472         else
1473         {
1474             break;
1475         }
1476     }
1477
1478     if (iDims == 1)
1479     {
1480         if (iNewDimSize == 0)
1481         {
1482             delete[] piNewDims;
1483             //free pArg content
1484             cleanIndexesArguments(_pArgs, &pArg);
1485             return new Sparse(0, 0);
1486         }
1487         else
1488         {
1489             //two cases, depends of original matrix/vector
1490             if ((*_pArgs)[0]->isColon() == false && m_iDims == 2 && m_piDims[0] == 1 && m_piDims[1] != 1)
1491             {
1492                 //special case for row vector
1493                 pOut = new Sparse(1, iNewDimSize, isComplex());
1494                 //in this case we have to care of 2nd dimension
1495                 //iNotEntire = 1;
1496             }
1497             else
1498             {
1499                 pOut = new Sparse(iNewDimSize, 1, isComplex());
1500             }
1501         }
1502     }
1503     else
1504     {
1505         pOut = new Sparse(piNewDims[0], piNewDims[0], isComplex());
1506     }
1507
1508     delete[] piNewDims;
1509     //find a way to copy existing data to new variable ...
1510     int iNewPos = 0;
1511     int* piIndexes = new int[iOrigDims];
1512     int* piViewDims = new int[iOrigDims];
1513     for (int i = 0 ; i < iOrigDims ; i++)
1514     {
1515         piViewDims[i] = getVarMaxDim(i, iOrigDims);
1516     }
1517
1518     for (int i = 0 ; i < getSize() ; i++)
1519     {
1520         bool bByPass = false;
1521         getIndexesWithDims(i, piIndexes, piViewDims, iOrigDims);
1522
1523         //check if piIndexes use removed indexes
1524         for (int j = 0 ; j < iNotEntireSize ; j++)
1525         {
1526             if ((piNotEntireIndex[j] - 1) == piIndexes[iNotEntire])
1527             {
1528                 //by pass this value
1529                 bByPass = true;
1530                 break;
1531             }
1532         }
1533
1534         if (bByPass == false)
1535         {
1536             //compute new index
1537             if (isComplex())
1538             {
1539                 pOut->set(iNewPos, getImg(i));
1540             }
1541             else
1542             {
1543                 pOut->set(iNewPos, get(i));
1544             }
1545             iNewPos++;
1546         }
1547     }
1548
1549     //free allocated data
1550     for (int i = 0 ; i < iDims ; i++)
1551     {
1552         if (pArg[i] != (*_pArgs)[i])
1553         {
1554             delete pArg[i];
1555         }
1556     }
1557
1558     delete[] piIndexes;
1559     delete[] piViewDims;
1560
1561     //free pArg content
1562     cleanIndexesArguments(_pArgs, &pArg);
1563
1564     return pOut;
1565 }
1566
1567 bool Sparse::append(int r, int c, types::Sparse SPARSE_CONST* src)
1568 {
1569     //        std::wcerr << L"to a sparse of size"<<getRows() << L","<<getCols() << L" should append @"<<r << L","<<c<< "a sparse:"<< src->toString(32,80)<<std::endl;
1570     if (src->isComplex())
1571     {
1572         toComplex();
1573     }
1574     if (isComplex())
1575     {
1576         if (src->isComplex())
1577         {
1578             doAppend(*(src->matrixCplx), r, c, *matrixCplx);
1579         }
1580         else
1581         {
1582             doAppend(*(src->matrixReal), r, c, *matrixCplx);
1583         }
1584     }
1585     else
1586     {
1587         doAppend(*(src->matrixReal), r, c, *matrixReal);
1588     }
1589
1590     finalize();
1591
1592     return true; // realloc is meaningless for sparse matrices
1593 }
1594
1595 /*
1596 * create a new Sparse of dims according to resSize and fill it from currentSparse (along coords)
1597 */
1598 InternalType* Sparse::extract(typed_list* _pArgs)
1599 {
1600     Sparse* pOut        = NULL;
1601     int iDims           = (int)_pArgs->size();
1602     typed_list pArg;
1603
1604     int* piMaxDim       = new int[iDims];
1605     int* piCountDim     = new int[iDims];
1606
1607     //evaluate each argument and replace by appropriate value and compute the count of combinations
1608     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1609     if (iSeqCount == 0)
1610     {
1611         //free pArg content
1612         cleanIndexesArguments(_pArgs, &pArg);
1613         if (_pArgs->size() == 0)
1614         {
1615             //a()
1616             delete[] piMaxDim;
1617             delete[] piCountDim;
1618             //free pArg content
1619             cleanIndexesArguments(_pArgs, &pArg);
1620             return this;
1621         }
1622         else
1623         {
1624             //a([])
1625             delete[] piMaxDim;
1626             delete[] piCountDim;
1627             //free pArg content
1628             cleanIndexesArguments(_pArgs, &pArg);
1629             return Double::Empty();
1630         }
1631     }
1632
1633     if (iDims < 2)
1634     {
1635         if (piMaxDim[0] <= getSize())
1636         {
1637             int iNewRows = 0;
1638             int iNewCols = 0;
1639
1640             if (getRows() == 1 && getCols() != 1 && (*_pArgs)[0]->isColon() == false)
1641             {
1642                 //special case for row vector
1643                 iNewRows = 1;
1644                 iNewCols = piCountDim[0];
1645             }
1646             else
1647             {
1648                 iNewRows = piCountDim[0];
1649                 iNewCols = 1;
1650             }
1651
1652             pOut = new Sparse(iNewRows, iNewCols, isComplex());
1653             double* pIdx = pArg[0]->getAs<Double>()->get();
1654             for (int i = 0 ; i < iSeqCount ; i++)
1655             {
1656                 if (pIdx[i] < 1)
1657                 {
1658                     delete pOut;
1659                     pOut = NULL;
1660                     break;
1661                 }
1662                 int iRowRead = static_cast<int>(pIdx[i] - 1) % getRows();
1663                 int iColRead = static_cast<int>(pIdx[i] - 1) / getRows();
1664
1665                 int iRowWrite = static_cast<int>(i) % iNewRows;
1666                 int iColWrite = static_cast<int>(i) / iNewRows;
1667                 if (isComplex())
1668                 {
1669                     std::complex<double> dbl = getImg(iRowRead, iColRead);
1670                     if (dbl.real() != 0 || dbl.imag() != 0)
1671                     {
1672                         //only non zero values
1673                         pOut->set(iRowWrite, iColWrite, dbl, false);
1674                     }
1675                 }
1676                 else
1677                 {
1678                     double dbl = get(iRowRead, iColRead);
1679                     if (dbl != 0)
1680                     {
1681                         //only non zero values
1682                         pOut->set(iRowWrite, iColWrite, dbl, false);
1683                     }
1684                 }
1685             }
1686         }
1687         else
1688         {
1689             delete[] piMaxDim;
1690             delete[] piCountDim;
1691             //free pArg content
1692             cleanIndexesArguments(_pArgs, &pArg);
1693             return NULL;
1694         }
1695     }
1696     else
1697     {
1698         if (piMaxDim[0] <= getRows() && piMaxDim[1] <= getCols())
1699         {
1700             double* pIdxRow = pArg[0]->getAs<Double>()->get();
1701             double* pIdxCol = pArg[1]->getAs<Double>()->get();
1702
1703             int iNewRows = pArg[0]->getAs<Double>()->getSize();
1704             int iNewCols = pArg[1]->getAs<Double>()->getSize();
1705
1706             pOut = new Sparse(iNewRows, iNewCols, isComplex());
1707
1708             int iPos = 0;
1709             for (int iRow = 0 ; iRow < iNewRows ; iRow++)
1710             {
1711                 for (int iCol = 0 ; iCol < iNewCols ; iCol++)
1712                 {
1713                     if ((pIdxRow[iRow] < 1) || (pIdxCol[iCol] < 1))
1714                     {
1715                         delete pOut;
1716                         pOut = NULL;
1717                         break;
1718                     }
1719                     if (isComplex())
1720                     {
1721                         std::complex<double> dbl = getImg((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
1722                         if (dbl.real() != 0 || dbl.imag() != 0)
1723                         {
1724                             //only non zero values
1725                             pOut->set(iRow, iCol, dbl, false);
1726                         }
1727                     }
1728                     else
1729                     {
1730                         double dbl = get((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
1731                         if (dbl != 0)
1732                         {
1733                             //only non zero values
1734                             pOut->set(iRow, iCol, dbl, false);
1735                         }
1736                     }
1737                     iPos++;
1738                 }
1739             }
1740         }
1741         else
1742         {
1743             delete[] piMaxDim;
1744             delete[] piCountDim;
1745             //free pArg content
1746             cleanIndexesArguments(_pArgs, &pArg);
1747             return NULL;
1748         }
1749     }
1750
1751     pOut->finalize();
1752
1753     delete[] piMaxDim;
1754     delete[] piCountDim;
1755     //free pArg content
1756     cleanIndexesArguments(_pArgs, &pArg);
1757
1758     return pOut;
1759 }
1760
1761 Sparse* Sparse::extract(int nbCoords, int SPARSE_CONST* coords, int SPARSE_CONST* maxCoords, int SPARSE_CONST* resSize, bool asVector) SPARSE_CONST
1762 {
1763     if ( (asVector && maxCoords[0] > getSize()) ||
1764     (asVector == false && maxCoords[0] > getRows()) ||
1765     (asVector == false && maxCoords[1] > getCols()))
1766     {
1767         return 0;
1768     }
1769
1770     bool const cplx(isComplex());
1771     Sparse * pSp (0);
1772     if (asVector)
1773     {
1774         pSp = (getRows() == 1) ?  new Sparse(1, resSize[0], cplx) : new Sparse(resSize[0], 1, cplx);
1775     }
1776     else
1777     {
1778         pSp = new Sparse(resSize[0], resSize[1], cplx);
1779     }
1780     //        std::cerr<<"extracted sparse:"<<pSp->getRows()<<", "<<pSp->getCols()<<"seqCount="<<nbCoords<<"maxDim="<<maxCoords[0] <<","<< maxCoords[1]<<std::endl;
1781     if (! (asVector
1782     ? copyToSparse(*this,  Coords<true>(coords, getRows()), nbCoords
1783     , *pSp, RowWiseFullIterator(pSp->getRows(), pSp->getCols()))
1784     : copyToSparse(*this,  Coords<false>(coords), nbCoords
1785     , *pSp, RowWiseFullIterator(pSp->getRows(), pSp->getCols()))))
1786     {
1787         delete pSp;
1788         pSp = NULL;
1789     }
1790     return pSp;
1791 }
1792 /*
1793 coords are Scilab 1-based
1794 extract std::make_pair(coords, asVector), rowIter
1795 */
1796 template<typename Src, typename SrcTraversal, typename Sz, typename DestTraversal>
1797 bool Sparse::copyToSparse(Src SPARSE_CONST& src, SrcTraversal srcTrav, Sz n, Sparse& sp, DestTraversal destTrav)
1798 {
1799     if (!(src.isComplex() || sp.isComplex()))
1800     {
1801         mycopy_n(makeMatrixIterator<double>(src, srcTrav), n
1802                  , makeMatrixIterator<double>(*sp.matrixReal, destTrav));
1803     }
1804     else
1805     {
1806         sp.toComplex();
1807         mycopy_n(makeMatrixIterator<std::complex<double> >(src, srcTrav), n
1808                  , makeMatrixIterator<std::complex<double> >(*sp.matrixCplx, destTrav));
1809     }
1810
1811     sp.finalize();
1812     return true;
1813 }
1814
1815 // GenericType because we might return a Double* for scalar operand
1816 Sparse* Sparse::add(Sparse const& o) const
1817 {
1818     RealSparse_t* realSp(0);
1819     CplxSparse_t* cplxSp(0);
1820     if (isComplex() == false && o.isComplex() == false)
1821     {
1822         //R + R -> R
1823         realSp = new RealSparse_t(*matrixReal + * (o.matrixReal));
1824     }
1825     else if (isComplex() == false && o.isComplex() == true)
1826     {
1827         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >() + * (o.matrixCplx));
1828     }
1829     else if (isComplex() == true && o.isComplex() == false)
1830     {
1831         cplxSp = new CplxSparse_t(*matrixCplx + o.matrixReal->cast<std::complex<double> >());
1832     }
1833     else if (isComplex() == true && o.isComplex() == true)
1834     {
1835         cplxSp = new CplxSparse_t(*matrixCplx + * (o.matrixCplx));
1836     }
1837
1838     return new Sparse(realSp, cplxSp);
1839 }
1840
1841 Sparse* Sparse::substract(Sparse const& o) const
1842 {
1843     RealSparse_t* realSp(0);
1844     CplxSparse_t* cplxSp(0);
1845     if (isComplex() == false && o.isComplex() == false)
1846     {
1847         //R - R -> R
1848         realSp = new RealSparse_t(*matrixReal - * (o.matrixReal));
1849     }
1850     else if (isComplex() == false && o.isComplex() == true)
1851     {
1852         //R - C -> C
1853         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >() - * (o.matrixCplx));
1854     }
1855     else if (isComplex() == true && o.isComplex() == false)
1856     {
1857         //C - R -> C
1858         cplxSp = new CplxSparse_t(*matrixCplx - o.matrixReal->cast<std::complex<double> >());
1859     }
1860     else if (isComplex() == true && o.isComplex() == true)
1861     {
1862         //C - C -> C
1863         cplxSp = new CplxSparse_t(*matrixCplx - * (o.matrixCplx));
1864     }
1865
1866     return new Sparse(realSp, cplxSp);
1867 }
1868
1869 Sparse* Sparse::multiply(double s) const
1870 {
1871     return new Sparse( isComplex() ? 0 : new RealSparse_t((*matrixReal)*s)
1872                        , isComplex() ? new CplxSparse_t((*matrixCplx)*s) : 0);
1873 }
1874
1875 Sparse* Sparse::multiply(std::complex<double> s) const
1876 {
1877     return new Sparse( 0
1878                        , isComplex() ? new CplxSparse_t((*matrixCplx) * s) : new CplxSparse_t((*matrixReal) * s));
1879 }
1880
1881 Sparse* Sparse::multiply(Sparse const& o) const
1882 {
1883     RealSparse_t* realSp(0);
1884     CplxSparse_t* cplxSp(0);
1885
1886     if (isComplex() == false && o.isComplex() == false)
1887     {
1888         realSp = new RealSparse_t(*matrixReal **(o.matrixReal));
1889     }
1890     else if (isComplex() == false && o.isComplex() == true)
1891     {
1892         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >() **(o.matrixCplx));
1893     }
1894     else if (isComplex() == true && o.isComplex() == false)
1895     {
1896         cplxSp = new CplxSparse_t(*matrixCplx * o.matrixReal->cast<std::complex<double> >());
1897     }
1898     else if (isComplex() == true && o.isComplex() == true)
1899     {
1900         cplxSp = new CplxSparse_t(*matrixCplx **(o.matrixCplx));
1901     }
1902
1903     return new Sparse(realSp, cplxSp);
1904 }
1905
1906 Sparse* Sparse::dotMultiply(Sparse SPARSE_CONST& o) const
1907 {
1908     RealSparse_t* realSp(0);
1909     CplxSparse_t* cplxSp(0);
1910     if (isComplex() == false && o.isComplex() == false)
1911     {
1912         realSp = new RealSparse_t(matrixReal->cwiseProduct(*(o.matrixReal)));
1913     }
1914     else if (isComplex() == false && o.isComplex() == true)
1915     {
1916         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >().cwiseProduct( *(o.matrixCplx)));
1917     }
1918     else if (isComplex() == true && o.isComplex() == false)
1919     {
1920         cplxSp = new CplxSparse_t(matrixCplx->cwiseProduct(o.matrixReal->cast<std::complex<double> >()));
1921     }
1922     else if (isComplex() == true && o.isComplex() == true)
1923     {
1924         cplxSp = new CplxSparse_t(matrixCplx->cwiseProduct(*(o.matrixCplx)));
1925     }
1926
1927     return new Sparse(realSp, cplxSp);
1928 }
1929
1930 Sparse* Sparse::dotDivide(Sparse SPARSE_CONST& o) const
1931 {
1932     RealSparse_t* realSp(0);
1933     CplxSparse_t* cplxSp(0);
1934     if (isComplex() == false && o.isComplex() == false)
1935     {
1936         realSp = new RealSparse_t(matrixReal->cwiseQuotient(*(o.matrixReal)));
1937     }
1938     else if (isComplex() == false && o.isComplex() == true)
1939     {
1940         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >().cwiseQuotient( *(o.matrixCplx)));
1941     }
1942     else if (isComplex() == true && o.isComplex() == false)
1943     {
1944         cplxSp = new CplxSparse_t(matrixCplx->cwiseQuotient(o.matrixReal->cast<std::complex<double> >()));
1945     }
1946     else if (isComplex() == true && o.isComplex() == true)
1947     {
1948         cplxSp = new CplxSparse_t(matrixCplx->cwiseQuotient(*(o.matrixCplx)));
1949     }
1950
1951     return new Sparse(realSp, cplxSp);
1952 }
1953
1954 int Sparse::newCholLLT(Sparse** _SpPermut, Sparse** _SpFactor) const
1955 {
1956     typedef Eigen::SparseMatrix<double, Eigen::ColMajor> RealSparseCol_t;
1957     RealSparseCol_t spColMajor = RealSparseCol_t((const RealSparse_t&) * matrixReal);
1958
1959     // Constructs and performs the LLT factorization of sparse
1960     Eigen::SimplicialLLT<RealSparseCol_t> pLLT(spColMajor);
1961     int iInfo = pLLT.info();
1962     if (iInfo != Eigen::Success)
1963     {
1964         *_SpFactor = NULL;
1965         *_SpPermut = NULL;
1966         return iInfo;
1967     }
1968
1969     // Get the lower matrix of factorization.
1970     // The new RealSparse_t will be setted in Sparse without copy.
1971     *_SpFactor = new Sparse(new RealSparse_t(pLLT.matrixL()), NULL);
1972
1973     // Get the permutation matrix.
1974     Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic, int> p = pLLT.permutationP();
1975     *_SpPermut = new Sparse(p.rows(), p.cols());
1976     for (int i = 0; i < p.rows(); i++)
1977     {
1978         (*_SpPermut)->set(i, p.indices()[i], 1, false);
1979     }
1980
1981     (*_SpPermut)->finalize();
1982
1983     return iInfo;
1984 }
1985
1986 bool Sparse::transpose(InternalType *& out)
1987 {
1988     out = new Sparse(matrixReal ? new RealSparse_t(matrixReal->transpose()) : 0, matrixCplx ? new CplxSparse_t(matrixCplx->transpose()) : 0);
1989     return true;
1990 }
1991
1992 bool Sparse::adjoint(InternalType *& out)
1993 {
1994     out = new Sparse(matrixReal ? new RealSparse_t(matrixReal->adjoint()) : 0, matrixCplx ? new CplxSparse_t(matrixCplx->adjoint()) : 0);
1995     return true;
1996 }
1997
1998 struct BoolCast
1999 {
2000     BoolCast(std::complex<double> const& c): b(c.real() || c.imag()) {}
2001     operator bool () const
2002     {
2003         return b;
2004     }
2005     operator double() const
2006     {
2007         return b ? 1. : 0.;
2008     }
2009     bool b;
2010 };
2011 Sparse* Sparse::newOnes() const
2012 {
2013     // result is never cplx
2014     return new Sparse( matrixReal
2015                        ? new RealSparse_t(matrixReal->cast<bool>().cast<double>())
2016                        : new RealSparse_t(matrixCplx->cast<BoolCast>().cast<double>())
2017                        , 0);
2018 }
2019
2020 struct RealCast
2021 {
2022     RealCast(std::complex<double> const& c): b(c.real()) {}
2023     operator bool () const
2024     {
2025         return b != 0;
2026     }
2027     operator double() const
2028     {
2029         return b;
2030     }
2031     double b;
2032 };
2033 Sparse* Sparse::newReal() const
2034 {
2035     return new Sparse( matrixReal
2036                        ? matrixReal
2037                        : new RealSparse_t(matrixCplx->cast<RealCast>().cast<double>())
2038                        , 0);
2039 }
2040
2041 std::size_t Sparse::nonZeros() const
2042 {
2043     if (isComplex())
2044     {
2045         return matrixCplx->nonZeros();
2046     }
2047     else
2048     {
2049         return matrixReal->nonZeros();
2050     }
2051 }
2052 std::size_t Sparse::nonZeros(std::size_t r) const
2053 {
2054     std::size_t res;
2055     if (matrixReal)
2056     {
2057         int* piIndex = matrixReal->outerIndexPtr();
2058         res = piIndex[r + 1] - piIndex[r];
2059     }
2060     else
2061     {
2062         int* piIndex = matrixCplx->outerIndexPtr();
2063         res = piIndex[r + 1] - piIndex[r];
2064     }
2065
2066     return res;
2067 }
2068
2069 int* Sparse::getNbItemByRow(int* _piNbItemByRows)
2070 {
2071     int* piNbItemByCols = new int[getRows() + 1];
2072     if (isComplex())
2073     {
2074         mycopy_n(matrixCplx->outerIndexPtr(), getRows() + 1, piNbItemByCols);
2075     }
2076     else
2077     {
2078         mycopy_n(matrixReal->outerIndexPtr(), getRows() + 1, piNbItemByCols);
2079     }
2080
2081     for (int i = 0 ; i < getRows() ; i++)
2082     {
2083         _piNbItemByRows[i] = piNbItemByCols[i + 1] - piNbItemByCols[i];
2084     }
2085
2086     delete[] piNbItemByCols;
2087     return _piNbItemByRows;
2088 }
2089
2090 int* Sparse::getColPos(int* _piColPos)
2091 {
2092     if (isComplex())
2093     {
2094         mycopy_n(matrixCplx->innerIndexPtr(), nonZeros(), _piColPos);
2095     }
2096     else
2097     {
2098         mycopy_n(matrixReal->innerIndexPtr(), nonZeros(), _piColPos);
2099     }
2100
2101     for (int i = 0; i < nonZeros(); i++)
2102     {
2103         _piColPos[i]++;
2104     }
2105
2106     return _piColPos;
2107 }
2108
2109 int* Sparse::getInnerPtr(int* count)
2110 {
2111     int* ret = nullptr;
2112     if (isComplex())
2113     {
2114         ret = matrixCplx->innerIndexPtr();
2115         *count = matrixCplx->innerSize();
2116     }
2117     else
2118     {
2119         ret = matrixReal->innerIndexPtr();
2120         *count = matrixReal->innerSize();
2121     }
2122
2123     return ret;
2124 }
2125
2126 int* Sparse::getOuterPtr(int* count)
2127 {
2128     int* ret = nullptr;
2129     if (isComplex())
2130     {
2131         ret = matrixCplx->outerIndexPtr();
2132         *count = matrixCplx->outerSize();
2133     }
2134     else
2135     {
2136         ret = matrixReal->outerIndexPtr();
2137         *count = matrixReal->outerSize();
2138     }
2139
2140     return ret;
2141 }
2142
2143 namespace
2144 {
2145 template<typename S> struct GetReal: std::unary_function<typename S::InnerIterator, double>
2146 {
2147     double operator()(typename S::InnerIterator it) const
2148     {
2149         return it.value();
2150     }
2151 };
2152 template<> struct GetReal< Eigen::SparseMatrix<std::complex<double >, Eigen::RowMajor > >
2153     : std::unary_function<Sparse::CplxSparse_t::InnerIterator, double>
2154 {
2155     double operator()( Sparse::CplxSparse_t::InnerIterator it) const
2156     {
2157         return it.value().real();
2158     }
2159 };
2160 template<typename S> struct GetImag: std::unary_function<typename S::InnerIterator, double>
2161 {
2162     double operator()(typename S::InnerIterator it) const
2163     {
2164         return it.value().imag();
2165     }
2166 };
2167 template<typename S> struct GetRow: std::unary_function<typename S::InnerIterator, int>
2168 {
2169     int operator()(typename S::InnerIterator it) const
2170     {
2171         return it.row() + 1;
2172     }
2173 };
2174 template<typename S> struct GetCol: std::unary_function<typename S::InnerIterator, int>
2175 {
2176     int operator()(typename S::InnerIterator it) const
2177     {
2178         return it.col() + 1;
2179     }
2180 };
2181
2182 template<typename S, typename Out, typename F> Out sparseTransform(S& s, Out o, F f)
2183 {
2184     for (std::size_t k(0); k < s.outerSize(); ++k)
2185     {
2186         for (typename S::InnerIterator it(s, (int)k); it; ++it, ++o)
2187         {
2188             *o = f(it);
2189         }
2190     }
2191     return o;
2192 }
2193 }
2194
2195 std::pair<double*, double*> Sparse::outputValues(double* outReal, double* outImag)const
2196 {
2197     return matrixReal
2198            ? std::make_pair(sparseTransform(*matrixReal, outReal, GetReal<RealSparse_t>()), outImag)
2199            : std::make_pair(sparseTransform(*matrixCplx, outReal, GetReal<CplxSparse_t>())
2200                             , sparseTransform(*matrixCplx, outImag, GetImag<CplxSparse_t>()));
2201 }
2202
2203 int* Sparse::outputRowCol(int* out)const
2204 {
2205     return matrixReal
2206            ? sparseTransform(*matrixReal, sparseTransform(*matrixReal, out, GetRow<RealSparse_t>()), GetCol<RealSparse_t>())
2207            : sparseTransform(*matrixCplx, sparseTransform(*matrixCplx, out, GetRow<CplxSparse_t>()), GetCol<CplxSparse_t>());
2208 }
2209 double* Sparse::outputCols(double* out) const
2210 {
2211     if (isComplex())
2212     {
2213         mycopy_n(matrixCplx->innerIndexPtr(), nonZeros(), out);
2214     }
2215     else
2216     {
2217         mycopy_n(matrixReal->innerIndexPtr(), nonZeros(), out);
2218     }
2219
2220     return std::transform(out, out, out, std::bind2nd(std::plus<double>(), 1));
2221
2222 }
2223
2224 void Sparse::opposite(void)
2225 {
2226     if (isComplex())
2227     {
2228         std::complex<double>* data = matrixCplx->valuePtr();
2229         std::transform(data, data + matrixCplx->nonZeros(), data, std::negate<std::complex<double> >());
2230     }
2231     else
2232     {
2233         double* data = matrixReal->valuePtr();
2234         std::transform(data, data + matrixReal->nonZeros(), data, std::negate<double>());
2235     }
2236 }
2237
2238 SparseBool* Sparse::newLessThan(Sparse const&o) const
2239 {
2240     return cwiseOp<std::less>(*this, o);
2241 }
2242
2243 SparseBool* Sparse::newGreaterThan(Sparse const&o) const
2244 {
2245     return cwiseOp<std::greater>(*this, o);
2246 }
2247
2248 SparseBool* Sparse::newNotEqualTo(Sparse const&o) const
2249 {
2250     return cwiseOp<std::not_equal_to>(*this, o);
2251 }
2252
2253 SparseBool* Sparse::newLessOrEqual(Sparse const&o) const
2254 {
2255     return cwiseOp<std::less_equal>(*this, o);
2256 }
2257
2258 SparseBool* Sparse::newGreaterOrEqual(Sparse const&o) const
2259 {
2260     return cwiseOp<std::greater_equal>(*this, o);
2261 }
2262
2263 SparseBool* Sparse::newEqualTo(Sparse &o)
2264 {
2265     int rowL = getRows();
2266     int colL = getCols();
2267
2268     int rowR = o.getRows();
2269     int colR = o.getCols();
2270     int row = std::max(rowL, rowR);
2271     int col = std::max(colL, colR);
2272
2273     //create a boolean sparse matrix with dims of sparses
2274     types::SparseBool* ret = new types::SparseBool(row, col);
2275     if (isScalar() && o.isScalar())
2276     {
2277         if (isComplex() || o.isComplex())
2278         {
2279             std::complex<double> l = getImg(0, 0);
2280             std::complex<double> r = o.getImg(0, 0);
2281             ret->set(0, 0, l == r, false);
2282         }
2283         else
2284         {
2285             double l = get(0, 0);
2286             double r = o.get(0, 0);
2287             ret->set(0, 0, l == r, false);
2288         }
2289     }
2290     else if (isScalar())
2291     {
2292         int nnzR = static_cast<int>(o.nonZeros());
2293         std::vector<int> rowcolR(nnzR * 2, 0);
2294         o.outputRowCol(rowcolR.data());
2295
2296         //compare all items of R with R[0]
2297         if (isComplex() || o.isComplex())
2298         {
2299             std::complex<double> l = getImg(0, 0);
2300             for (int i = 0; i < nnzR; ++i)
2301             {
2302                 std::complex<double> r = o.getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2303                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2304             }
2305         }
2306         else
2307         {
2308             double l = get(0, 0);
2309             for (int i = 0; i < nnzR; ++i)
2310             {
2311                 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2312                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2313             }
2314         }
2315     }
2316     else if (o.isScalar())
2317     {
2318         int nnzL = static_cast<int>(nonZeros());
2319         std::vector<int> rowcolL(nnzL * 2, 0);
2320         outputRowCol(rowcolL.data());
2321
2322         if (isComplex() || o.isComplex())
2323         {
2324             std::complex<double> r = o.getImg(0, 0);
2325             for (int i = 0; i < nnzL; ++i)
2326             {
2327                 std::complex<double> l = getImg(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2328                 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
2329             }
2330         }
2331         else
2332         {
2333             double r = get(0, 0);
2334             for (int i = 0; i < nnzL; ++i)
2335             {
2336                 double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2337                 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
2338             }
2339         }
2340     }
2341     else
2342     {
2343         int nnzR = static_cast<int>(o.nonZeros());
2344         std::vector<int> rowcolR(nnzR * 2, 0);
2345         o.outputRowCol(rowcolR.data());
2346         int nnzL = static_cast<int>(nonZeros());
2347         std::vector<int> rowcolL(nnzL * 2, 0);
2348         outputRowCol(rowcolL.data());
2349         //set all values to %t
2350         ret->setTrue(false);
2351         //set %f in each pL values
2352         for (int i = 0; i < nnzL; ++i)
2353         {
2354             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, false, false);
2355         }
2356         ret->finalize();
2357
2358         //set _pR[i] == _pL[i] for each _pR values
2359         if (isComplex() || o.isComplex())
2360         {
2361             for (int i = 0; i < nnzR; ++i)
2362             {
2363                 //get l and r following non zeros value of R
2364                 std::complex<double> l = getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2365                 std::complex<double> r = o.getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2366                 //set value following non zeros value of R
2367                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2368             }
2369         }
2370         else
2371         {
2372             for (int i = 0; i < nnzR; ++i)
2373             {
2374                 //get l and r following non zeros value of R
2375                 double l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2376                 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2377                 //set value following non zeros value of R
2378                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2379             }
2380         }
2381     }
2382
2383     ret->finalize();
2384     return ret;
2385 }
2386
2387 bool Sparse::reshape(int* _piDims, int _iDims)
2388 {
2389     bool bOk = false;
2390     int iCols = 1;
2391
2392     if (_iDims == 2)
2393     {
2394         iCols = _piDims[1];
2395     }
2396
2397     if (_iDims <= 2)
2398     {
2399         bOk = reshape(_piDims[0], iCols);
2400     }
2401
2402     return bOk;
2403 }
2404
2405 bool Sparse::reshape(int _iNewRows, int _iNewCols)
2406 {
2407     if (_iNewRows * _iNewCols != getRows() * getCols())
2408     {
2409         return false;
2410     }
2411
2412     bool res = false;
2413     try
2414     {
2415         if (matrixReal)
2416         {
2417             //item count
2418             size_t iNonZeros = nonZeros();
2419             RealSparse_t *newReal = new RealSparse_t(_iNewRows, _iNewCols);
2420             newReal->reserve((int)iNonZeros);
2421
2422             //coords
2423             int* pRows = new int[iNonZeros * 2];
2424             outputRowCol(pRows);
2425             int* pCols = pRows + iNonZeros;
2426
2427             //values
2428             double* pNonZeroR = new double[iNonZeros];
2429             double* pNonZeroI = new double[iNonZeros];
2430             outputValues(pNonZeroR, pNonZeroI);
2431
2432             typedef Eigen::Triplet<double> triplet;
2433             std::vector<triplet> tripletList;
2434
2435             for (size_t i = 0 ; i < iNonZeros ; i++)
2436             {
2437                 int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
2438                 tripletList.push_back(triplet((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), pNonZeroR[i]));
2439             }
2440
2441             newReal->setFromTriplets(tripletList.begin(), tripletList.end());
2442
2443             delete matrixReal;
2444             matrixReal = newReal;
2445             delete[] pRows;
2446             delete[] pNonZeroR;
2447             delete[] pNonZeroI;
2448         }
2449         else
2450         {
2451             //item count
2452             size_t iNonZeros = nonZeros();
2453             CplxSparse_t *newCplx = new CplxSparse_t(_iNewRows, _iNewCols);
2454             newCplx->reserve((int)iNonZeros);
2455
2456             //coords
2457             int* pRows = new int[iNonZeros * 2];
2458             outputRowCol(pRows);
2459             int* pCols = pRows + iNonZeros;
2460
2461             //values
2462             double* pNonZeroR = new double[iNonZeros];
2463             double* pNonZeroI = new double[iNonZeros];
2464             outputValues(pNonZeroR, pNonZeroI);
2465
2466             typedef Eigen::Triplet<std::complex<double> > triplet;
2467             std::vector<triplet> tripletList;
2468
2469             for (size_t i = 0 ; i < iNonZeros ; i++)
2470             {
2471                 int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
2472                 tripletList.push_back(triplet((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), std::complex<double>(pNonZeroR[i], pNonZeroI[i])));
2473             }
2474
2475             newCplx->setFromTriplets(tripletList.begin(), tripletList.end());
2476
2477             delete matrixCplx;
2478             matrixCplx = newCplx;
2479             delete[] pRows;
2480             delete[] pNonZeroR;
2481             delete[] pNonZeroI;
2482         }
2483
2484         m_iRows = _iNewRows;
2485         m_iCols = _iNewCols;
2486         m_iSize = _iNewRows * _iNewCols;
2487
2488         m_iDims = 2;
2489         m_piDims[0] = m_iRows;
2490         m_piDims[1] = m_iCols;
2491
2492         finalize();
2493
2494         res = true;
2495     }
2496     catch (...)
2497     {
2498         res = false;
2499     }
2500     return res;
2501 }
2502
2503 //    SparseBool* SparseBool::new
2504
2505 SparseBool::SparseBool(Bool SPARSE_CONST& src)
2506 {
2507     //compute idx
2508     int size = src.getSize();
2509     int row = src.getRows();
2510     Double* idx = new Double(src.getSize(), 2);
2511     double* p = idx->get();
2512     for (int i = 0; i < size; ++i)
2513     {
2514         p[i] = (double)(i % row) + 1;
2515         p[i + size] = (double)(i / row) + 1;
2516     }
2517     create2(src.getRows(), src.getCols(), src, *idx);
2518     idx->killMe();
2519 #ifndef NDEBUG
2520     Inspector::addItem(this);
2521 #endif
2522 }
2523 /* @param src : Bool matrix to copy into a new sparse matrix
2524 @param idx : Double matrix to use as indexes to get values from the src
2525 **/
2526 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
2527 {
2528     int idxrow = idx.getRows();
2529     int rows = static_cast<int>(*std::max_element(idx.get(), idx.get() + idxrow));
2530     int cols = static_cast<int>(*std::max_element(idx.get() + idxrow, idx.get() + idxrow * 2));
2531     create2(rows, cols, src, idx);
2532 #ifndef NDEBUG
2533     Inspector::addItem(this);
2534 #endif
2535 }
2536
2537 /* @param src : Bool matrix to copy into a new sparse matrix
2538 @param idx : Double matrix to use as indexes to get values from the src
2539 @param dims : Double matrix containing the dimensions of the new matrix
2540 **/
2541 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx, Double SPARSE_CONST& dims)
2542 {
2543     create2(static_cast<int>(dims.get(0)), static_cast<int>(dims.get(1)), src, idx);
2544 #ifndef NDEBUG
2545     Inspector::addItem(this);
2546 #endif
2547 }
2548
2549 SparseBool::SparseBool(int _iRows, int _iCols) : matrixBool(new BoolSparse_t(_iRows, _iCols))
2550 {
2551     m_iRows = _iRows;
2552     m_iCols = _iCols;
2553     m_iSize = _iRows * _iCols;
2554     m_iDims = 2;
2555     m_piDims[0] = _iRows;
2556     m_piDims[1] = _iCols;
2557 #ifndef NDEBUG
2558     Inspector::addItem(this);
2559 #endif
2560 }
2561
2562 SparseBool::SparseBool(SparseBool const& src) : matrixBool(new BoolSparse_t(*src.matrixBool))
2563 {
2564     m_iDims = 2;
2565     m_iRows = const_cast<SparseBool*>(&src)->getRows();
2566     m_iCols = const_cast<SparseBool*>(&src)->getCols();
2567     m_iSize = m_iRows * m_iCols;
2568     m_piDims[0] = m_iRows;
2569     m_piDims[1] = m_iCols;
2570 #ifndef NDEBUG
2571     Inspector::addItem(this);
2572 #endif
2573 }
2574
2575 SparseBool::SparseBool(BoolSparse_t* src) : matrixBool(src)
2576 {
2577     m_iRows = src->rows();
2578     m_iCols = src->cols();
2579     m_iSize = m_iRows * m_iCols;
2580     m_iDims = 2;
2581     m_piDims[0] = m_iRows;
2582     m_piDims[1] = m_iCols;
2583 #ifndef NDEBUG
2584     Inspector::addItem(this);
2585 #endif
2586 }
2587
2588 SparseBool::SparseBool(int rows, int cols, int trues, int* inner, int* outer)
2589 {
2590     int* out = nullptr;
2591     int* in = nullptr;
2592
2593     matrixBool = new BoolSparse_t(rows, cols);
2594     matrixBool->reserve((int)trues);
2595     out = matrixBool->outerIndexPtr();
2596     in = matrixBool->innerIndexPtr();
2597
2598     //update outerIndexPtr
2599     memcpy(out, outer, sizeof(int) * (rows + 1));
2600     //update innerIndexPtr
2601     memcpy(in, inner, sizeof(int) * trues);
2602
2603     bool* data = matrixBool->valuePtr();
2604     for (int i = 0; i < trues; ++i)
2605     {
2606         data[i] = true;
2607     }
2608
2609     m_iCols = cols;
2610     m_iRows = rows;
2611     m_iSize = cols * rows;
2612     m_iDims = 2;
2613     m_piDims[0] = m_iRows;
2614     m_piDims[1] = m_iCols;
2615
2616     matrixBool->resizeNonZeros(trues);
2617 }
2618
2619 void SparseBool::create2(int rows, int cols, Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
2620 {
2621     int nnz = src.getSize();
2622     double* i = idx.get();
2623     double* j = i + idx.getRows();
2624     int* val = src.get();
2625
2626     typedef Eigen::Triplet<bool> T;
2627     std::vector<T> tripletList;
2628     tripletList.reserve((int)nnz);
2629
2630     for (int k = 0; k < nnz; ++k)
2631     {
2632         tripletList.push_back(T(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, val[k] == 1));
2633     }
2634
2635     matrixBool = new BoolSparse_t(rows, cols);
2636     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
2637
2638     m_iRows = matrixBool->rows();
2639     m_iCols = matrixBool->cols();
2640     m_iSize = cols * rows;
2641     m_iDims = 2;
2642     m_piDims[0] = m_iRows;
2643     m_piDims[1] = m_iCols;
2644     finalize();
2645 }
2646
2647 SparseBool::~SparseBool()
2648 {
2649     delete matrixBool;
2650 #ifndef NDEBUG
2651     Inspector::removeItem(this);
2652 #endif
2653 }
2654
2655 bool SparseBool::toString(std::wostringstream& ostr) const
2656 {
2657     ostr << ::toString(*matrixBool, 0);
2658     return true;
2659 }
2660
2661 void SparseBool::whoAmI() SPARSE_CONST
2662 {
2663     std::cout << "types::SparseBool";
2664 }
2665
2666 SparseBool* SparseBool::clone(void) const
2667 {
2668     return new SparseBool(*this);
2669 }
2670
2671 bool SparseBool::resize(int _iNewRows, int _iNewCols)
2672 {
2673     if (_iNewRows <= getRows() && _iNewCols <= getCols())
2674     {
2675         //nothing to do: hence we do NOT fail
2676         return true;
2677     }
2678
2679     bool res = false;
2680     try
2681     {
2682         //item count
2683         size_t iNonZeros = nbTrue();
2684
2685         BoolSparse_t *newBool = new BoolSparse_t(_iNewRows, _iNewCols);
2686         newBool->reserve((int)iNonZeros);
2687
2688         //coords
2689         int* pRows = new int[iNonZeros * 2];
2690         outputRowCol(pRows);
2691         int* pCols = pRows + iNonZeros;
2692
2693         typedef Eigen::Triplet<bool> triplet;
2694         std::vector<triplet> tripletList;
2695
2696         for (size_t i = 0 ; i < iNonZeros ; i++)
2697         {
2698             tripletList.push_back(triplet((int)pRows[i] - 1, (int)pCols[i] - 1, true));
2699         }
2700
2701         newBool->setFromTriplets(tripletList.begin(), tripletList.end());
2702
2703         delete matrixBool;
2704         matrixBool = newBool;
2705         delete[] pRows;
2706
2707         m_iRows = _iNewRows;
2708         m_iCols = _iNewCols;
2709         m_iSize = _iNewRows * _iNewCols;
2710         m_piDims[0] = m_iRows;
2711         m_piDims[1] = m_iCols;
2712
2713         res = true;
2714
2715     }
2716     catch (...)
2717     {
2718         res = false;
2719     }
2720     return res;
2721 }
2722
2723 SparseBool* SparseBool::insert(typed_list* _pArgs, SparseBool* _pSource)
2724 {
2725     bool bNeedToResize  = false;
2726     int iDims           = (int)_pArgs->size();
2727     if (iDims > 2)
2728     {
2729         //sparse are only in 2 dims
2730         return NULL;
2731     }
2732
2733     typed_list pArg;
2734
2735     int piMaxDim[2];
2736     int piCountDim[2];
2737
2738     //on case of resize
2739     int iNewRows    = 0;
2740     int iNewCols    = 0;
2741
2742     //evaluate each argument and replace by appropriate value and compute the count of combinations
2743     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
2744     if (iSeqCount == 0)
2745     {
2746         //free pArg content
2747         cleanIndexesArguments(_pArgs, &pArg);
2748         return this;
2749     }
2750
2751     if (iDims < 2)
2752     {
2753         //see as vector
2754         if (getRows() == 1 || getCols() == 1)
2755         {
2756             //vector or scalar
2757             if (getSize() < piMaxDim[0])
2758             {
2759                 bNeedToResize = true;
2760
2761                 //need to enlarge sparse dimensions
2762                 if (getCols() == 1 || getSize() == 0)
2763                 {
2764                     //column vector
2765                     iNewRows    = piMaxDim[0];
2766                     iNewCols    = 1;
2767                 }
2768                 else if (getRows() == 1)
2769                 {
2770                     //row vector
2771                     iNewRows    = 1;
2772                     iNewCols    = piMaxDim[0];
2773                 }
2774             }
2775         }
2776         else if (getSize() < piMaxDim[0])
2777         {
2778             //free pArg content
2779             cleanIndexesArguments(_pArgs, &pArg);
2780             //out of range
2781             return NULL;
2782         }
2783     }
2784     else
2785     {
2786         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
2787         {
2788             bNeedToResize = true;
2789             iNewRows = std::max(getRows(), piMaxDim[0]);
2790             iNewCols = std::max(getCols(), piMaxDim[1]);
2791         }
2792     }
2793
2794     //check number of insertion
2795     if (_pSource->isScalar() == false && _pSource->getSize() != iSeqCount)
2796     {
2797         //free pArg content
2798         cleanIndexesArguments(_pArgs, &pArg);
2799         return NULL;
2800     }
2801
2802     //now you are sure to be able to insert values
2803     if (bNeedToResize)
2804     {
2805         if (resize(iNewRows, iNewCols) == false)
2806         {
2807             //free pArg content
2808             cleanIndexesArguments(_pArgs, &pArg);
2809             return NULL;
2810         }
2811     }
2812
2813     if (iDims == 1)
2814     {
2815         double* pIdx = pArg[0]->getAs<Double>()->get();
2816         for (int i = 0 ; i < iSeqCount ; i++)
2817         {
2818             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
2819             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
2820
2821             if (_pSource->isScalar())
2822             {
2823                 set(iRow, iCol, _pSource->get(0, 0), false);
2824             }
2825             else
2826             {
2827                 int iRowOrig = i % _pSource->getRows();
2828                 int iColOrig = i / _pSource->getRows();
2829                 set(iRow, iCol, _pSource->get(iRowOrig, iColOrig), false);
2830             }
2831         }
2832     }
2833     else
2834     {
2835         double* pIdxRow = pArg[0]->getAs<Double>()->get();
2836         int iRowSize    = pArg[0]->getAs<Double>()->getSize();
2837         double* pIdxCol = pArg[1]->getAs<Double>()->get();
2838
2839         for (int i = 0 ; i < iSeqCount ; i++)
2840         {
2841             if (_pSource->isScalar())
2842             {
2843                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(0, 0), false);
2844             }
2845             else
2846             {
2847                 int iRowOrig = i % _pSource->getRows();
2848                 int iColOrig = i / _pSource->getRows();
2849                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(iRowOrig, iColOrig), false);
2850             }
2851         }
2852     }
2853
2854     finalize();
2855
2856     //free pArg content
2857     cleanIndexesArguments(_pArgs, &pArg);
2858
2859     return this;
2860 }
2861
2862 SparseBool* SparseBool::insert(typed_list* _pArgs, InternalType* _pSource)
2863 {
2864     bool bNeedToResize  = false;
2865     int iDims           = (int)_pArgs->size();
2866     if (iDims > 2)
2867     {
2868         //sparse are only in 2 dims
2869         return NULL;
2870     }
2871
2872     typed_list pArg;
2873
2874     int piMaxDim[2];
2875     int piCountDim[2];
2876
2877     //on case of resize
2878     int iNewRows    = 0;
2879     int iNewCols    = 0;
2880     Bool* pSource = _pSource->getAs<Bool>();
2881
2882     //evaluate each argument and replace by appropriate value and compute the count of combinations
2883     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
2884     if (iSeqCount == 0)
2885     {
2886         //free pArg content
2887         cleanIndexesArguments(_pArgs, &pArg);
2888         return this;
2889     }
2890
2891     if (iDims < 2)
2892     {
2893         //see as vector
2894         if (getRows() == 1 || getCols() == 1)
2895         {
2896             //vector or scalar
2897             bNeedToResize = true;
2898             if (getSize() < piMaxDim[0])
2899             {
2900                 //need to enlarge sparse dimensions
2901                 if (getCols() == 1 || getSize() == 0)
2902                 {
2903                     //column vector
2904                     iNewRows    = piMaxDim[0];
2905                     iNewCols    = 1;
2906                 }
2907                 else if (getRows() == 1)
2908                 {
2909                     //row vector
2910                     iNewRows    = 1;
2911                     iNewCols    = piMaxDim[0];
2912                 }
2913             }
2914         }
2915         else if (getSize() < piMaxDim[0])
2916         {
2917             //free pArg content
2918             cleanIndexesArguments(_pArgs, &pArg);
2919             //out of range
2920             return NULL;
2921         }
2922     }
2923     else
2924     {
2925         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
2926         {
2927             bNeedToResize = true;
2928             iNewRows = std::max(getRows(), piMaxDim[0]);
2929             iNewCols = std::max(getCols(), piMaxDim[1]);
2930         }
2931     }
2932
2933     //check number of insertion
2934     if (pSource->isScalar() == false && pSource->getSize() != iSeqCount)
2935     {
2936         //free pArg content
2937         cleanIndexesArguments(_pArgs, &pArg);
2938         return NULL;
2939     }
2940
2941     //now you are sure to be able to insert values
2942     if (bNeedToResize)
2943     {
2944         if (resize(iNewRows, iNewCols) == false)
2945         {
2946             //free pArg content
2947             cleanIndexesArguments(_pArgs, &pArg);
2948             return NULL;
2949         }
2950     }
2951
2952     if (iDims == 1)
2953     {
2954         double* pIdx = pArg[0]->getAs<Double>()->get();
2955         for (int i = 0 ; i < iSeqCount ; i++)
2956         {
2957             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
2958             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
2959             if (pSource->isScalar())
2960             {
2961                 set(iRow, iCol, pSource->get(0) != 0, false);
2962             }
2963             else
2964             {
2965                 set(iRow, iCol, pSource->get(i) != 0, false);
2966             }
2967         }
2968     }
2969     else
2970     {
2971         double* pIdxRow = pArg[0]->getAs<Double>()->get();
2972         int iRowSize    = pArg[0]->getAs<Double>()->getSize();
2973         double* pIdxCol = pArg[1]->getAs<Double>()->get();
2974
2975         for (int i = 0 ; i < iSeqCount ; i++)
2976         {
2977             if (pSource->isScalar())
2978             {
2979                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(0) != 0, false);
2980             }
2981             else
2982             {
2983                 int iRowOrig = i % pSource->getRows();
2984                 int iColOrig = i / pSource->getRows();
2985
2986                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(iRowOrig, iColOrig) != 0, false);
2987             }
2988         }
2989     }
2990
2991     finalize();
2992
2993     //free pArg content
2994     cleanIndexesArguments(_pArgs, &pArg);
2995     return this;
2996 }
2997
2998 SparseBool* SparseBool::remove(typed_list* _pArgs)
2999 {
3000     SparseBool* pOut = NULL;
3001     int iDims = (int)_pArgs->size();
3002     if (iDims > 2)
3003     {
3004         //sparse are only in 2 dims
3005         return NULL;
3006     }
3007
3008     typed_list pArg;
3009
3010     int piMaxDim[2];
3011     int piCountDim[2];
3012
3013     //evaluate each argument and replace by appropriate value and compute the count of combinations
3014     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3015     if (iSeqCount == 0)
3016     {
3017         //free pArg content
3018         cleanIndexesArguments(_pArgs, &pArg);
3019         return this;
3020     }
3021
3022     bool* pbFull = new bool[iDims];
3023     //coord must represent all values on a dimension
3024     for (int i = 0 ; i < iDims ; i++)
3025     {
3026         pbFull[i]       = false;
3027         int iDimToCheck = getVarMaxDim(i, iDims);
3028         int iIndexSize  = pArg[i]->getAs<GenericType>()->getSize();
3029
3030         //we can have index more than once
3031         if (iIndexSize >= iDimToCheck)
3032         {
3033             //size is good, now check datas
3034             double* pIndexes = getDoubleArrayFromDouble(pArg[i]);
3035             for (int j = 0 ; j < iDimToCheck ; j++)
3036             {
3037                 bool bFind = false;
3038                 for (int k = 0 ; k < iIndexSize ; k++)
3039                 {
3040                     if ((int)pIndexes[k] == j + 1)
3041                     {
3042                         bFind = true;
3043                         break;
3044                     }
3045                 }
3046                 pbFull[i]  = bFind;
3047             }
3048         }
3049     }
3050
3051     //only one dims can be not full/entire
3052     bool bNotEntire = false;
3053     int iNotEntire  = 0;
3054     bool bTooMuchNotEntire = false;
3055     for (int i = 0 ; i < iDims ; i++)
3056     {
3057         if (pbFull[i] == false)
3058         {
3059             if (bNotEntire == false)
3060             {
3061                 bNotEntire = true;
3062                 iNotEntire = i;
3063             }
3064             else
3065             {
3066                 bTooMuchNotEntire = true;
3067                 break;
3068             }
3069         }
3070     }
3071
3072     if (bTooMuchNotEntire == true)
3073     {
3074         //free pArg content
3075         cleanIndexesArguments(_pArgs, &pArg);
3076         return NULL;
3077     }
3078
3079     delete[] pbFull;
3080
3081     //find index to keep
3082     int iNotEntireSize          = pArg[iNotEntire]->getAs<GenericType>()->getSize();
3083     double* piNotEntireIndex    = getDoubleArrayFromDouble(pArg[iNotEntire]);
3084     int iKeepSize               = getVarMaxDim(iNotEntire, iDims);
3085     bool* pbKeep                = new bool[iKeepSize];
3086
3087     //fill pbKeep with true value
3088     for (int i = 0 ; i < iKeepSize ; i++)
3089     {
3090         pbKeep[i] = true;
3091     }
3092
3093     for (int i = 0 ; i < iNotEntireSize ; i++)
3094     {
3095         int idx = (int)piNotEntireIndex[i] - 1;
3096
3097         //don't care of value out of bounds
3098         if (idx < iKeepSize)
3099         {
3100             pbKeep[idx] = false;
3101         }
3102     }
3103
3104     int iNewDimSize = 0;
3105     for (int i = 0 ; i < iKeepSize ; i++)
3106     {
3107         if (pbKeep[i] == true)
3108         {
3109             iNewDimSize++;
3110         }
3111     }
3112     delete[] pbKeep;
3113
3114     int* piNewDims = new int[iDims];
3115     for (int i = 0 ; i < iDims ; i++)
3116     {
3117         if (i == iNotEntire)
3118         {
3119             piNewDims[i] = iNewDimSize;
3120         }
3121         else
3122         {
3123             piNewDims[i] = getVarMaxDim(i, iDims);
3124         }
3125     }
3126
3127     //remove last dimension if are == 1
3128     int iOrigDims = iDims;
3129     for (int i = (iDims - 1) ; i >= 2 ; i--)
3130     {
3131         if (piNewDims[i] == 1)
3132         {
3133             iDims--;
3134         }
3135         else
3136         {
3137             break;
3138         }
3139     }
3140
3141     if (iDims == 1)
3142     {
3143         if (iNewDimSize == 0)
3144         {
3145             //free pArg content
3146             cleanIndexesArguments(_pArgs, &pArg);
3147             return new SparseBool(0, 0);
3148         }
3149         else
3150         {
3151             //two cases, depends of original matrix/vector
3152             if ((*_pArgs)[0]->isColon() == false && m_iDims == 2 && m_piDims[0] == 1 && m_piDims[1] != 1)
3153             {
3154                 //special case for row vector
3155                 pOut = new SparseBool(1, iNewDimSize);
3156                 //in this case we have to care of 2nd dimension
3157                 //iNotEntire = 1;
3158             }
3159             else
3160             {
3161                 pOut = new SparseBool(iNewDimSize, 1);
3162             }
3163         }
3164     }
3165     else
3166     {
3167         pOut = new SparseBool(piNewDims[0], piNewDims[0]);
3168     }
3169
3170     delete[] piNewDims;
3171     //find a way to copy existing data to new variable ...
3172     int iNewPos = 0;
3173     int* piIndexes = new int[iOrigDims];
3174     int* piViewDims = new int[iOrigDims];
3175     for (int i = 0 ; i < iOrigDims ; i++)
3176     {
3177         piViewDims[i] = getVarMaxDim(i, iOrigDims);
3178     }
3179
3180     for (int i = 0 ; i < getSize() ; i++)
3181     {
3182         bool bByPass = false;
3183         getIndexesWithDims(i, piIndexes, piViewDims, iOrigDims);
3184
3185         //check if piIndexes use removed indexes
3186         for (int j = 0 ; j < iNotEntireSize ; j++)
3187         {
3188             if ((piNotEntireIndex[j] - 1) == piIndexes[iNotEntire])
3189             {
3190                 //by pass this value
3191                 bByPass = true;
3192                 break;
3193             }
3194         }
3195
3196         if (bByPass == false)
3197         {
3198             //compute new index
3199             pOut->set(iNewPos, get(i));
3200             iNewPos++;
3201         }
3202     }
3203
3204     //free allocated data
3205     for (int i = 0 ; i < iDims ; i++)
3206     {
3207         if (pArg[i] != (*_pArgs)[i])
3208         {
3209             delete pArg[i];
3210         }
3211     }
3212
3213     delete[] piIndexes;
3214     delete[] piViewDims;
3215
3216     //free pArg content
3217     cleanIndexesArguments(_pArgs, &pArg);
3218
3219     return pOut;
3220 }
3221
3222 bool SparseBool::append(int r, int c, SparseBool SPARSE_CONST* src)
3223 {
3224     doAppend(*src, r, c, *matrixBool);
3225     finalize();
3226     return true;
3227 }
3228
3229 InternalType* SparseBool::insertNew(typed_list* _pArgs, InternalType* _pSource)
3230 {
3231     typed_list pArg;
3232     InternalType *pOut  = NULL;
3233     SparseBool* pSource = _pSource->getAs<SparseBool>();
3234
3235     int iDims           = (int)_pArgs->size();
3236     int* piMaxDim       = new int[iDims];
3237     int* piCountDim     = new int[iDims];
3238     bool bUndefine      = false;
3239
3240     //evaluate each argument and replace by appropriate value and compute the count of combinations
3241     int iSeqCount = checkIndexesArguments(NULL, _pArgs, &pArg, piMaxDim, piCountDim);
3242     if (iSeqCount == 0)
3243     {
3244         //free pArg content
3245         cleanIndexesArguments(_pArgs, &pArg);
3246         return createEmptyDouble();
3247     }
3248
3249     if (iSeqCount < 0)
3250     {
3251         iSeqCount = -iSeqCount;
3252         bUndefine = true;
3253     }
3254
3255     if (bUndefine)
3256     {
3257         //manage : and $ in creation by insertion
3258         int iSource         = 0;
3259         int *piSourceDims   = pSource->getDimsArray();
3260
3261         for (int i = 0 ; i < iDims ; i++)
3262         {
3263             if (pArg[i] == NULL)
3264             {
3265                 //undefine value
3266                 if (pSource->isScalar())
3267                 {
3268                     piMaxDim[i]     = 1;
3269                     piCountDim[i]   = 1;
3270                 }
3271                 else
3272                 {
3273                     piMaxDim[i]     = piSourceDims[iSource];
3274                     piCountDim[i]   = piSourceDims[iSource];
3275                 }
3276                 iSource++;
3277                 //replace pArg value by the new one
3278                 pArg[i] = createDoubleVector(piMaxDim[i]);
3279             }
3280             //else
3281             //{
3282             //    piMaxDim[i] = piCountDim[i];
3283             //}
3284         }
3285     }
3286
3287     //remove last dimension at size 1
3288     //remove last dimension if are == 1
3289     for (int i = (iDims - 1) ; i >= 2 ; i--)
3290     {
3291         if (piMaxDim[i] == 1)
3292         {
3293             iDims--;
3294             pArg.pop_back();
3295         }
3296         else
3297         {
3298             break;
3299         }
3300     }
3301
3302     if (checkArgValidity(pArg) == false)
3303     {
3304         //free pArg content
3305         cleanIndexesArguments(_pArgs, &pArg);
3306         //contain bad index, like <= 0, ...
3307         return NULL;
3308     }
3309
3310     if (iDims == 1)
3311     {
3312         if (pSource->getCols() == 1)
3313         {
3314             pOut = new SparseBool(piCountDim[0], 1);
3315         }
3316         else
3317         {
3318             //rows == 1
3319             pOut = new SparseBool(1, piCountDim[0]);
3320         }
3321     }
3322     else
3323     {
3324         pOut = new SparseBool(piMaxDim[0], piMaxDim[1]);
3325     }
3326
3327     //fill with null item
3328     SparseBool* pSpOut = pOut->getAs<SparseBool>();
3329
3330     //insert values in new matrix
3331     InternalType* pOut2 = pSpOut->insert(&pArg, pSource);
3332     if (pOut != pOut2)
3333     {
3334         delete pOut;
3335     }
3336
3337     //free pArg content
3338     cleanIndexesArguments(_pArgs, &pArg);
3339
3340     return pOut2;
3341 }
3342
3343 SparseBool* SparseBool::extract(int nbCoords, int SPARSE_CONST* coords, int SPARSE_CONST* maxCoords, int SPARSE_CONST* resSize, bool asVector) SPARSE_CONST
3344 {
3345     if ( (asVector && maxCoords[0] > getSize()) ||
3346     (asVector == false && maxCoords[0] > getRows()) ||
3347     (asVector == false && maxCoords[1] > getCols()))
3348     {
3349         return 0;
3350     }
3351
3352     SparseBool * pSp (0);
3353     if (asVector)
3354     {
3355         pSp = (getRows() == 1) ?  new SparseBool(1, resSize[0]) : new SparseBool(resSize[0], 1);
3356         mycopy_n(makeMatrixIterator<bool>(*this,  Coords<true>(coords, getRows())), nbCoords
3357         , makeMatrixIterator<bool>(*(pSp->matrixBool), RowWiseFullIterator(pSp->getRows(), pSp->getCols())));
3358     }
3359     else
3360     {
3361         pSp = new SparseBool(resSize[0], resSize[1]);
3362         mycopy_n(makeMatrixIterator<bool>(*this,  Coords<false>(coords, getRows())), nbCoords
3363         , makeMatrixIterator<bool>(*(pSp->matrixBool), RowWiseFullIterator(pSp->getRows(), pSp->getCols())));
3364
3365     }
3366     return pSp;
3367 }
3368
3369 /*
3370 * create a new SparseBool of dims according to resSize and fill it from currentSparseBool (along coords)
3371 */
3372 InternalType* SparseBool::extract(typed_list* _pArgs)
3373 {
3374     SparseBool* pOut    = NULL;
3375     int iDims           = (int)_pArgs->size();
3376     typed_list pArg;
3377
3378     int* piMaxDim       = new int[iDims];
3379     int* piCountDim     = new int[iDims];
3380
3381     //evaluate each argument and replace by appropriate value and compute the count of combinations
3382     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3383     if (iSeqCount == 0)
3384     {
3385         //free pArg content
3386         cleanIndexesArguments(_pArgs, &pArg);
3387         if (_pArgs->size() == 0)
3388         {
3389             delete[] piMaxDim;
3390             delete[] piCountDim;
3391             //a()
3392             return this;
3393         }
3394         else
3395         {
3396             delete[] piMaxDim;
3397             delete[] piCountDim;
3398             //a([])
3399             return Double::Empty();
3400         }
3401     }
3402
3403     if (iDims < 2)
3404     {
3405         // Check that we stay inside the input size.
3406         if (piMaxDim[0] <= getSize())
3407         {
3408             int iNewRows = 0;
3409             int iNewCols = 0;
3410
3411             if (getRows() == 1 && getCols() != 1 && (*_pArgs)[0]->isColon() == false)
3412             {
3413                 //special case for row vector
3414                 iNewRows = 1;
3415                 iNewCols = piCountDim[0];
3416             }
3417             else
3418             {
3419                 iNewRows = piCountDim[0];
3420                 iNewCols = 1;
3421             }
3422
3423             pOut = new SparseBool(iNewRows, iNewCols);
3424             double* pIdx = pArg[0]->getAs<Double>()->get();
3425             // Write in output all elements extract from input.
3426             for (int i = 0 ; i < iSeqCount ; i++)
3427             {
3428                 if (pIdx[i] < 1)
3429                 {
3430                     delete pOut;
3431                     pOut = NULL;
3432                     break;
3433                 }
3434                 int iRowRead = static_cast<int>(pIdx[i] - 1) % getRows();
3435                 int iColRead = static_cast<int>(pIdx[i] - 1) / getRows();
3436
3437                 int iRowWrite = static_cast<int>(i) % iNewRows;
3438                 int iColWrite = static_cast<int>(i) / iNewRows;
3439
3440                 bool bValue = get(iRowRead, iColRead);
3441                 if (bValue)
3442                 {
3443                     //only non zero values
3444                     pOut->set(iRowWrite, iColWrite, true, false);
3445                 }
3446             }
3447         }
3448         else
3449         {
3450             delete[] piMaxDim;
3451             delete[] piCountDim;
3452             //free pArg content
3453             cleanIndexesArguments(_pArgs, &pArg);
3454             return NULL;
3455         }
3456     }
3457     else
3458     {
3459         // Check that we stay inside the input size.
3460         if (piMaxDim[0] <= getRows() && piMaxDim[1] <= getCols())
3461         {
3462             double* pIdxRow = pArg[0]->getAs<Double>()->get();
3463             double* pIdxCol = pArg[1]->getAs<Double>()->get();
3464
3465             int iNewRows = pArg[0]->getAs<Double>()->getSize();
3466             int iNewCols = pArg[1]->getAs<Double>()->getSize();
3467
3468             pOut = new SparseBool(iNewRows, iNewCols);
3469
3470             int iPos = 0;
3471             // Write in output all elements extract from input.
3472             for (int iRow = 0 ; iRow < iNewRows ; iRow++)
3473             {
3474                 for (int iCol = 0 ; iCol < iNewCols ; iCol++)
3475                 {
3476                     if ((pIdxRow[iRow] < 1) || (pIdxCol[iCol] < 1))
3477                     {
3478                         delete pOut;
3479                         pOut = NULL;
3480                         break;
3481                     }
3482                     bool bValue = get((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
3483                     if (bValue)
3484                     {
3485                         //only non zero values
3486                         pOut->set(iRow, iCol, true, false);
3487                     }
3488                     iPos++;
3489                 }
3490             }
3491         }
3492         else
3493         {
3494             delete[] piMaxDim;
3495             delete[] piCountDim;
3496             //free pArg content
3497             cleanIndexesArguments(_pArgs, &pArg);
3498             return NULL;
3499         }
3500     }
3501
3502     finalize();
3503
3504     delete[] piMaxDim;
3505     delete[] piCountDim;
3506     //free pArg content
3507     cleanIndexesArguments(_pArgs, &pArg);
3508
3509     return pOut;
3510 }
3511
3512 std::size_t SparseBool::nbTrue() const
3513 {
3514     return  matrixBool->nonZeros() ;
3515 }
3516 std::size_t SparseBool::nbTrue(std::size_t r) const
3517 {
3518     int* piIndex = matrixBool->outerIndexPtr();
3519     return piIndex[r + 1] - piIndex[r];
3520 }
3521
3522
3523 void SparseBool::setTrue(bool finalize)
3524 {
3525     int rows = getRows();
3526     int cols = getCols();
3527
3528     typedef Eigen::Triplet<bool> triplet;
3529     std::vector<triplet> tripletList;
3530
3531     for (int i = 0; i < rows; ++i)
3532     {
3533         for (int j = 0; j < cols; ++j)
3534         {
3535             tripletList.push_back(triplet(i, j, true));
3536         }
3537     }
3538
3539     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
3540
3541     if (finalize)
3542     {
3543         matrixBool->finalize();
3544     }
3545 }
3546
3547 void SparseBool::setFalse(bool finalize)
3548 {
3549     int rows = getRows();
3550     int cols = getCols();
3551
3552     typedef Eigen::Triplet<bool> triplet;
3553     std::vector<triplet> tripletList;
3554
3555     for (int i = 0; i < rows; ++i)
3556     {
3557         for (int j = 0; j < cols; ++j)
3558         {
3559             tripletList.push_back(triplet(i, j, false));
3560         }
3561     }
3562
3563     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
3564
3565     if (finalize)
3566     {
3567         matrixBool->finalize();
3568     }
3569 }
3570
3571 int* SparseBool::getNbItemByRow(int* _piNbItemByRows)
3572 {
3573     int* piNbItemByRows = new int[getRows() + 1];
3574     mycopy_n(matrixBool->outerIndexPtr(), getRows() + 1, piNbItemByRows);
3575
3576     for (int i = 0 ; i < getRows() ; i++)
3577     {
3578         _piNbItemByRows[i] = piNbItemByRows[i + 1] - piNbItemByRows[i];
3579     }
3580
3581     delete[] piNbItemByRows;
3582     return _piNbItemByRows;
3583 }
3584
3585 int* SparseBool::getColPos(int* _piColPos)
3586 {
3587     mycopy_n(matrixBool->innerIndexPtr(), nbTrue(), _piColPos);
3588     for (int i = 0; i < nbTrue(); i++)
3589     {
3590         _piColPos[i]++;
3591     }
3592
3593     return _piColPos;
3594 }
3595
3596 int* SparseBool::outputRowCol(int* out)const
3597 {
3598     return sparseTransform(*matrixBool, sparseTransform(*matrixBool, out, GetRow<BoolSparse_t>()), GetCol<BoolSparse_t>());
3599 }
3600
3601 int* SparseBool::getInnerPtr(int* count)
3602 {
3603     *count = matrixBool->innerSize();
3604     return matrixBool->innerIndexPtr();
3605 }
3606
3607 int* SparseBool::getOuterPtr(int* count)
3608 {
3609     *count = matrixBool->outerSize();
3610     return matrixBool->outerIndexPtr();
3611 }
3612
3613
3614 bool SparseBool::operator==(const InternalType& it) SPARSE_CONST
3615 {
3616     SparseBool* otherSparse = const_cast<SparseBool*>(dynamic_cast<SparseBool const*>(&it));/* types::GenericType is not const-correct :( */
3617     return (otherSparse
3618     && (otherSparse->getRows() == getRows())
3619     && (otherSparse->getCols() == getCols())
3620     && equal(*matrixBool, *otherSparse->matrixBool));
3621 }
3622
3623 bool SparseBool::operator!=(const InternalType& it) SPARSE_CONST
3624 {
3625     return !(*this == it);
3626 }
3627
3628 void SparseBool::finalize()
3629 {
3630     matrixBool->prune(&keepForSparse<bool>);
3631     matrixBool->finalize();
3632 }
3633
3634 bool SparseBool::get(int r, int c) SPARSE_CONST
3635 {
3636     return matrixBool->coeff(r, c);
3637 }
3638
3639 bool SparseBool::set(int _iRows, int _iCols, bool _bVal, bool _bFinalize) SPARSE_CONST
3640 {
3641     matrixBool->coeffRef(_iRows, _iCols) = _bVal;
3642
3643     if (_bFinalize)
3644     {
3645         finalize();
3646     }
3647
3648     return true;
3649 }
3650
3651 void SparseBool::fill(Bool& dest, int r, int c) SPARSE_CONST
3652 {
3653     mycopy_n(makeMatrixIterator<bool >(*matrixBool, RowWiseFullIterator(getRows(), getCols())), getSize()
3654     , makeMatrixIterator<bool >(dest, RowWiseFullIterator(dest.getRows(), dest.getCols(), r, c)));
3655 }
3656
3657 Sparse* SparseBool::newOnes() const
3658 {
3659     return new Sparse(new types::Sparse::RealSparse_t(matrixBool->cast<double>()), 0);
3660 }
3661
3662 SparseBool* SparseBool::newNotEqualTo(SparseBool const&o) const
3663 {
3664     return cwiseOp<std::not_equal_to>(*this, o);
3665 }
3666
3667 SparseBool* SparseBool::newEqualTo(SparseBool& o)
3668 {
3669     int rowL = getRows();
3670     int colL = getCols();
3671
3672     int rowR = o.getRows();
3673     int colR = o.getCols();
3674     int row = std::max(rowL, rowR);
3675     int col = std::max(colL, colR);
3676
3677     //create a boolean sparse matrix with dims of sparses
3678     types::SparseBool* ret = new types::SparseBool(row, col);
3679
3680     if (isScalar() && o.isScalar())
3681     {
3682         bool l = get(0, 0);
3683         bool r = o.get(0, 0);
3684         ret->set(0, 0, l == r, false);
3685     }
3686     else if (isScalar())
3687     {
3688         int nnzR = static_cast<int>(o.nbTrue());
3689         std::vector<int> rowcolR(nnzR * 2, 0);
3690         o.outputRowCol(rowcolR.data());
3691
3692         //compare all items of R with R[0]
3693         bool l = get(0, 0);
3694         for (int i = 0; i < nnzR; ++i)
3695         {
3696             bool r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
3697             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
3698         }
3699     }
3700     else if (o.isScalar())
3701     {
3702         int nnzL = static_cast<int>(nbTrue());
3703         std::vector<int> rowcolL(nnzL * 2, 0);
3704         outputRowCol(rowcolL.data());
3705
3706         bool r = get(0, 0);
3707         for (int i = 0; i < nnzL; ++i)
3708         {
3709             bool l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
3710             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
3711         }
3712     }
3713     else
3714     {
3715         int nnzR = static_cast<int>(o.nbTrue());
3716         std::vector<int> rowcolR(nnzR * 2, 0);
3717         o.outputRowCol(rowcolR.data());
3718         int nnzL = static_cast<int>(nbTrue());
3719         std::vector<int> rowcolL(nnzL * 2, 0);
3720         outputRowCol(rowcolL.data());
3721         //set all values to %t
3722         ret->setTrue(false);
3723         //set %f in each pL values
3724         for (int i = 0; i < nnzL; ++i)
3725         {
3726             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, false, false);
3727         }
3728         ret->finalize();
3729
3730         //set _pR[i] == _pL[i] for each _pR values
3731         for (int i = 0; i < nnzR; ++i)
3732         {
3733             //get l and r following non zeros value of R
3734             bool l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
3735             bool r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
3736             //set value following non zeros value of R
3737             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
3738         }
3739     }
3740
3741     ret->finalize();
3742     return ret;
3743 }
3744
3745 SparseBool* SparseBool::newLogicalOr(SparseBool const&o) const
3746 {
3747     return cwiseOp<std::logical_or>(*this, o);
3748 }
3749
3750 SparseBool* SparseBool::newLogicalAnd(SparseBool const&o) const
3751 {
3752     return cwiseOp<std::logical_and>(*this, o);
3753 }
3754
3755 bool SparseBool::reshape(int* _piDims, int _iDims)
3756 {
3757     bool bOk = false;
3758     int iCols = 1;
3759
3760     if (_iDims == 2)
3761     {
3762         iCols = _piDims[1];
3763     }
3764
3765     if (_iDims <= 2)
3766     {
3767         bOk = reshape(_piDims[0], iCols);
3768     }
3769
3770     return bOk;
3771 }
3772
3773 bool SparseBool::reshape(int _iNewRows, int _iNewCols)
3774 {
3775     if (_iNewRows * _iNewCols != getRows() * getCols())
3776     {
3777         return false;
3778     }
3779
3780     bool res = false;
3781     try
3782     {
3783         //item count
3784         size_t iNonZeros = matrixBool->nonZeros();
3785         BoolSparse_t *newBool = new BoolSparse_t(_iNewRows, _iNewCols);
3786         newBool->reserve((int)iNonZeros);
3787
3788         //coords
3789         int* pRows = new int[iNonZeros * 2];
3790         outputRowCol(pRows);
3791         int* pCols = pRows + iNonZeros;
3792
3793         typedef Eigen::Triplet<bool> triplet;
3794         std::vector<triplet> tripletList;
3795
3796         for (size_t i = 0 ; i < iNonZeros ; i++)
3797         {
3798             int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
3799             tripletList.push_back(triplet((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), true));
3800         }
3801
3802         newBool->setFromTriplets(tripletList.begin(), tripletList.end());
3803
3804         delete matrixBool;
3805         matrixBool = newBool;
3806         delete[] pRows;
3807
3808         m_iRows = _iNewRows;
3809         m_iCols = _iNewCols;
3810         m_iSize = _iNewRows * _iNewCols;
3811
3812         m_iDims = 2;
3813         m_piDims[0] = m_iRows;
3814         m_piDims[1] = m_iCols;
3815
3816         finalize();
3817
3818         res = true;
3819     }
3820     catch (...)
3821     {
3822         res = false;
3823     }
3824     return res;
3825 }
3826
3827 bool SparseBool::transpose(InternalType *& out)
3828 {
3829     out = new SparseBool(new BoolSparse_t(matrixBool->transpose()));
3830     return true;
3831 }
3832
3833 template<typename T>
3834 void neg(const int r, const int c, const T * const in, Eigen::SparseMatrix<bool, 1> * const out)
3835 {
3836     for (int i = 0; i < r; i++)
3837     {
3838         for (int j = 0; j < c; j++)
3839         {
3840             out->coeffRef(i, j) = !in->coeff(i, j);
3841         }
3842     }
3843
3844     out->prune(&keepForSparse<bool>);
3845     out->finalize();
3846 }
3847
3848
3849 }