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