fix comparison between sparse and sparse or sparsebool and sparsebool ( intel compile...
[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
2276     std::cout << "row " << row << std::endl << "col " << col << std::endl;
2277
2278     if (isScalar() && o.isScalar())
2279     {
2280         if (isComplex() || o.isComplex())
2281         {
2282             std::complex<double> l = getImg(0, 0);
2283             std::complex<double> r = o.getImg(0, 0);
2284             ret->set(0, 0, l == r, false);
2285         }
2286         else
2287         {
2288             double l = get(0, 0);
2289             double r = o.get(0, 0);
2290             ret->set(0, 0, l == r, false);
2291         }
2292     }
2293     else if (isScalar())
2294     {
2295         int nnzR = static_cast<int>(o.nonZeros());
2296         std::vector<int> rowcolR(nnzR * 2, 0);
2297         o.outputRowCol(rowcolR.data());
2298
2299         //compare all items of R with R[0]
2300         if (isComplex() || o.isComplex())
2301         {
2302             std::complex<double> l = getImg(0, 0);
2303             for (int i = 0; i < nnzR; ++i)
2304             {
2305                 std::complex<double> r = o.getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2306                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2307             }
2308         }
2309         else
2310         {
2311             double l = get(0, 0);
2312             for (int i = 0; i < nnzR; ++i)
2313             {
2314                 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2315                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2316             }
2317         }
2318     }
2319     else if (o.isScalar())
2320     {
2321         int nnzL = static_cast<int>(nonZeros());
2322         std::vector<int> rowcolL(nnzL * 2, 0);
2323         outputRowCol(rowcolL.data());
2324
2325         if (isComplex() || o.isComplex())
2326         {
2327             std::complex<double> r = o.getImg(0, 0);
2328             for (int i = 0; i < nnzL; ++i)
2329             {
2330                 std::complex<double> l = getImg(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2331                 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
2332             }
2333         }
2334         else
2335         {
2336             double r = get(0, 0);
2337             for (int i = 0; i < nnzL; ++i)
2338             {
2339                 double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2340                 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
2341             }
2342         }
2343     }
2344     else
2345     {
2346         int nnzR = static_cast<int>(o.nonZeros());
2347         std::vector<int> rowcolR(nnzR * 2, 0);
2348         o.outputRowCol(rowcolR.data());
2349         int nnzL = static_cast<int>(nonZeros());
2350         std::vector<int> rowcolL(nnzL * 2, 0);
2351         outputRowCol(rowcolL.data());
2352         //set all values to %t
2353         ret->setTrue(false);
2354         //set %f in each pL values
2355         for (int i = 0; i < nnzL; ++i)
2356         {
2357             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, false, false);
2358         }
2359         ret->finalize();
2360
2361         //set _pR[i] == _pL[i] for each _pR values
2362         if (isComplex() || o.isComplex())
2363         {
2364             for (int i = 0; i < nnzR; ++i)
2365             {
2366                 //get l and r following non zeros value of R
2367                 std::complex<double> l = getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2368                 std::complex<double> r = o.getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2369                 //set value following non zeros value of R
2370                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2371             }
2372         }
2373         else
2374         {
2375             for (int i = 0; i < nnzR; ++i)
2376             {
2377                 //get l and r following non zeros value of R
2378                 double l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2379                 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2380                 //set value following non zeros value of R
2381                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2382             }
2383         }
2384     }
2385
2386     ret->finalize();
2387     return ret;
2388 }
2389
2390 bool Sparse::reshape(int* _piDims, int _iDims)
2391 {
2392     bool bOk = false;
2393     int iCols = 1;
2394
2395     if (_iDims == 2)
2396     {
2397         iCols = _piDims[1];
2398     }
2399
2400     if (_iDims <= 2)
2401     {
2402         bOk = reshape(_piDims[0], iCols);
2403     }
2404
2405     return bOk;
2406 }
2407
2408 bool Sparse::reshape(int _iNewRows, int _iNewCols)
2409 {
2410     if (_iNewRows * _iNewCols != getRows() * getCols())
2411     {
2412         return false;
2413     }
2414
2415     bool res = false;
2416     try
2417     {
2418         if (matrixReal)
2419         {
2420             //item count
2421             size_t iNonZeros = nonZeros();
2422             RealSparse_t *newReal = new RealSparse_t(_iNewRows, _iNewCols);
2423             newReal->reserve((int)iNonZeros);
2424
2425             //coords
2426             int* pRows = new int[iNonZeros * 2];
2427             outputRowCol(pRows);
2428             int* pCols = pRows + iNonZeros;
2429
2430             //values
2431             double* pNonZeroR = new double[iNonZeros];
2432             double* pNonZeroI = new double[iNonZeros];
2433             outputValues(pNonZeroR, pNonZeroI);
2434
2435             typedef Eigen::Triplet<double> triplet;
2436             std::vector<triplet> tripletList;
2437
2438             for (size_t i = 0 ; i < iNonZeros ; i++)
2439             {
2440                 int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
2441                 tripletList.push_back(triplet((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), pNonZeroR[i]));
2442             }
2443
2444             newReal->setFromTriplets(tripletList.begin(), tripletList.end());
2445
2446             delete matrixReal;
2447             matrixReal = newReal;
2448             delete[] pRows;
2449             delete[] pNonZeroR;
2450             delete[] pNonZeroI;
2451         }
2452         else
2453         {
2454             //item count
2455             size_t iNonZeros = nonZeros();
2456             CplxSparse_t *newCplx = new CplxSparse_t(_iNewRows, _iNewCols);
2457             newCplx->reserve((int)iNonZeros);
2458
2459             //coords
2460             int* pRows = new int[iNonZeros * 2];
2461             outputRowCol(pRows);
2462             int* pCols = pRows + iNonZeros;
2463
2464             //values
2465             double* pNonZeroR = new double[iNonZeros];
2466             double* pNonZeroI = new double[iNonZeros];
2467             outputValues(pNonZeroR, pNonZeroI);
2468
2469             typedef Eigen::Triplet<std::complex<double> > triplet;
2470             std::vector<triplet> tripletList;
2471
2472             for (size_t i = 0 ; i < iNonZeros ; i++)
2473             {
2474                 int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
2475                 tripletList.push_back(triplet((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), std::complex<double>(pNonZeroR[i], pNonZeroI[i])));
2476             }
2477
2478             newCplx->setFromTriplets(tripletList.begin(), tripletList.end());
2479
2480             delete matrixCplx;
2481             matrixCplx = newCplx;
2482             delete[] pRows;
2483             delete[] pNonZeroR;
2484             delete[] pNonZeroI;
2485         }
2486
2487         m_iRows = _iNewRows;
2488         m_iCols = _iNewCols;
2489         m_iSize = _iNewRows * _iNewCols;
2490
2491         m_iDims = 2;
2492         m_piDims[0] = m_iRows;
2493         m_piDims[1] = m_iCols;
2494
2495         finalize();
2496
2497         res = true;
2498     }
2499     catch (...)
2500     {
2501         res = false;
2502     }
2503     return res;
2504 }
2505
2506 //    SparseBool* SparseBool::new
2507
2508 SparseBool::SparseBool(Bool SPARSE_CONST& src)
2509 {
2510     //compute idx
2511     int size = src.getSize();
2512     int row = src.getRows();
2513     Double* idx = new Double(src.getSize(), 2);
2514     double* p = idx->get();
2515     for (int i = 0; i < size; ++i)
2516     {
2517         p[i] = (double)(i % row) + 1;
2518         p[i + size] = (double)(i / row) + 1;
2519     }
2520     create2(src.getRows(), src.getCols(), src, *idx);
2521     idx->killMe();
2522 #ifndef NDEBUG
2523     Inspector::addItem(this);
2524 #endif
2525 }
2526 /* @param src : Bool matrix to copy into a new sparse matrix
2527 @param idx : Double matrix to use as indexes to get values from the src
2528 **/
2529 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
2530 {
2531     int idxrow = idx.getRows();
2532     int rows = static_cast<int>(*std::max_element(idx.get(), idx.get() + idxrow));
2533     int cols = static_cast<int>(*std::max_element(idx.get() + idxrow, idx.get() + idxrow * 2));
2534     create2(rows, cols, src, idx);
2535 #ifndef NDEBUG
2536     Inspector::addItem(this);
2537 #endif
2538 }
2539
2540 /* @param src : Bool matrix to copy into a new sparse matrix
2541 @param idx : Double matrix to use as indexes to get values from the src
2542 @param dims : Double matrix containing the dimensions of the new matrix
2543 **/
2544 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx, Double SPARSE_CONST& dims)
2545 {
2546     create2(static_cast<int>(dims.get(0)), static_cast<int>(dims.get(1)), src, idx);
2547 #ifndef NDEBUG
2548     Inspector::addItem(this);
2549 #endif
2550 }
2551
2552 SparseBool::SparseBool(int _iRows, int _iCols) : matrixBool(new BoolSparse_t(_iRows, _iCols))
2553 {
2554     m_iRows = _iRows;
2555     m_iCols = _iCols;
2556     m_iSize = _iRows * _iCols;
2557     m_iDims = 2;
2558     m_piDims[0] = _iRows;
2559     m_piDims[1] = _iCols;
2560 #ifndef NDEBUG
2561     Inspector::addItem(this);
2562 #endif
2563 }
2564
2565 SparseBool::SparseBool(SparseBool const& src) : matrixBool(new BoolSparse_t(*src.matrixBool))
2566 {
2567     m_iDims = 2;
2568     m_iRows = const_cast<SparseBool*>(&src)->getRows();
2569     m_iCols = const_cast<SparseBool*>(&src)->getCols();
2570     m_iSize = m_iRows * m_iCols;
2571     m_piDims[0] = m_iRows;
2572     m_piDims[1] = m_iCols;
2573 #ifndef NDEBUG
2574     Inspector::addItem(this);
2575 #endif
2576 }
2577
2578 SparseBool::SparseBool(BoolSparse_t* src) : matrixBool(src)
2579 {
2580     m_iRows = src->rows();
2581     m_iCols = src->cols();
2582     m_iSize = m_iRows * m_iCols;
2583     m_iDims = 2;
2584     m_piDims[0] = m_iRows;
2585     m_piDims[1] = m_iCols;
2586 #ifndef NDEBUG
2587     Inspector::addItem(this);
2588 #endif
2589 }
2590
2591 SparseBool::SparseBool(int rows, int cols, int trues, int* inner, int* outer)
2592 {
2593     int* out = nullptr;
2594     int* in = nullptr;
2595
2596     matrixBool = new BoolSparse_t(rows, cols);
2597     matrixBool->reserve((int)trues);
2598     out = matrixBool->outerIndexPtr();
2599     in = matrixBool->innerIndexPtr();
2600
2601     //update outerIndexPtr
2602     memcpy(out, outer, sizeof(int) * (rows + 1));
2603     //update innerIndexPtr
2604     memcpy(in, inner, sizeof(int) * trues);
2605
2606     bool* data = matrixBool->valuePtr();
2607     for (int i = 0; i < trues; ++i)
2608     {
2609         data[i] = true;
2610     }
2611
2612     m_iCols = cols;
2613     m_iRows = rows;
2614     m_iSize = cols * rows;
2615     m_iDims = 2;
2616     m_piDims[0] = m_iRows;
2617     m_piDims[1] = m_iCols;
2618
2619     matrixBool->resizeNonZeros(trues);
2620 }
2621
2622 void SparseBool::create2(int rows, int cols, Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
2623 {
2624     int nnz = src.getSize();
2625     double* i = idx.get();
2626     double* j = i + idx.getRows();
2627     int* val = src.get();
2628
2629     typedef Eigen::Triplet<bool> T;
2630     std::vector<T> tripletList;
2631     tripletList.reserve((int)nnz);
2632
2633     for (int k = 0; k < nnz; ++k)
2634     {
2635         tripletList.push_back(T(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, val[k] == 1));
2636     }
2637
2638     matrixBool = new BoolSparse_t(rows, cols);
2639     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
2640
2641     m_iRows = matrixBool->rows();
2642     m_iCols = matrixBool->cols();
2643     m_iSize = cols * rows;
2644     m_iDims = 2;
2645     m_piDims[0] = m_iRows;
2646     m_piDims[1] = m_iCols;
2647     finalize();
2648 }
2649
2650 SparseBool::~SparseBool()
2651 {
2652     delete matrixBool;
2653 #ifndef NDEBUG
2654     Inspector::removeItem(this);
2655 #endif
2656 }
2657
2658 bool SparseBool::toString(std::wostringstream& ostr) const
2659 {
2660     ostr << ::toString(*matrixBool, 0);
2661     return true;
2662 }
2663
2664 void SparseBool::whoAmI() SPARSE_CONST
2665 {
2666     std::cout << "types::SparseBool";
2667 }
2668
2669 SparseBool* SparseBool::clone(void) const
2670 {
2671     return new SparseBool(*this);
2672 }
2673
2674 bool SparseBool::resize(int _iNewRows, int _iNewCols)
2675 {
2676     if (_iNewRows <= getRows() && _iNewCols <= getCols())
2677     {
2678         //nothing to do: hence we do NOT fail
2679         return true;
2680     }
2681
2682     bool res = false;
2683     try
2684     {
2685         //item count
2686         size_t iNonZeros = nbTrue();
2687
2688         BoolSparse_t *newBool = new BoolSparse_t(_iNewRows, _iNewCols);
2689         newBool->reserve((int)iNonZeros);
2690
2691         //coords
2692         int* pRows = new int[iNonZeros * 2];
2693         outputRowCol(pRows);
2694         int* pCols = pRows + iNonZeros;
2695
2696         typedef Eigen::Triplet<bool> triplet;
2697         std::vector<triplet> tripletList;
2698
2699         for (size_t i = 0 ; i < iNonZeros ; i++)
2700         {
2701             tripletList.push_back(triplet((int)pRows[i] - 1, (int)pCols[i] - 1, true));
2702         }
2703
2704         newBool->setFromTriplets(tripletList.begin(), tripletList.end());
2705
2706         delete matrixBool;
2707         matrixBool = newBool;
2708         delete[] pRows;
2709
2710         m_iRows = _iNewRows;
2711         m_iCols = _iNewCols;
2712         m_iSize = _iNewRows * _iNewCols;
2713         m_piDims[0] = m_iRows;
2714         m_piDims[1] = m_iCols;
2715
2716         res = true;
2717
2718     }
2719     catch (...)
2720     {
2721         res = false;
2722     }
2723     return res;
2724 }
2725
2726 SparseBool* SparseBool::insert(typed_list* _pArgs, SparseBool* _pSource)
2727 {
2728     bool bNeedToResize  = false;
2729     int iDims           = (int)_pArgs->size();
2730     if (iDims > 2)
2731     {
2732         //sparse are only in 2 dims
2733         return NULL;
2734     }
2735
2736     typed_list pArg;
2737
2738     int piMaxDim[2];
2739     int piCountDim[2];
2740
2741     //on case of resize
2742     int iNewRows    = 0;
2743     int iNewCols    = 0;
2744
2745     //evaluate each argument and replace by appropriate value and compute the count of combinations
2746     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
2747     if (iSeqCount == 0)
2748     {
2749         //free pArg content
2750         cleanIndexesArguments(_pArgs, &pArg);
2751         return this;
2752     }
2753
2754     if (iDims < 2)
2755     {
2756         //see as vector
2757         if (getRows() == 1 || getCols() == 1)
2758         {
2759             //vector or scalar
2760             if (getSize() < piMaxDim[0])
2761             {
2762                 bNeedToResize = true;
2763
2764                 //need to enlarge sparse dimensions
2765                 if (getCols() == 1 || getSize() == 0)
2766                 {
2767                     //column vector
2768                     iNewRows    = piMaxDim[0];
2769                     iNewCols    = 1;
2770                 }
2771                 else if (getRows() == 1)
2772                 {
2773                     //row vector
2774                     iNewRows    = 1;
2775                     iNewCols    = piMaxDim[0];
2776                 }
2777             }
2778         }
2779         else if (getSize() < piMaxDim[0])
2780         {
2781             //free pArg content
2782             cleanIndexesArguments(_pArgs, &pArg);
2783             //out of range
2784             return NULL;
2785         }
2786     }
2787     else
2788     {
2789         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
2790         {
2791             bNeedToResize = true;
2792             iNewRows = std::max(getRows(), piMaxDim[0]);
2793             iNewCols = std::max(getCols(), piMaxDim[1]);
2794         }
2795     }
2796
2797     //check number of insertion
2798     if (_pSource->isScalar() == false && _pSource->getSize() != iSeqCount)
2799     {
2800         //free pArg content
2801         cleanIndexesArguments(_pArgs, &pArg);
2802         return NULL;
2803     }
2804
2805     //now you are sure to be able to insert values
2806     if (bNeedToResize)
2807     {
2808         if (resize(iNewRows, iNewCols) == false)
2809         {
2810             //free pArg content
2811             cleanIndexesArguments(_pArgs, &pArg);
2812             return NULL;
2813         }
2814     }
2815
2816     if (iDims == 1)
2817     {
2818         double* pIdx = pArg[0]->getAs<Double>()->get();
2819         for (int i = 0 ; i < iSeqCount ; i++)
2820         {
2821             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
2822             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
2823
2824             if (_pSource->isScalar())
2825             {
2826                 set(iRow, iCol, _pSource->get(0, 0), false);
2827             }
2828             else
2829             {
2830                 int iRowOrig = i % _pSource->getRows();
2831                 int iColOrig = i / _pSource->getRows();
2832                 set(iRow, iCol, _pSource->get(iRowOrig, iColOrig), false);
2833             }
2834         }
2835     }
2836     else
2837     {
2838         double* pIdxRow = pArg[0]->getAs<Double>()->get();
2839         int iRowSize    = pArg[0]->getAs<Double>()->getSize();
2840         double* pIdxCol = pArg[1]->getAs<Double>()->get();
2841
2842         for (int i = 0 ; i < iSeqCount ; i++)
2843         {
2844             if (_pSource->isScalar())
2845             {
2846                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(0, 0), false);
2847             }
2848             else
2849             {
2850                 int iRowOrig = i % _pSource->getRows();
2851                 int iColOrig = i / _pSource->getRows();
2852                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(iRowOrig, iColOrig), false);
2853             }
2854         }
2855     }
2856
2857     finalize();
2858
2859     //free pArg content
2860     cleanIndexesArguments(_pArgs, &pArg);
2861
2862     return this;
2863 }
2864
2865 SparseBool* SparseBool::insert(typed_list* _pArgs, InternalType* _pSource)
2866 {
2867     bool bNeedToResize  = false;
2868     int iDims           = (int)_pArgs->size();
2869     if (iDims > 2)
2870     {
2871         //sparse are only in 2 dims
2872         return NULL;
2873     }
2874
2875     typed_list pArg;
2876
2877     int piMaxDim[2];
2878     int piCountDim[2];
2879
2880     //on case of resize
2881     int iNewRows    = 0;
2882     int iNewCols    = 0;
2883     Bool* pSource = _pSource->getAs<Bool>();
2884
2885     //evaluate each argument and replace by appropriate value and compute the count of combinations
2886     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
2887     if (iSeqCount == 0)
2888     {
2889         //free pArg content
2890         cleanIndexesArguments(_pArgs, &pArg);
2891         return this;
2892     }
2893
2894     if (iDims < 2)
2895     {
2896         //see as vector
2897         if (getRows() == 1 || getCols() == 1)
2898         {
2899             //vector or scalar
2900             bNeedToResize = true;
2901             if (getSize() < piMaxDim[0])
2902             {
2903                 //need to enlarge sparse dimensions
2904                 if (getCols() == 1 || getSize() == 0)
2905                 {
2906                     //column vector
2907                     iNewRows    = piMaxDim[0];
2908                     iNewCols    = 1;
2909                 }
2910                 else if (getRows() == 1)
2911                 {
2912                     //row vector
2913                     iNewRows    = 1;
2914                     iNewCols    = piMaxDim[0];
2915                 }
2916             }
2917         }
2918         else if (getSize() < piMaxDim[0])
2919         {
2920             //free pArg content
2921             cleanIndexesArguments(_pArgs, &pArg);
2922             //out of range
2923             return NULL;
2924         }
2925     }
2926     else
2927     {
2928         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
2929         {
2930             bNeedToResize = true;
2931             iNewRows = std::max(getRows(), piMaxDim[0]);
2932             iNewCols = std::max(getCols(), piMaxDim[1]);
2933         }
2934     }
2935
2936     //check number of insertion
2937     if (pSource->isScalar() == false && pSource->getSize() != iSeqCount)
2938     {
2939         //free pArg content
2940         cleanIndexesArguments(_pArgs, &pArg);
2941         return NULL;
2942     }
2943
2944     //now you are sure to be able to insert values
2945     if (bNeedToResize)
2946     {
2947         if (resize(iNewRows, iNewCols) == false)
2948         {
2949             //free pArg content
2950             cleanIndexesArguments(_pArgs, &pArg);
2951             return NULL;
2952         }
2953     }
2954
2955     if (iDims == 1)
2956     {
2957         double* pIdx = pArg[0]->getAs<Double>()->get();
2958         for (int i = 0 ; i < iSeqCount ; i++)
2959         {
2960             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
2961             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
2962             if (pSource->isScalar())
2963             {
2964                 set(iRow, iCol, pSource->get(0) != 0, false);
2965             }
2966             else
2967             {
2968                 set(iRow, iCol, pSource->get(i) != 0, false);
2969             }
2970         }
2971     }
2972     else
2973     {
2974         double* pIdxRow = pArg[0]->getAs<Double>()->get();
2975         int iRowSize    = pArg[0]->getAs<Double>()->getSize();
2976         double* pIdxCol = pArg[1]->getAs<Double>()->get();
2977
2978         for (int i = 0 ; i < iSeqCount ; i++)
2979         {
2980             if (pSource->isScalar())
2981             {
2982                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(0) != 0, false);
2983             }
2984             else
2985             {
2986                 int iRowOrig = i % pSource->getRows();
2987                 int iColOrig = i / pSource->getRows();
2988
2989                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(iRowOrig, iColOrig) != 0, false);
2990             }
2991         }
2992     }
2993
2994     finalize();
2995
2996     //free pArg content
2997     cleanIndexesArguments(_pArgs, &pArg);
2998     return this;
2999 }
3000
3001 SparseBool* SparseBool::remove(typed_list* _pArgs)
3002 {
3003     SparseBool* pOut = NULL;
3004     int iDims = (int)_pArgs->size();
3005     if (iDims > 2)
3006     {
3007         //sparse are only in 2 dims
3008         return NULL;
3009     }
3010
3011     typed_list pArg;
3012
3013     int piMaxDim[2];
3014     int piCountDim[2];
3015
3016     //evaluate each argument and replace by appropriate value and compute the count of combinations
3017     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3018     if (iSeqCount == 0)
3019     {
3020         //free pArg content
3021         cleanIndexesArguments(_pArgs, &pArg);
3022         return this;
3023     }
3024
3025     bool* pbFull = new bool[iDims];
3026     //coord must represent all values on a dimension
3027     for (int i = 0 ; i < iDims ; i++)
3028     {
3029         pbFull[i]       = false;
3030         int iDimToCheck = getVarMaxDim(i, iDims);
3031         int iIndexSize  = pArg[i]->getAs<GenericType>()->getSize();
3032
3033         //we can have index more than once
3034         if (iIndexSize >= iDimToCheck)
3035         {
3036             //size is good, now check datas
3037             double* pIndexes = getDoubleArrayFromDouble(pArg[i]);
3038             for (int j = 0 ; j < iDimToCheck ; j++)
3039             {
3040                 bool bFind = false;
3041                 for (int k = 0 ; k < iIndexSize ; k++)
3042                 {
3043                     if ((int)pIndexes[k] == j + 1)
3044                     {
3045                         bFind = true;
3046                         break;
3047                     }
3048                 }
3049                 pbFull[i]  = bFind;
3050             }
3051         }
3052     }
3053
3054     //only one dims can be not full/entire
3055     bool bNotEntire = false;
3056     int iNotEntire  = 0;
3057     bool bTooMuchNotEntire = false;
3058     for (int i = 0 ; i < iDims ; i++)
3059     {
3060         if (pbFull[i] == false)
3061         {
3062             if (bNotEntire == false)
3063             {
3064                 bNotEntire = true;
3065                 iNotEntire = i;
3066             }
3067             else
3068             {
3069                 bTooMuchNotEntire = true;
3070                 break;
3071             }
3072         }
3073     }
3074
3075     if (bTooMuchNotEntire == true)
3076     {
3077         //free pArg content
3078         cleanIndexesArguments(_pArgs, &pArg);
3079         return NULL;
3080     }
3081
3082     delete[] pbFull;
3083
3084     //find index to keep
3085     int iNotEntireSize          = pArg[iNotEntire]->getAs<GenericType>()->getSize();
3086     double* piNotEntireIndex    = getDoubleArrayFromDouble(pArg[iNotEntire]);
3087     int iKeepSize               = getVarMaxDim(iNotEntire, iDims);
3088     bool* pbKeep                = new bool[iKeepSize];
3089
3090     //fill pbKeep with true value
3091     for (int i = 0 ; i < iKeepSize ; i++)
3092     {
3093         pbKeep[i] = true;
3094     }
3095
3096     for (int i = 0 ; i < iNotEntireSize ; i++)
3097     {
3098         int idx = (int)piNotEntireIndex[i] - 1;
3099
3100         //don't care of value out of bounds
3101         if (idx < iKeepSize)
3102         {
3103             pbKeep[idx] = false;
3104         }
3105     }
3106
3107     int iNewDimSize = 0;
3108     for (int i = 0 ; i < iKeepSize ; i++)
3109     {
3110         if (pbKeep[i] == true)
3111         {
3112             iNewDimSize++;
3113         }
3114     }
3115     delete[] pbKeep;
3116
3117     int* piNewDims = new int[iDims];
3118     for (int i = 0 ; i < iDims ; i++)
3119     {
3120         if (i == iNotEntire)
3121         {
3122             piNewDims[i] = iNewDimSize;
3123         }
3124         else
3125         {
3126             piNewDims[i] = getVarMaxDim(i, iDims);
3127         }
3128     }
3129
3130     //remove last dimension if are == 1
3131     int iOrigDims = iDims;
3132     for (int i = (iDims - 1) ; i >= 2 ; i--)
3133     {
3134         if (piNewDims[i] == 1)
3135         {
3136             iDims--;
3137         }
3138         else
3139         {
3140             break;
3141         }
3142     }
3143
3144     if (iDims == 1)
3145     {
3146         if (iNewDimSize == 0)
3147         {
3148             //free pArg content
3149             cleanIndexesArguments(_pArgs, &pArg);
3150             return new SparseBool(0, 0);
3151         }
3152         else
3153         {
3154             //two cases, depends of original matrix/vector
3155             if ((*_pArgs)[0]->isColon() == false && m_iDims == 2 && m_piDims[0] == 1 && m_piDims[1] != 1)
3156             {
3157                 //special case for row vector
3158                 pOut = new SparseBool(1, iNewDimSize);
3159                 //in this case we have to care of 2nd dimension
3160                 //iNotEntire = 1;
3161             }
3162             else
3163             {
3164                 pOut = new SparseBool(iNewDimSize, 1);
3165             }
3166         }
3167     }
3168     else
3169     {
3170         pOut = new SparseBool(piNewDims[0], piNewDims[0]);
3171     }
3172
3173     delete[] piNewDims;
3174     //find a way to copy existing data to new variable ...
3175     int iNewPos = 0;
3176     int* piIndexes = new int[iOrigDims];
3177     int* piViewDims = new int[iOrigDims];
3178     for (int i = 0 ; i < iOrigDims ; i++)
3179     {
3180         piViewDims[i] = getVarMaxDim(i, iOrigDims);
3181     }
3182
3183     for (int i = 0 ; i < getSize() ; i++)
3184     {
3185         bool bByPass = false;
3186         getIndexesWithDims(i, piIndexes, piViewDims, iOrigDims);
3187
3188         //check if piIndexes use removed indexes
3189         for (int j = 0 ; j < iNotEntireSize ; j++)
3190         {
3191             if ((piNotEntireIndex[j] - 1) == piIndexes[iNotEntire])
3192             {
3193                 //by pass this value
3194                 bByPass = true;
3195                 break;
3196             }
3197         }
3198
3199         if (bByPass == false)
3200         {
3201             //compute new index
3202             pOut->set(iNewPos, get(i));
3203             iNewPos++;
3204         }
3205     }
3206
3207     //free allocated data
3208     for (int i = 0 ; i < iDims ; i++)
3209     {
3210         if (pArg[i] != (*_pArgs)[i])
3211         {
3212             delete pArg[i];
3213         }
3214     }
3215
3216     delete[] piIndexes;
3217     delete[] piViewDims;
3218
3219     //free pArg content
3220     cleanIndexesArguments(_pArgs, &pArg);
3221
3222     return pOut;
3223 }
3224
3225 bool SparseBool::append(int r, int c, SparseBool SPARSE_CONST* src)
3226 {
3227     doAppend(*src, r, c, *matrixBool);
3228     finalize();
3229     return true;
3230 }
3231
3232 InternalType* SparseBool::insertNew(typed_list* _pArgs, InternalType* _pSource)
3233 {
3234     typed_list pArg;
3235     InternalType *pOut  = NULL;
3236     SparseBool* pSource = _pSource->getAs<SparseBool>();
3237
3238     int iDims           = (int)_pArgs->size();
3239     int* piMaxDim       = new int[iDims];
3240     int* piCountDim     = new int[iDims];
3241     bool bUndefine      = false;
3242
3243     //evaluate each argument and replace by appropriate value and compute the count of combinations
3244     int iSeqCount = checkIndexesArguments(NULL, _pArgs, &pArg, piMaxDim, piCountDim);
3245     if (iSeqCount == 0)
3246     {
3247         //free pArg content
3248         cleanIndexesArguments(_pArgs, &pArg);
3249         return createEmptyDouble();
3250     }
3251
3252     if (iSeqCount < 0)
3253     {
3254         iSeqCount = -iSeqCount;
3255         bUndefine = true;
3256     }
3257
3258     if (bUndefine)
3259     {
3260         //manage : and $ in creation by insertion
3261         int iSource         = 0;
3262         int *piSourceDims   = pSource->getDimsArray();
3263
3264         for (int i = 0 ; i < iDims ; i++)
3265         {
3266             if (pArg[i] == NULL)
3267             {
3268                 //undefine value
3269                 if (pSource->isScalar())
3270                 {
3271                     piMaxDim[i]     = 1;
3272                     piCountDim[i]   = 1;
3273                 }
3274                 else
3275                 {
3276                     piMaxDim[i]     = piSourceDims[iSource];
3277                     piCountDim[i]   = piSourceDims[iSource];
3278                 }
3279                 iSource++;
3280                 //replace pArg value by the new one
3281                 pArg[i] = createDoubleVector(piMaxDim[i]);
3282             }
3283             //else
3284             //{
3285             //    piMaxDim[i] = piCountDim[i];
3286             //}
3287         }
3288     }
3289
3290     //remove last dimension at size 1
3291     //remove last dimension if are == 1
3292     for (int i = (iDims - 1) ; i >= 2 ; i--)
3293     {
3294         if (piMaxDim[i] == 1)
3295         {
3296             iDims--;
3297             pArg.pop_back();
3298         }
3299         else
3300         {
3301             break;
3302         }
3303     }
3304
3305     if (checkArgValidity(pArg) == false)
3306     {
3307         //free pArg content
3308         cleanIndexesArguments(_pArgs, &pArg);
3309         //contain bad index, like <= 0, ...
3310         return NULL;
3311     }
3312
3313     if (iDims == 1)
3314     {
3315         if (pSource->getCols() == 1)
3316         {
3317             pOut = new SparseBool(piCountDim[0], 1);
3318         }
3319         else
3320         {
3321             //rows == 1
3322             pOut = new SparseBool(1, piCountDim[0]);
3323         }
3324     }
3325     else
3326     {
3327         pOut = new SparseBool(piMaxDim[0], piMaxDim[1]);
3328     }
3329
3330     //fill with null item
3331     SparseBool* pSpOut = pOut->getAs<SparseBool>();
3332
3333     //insert values in new matrix
3334     InternalType* pOut2 = pSpOut->insert(&pArg, pSource);
3335     if (pOut != pOut2)
3336     {
3337         delete pOut;
3338     }
3339
3340     //free pArg content
3341     cleanIndexesArguments(_pArgs, &pArg);
3342
3343     return pOut2;
3344 }
3345
3346 SparseBool* SparseBool::extract(int nbCoords, int SPARSE_CONST* coords, int SPARSE_CONST* maxCoords, int SPARSE_CONST* resSize, bool asVector) SPARSE_CONST
3347 {
3348     if ( (asVector && maxCoords[0] > getSize()) ||
3349     (asVector == false && maxCoords[0] > getRows()) ||
3350     (asVector == false && maxCoords[1] > getCols()))
3351     {
3352         return 0;
3353     }
3354
3355     SparseBool * pSp (0);
3356     if (asVector)
3357     {
3358         pSp = (getRows() == 1) ?  new SparseBool(1, resSize[0]) : new SparseBool(resSize[0], 1);
3359         mycopy_n(makeMatrixIterator<bool>(*this,  Coords<true>(coords, getRows())), nbCoords
3360         , makeMatrixIterator<bool>(*(pSp->matrixBool), RowWiseFullIterator(pSp->getRows(), pSp->getCols())));
3361     }
3362     else
3363     {
3364         pSp = new SparseBool(resSize[0], resSize[1]);
3365         mycopy_n(makeMatrixIterator<bool>(*this,  Coords<false>(coords, getRows())), nbCoords
3366         , makeMatrixIterator<bool>(*(pSp->matrixBool), RowWiseFullIterator(pSp->getRows(), pSp->getCols())));
3367
3368     }
3369     return pSp;
3370 }
3371
3372 /*
3373 * create a new SparseBool of dims according to resSize and fill it from currentSparseBool (along coords)
3374 */
3375 InternalType* SparseBool::extract(typed_list* _pArgs)
3376 {
3377     SparseBool* pOut    = NULL;
3378     int iDims           = (int)_pArgs->size();
3379     typed_list pArg;
3380
3381     int* piMaxDim       = new int[iDims];
3382     int* piCountDim     = new int[iDims];
3383
3384     //evaluate each argument and replace by appropriate value and compute the count of combinations
3385     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3386     if (iSeqCount == 0)
3387     {
3388         //free pArg content
3389         cleanIndexesArguments(_pArgs, &pArg);
3390         if (_pArgs->size() == 0)
3391         {
3392             delete[] piMaxDim;
3393             delete[] piCountDim;
3394             //a()
3395             return this;
3396         }
3397         else
3398         {
3399             delete[] piMaxDim;
3400             delete[] piCountDim;
3401             //a([])
3402             return Double::Empty();
3403         }
3404     }
3405
3406     if (iDims < 2)
3407     {
3408         // Check that we stay inside the input size.
3409         if (piMaxDim[0] <= getSize())
3410         {
3411             int iNewRows = 0;
3412             int iNewCols = 0;
3413
3414             if (getRows() == 1 && getCols() != 1 && (*_pArgs)[0]->isColon() == false)
3415             {
3416                 //special case for row vector
3417                 iNewRows = 1;
3418                 iNewCols = piCountDim[0];
3419             }
3420             else
3421             {
3422                 iNewRows = piCountDim[0];
3423                 iNewCols = 1;
3424             }
3425
3426             pOut = new SparseBool(iNewRows, iNewCols);
3427             double* pIdx = pArg[0]->getAs<Double>()->get();
3428             // Write in output all elements extract from input.
3429             for (int i = 0 ; i < iSeqCount ; i++)
3430             {
3431                 if (pIdx[i] < 1)
3432                 {
3433                     delete pOut;
3434                     pOut = NULL;
3435                     break;
3436                 }
3437                 int iRowRead = static_cast<int>(pIdx[i] - 1) % getRows();
3438                 int iColRead = static_cast<int>(pIdx[i] - 1) / getRows();
3439
3440                 int iRowWrite = static_cast<int>(i) % iNewRows;
3441                 int iColWrite = static_cast<int>(i) / iNewRows;
3442
3443                 bool bValue = get(iRowRead, iColRead);
3444                 if (bValue)
3445                 {
3446                     //only non zero values
3447                     pOut->set(iRowWrite, iColWrite, true, false);
3448                 }
3449             }
3450         }
3451         else
3452         {
3453             delete[] piMaxDim;
3454             delete[] piCountDim;
3455             //free pArg content
3456             cleanIndexesArguments(_pArgs, &pArg);
3457             return NULL;
3458         }
3459     }
3460     else
3461     {
3462         // Check that we stay inside the input size.
3463         if (piMaxDim[0] <= getRows() && piMaxDim[1] <= getCols())
3464         {
3465             double* pIdxRow = pArg[0]->getAs<Double>()->get();
3466             double* pIdxCol = pArg[1]->getAs<Double>()->get();
3467
3468             int iNewRows = pArg[0]->getAs<Double>()->getSize();
3469             int iNewCols = pArg[1]->getAs<Double>()->getSize();
3470
3471             pOut = new SparseBool(iNewRows, iNewCols);
3472
3473             int iPos = 0;
3474             // Write in output all elements extract from input.
3475             for (int iRow = 0 ; iRow < iNewRows ; iRow++)
3476             {
3477                 for (int iCol = 0 ; iCol < iNewCols ; iCol++)
3478                 {
3479                     if ((pIdxRow[iRow] < 1) || (pIdxCol[iCol] < 1))
3480                     {
3481                         delete pOut;
3482                         pOut = NULL;
3483                         break;
3484                     }
3485                     bool bValue = get((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
3486                     if (bValue)
3487                     {
3488                         //only non zero values
3489                         pOut->set(iRow, iCol, true, false);
3490                     }
3491                     iPos++;
3492                 }
3493             }
3494         }
3495         else
3496         {
3497             delete[] piMaxDim;
3498             delete[] piCountDim;
3499             //free pArg content
3500             cleanIndexesArguments(_pArgs, &pArg);
3501             return NULL;
3502         }
3503     }
3504
3505     finalize();
3506
3507     delete[] piMaxDim;
3508     delete[] piCountDim;
3509     //free pArg content
3510     cleanIndexesArguments(_pArgs, &pArg);
3511
3512     return pOut;
3513 }
3514
3515 std::size_t SparseBool::nbTrue() const
3516 {
3517     return  matrixBool->nonZeros() ;
3518 }
3519 std::size_t SparseBool::nbTrue(std::size_t r) const
3520 {
3521     int* piIndex = matrixBool->outerIndexPtr();
3522     return piIndex[r + 1] - piIndex[r];
3523 }
3524
3525
3526 void SparseBool::setTrue(bool finalize)
3527 {
3528     int rows = getRows();
3529     int cols = getCols();
3530
3531     typedef Eigen::Triplet<bool> triplet;
3532     std::vector<triplet> tripletList;
3533
3534     for (int i = 0; i < rows; ++i)
3535     {
3536         for (int j = 0; j < cols; ++j)
3537         {
3538             tripletList.push_back(triplet(i, j, true));
3539         }
3540     }
3541
3542     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
3543
3544     if (finalize)
3545     {
3546         matrixBool->finalize();
3547     }
3548 }
3549
3550 void SparseBool::setFalse(bool finalize)
3551 {
3552     int rows = getRows();
3553     int cols = getCols();
3554
3555     typedef Eigen::Triplet<bool> triplet;
3556     std::vector<triplet> tripletList;
3557
3558     for (int i = 0; i < rows; ++i)
3559     {
3560         for (int j = 0; j < cols; ++j)
3561         {
3562             tripletList.push_back(triplet(i, j, false));
3563         }
3564     }
3565
3566     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
3567
3568     if (finalize)
3569     {
3570         matrixBool->finalize();
3571     }
3572 }
3573
3574 int* SparseBool::getNbItemByRow(int* _piNbItemByRows)
3575 {
3576     int* piNbItemByRows = new int[getRows() + 1];
3577     mycopy_n(matrixBool->outerIndexPtr(), getRows() + 1, piNbItemByRows);
3578
3579     for (int i = 0 ; i < getRows() ; i++)
3580     {
3581         _piNbItemByRows[i] = piNbItemByRows[i + 1] - piNbItemByRows[i];
3582     }
3583
3584     delete[] piNbItemByRows;
3585     return _piNbItemByRows;
3586 }
3587
3588 int* SparseBool::getColPos(int* _piColPos)
3589 {
3590     mycopy_n(matrixBool->innerIndexPtr(), nbTrue(), _piColPos);
3591     for (int i = 0; i < nbTrue(); i++)
3592     {
3593         _piColPos[i]++;
3594     }
3595
3596     return _piColPos;
3597 }
3598
3599 int* SparseBool::outputRowCol(int* out)const
3600 {
3601     return sparseTransform(*matrixBool, sparseTransform(*matrixBool, out, GetRow<BoolSparse_t>()), GetCol<BoolSparse_t>());
3602 }
3603
3604 int* SparseBool::getInnerPtr(int* count)
3605 {
3606     *count = matrixBool->innerSize();
3607     return matrixBool->innerIndexPtr();
3608 }
3609
3610 int* SparseBool::getOuterPtr(int* count)
3611 {
3612     *count = matrixBool->outerSize();
3613     return matrixBool->outerIndexPtr();
3614 }
3615
3616
3617 bool SparseBool::operator==(const InternalType& it) SPARSE_CONST
3618 {
3619     SparseBool* otherSparse = const_cast<SparseBool*>(dynamic_cast<SparseBool const*>(&it));/* types::GenericType is not const-correct :( */
3620     return (otherSparse
3621     && (otherSparse->getRows() == getRows())
3622     && (otherSparse->getCols() == getCols())
3623     && equal(*matrixBool, *otherSparse->matrixBool));
3624 }
3625
3626 bool SparseBool::operator!=(const InternalType& it) SPARSE_CONST
3627 {
3628     return !(*this == it);
3629 }
3630
3631 void SparseBool::finalize()
3632 {
3633     matrixBool->prune(&keepForSparse<bool>);
3634     matrixBool->finalize();
3635 }
3636
3637 bool SparseBool::get(int r, int c) SPARSE_CONST
3638 {
3639     return matrixBool->coeff(r, c);
3640 }
3641
3642 bool SparseBool::set(int _iRows, int _iCols, bool _bVal, bool _bFinalize) SPARSE_CONST
3643 {
3644     matrixBool->coeffRef(_iRows, _iCols) = _bVal;
3645
3646     if (_bFinalize)
3647     {
3648         finalize();
3649     }
3650
3651     return true;
3652 }
3653
3654 void SparseBool::fill(Bool& dest, int r, int c) SPARSE_CONST
3655 {
3656     mycopy_n(makeMatrixIterator<bool >(*matrixBool, RowWiseFullIterator(getRows(), getCols())), getSize()
3657     , makeMatrixIterator<bool >(dest, RowWiseFullIterator(dest.getRows(), dest.getCols(), r, c)));
3658 }
3659
3660 Sparse* SparseBool::newOnes() const
3661 {
3662     return new Sparse(new types::Sparse::RealSparse_t(matrixBool->cast<double>()), 0);
3663 }
3664
3665 SparseBool* SparseBool::newNotEqualTo(SparseBool const&o) const
3666 {
3667     return cwiseOp<std::not_equal_to>(*this, o);
3668 }
3669
3670 SparseBool* SparseBool::newEqualTo(SparseBool& o)
3671 {
3672     int rowL = getRows();
3673     int colL = getCols();
3674
3675     int rowR = o.getRows();
3676     int colR = o.getCols();
3677     int row = std::max(rowL, rowR);
3678     int col = std::max(colL, colR);
3679
3680     //create a boolean sparse matrix with dims of sparses
3681     types::SparseBool* ret = new types::SparseBool(row, col);
3682
3683     if (isScalar() && o.isScalar())
3684     {
3685         bool l = get(0, 0);
3686         bool r = o.get(0, 0);
3687         ret->set(0, 0, l == r, false);
3688     }
3689     else if (isScalar())
3690     {
3691         int nnzR = static_cast<int>(o.nbTrue());
3692         std::vector<int> rowcolR(nnzR * 2, 0);
3693         o.outputRowCol(rowcolR.data());
3694
3695         //compare all items of R with R[0]
3696         bool l = get(0, 0);
3697         for (int i = 0; i < nnzR; ++i)
3698         {
3699             bool r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
3700             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
3701         }
3702     }
3703     else if (o.isScalar())
3704     {
3705         int nnzL = static_cast<int>(nbTrue());
3706         std::vector<int> rowcolL(nnzL * 2, 0);
3707         outputRowCol(rowcolL.data());
3708
3709         bool r = get(0, 0);
3710         for (int i = 0; i < nnzL; ++i)
3711         {
3712             bool l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
3713             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
3714         }
3715     }
3716     else
3717     {
3718         int nnzR = static_cast<int>(o.nbTrue());
3719         std::vector<int> rowcolR(nnzR * 2, 0);
3720         o.outputRowCol(rowcolR.data());
3721         int nnzL = static_cast<int>(nbTrue());
3722         std::vector<int> rowcolL(nnzL * 2, 0);
3723         outputRowCol(rowcolL.data());
3724         //set all values to %t
3725         ret->setTrue(false);
3726         //set %f in each pL values
3727         for (int i = 0; i < nnzL; ++i)
3728         {
3729             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, false, false);
3730         }
3731         ret->finalize();
3732
3733         //set _pR[i] == _pL[i] for each _pR values
3734         for (int i = 0; i < nnzR; ++i)
3735         {
3736             //get l and r following non zeros value of R
3737             bool l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
3738             bool r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
3739             //set value following non zeros value of R
3740             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
3741         }
3742     }
3743
3744     ret->finalize();
3745     return ret;
3746 }
3747
3748 SparseBool* SparseBool::newLogicalOr(SparseBool const&o) const
3749 {
3750     return cwiseOp<std::logical_or>(*this, o);
3751 }
3752
3753 SparseBool* SparseBool::newLogicalAnd(SparseBool const&o) const
3754 {
3755     return cwiseOp<std::logical_and>(*this, o);
3756 }
3757
3758 bool SparseBool::reshape(int* _piDims, int _iDims)
3759 {
3760     bool bOk = false;
3761     int iCols = 1;
3762
3763     if (_iDims == 2)
3764     {
3765         iCols = _piDims[1];
3766     }
3767
3768     if (_iDims <= 2)
3769     {
3770         bOk = reshape(_piDims[0], iCols);
3771     }
3772
3773     return bOk;
3774 }
3775
3776 bool SparseBool::reshape(int _iNewRows, int _iNewCols)
3777 {
3778     if (_iNewRows * _iNewCols != getRows() * getCols())
3779     {
3780         return false;
3781     }
3782
3783     bool res = false;
3784     try
3785     {
3786         //item count
3787         size_t iNonZeros = matrixBool->nonZeros();
3788         BoolSparse_t *newBool = new BoolSparse_t(_iNewRows, _iNewCols);
3789         newBool->reserve((int)iNonZeros);
3790
3791         //coords
3792         int* pRows = new int[iNonZeros * 2];
3793         outputRowCol(pRows);
3794         int* pCols = pRows + iNonZeros;
3795
3796         typedef Eigen::Triplet<bool> triplet;
3797         std::vector<triplet> tripletList;
3798
3799         for (size_t i = 0 ; i < iNonZeros ; i++)
3800         {
3801             int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
3802             tripletList.push_back(triplet((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), true));
3803         }
3804
3805         newBool->setFromTriplets(tripletList.begin(), tripletList.end());
3806
3807         delete matrixBool;
3808         matrixBool = newBool;
3809         delete[] pRows;
3810
3811         m_iRows = _iNewRows;
3812         m_iCols = _iNewCols;
3813         m_iSize = _iNewRows * _iNewCols;
3814
3815         m_iDims = 2;
3816         m_piDims[0] = m_iRows;
3817         m_piDims[1] = m_iCols;
3818
3819         finalize();
3820
3821         res = true;
3822     }
3823     catch (...)
3824     {
3825         res = false;
3826     }
3827     return res;
3828 }
3829
3830 bool SparseBool::transpose(InternalType *& out)
3831 {
3832     out = new SparseBool(new BoolSparse_t(matrixBool->transpose()));
3833     return true;
3834 }
3835
3836 template<typename T>
3837 void neg(const int r, const int c, const T * const in, Eigen::SparseMatrix<bool, 1> * const out)
3838 {
3839     for (int i = 0; i < r; i++)
3840     {
3841         for (int j = 0; j < c; j++)
3842         {
3843             out->coeffRef(i, j) = !in->coeff(i, j);
3844         }
3845     }
3846
3847     out->prune(&keepForSparse<bool>);
3848     out->finalize();
3849 }
3850
3851
3852 }