Fix some Coverity issues after 33b80b 50fe1d 5f73e9 and c473f9
[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(int *_piSize, int* _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 (iDims == 1)
1864     {
1865         if (iNewDimSize == 0)
1866         {
1867             delete[] piNewDims;
1868             //free pArg content
1869             cleanIndexesArguments(_pArgs, &pArg);
1870             return new Sparse(0, 0);
1871         }
1872         else
1873         {
1874             //two cases, depends of original matrix/vector
1875             if ((*_pArgs)[0]->isColon() == false && m_iDims == 2 && m_piDims[0] == 1 && m_piDims[1] != 1)
1876             {
1877                 //special case for row vector
1878                 pOut = new Sparse(1, iNewDimSize, isComplex());
1879                 //in this case we have to care of 2nd dimension
1880                 //iNotEntire = 1;
1881             }
1882             else
1883             {
1884                 pOut = new Sparse(iNewDimSize, 1, isComplex());
1885             }
1886         }
1887     }
1888     else
1889     {
1890         pOut = new Sparse(piNewDims[0], piNewDims[0], isComplex());
1891     }
1892
1893     delete[] piNewDims;
1894     //find a way to copy existing data to new variable ...
1895     int iNewPos = 0;
1896     int* piIndexes = new int[iOrigDims];
1897     int* piViewDims = new int[iOrigDims];
1898     for (int i = 0; i < iOrigDims; i++)
1899     {
1900         piViewDims[i] = getVarMaxDim(i, iOrigDims);
1901     }
1902
1903     for (int i = 0; i < getSize(); i++)
1904     {
1905         bool bByPass = false;
1906         getIndexesWithDims(i, piIndexes, piViewDims, iOrigDims);
1907
1908         //check if piIndexes use removed indexes
1909         for (int j = 0; j < iNotEntireSize; j++)
1910         {
1911             if ((piNotEntireIndex[j] - 1) == piIndexes[iNotEntire])
1912             {
1913                 //by pass this value
1914                 bByPass = true;
1915                 break;
1916             }
1917         }
1918
1919         if (bByPass == false)
1920         {
1921             //compute new index
1922             if (isComplex())
1923             {
1924                 pOut->set(iNewPos, getImg(i));
1925             }
1926             else
1927             {
1928                 pOut->set(iNewPos, get(i));
1929             }
1930             iNewPos++;
1931         }
1932     }
1933
1934     delete[] piIndexes;
1935     delete[] piViewDims;
1936
1937     //free pArg content
1938     cleanIndexesArguments(_pArgs, &pArg);
1939
1940     return pOut;
1941 }
1942
1943 Sparse* Sparse::append(int r, int c, types::Sparse SPARSE_CONST* src)
1944 {
1945     Sparse* pIT = checkRef(this, &Sparse::append, r, c, src);
1946     if (pIT != this)
1947     {
1948         return pIT;
1949     }
1950
1951     //        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;
1952     if (src->isComplex())
1953     {
1954         toComplex();
1955     }
1956     if (isComplex())
1957     {
1958         if (src->isComplex())
1959         {
1960             doAppend(*(src->matrixCplx), r, c, *matrixCplx);
1961         }
1962         else
1963         {
1964             doAppend(*(src->matrixReal), r, c, *matrixCplx);
1965         }
1966     }
1967     else
1968     {
1969         doAppend(*(src->matrixReal), r, c, *matrixReal);
1970     }
1971
1972     finalize();
1973
1974     return this; // realloc is meaningless for sparse matrices
1975 }
1976
1977 /*
1978 * create a new Sparse of dims according to resSize and fill it from currentSparse (along coords)
1979 */
1980 GenericType* Sparse::extract(typed_list* _pArgs)
1981 {
1982     Sparse* pOut = NULL;
1983     int iDims = (int)_pArgs->size();
1984     typed_list pArg;
1985
1986     int* piMaxDim = new int[iDims];
1987     int* piCountDim = new int[iDims];
1988
1989     //evaluate each argument and replace by appropriate value and compute the count of combinations
1990     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1991     if (iSeqCount == 0)
1992     {
1993         //free pArg content
1994         cleanIndexesArguments(_pArgs, &pArg);
1995         if (_pArgs->size() == 0)
1996         {
1997             //a()
1998             delete[] piMaxDim;
1999             delete[] piCountDim;
2000             //free pArg content
2001             cleanIndexesArguments(_pArgs, &pArg);
2002             return this;
2003         }
2004         else
2005         {
2006             //a([])
2007             delete[] piMaxDim;
2008             delete[] piCountDim;
2009             //free pArg content
2010             cleanIndexesArguments(_pArgs, &pArg);
2011             return Double::Empty();
2012         }
2013     }
2014
2015     if (iDims < 2)
2016     {
2017         if (piMaxDim[0] <= getSize())
2018         {
2019             int iNewRows = 0;
2020             int iNewCols = 0;
2021
2022             if (getRows() == 1 && getCols() != 1 && (*_pArgs)[0]->isColon() == false)
2023             {
2024                 //special case for row vector
2025                 iNewRows = 1;
2026                 iNewCols = piCountDim[0];
2027             }
2028             else
2029             {
2030                 iNewRows = piCountDim[0];
2031                 iNewCols = 1;
2032             }
2033
2034             double* pIdx = pArg[0]->getAs<Double>()->get();
2035             if (isComplex())
2036             {
2037                 CplxSparse_t* cplx = new CplxSparse_t(iNewRows, iNewCols);
2038                 std::vector<CplxTriplet_t> tripletList;
2039                 int row = getRows();
2040                 for (int i = 0; i < iSeqCount; i++)
2041                 {
2042                     int iRowRead = static_cast<int>(pIdx[i] - 1) % row;
2043                     int iColRead = static_cast<int>(pIdx[i] - 1) / row;
2044                     int iRowWrite = i % iNewRows;
2045                     int iColWrite = i / iNewRows;
2046
2047                     std::complex<double> dbl = getImg(iRowRead, iColRead);
2048                     if (dbl != 0.)
2049                     {
2050                         //only non zero values
2051                         tripletList.emplace_back(iRowWrite, iColWrite, dbl);
2052                     }
2053                 }
2054
2055                 cplx->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<std::complex<double>>());
2056                 pOut = new Sparse(nullptr, cplx);
2057             }
2058             else
2059             {
2060                 RealSparse_t* real = new RealSparse_t(iNewRows, iNewCols);
2061                 std::vector<RealTriplet_t> tripletList;
2062                 int row = getRows();
2063                 for (int i = 0; i < iSeqCount; i++)
2064                 {
2065                     int iRowRead = static_cast<int>(pIdx[i] - 1) % row;
2066                     int iColRead = static_cast<int>(pIdx[i] - 1) / row;
2067                     int iRowWrite = i % iNewRows;
2068                     int iColWrite = i / iNewRows;
2069
2070                     double dbl = get(iRowRead, iColRead);
2071                     if (dbl != 0)
2072                     {
2073                         //only non zero values
2074                         tripletList.emplace_back(iRowWrite, iColWrite, dbl);
2075                     }
2076                 }
2077
2078                 real->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<double>());
2079                 pOut = new Sparse(real, nullptr);
2080             }
2081         }
2082         else
2083         {
2084             delete[] piMaxDim;
2085             delete[] piCountDim;
2086             //free pArg content
2087             cleanIndexesArguments(_pArgs, &pArg);
2088             return NULL;
2089         }
2090     }
2091     else
2092     {
2093         if (piMaxDim[0] <= getRows() && piMaxDim[1] <= getCols())
2094         {
2095             double* pIdxRow = pArg[0]->getAs<Double>()->get();
2096             double* pIdxCol = pArg[1]->getAs<Double>()->get();
2097
2098             int iNewRows = pArg[0]->getAs<Double>()->getSize();
2099             int iNewCols = pArg[1]->getAs<Double>()->getSize();
2100
2101             if (isComplex())
2102             {
2103                 CplxSparse_t* cplx = new CplxSparse_t(iNewRows, iNewCols);
2104                 std::vector<CplxTriplet_t> tripletList;
2105                 for (int iRow = 0; iRow < iNewRows; iRow++)
2106                 {
2107                     for (int iCol = 0; iCol < iNewCols; iCol++)
2108                     {
2109                         std::complex<double> dbl = getImg((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
2110                         if (dbl != 0.)
2111                         {
2112                             //only non zero values
2113                             tripletList.emplace_back(iRow, iCol, dbl);
2114                         }
2115                     }
2116                 }
2117
2118                 cplx->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<std::complex<double>>());
2119                 pOut = new Sparse(nullptr, cplx);
2120             }
2121             else
2122             {
2123                 RealSparse_t* real = new RealSparse_t(iNewRows, iNewCols);
2124                 std::vector<RealTriplet_t> tripletList;
2125                 for (int iRow = 0; iRow < iNewRows; iRow++)
2126                 {
2127                     for (int iCol = 0; iCol < iNewCols; iCol++)
2128                     {
2129                         double dbl = get((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
2130                         if (dbl != 0.)
2131                         {
2132                             //only non zero values
2133                             tripletList.emplace_back(iRow, iCol, dbl);
2134                         }
2135                     }
2136                 }
2137
2138                 real->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<double>());
2139                 pOut = new Sparse(real, nullptr);
2140             }
2141         }
2142         else
2143         {
2144             delete[] piMaxDim;
2145             delete[] piCountDim;
2146             //free pArg content
2147             cleanIndexesArguments(_pArgs, &pArg);
2148             return NULL;
2149         }
2150     }
2151
2152     pOut->finalize();
2153
2154     delete[] piMaxDim;
2155     delete[] piCountDim;
2156     //free pArg content
2157     cleanIndexesArguments(_pArgs, &pArg);
2158
2159     return pOut;
2160 }
2161
2162 Sparse* Sparse::extract(int nbCoords, int SPARSE_CONST* coords, int SPARSE_CONST* maxCoords, int SPARSE_CONST* resSize, bool asVector) SPARSE_CONST
2163 {
2164     if ((asVector && maxCoords[0] > getSize()) ||
2165             (asVector == false && maxCoords[0] > getRows()) ||
2166             (asVector == false && maxCoords[1] > getCols()))
2167     {
2168         return 0;
2169     }
2170
2171     bool const cplx(isComplex());
2172     Sparse * pSp(0);
2173     if (asVector)
2174     {
2175         pSp = (getRows() == 1) ? new Sparse(1, resSize[0], cplx) : new Sparse(resSize[0], 1, cplx);
2176     }
2177     else
2178     {
2179         pSp = new Sparse(resSize[0], resSize[1], cplx);
2180     }
2181     //        std::cerr<<"extracted sparse:"<<pSp->getRows()<<", "<<pSp->getCols()<<"seqCount="<<nbCoords<<"maxDim="<<maxCoords[0] <<","<< maxCoords[1]<<std::endl;
2182     if (!(asVector
2183             ? copyToSparse(*this, Coords<true>(coords, getRows()), nbCoords
2184                            , *pSp, RowWiseFullIterator(pSp->getRows(), pSp->getCols()))
2185             : copyToSparse(*this, Coords<false>(coords), nbCoords
2186                            , *pSp, RowWiseFullIterator(pSp->getRows(), pSp->getCols()))))
2187     {
2188         delete pSp;
2189         pSp = NULL;
2190     }
2191     return pSp;
2192 }
2193
2194 bool Sparse::invoke(typed_list & in, optional_list & /*opt*/, int /*_iRetCount*/, typed_list & out, const ast::Exp & e)
2195 {
2196     if (in.size() == 0)
2197     {
2198         out.push_back(this);
2199     }
2200     else
2201     {
2202         InternalType * _out = extract(&in);
2203         if (!_out)
2204         {
2205             std::wostringstream os;
2206             os << _W("Invalid index.\n");
2207             throw ast::InternalError(os.str(), 999, e.getLocation());
2208         }
2209         out.push_back(_out);
2210     }
2211
2212     return true;
2213 }
2214
2215
2216 bool Sparse::isInvokable() const
2217 {
2218     return true;
2219 }
2220
2221 bool Sparse::hasInvokeOption() const
2222 {
2223     return false;
2224 }
2225
2226 int Sparse::getInvokeNbIn()
2227 {
2228     return -1;
2229 }
2230
2231 int Sparse::getInvokeNbOut()
2232 {
2233     return 1;
2234 }
2235
2236 /*
2237 coords are Scilab 1-based
2238 extract std::make_pair(coords, asVector), rowIter
2239 */
2240 template<typename Src, typename SrcTraversal, typename Sz, typename DestTraversal>
2241 bool Sparse::copyToSparse(Src SPARSE_CONST& src, SrcTraversal srcTrav, Sz n, Sparse& sp, DestTraversal destTrav)
2242 {
2243     if (!(src.isComplex() || sp.isComplex()))
2244     {
2245         mycopy_n(makeMatrixIterator<double>(src, srcTrav), n
2246                  , makeMatrixIterator<double>(*sp.matrixReal, destTrav));
2247     }
2248     else
2249     {
2250         sp.toComplex();
2251         mycopy_n(makeMatrixIterator<std::complex<double> >(src, srcTrav), n
2252                  , makeMatrixIterator<std::complex<double> >(*sp.matrixCplx, destTrav));
2253     }
2254
2255     sp.finalize();
2256     return true;
2257 }
2258
2259 // GenericType because we might return a Double* for scalar operand
2260 Sparse* Sparse::add(Sparse const& o) const
2261 {
2262     RealSparse_t* realSp(0);
2263     CplxSparse_t* cplxSp(0);
2264     if (isComplex() == false && o.isComplex() == false)
2265     {
2266         //R + R -> R
2267         realSp = new RealSparse_t(*matrixReal + * (o.matrixReal));
2268     }
2269     else if (isComplex() == false && o.isComplex() == true)
2270     {
2271         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >() + * (o.matrixCplx));
2272     }
2273     else if (isComplex() == true && o.isComplex() == false)
2274     {
2275         cplxSp = new CplxSparse_t(*matrixCplx + o.matrixReal->cast<std::complex<double> >());
2276     }
2277     else if (isComplex() == true && o.isComplex() == true)
2278     {
2279         cplxSp = new CplxSparse_t(*matrixCplx + * (o.matrixCplx));
2280     }
2281
2282     return new Sparse(realSp, cplxSp);
2283 }
2284
2285 Sparse* Sparse::substract(Sparse const& o) const
2286 {
2287     RealSparse_t* realSp(0);
2288     CplxSparse_t* cplxSp(0);
2289     if (isComplex() == false && o.isComplex() == false)
2290     {
2291         //R - R -> R
2292         realSp = new RealSparse_t(*matrixReal - * (o.matrixReal));
2293     }
2294     else if (isComplex() == false && o.isComplex() == true)
2295     {
2296         //R - C -> C
2297         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >() - * (o.matrixCplx));
2298     }
2299     else if (isComplex() == true && o.isComplex() == false)
2300     {
2301         //C - R -> C
2302         cplxSp = new CplxSparse_t(*matrixCplx - o.matrixReal->cast<std::complex<double> >());
2303     }
2304     else if (isComplex() == true && o.isComplex() == true)
2305     {
2306         //C - C -> C
2307         cplxSp = new CplxSparse_t(*matrixCplx - * (o.matrixCplx));
2308     }
2309
2310     return new Sparse(realSp, cplxSp);
2311 }
2312
2313 Sparse* Sparse::multiply(double s) const
2314 {
2315     return new Sparse(isComplex() ? 0 : new RealSparse_t((*matrixReal)*s)
2316                       , isComplex() ? new CplxSparse_t((*matrixCplx)*s) : 0);
2317 }
2318
2319 Sparse* Sparse::multiply(std::complex<double> s) const
2320 {
2321     return new Sparse(0
2322                       , isComplex() ? new CplxSparse_t((*matrixCplx) * s) : new CplxSparse_t((*matrixReal) * s));
2323 }
2324
2325 Sparse* Sparse::multiply(Sparse const& o) const
2326 {
2327     RealSparse_t* realSp(0);
2328     CplxSparse_t* cplxSp(0);
2329
2330     if (isComplex() == false && o.isComplex() == false)
2331     {
2332         realSp = new RealSparse_t(*matrixReal **(o.matrixReal));
2333     }
2334     else if (isComplex() == false && o.isComplex() == true)
2335     {
2336         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >() **(o.matrixCplx));
2337     }
2338     else if (isComplex() == true && o.isComplex() == false)
2339     {
2340         cplxSp = new CplxSparse_t(*matrixCplx * o.matrixReal->cast<std::complex<double> >());
2341     }
2342     else if (isComplex() == true && o.isComplex() == true)
2343     {
2344         cplxSp = new CplxSparse_t(*matrixCplx **(o.matrixCplx));
2345     }
2346
2347     return new Sparse(realSp, cplxSp);
2348 }
2349
2350 Sparse* Sparse::dotMultiply(Sparse SPARSE_CONST& o) const
2351 {
2352     RealSparse_t* realSp(0);
2353     CplxSparse_t* cplxSp(0);
2354     if (isComplex() == false && o.isComplex() == false)
2355     {
2356         realSp = new RealSparse_t(matrixReal->cwiseProduct(*(o.matrixReal)));
2357     }
2358     else if (isComplex() == false && o.isComplex() == true)
2359     {
2360         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >().cwiseProduct(*(o.matrixCplx)));
2361     }
2362     else if (isComplex() == true && o.isComplex() == false)
2363     {
2364         cplxSp = new CplxSparse_t(matrixCplx->cwiseProduct(o.matrixReal->cast<std::complex<double> >()));
2365     }
2366     else if (isComplex() == true && o.isComplex() == true)
2367     {
2368         cplxSp = new CplxSparse_t(matrixCplx->cwiseProduct(*(o.matrixCplx)));
2369     }
2370
2371     return new Sparse(realSp, cplxSp);
2372 }
2373
2374 Sparse* Sparse::dotDivide(Sparse SPARSE_CONST& o) const
2375 {
2376     RealSparse_t* realSp(0);
2377     CplxSparse_t* cplxSp(0);
2378     if (isComplex() == false && o.isComplex() == false)
2379     {
2380         realSp = new RealSparse_t(matrixReal->cwiseQuotient(*(o.matrixReal)));
2381     }
2382     else if (isComplex() == false && o.isComplex() == true)
2383     {
2384         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >().cwiseQuotient(*(o.matrixCplx)));
2385     }
2386     else if (isComplex() == true && o.isComplex() == false)
2387     {
2388         cplxSp = new CplxSparse_t(matrixCplx->cwiseQuotient(o.matrixReal->cast<std::complex<double> >()));
2389     }
2390     else if (isComplex() == true && o.isComplex() == true)
2391     {
2392         cplxSp = new CplxSparse_t(matrixCplx->cwiseQuotient(*(o.matrixCplx)));
2393     }
2394
2395     return new Sparse(realSp, cplxSp);
2396 }
2397
2398 int Sparse::newCholLLT(Sparse** _SpPermut, Sparse** _SpFactor) const
2399 {
2400     typedef Eigen::SparseMatrix<double, Eigen::ColMajor> RealSparseCol_t;
2401     RealSparseCol_t spColMajor = RealSparseCol_t((const RealSparse_t&) * matrixReal);
2402
2403     // Constructs and performs the LLT factorization of sparse
2404     Eigen::SimplicialLLT<RealSparseCol_t> pLLT(spColMajor);
2405     int iInfo = pLLT.info();
2406     if (iInfo != Eigen::Success)
2407     {
2408         *_SpFactor = NULL;
2409         *_SpPermut = NULL;
2410         return iInfo;
2411     }
2412
2413     // Get the lower matrix of factorization.
2414     // The new RealSparse_t will be set in Sparse without copy.
2415     *_SpFactor = new Sparse(new RealSparse_t(pLLT.matrixL()), NULL);
2416
2417     // Get the permutation matrix.
2418     Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic, int> p = pLLT.permutationP();
2419     *_SpPermut = new Sparse(static_cast<int>(p.rows()), static_cast<int>(p.cols()));
2420     for (int i = 0; i < p.rows(); i++)
2421     {
2422         (*_SpPermut)->set(i, p.indices()[i], 1, false);
2423     }
2424
2425     (*_SpPermut)->finalize();
2426
2427     return iInfo;
2428 }
2429
2430 bool Sparse::transpose(InternalType *& out)
2431 {
2432     out = new Sparse(matrixReal ? new RealSparse_t(matrixReal->transpose()) : 0, matrixCplx ? new CplxSparse_t(matrixCplx->transpose()) : 0);
2433     return true;
2434 }
2435
2436 bool Sparse::adjoint(InternalType *& out)
2437 {
2438     out = new Sparse(matrixReal ? new RealSparse_t(matrixReal->adjoint()) : 0, matrixCplx ? new CplxSparse_t(matrixCplx->adjoint()) : 0);
2439     return true;
2440 }
2441
2442 struct BoolCast
2443 {
2444     BoolCast(std::complex<double> const& c) : b(c.real() || c.imag()) {}
2445     operator bool() const
2446     {
2447         return b;
2448     }
2449     operator double() const
2450     {
2451         return b ? 1. : 0.;
2452     }
2453     bool b;
2454 };
2455 Sparse* Sparse::newOnes() const
2456 {
2457     // result is never cplx
2458     return new Sparse(matrixReal
2459                       ? new RealSparse_t(matrixReal->cast<bool>().cast<double>())
2460                       : new RealSparse_t(matrixCplx->cast<BoolCast>().cast<double>())
2461                       , 0);
2462 }
2463
2464 struct RealCast
2465 {
2466     RealCast(std::complex<double> const& c) : b(c.real()) {}
2467     operator bool() const
2468     {
2469         return b != 0;
2470     }
2471     operator double() const
2472     {
2473         return b;
2474     }
2475     double b;
2476 };
2477 Sparse* Sparse::newReal() const
2478 {
2479     return new Sparse(matrixReal
2480                       ? matrixReal
2481                       : new RealSparse_t(matrixCplx->cast<RealCast>().cast<double>())
2482                       , 0);
2483 }
2484
2485 std::size_t Sparse::nonZeros() const
2486 {
2487     if (isComplex())
2488     {
2489         return matrixCplx->nonZeros();
2490     }
2491     else
2492     {
2493         return matrixReal->nonZeros();
2494     }
2495 }
2496 std::size_t Sparse::nonZeros(std::size_t r) const
2497 {
2498     std::size_t res;
2499     if (matrixReal)
2500     {
2501         int* piIndex = matrixReal->outerIndexPtr();
2502         res = piIndex[r + 1] - piIndex[r];
2503     }
2504     else
2505     {
2506         int* piIndex = matrixCplx->outerIndexPtr();
2507         res = piIndex[r + 1] - piIndex[r];
2508     }
2509
2510     return res;
2511 }
2512
2513 int* Sparse::getNbItemByRow(int* _piNbItemByRows)
2514 {
2515     int* piNbItemByCols = new int[getRows() + 1];
2516     if (isComplex())
2517     {
2518         mycopy_n(matrixCplx->outerIndexPtr(), getRows() + 1, piNbItemByCols);
2519     }
2520     else
2521     {
2522         mycopy_n(matrixReal->outerIndexPtr(), getRows() + 1, piNbItemByCols);
2523     }
2524
2525     for (int i = 0; i < getRows(); i++)
2526     {
2527         _piNbItemByRows[i] = piNbItemByCols[i + 1] - piNbItemByCols[i];
2528     }
2529
2530     delete[] piNbItemByCols;
2531     return _piNbItemByRows;
2532 }
2533
2534 int* Sparse::getColPos(int* _piColPos)
2535 {
2536     if (isComplex())
2537     {
2538         mycopy_n(matrixCplx->innerIndexPtr(), nonZeros(), _piColPos);
2539     }
2540     else
2541     {
2542         mycopy_n(matrixReal->innerIndexPtr(), nonZeros(), _piColPos);
2543     }
2544
2545     for (size_t i = 0; i < nonZeros(); i++)
2546     {
2547         _piColPos[i]++;
2548     }
2549
2550     return _piColPos;
2551 }
2552
2553 int* Sparse::getInnerPtr(int* count)
2554 {
2555     int* ret = nullptr;
2556     if (isComplex())
2557     {
2558         ret = matrixCplx->innerIndexPtr();
2559         *count = static_cast<int>(matrixCplx->innerSize());
2560     }
2561     else
2562     {
2563         ret = matrixReal->innerIndexPtr();
2564         *count = static_cast<int>(matrixReal->innerSize());
2565     }
2566
2567     return ret;
2568 }
2569
2570 int* Sparse::getOuterPtr(int* count)
2571 {
2572     int* ret = nullptr;
2573     if (isComplex())
2574     {
2575         ret = matrixCplx->outerIndexPtr();
2576         *count = static_cast<int>(matrixCplx->outerSize());
2577     }
2578     else
2579     {
2580         ret = matrixReal->outerIndexPtr();
2581         *count = static_cast<int>(matrixReal->outerSize());
2582     }
2583
2584     return ret;
2585 }
2586
2587 namespace
2588 {
2589 template<typename S> struct GetReal : std::unary_function<typename S::InnerIterator, double>
2590 {
2591     double operator()(typename S::InnerIterator it) const
2592     {
2593         return it.value();
2594     }
2595 };
2596 template<> struct GetReal< Eigen::SparseMatrix<std::complex<double >, Eigen::RowMajor > >
2597     : std::unary_function<Sparse::CplxSparse_t::InnerIterator, double>
2598 {
2599     double operator()(Sparse::CplxSparse_t::InnerIterator it) const
2600     {
2601         return it.value().real();
2602     }
2603 };
2604 template<typename S> struct GetImag : std::unary_function<typename S::InnerIterator, double>
2605 {
2606     double operator()(typename S::InnerIterator it) const
2607     {
2608         return it.value().imag();
2609     }
2610 };
2611 template<typename S> struct GetRow : std::unary_function<typename S::InnerIterator, int>
2612 {
2613     int operator()(typename S::InnerIterator it) const
2614     {
2615         return static_cast<int>(it.row() + 1);
2616     }
2617 };
2618 template<typename S> struct GetCol : std::unary_function<typename S::InnerIterator, int>
2619 {
2620     int operator()(typename S::InnerIterator it) const
2621     {
2622         return static_cast<int>(it.col() + 1);
2623     }
2624 };
2625
2626 template<typename S, typename Out, typename F> Out sparseTransform(S& s, Out o, F f)
2627 {
2628     int size = static_cast<int>(s.outerSize());
2629     for (std::size_t k(0); k < size; ++k)
2630     {
2631         for (typename S::InnerIterator it(s, (int)k); it; ++it, ++o)
2632         {
2633             *o = f(it);
2634         }
2635     }
2636     return o;
2637 }
2638 }
2639
2640 std::pair<double*, double*> Sparse::outputValues(double* outReal, double* outImag)const
2641 {
2642     return matrixReal
2643            ? std::make_pair(sparseTransform(*matrixReal, outReal, GetReal<RealSparse_t>()), outImag)
2644            : std::make_pair(sparseTransform(*matrixCplx, outReal, GetReal<CplxSparse_t>())
2645                             , sparseTransform(*matrixCplx, outImag, GetImag<CplxSparse_t>()));
2646 }
2647
2648 int* Sparse::outputRowCol(int* out)const
2649 {
2650     return matrixReal
2651            ? sparseTransform(*matrixReal, sparseTransform(*matrixReal, out, GetRow<RealSparse_t>()), GetCol<RealSparse_t>())
2652            : sparseTransform(*matrixCplx, sparseTransform(*matrixCplx, out, GetRow<CplxSparse_t>()), GetCol<CplxSparse_t>());
2653 }
2654 double* Sparse::outputCols(double* out) const
2655 {
2656     if (isComplex())
2657     {
2658         mycopy_n(matrixCplx->innerIndexPtr(), nonZeros(), out);
2659     }
2660     else
2661     {
2662         mycopy_n(matrixReal->innerIndexPtr(), nonZeros(), out);
2663     }
2664
2665     return std::transform(out, out, out, std::bind2nd(std::plus<double>(), 1));
2666
2667 }
2668
2669 void Sparse::opposite(void)
2670 {
2671     if (isComplex())
2672     {
2673         std::complex<double>* data = matrixCplx->valuePtr();
2674         std::transform(data, data + matrixCplx->nonZeros(), data, std::negate<std::complex<double> >());
2675     }
2676     else
2677     {
2678         double* data = matrixReal->valuePtr();
2679         std::transform(data, data + matrixReal->nonZeros(), data, std::negate<double>());
2680     }
2681 }
2682
2683 SparseBool* Sparse::newLessThan(Sparse &o)
2684 {
2685     //only real values !
2686
2687     //return cwiseOp<std::less>(*this, o);
2688     int rowL = getRows();
2689     int colL = getCols();
2690
2691     int rowR = o.getRows();
2692     int colR = o.getCols();
2693     int row = std::max(rowL, rowR);
2694     int col = std::max(colL, colR);
2695
2696     //create a boolean sparse matrix with dims of sparses
2697     types::SparseBool* ret = new types::SparseBool(row, col);
2698     if (isScalar() && o.isScalar())
2699     {
2700         double l = get(0, 0);
2701         double r = o.get(0, 0);
2702         ret->set(0, 0, l < r, false);
2703     }
2704     else if (isScalar())
2705     {
2706         int nnzR = static_cast<int>(o.nonZeros());
2707         std::vector<int> rowcolR(nnzR * 2, 0);
2708         o.outputRowCol(rowcolR.data());
2709
2710         //compare all items of R with R[0]
2711         double l = get(0, 0);
2712         if (l < 0)
2713         {
2714             //set true
2715             ret->setTrue(false);
2716         }
2717
2718         for (int i = 0; i < nnzR; ++i)
2719         {
2720             double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2721             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l < r, false);
2722         }
2723     }
2724     else if (o.isScalar())
2725     {
2726         int nnzL = static_cast<int>(nonZeros());
2727         std::vector<int> rowcolL(nnzL * 2, 0);
2728         outputRowCol(rowcolL.data());
2729
2730         double r = o.get(0, 0);
2731         if (r > 0)
2732         {
2733             ret->setTrue(true);
2734         }
2735
2736         for (int i = 0; i < nnzL; ++i)
2737         {
2738             double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2739             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l < r, false);
2740         }
2741     }
2742     else
2743     {
2744         int nnzR = static_cast<int>(o.nonZeros());
2745         std::vector<int> rowcolR(nnzR * 2, 0);
2746         o.outputRowCol(rowcolR.data());
2747         int nnzL = static_cast<int>(nonZeros());
2748         std::vector<int> rowcolL(nnzL * 2, 0);
2749         outputRowCol(rowcolL.data());
2750         //set all values to %t
2751         ret->setFalse(false);
2752
2753         for (int i = 0; i < nnzL; ++i)
2754         {
2755             double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2756             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l < 0, false);
2757         }
2758         ret->finalize();
2759
2760         //set _pR[i] == _pL[i] for each _pR values
2761         for (int i = 0; i < nnzR; ++i)
2762         {
2763             //get l and r following non zeros value of R
2764             double l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2765             double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2766             //set value following non zeros value of R
2767             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l < r, false);
2768         }
2769     }
2770
2771     ret->finalize();
2772     return ret;
2773 }
2774
2775 SparseBool* Sparse::newNotEqualTo(Sparse const&o) const
2776 {
2777     return cwiseOp<std::not_equal_to>(*this, o);
2778 }
2779
2780 SparseBool* Sparse::newLessOrEqual(Sparse &o)
2781 {
2782     //only real values !
2783
2784     //return cwiseOp<std::less>(*this, o);
2785     int rowL = getRows();
2786     int colL = getCols();
2787
2788     int rowR = o.getRows();
2789     int colR = o.getCols();
2790     int row = std::max(rowL, rowR);
2791     int col = std::max(colL, colR);
2792
2793     //create a boolean sparse matrix with dims of sparses
2794     types::SparseBool* ret = new types::SparseBool(row, col);
2795     if (isScalar() && o.isScalar())
2796     {
2797         double l = get(0, 0);
2798         double r = o.get(0, 0);
2799         ret->set(0, 0, l <= r, false);
2800     }
2801     else if (isScalar())
2802     {
2803         int nnzR = static_cast<int>(o.nonZeros());
2804         std::vector<int> rowcolR(nnzR * 2, 0);
2805         o.outputRowCol(rowcolR.data());
2806
2807         //compare all items of R with R[0]
2808         double l = get(0, 0);
2809         if (l <= 0)
2810         {
2811             //set true
2812             ret->setTrue(false);
2813         }
2814
2815         for (int i = 0; i < nnzR; ++i)
2816         {
2817             double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2818             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l <= r, false);
2819         }
2820     }
2821     else if (o.isScalar())
2822     {
2823         int nnzL = static_cast<int>(nonZeros());
2824         std::vector<int> rowcolL(nnzL * 2, 0);
2825         outputRowCol(rowcolL.data());
2826
2827         double r = o.get(0, 0);
2828         if (r > 0)
2829         {
2830             ret->setTrue(true);
2831         }
2832
2833         for (int i = 0; i < nnzL; ++i)
2834         {
2835             double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2836             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l <= r, false);
2837         }
2838     }
2839     else
2840     {
2841         int nnzR = static_cast<int>(o.nonZeros());
2842         std::vector<int> rowcolR(nnzR * 2, 0);
2843         o.outputRowCol(rowcolR.data());
2844         int nnzL = static_cast<int>(nonZeros());
2845         std::vector<int> rowcolL(nnzL * 2, 0);
2846         outputRowCol(rowcolL.data());
2847         //set all values to %t
2848         ret->setTrue(false);
2849
2850         for (int i = 0; i < nnzL; ++i)
2851         {
2852             double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2853             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l <= 0, false);
2854         }
2855         ret->finalize();
2856
2857         //set _pR[i] == _pL[i] for each _pR values
2858         for (int i = 0; i < nnzR; ++i)
2859         {
2860             //get l and r following non zeros value of R
2861             double l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2862             double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2863             //set value following non zeros value of R
2864             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l <= r, false);
2865         }
2866     }
2867
2868     ret->finalize();
2869     return ret;
2870 }
2871
2872 SparseBool* Sparse::newEqualTo(Sparse &o)
2873 {
2874     int rowL = getRows();
2875     int colL = getCols();
2876
2877     int rowR = o.getRows();
2878     int colR = o.getCols();
2879     int row = std::max(rowL, rowR);
2880     int col = std::max(colL, colR);
2881
2882     //create a boolean sparse matrix with dims of sparses
2883     types::SparseBool* ret = new types::SparseBool(row, col);
2884     if (isScalar() && o.isScalar())
2885     {
2886         if (isComplex() || o.isComplex())
2887         {
2888             std::complex<double> l = getImg(0, 0);
2889             std::complex<double> r = o.getImg(0, 0);
2890             ret->set(0, 0, l == r, false);
2891         }
2892         else
2893         {
2894             double l = get(0, 0);
2895             double r = o.get(0, 0);
2896             ret->set(0, 0, l == r, false);
2897         }
2898     }
2899     else if (isScalar())
2900     {
2901         int nnzR = static_cast<int>(o.nonZeros());
2902         std::vector<int> rowcolR(nnzR * 2, 0);
2903         o.outputRowCol(rowcolR.data());
2904
2905         //compare all items of R with R[0]
2906         if (isComplex() || o.isComplex())
2907         {
2908             std::complex<double> l = getImg(0, 0);
2909             for (int i = 0; i < nnzR; ++i)
2910             {
2911                 std::complex<double> r = o.getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2912                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2913             }
2914         }
2915         else
2916         {
2917             double l = get(0, 0);
2918             for (int i = 0; i < nnzR; ++i)
2919             {
2920                 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2921                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2922             }
2923         }
2924     }
2925     else if (o.isScalar())
2926     {
2927         int nnzL = static_cast<int>(nonZeros());
2928         std::vector<int> rowcolL(nnzL * 2, 0);
2929         outputRowCol(rowcolL.data());
2930
2931         if (isComplex() || o.isComplex())
2932         {
2933             std::complex<double> r = o.getImg(0, 0);
2934             for (int i = 0; i < nnzL; ++i)
2935             {
2936                 std::complex<double> l = getImg(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2937                 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
2938             }
2939         }
2940         else
2941         {
2942             double r = get(0, 0);
2943             for (int i = 0; i < nnzL; ++i)
2944             {
2945                 double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2946                 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
2947             }
2948         }
2949     }
2950     else
2951     {
2952         int nnzR = static_cast<int>(o.nonZeros());
2953         std::vector<int> rowcolR(nnzR * 2, 0);
2954         o.outputRowCol(rowcolR.data());
2955         int nnzL = static_cast<int>(nonZeros());
2956         std::vector<int> rowcolL(nnzL * 2, 0);
2957         outputRowCol(rowcolL.data());
2958         //set all values to %t
2959         ret->setTrue(false);
2960         //set %f in each pL values
2961         for (int i = 0; i < nnzL; ++i)
2962         {
2963             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, false, false);
2964         }
2965         ret->finalize();
2966
2967         //set _pR[i] == _pL[i] for each _pR values
2968         if (isComplex() || o.isComplex())
2969         {
2970             for (int i = 0; i < nnzR; ++i)
2971             {
2972                 //get l and r following non zeros value of R
2973                 std::complex<double> l = getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2974                 std::complex<double> r = o.getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2975                 //set value following non zeros value of R
2976                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2977             }
2978         }
2979         else
2980         {
2981             for (int i = 0; i < nnzR; ++i)
2982             {
2983                 //get l and r following non zeros value of R
2984                 double l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2985                 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2986                 //set value following non zeros value of R
2987                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2988             }
2989         }
2990     }
2991
2992     ret->finalize();
2993     return ret;
2994 }
2995
2996 Sparse* Sparse::reshape(int* _piDims, int _iDims)
2997 {
2998     Sparse* pSp = NULL;
2999     int iCols = 1;
3000
3001     if (_iDims == 2)
3002     {
3003         iCols = _piDims[1];
3004     }
3005
3006     if (_iDims <= 2)
3007     {
3008         pSp = reshape(_piDims[0], iCols);
3009     }
3010
3011     return pSp;
3012 }
3013
3014 Sparse* Sparse::reshape(int _iNewRows, int _iNewCols)
3015 {
3016     typedef Sparse* (Sparse::*reshape_t)(int, int);
3017     Sparse* pIT = checkRef(this, (reshape_t)&Sparse::reshape, _iNewRows, _iNewCols);
3018     if (pIT != this)
3019     {
3020         return pIT;
3021     }
3022
3023     if (_iNewRows * _iNewCols != getRows() * getCols())
3024     {
3025         return NULL;
3026     }
3027
3028     Sparse* res = NULL;
3029     try
3030     {
3031         if (matrixReal)
3032         {
3033             //item count
3034             size_t iNonZeros = nonZeros();
3035             RealSparse_t *newReal = new RealSparse_t(_iNewRows, _iNewCols);
3036             newReal->reserve((int)iNonZeros);
3037
3038             //coords
3039             int* pRows = new int[iNonZeros * 2];
3040             outputRowCol(pRows);
3041             int* pCols = pRows + iNonZeros;
3042
3043             //values
3044             double* pNonZeroR = new double[iNonZeros];
3045             double* pNonZeroI = new double[iNonZeros];
3046             outputValues(pNonZeroR, pNonZeroI);
3047
3048             std::vector<RealTriplet_t> tripletList;
3049             for (size_t i = 0; i < iNonZeros; i++)
3050             {
3051                 int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
3052                 tripletList.emplace_back((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), pNonZeroR[i]);
3053             }
3054
3055             newReal->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<double>());
3056
3057             delete matrixReal;
3058             matrixReal = newReal;
3059             delete[] pRows;
3060             delete[] pNonZeroR;
3061             delete[] pNonZeroI;
3062         }
3063         else
3064         {
3065             //item count
3066             size_t iNonZeros = nonZeros();
3067             CplxSparse_t *newCplx = new CplxSparse_t(_iNewRows, _iNewCols);
3068             newCplx->reserve((int)iNonZeros);
3069
3070             //coords
3071             int* pRows = new int[iNonZeros * 2];
3072             outputRowCol(pRows);
3073             int* pCols = pRows + iNonZeros;
3074
3075             //values
3076             double* pNonZeroR = new double[iNonZeros];
3077             double* pNonZeroI = new double[iNonZeros];
3078             outputValues(pNonZeroR, pNonZeroI);
3079
3080             std::vector<CplxTriplet_t> tripletList;
3081
3082             for (size_t i = 0; i < iNonZeros; i++)
3083             {
3084                 int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
3085                 tripletList.emplace_back((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), std::complex<double>(pNonZeroR[i], pNonZeroI[i]));
3086             }
3087
3088             newCplx->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<std::complex<double>>());
3089
3090             delete matrixCplx;
3091             matrixCplx = newCplx;
3092             delete[] pRows;
3093             delete[] pNonZeroR;
3094             delete[] pNonZeroI;
3095         }
3096
3097         m_iRows = _iNewRows;
3098         m_iCols = _iNewCols;
3099         m_iSize = _iNewRows * _iNewCols;
3100
3101         m_iDims = 2;
3102         m_piDims[0] = m_iRows;
3103         m_piDims[1] = m_iCols;
3104
3105         finalize();
3106
3107         res = this;
3108     }
3109     catch (...)
3110     {
3111         res = NULL;
3112     }
3113     return res;
3114 }
3115
3116 //    SparseBool* SparseBool::new
3117
3118 SparseBool::SparseBool(Bool SPARSE_CONST& src)
3119 {
3120     //compute idx
3121     int size = src.getSize();
3122     int row = src.getRows();
3123     Double* idx = new Double(src.getSize(), 2);
3124     double* p = idx->get();
3125     for (int i = 0; i < size; ++i)
3126     {
3127         p[i] = (double)(i % row) + 1;
3128         p[i + size] = (double)(i / row) + 1;
3129     }
3130     create2(src.getRows(), src.getCols(), src, *idx);
3131     idx->killMe();
3132 #ifndef NDEBUG
3133     Inspector::addItem(this);
3134 #endif
3135 }
3136 /* @param src : Bool matrix to copy into a new sparse matrix
3137 @param idx : Double matrix to use as indexes to get values from the src
3138 **/
3139 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
3140 {
3141     int idxrow = idx.getRows();
3142     int rows = static_cast<int>(*std::max_element(idx.get(), idx.get() + idxrow));
3143     int cols = static_cast<int>(*std::max_element(idx.get() + idxrow, idx.get() + idxrow * 2));
3144     create2(rows, cols, src, idx);
3145 #ifndef NDEBUG
3146     Inspector::addItem(this);
3147 #endif
3148 }
3149
3150 /* @param src : Bool matrix to copy into a new sparse matrix
3151 @param idx : Double matrix to use as indexes to get values from the src
3152 @param dims : Double matrix containing the dimensions of the new matrix
3153 **/
3154 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx, Double SPARSE_CONST& dims)
3155 {
3156     create2(static_cast<int>(dims.get(0)), static_cast<int>(dims.get(1)), src, idx);
3157 #ifndef NDEBUG
3158     Inspector::addItem(this);
3159 #endif
3160 }
3161
3162 SparseBool::SparseBool(int _iRows, int _iCols) : matrixBool(new BoolSparse_t(_iRows, _iCols))
3163 {
3164     m_iRows = _iRows;
3165     m_iCols = _iCols;
3166     m_iSize = _iRows * _iCols;
3167     m_iDims = 2;
3168     m_piDims[0] = _iRows;
3169     m_piDims[1] = _iCols;
3170 #ifndef NDEBUG
3171     Inspector::addItem(this);
3172 #endif
3173 }
3174
3175 SparseBool::SparseBool(SparseBool const& src) : matrixBool(new BoolSparse_t(*src.matrixBool))
3176 {
3177     m_iDims = 2;
3178     m_iRows = const_cast<SparseBool*>(&src)->getRows();
3179     m_iCols = const_cast<SparseBool*>(&src)->getCols();
3180     m_iSize = m_iRows * m_iCols;
3181     m_piDims[0] = m_iRows;
3182     m_piDims[1] = m_iCols;
3183 #ifndef NDEBUG
3184     Inspector::addItem(this);
3185 #endif
3186 }
3187
3188 SparseBool::SparseBool(BoolSparse_t* src) : matrixBool(src)
3189 {
3190     m_iRows = static_cast<int>(src->rows());
3191     m_iCols = static_cast<int>(src->cols());
3192     m_iSize = m_iRows * m_iCols;
3193     m_iDims = 2;
3194     m_piDims[0] = m_iRows;
3195     m_piDims[1] = m_iCols;
3196 #ifndef NDEBUG
3197     Inspector::addItem(this);
3198 #endif
3199 }
3200
3201 SparseBool::SparseBool(int rows, int cols, int trues, int* inner, int* outer)
3202 {
3203     int* out = nullptr;
3204     int* in = nullptr;
3205
3206     matrixBool = new BoolSparse_t(rows, cols);
3207     matrixBool->reserve((int)trues);
3208     out = matrixBool->outerIndexPtr();
3209     in = matrixBool->innerIndexPtr();
3210
3211     //update outerIndexPtr
3212     memcpy(out, outer, sizeof(int) * (rows + 1));
3213     //update innerIndexPtr
3214     memcpy(in, inner, sizeof(int) * trues);
3215
3216     bool* data = matrixBool->valuePtr();
3217     for (int i = 0; i < trues; ++i)
3218     {
3219         data[i] = true;
3220     }
3221
3222     m_iCols = cols;
3223     m_iRows = rows;
3224     m_iSize = cols * rows;
3225     m_iDims = 2;
3226     m_piDims[0] = m_iRows;
3227     m_piDims[1] = m_iCols;
3228
3229     matrixBool->resizeNonZeros(trues);
3230 }
3231
3232 bool SparseBool::getMemory(int *_piSize, int* _piSizePlusType)
3233 {
3234     *_piSize = nbTrue() * sizeof(bool);
3235     *_piSizePlusType = *_piSize + sizeof(*this);
3236     return true;
3237 }
3238
3239 void SparseBool::create2(int rows, int cols, Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
3240 {
3241     int nnz = src.getSize();
3242     double* i = idx.get();
3243     double* j = i + idx.getRows();
3244     int* val = src.get();
3245
3246     std::vector<BoolTriplet_t> tripletList;
3247     tripletList.reserve((int)nnz);
3248
3249     for (int k = 0; k < nnz; ++k)
3250     {
3251         tripletList.emplace_back(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, val[k] == 1);
3252     }
3253
3254     matrixBool = new BoolSparse_t(rows, cols);
3255     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
3256
3257     m_iRows = static_cast<int>(matrixBool->rows());
3258     m_iCols = static_cast<int>(matrixBool->cols());
3259     m_iSize = cols * rows;
3260     m_iDims = 2;
3261     m_piDims[0] = m_iRows;
3262     m_piDims[1] = m_iCols;
3263     finalize();
3264 }
3265
3266 SparseBool::~SparseBool()
3267 {
3268     delete matrixBool;
3269 #ifndef NDEBUG
3270     Inspector::removeItem(this);
3271 #endif
3272 }
3273
3274 bool SparseBool::toString(std::wostringstream& ostr)
3275 {
3276     ostr << ::toString(*matrixBool, 0);
3277     return true;
3278 }
3279
3280 void SparseBool::whoAmI() SPARSE_CONST
3281 {
3282     std::cout << "types::SparseBool";
3283 }
3284
3285 SparseBool* SparseBool::clone(void)
3286 {
3287     return new SparseBool(*this);
3288 }
3289
3290 SparseBool* SparseBool::resize(int _iNewRows, int _iNewCols)
3291 {
3292     typedef SparseBool* (SparseBool::*resize_t)(int, int);
3293     SparseBool* pIT = checkRef(this, (resize_t)&SparseBool::resize, _iNewRows, _iNewCols);
3294     if (pIT != this)
3295     {
3296         return pIT;
3297     }
3298
3299     if (_iNewRows <= getRows() && _iNewCols <= getCols())
3300     {
3301         //nothing to do: hence we do NOT fail
3302         return this;
3303     }
3304
3305     SparseBool* res = NULL;
3306     try
3307     {
3308         //item count
3309         size_t iNonZeros = nbTrue();
3310
3311         BoolSparse_t *newBool = new BoolSparse_t(_iNewRows, _iNewCols);
3312         newBool->reserve((int)iNonZeros);
3313
3314         //coords
3315         int* pRows = new int[iNonZeros * 2];
3316         outputRowCol(pRows);
3317         int* pCols = pRows + iNonZeros;
3318
3319         std::vector<BoolTriplet_t> tripletList;
3320
3321         for (size_t i = 0; i < iNonZeros; i++)
3322         {
3323             tripletList.emplace_back((int)pRows[i] - 1, (int)pCols[i] - 1, true);
3324         }
3325
3326         newBool->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<bool>());
3327
3328         delete matrixBool;
3329         matrixBool = newBool;
3330         delete[] pRows;
3331
3332         m_iRows = _iNewRows;
3333         m_iCols = _iNewCols;
3334         m_iSize = _iNewRows * _iNewCols;
3335         m_piDims[0] = m_iRows;
3336         m_piDims[1] = m_iCols;
3337
3338         res = this;
3339     }
3340     catch (...)
3341     {
3342         res = NULL;
3343     }
3344     return res;
3345 }
3346
3347 SparseBool* SparseBool::insert(typed_list* _pArgs, SparseBool* _pSource)
3348 {
3349     bool bNeedToResize = false;
3350     int iDims = (int)_pArgs->size();
3351     if (iDims > 2)
3352     {
3353         //sparse are only in 2 dims
3354         return NULL;
3355     }
3356
3357     typed_list pArg;
3358
3359     int piMaxDim[2];
3360     int piCountDim[2];
3361
3362     //on case of resize
3363     int iNewRows = 0;
3364     int iNewCols = 0;
3365
3366     //evaluate each argument and replace by appropriate value and compute the count of combinations
3367     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3368     if (iSeqCount == 0)
3369     {
3370         //free pArg content
3371         cleanIndexesArguments(_pArgs, &pArg);
3372         return this;
3373     }
3374
3375     if (iDims < 2)
3376     {
3377         //see as vector
3378         if (getRows() == 1 || getCols() == 1)
3379         {
3380             //vector or scalar
3381             if (getRows() * getCols() < piMaxDim[0])
3382             {
3383                 bNeedToResize = true;
3384
3385                 //need to enlarge sparse dimensions
3386                 if (getCols() == 1 || getSize() == 0)
3387                 {
3388                     //column vector
3389                     iNewRows = piMaxDim[0];
3390                     iNewCols = 1;
3391                 }
3392                 else if (getRows() == 1)
3393                 {
3394                     //row vector
3395                     iNewRows = 1;
3396                     iNewCols = piMaxDim[0];
3397                 }
3398             }
3399         }
3400         else if (getRows() * getCols() < piMaxDim[0])
3401         {
3402             //free pArg content
3403             cleanIndexesArguments(_pArgs, &pArg);
3404             //out of range
3405             return NULL;
3406         }
3407     }
3408     else
3409     {
3410         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
3411         {
3412             bNeedToResize = true;
3413             iNewRows = std::max(getRows(), piMaxDim[0]);
3414             iNewCols = std::max(getCols(), piMaxDim[1]);
3415         }
3416     }
3417
3418     //check number of insertion
3419     if (_pSource->isScalar() == false && _pSource->getSize() != iSeqCount)
3420     {
3421         //free pArg content
3422         cleanIndexesArguments(_pArgs, &pArg);
3423         return NULL;
3424     }
3425
3426     //now you are sure to be able to insert values
3427     if (bNeedToResize)
3428     {
3429         if (resize(iNewRows, iNewCols) == NULL)
3430         {
3431             //free pArg content
3432             cleanIndexesArguments(_pArgs, &pArg);
3433             return NULL;
3434         }
3435     }
3436
3437     if (iDims == 1)
3438     {
3439         double* pIdx = pArg[0]->getAs<Double>()->get();
3440         for (int i = 0; i < iSeqCount; i++)
3441         {
3442             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
3443             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
3444
3445             if (_pSource->isScalar())
3446             {
3447                 set(iRow, iCol, _pSource->get(0, 0), false);
3448             }
3449             else
3450             {
3451                 int iRowOrig = i % _pSource->getRows();
3452                 int iColOrig = i / _pSource->getRows();
3453                 set(iRow, iCol, _pSource->get(iRowOrig, iColOrig), false);
3454             }
3455         }
3456     }
3457     else
3458     {
3459         double* pIdxRow = pArg[0]->getAs<Double>()->get();
3460         int iRowSize = pArg[0]->getAs<Double>()->getSize();
3461         double* pIdxCol = pArg[1]->getAs<Double>()->get();
3462
3463         for (int i = 0; i < iSeqCount; i++)
3464         {
3465             if (_pSource->isScalar())
3466             {
3467                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(0, 0), false);
3468             }
3469             else
3470             {
3471                 int iRowOrig = i % _pSource->getRows();
3472                 int iColOrig = i / _pSource->getRows();
3473                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(iRowOrig, iColOrig), false);
3474             }
3475         }
3476     }
3477
3478     finalize();
3479
3480     //free pArg content
3481     cleanIndexesArguments(_pArgs, &pArg);
3482
3483     return this;
3484 }
3485
3486 SparseBool* SparseBool::insert(typed_list* _pArgs, InternalType* _pSource)
3487 {
3488     typedef SparseBool* (SparseBool::*insert_t)(typed_list*, InternalType*);
3489     SparseBool* pIT = checkRef(this, (insert_t)&SparseBool::insert, _pArgs, _pSource);
3490
3491
3492     if (pIT != this)
3493     {
3494         return pIT;
3495     }
3496
3497     if (_pSource->isSparseBool())
3498     {
3499         return insert(_pArgs, _pSource->getAs<SparseBool>());
3500     }
3501
3502     bool bNeedToResize  = false;
3503     int iDims           = (int)_pArgs->size();
3504     if (iDims > 2)
3505     {
3506         //sparse are only in 2 dims
3507         return NULL;
3508     }
3509
3510     typed_list pArg;
3511
3512     int piMaxDim[2];
3513     int piCountDim[2];
3514
3515     //on case of resize
3516     int iNewRows = 0;
3517     int iNewCols = 0;
3518     Bool* pSource = _pSource->getAs<Bool>();
3519
3520     //evaluate each argument and replace by appropriate value and compute the count of combinations
3521     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3522     if (iSeqCount == 0)
3523     {
3524         //free pArg content
3525         cleanIndexesArguments(_pArgs, &pArg);
3526         return this;
3527     }
3528
3529     if (iDims < 2)
3530     {
3531         //see as vector
3532         if (getRows() == 1 || getCols() == 1)
3533         {
3534             //vector or scalar
3535             bNeedToResize = true;
3536             if (getRows() * getCols() < piMaxDim[0])
3537             {
3538                 //need to enlarge sparse dimensions
3539                 if (getCols() == 1 || getRows() * getCols() == 0)
3540                 {
3541                     //column vector
3542                     iNewRows = piMaxDim[0];
3543                     iNewCols = 1;
3544                 }
3545                 else if (getRows() == 1)
3546                 {
3547                     //row vector
3548                     iNewRows = 1;
3549                     iNewCols = piMaxDim[0];
3550                 }
3551             }
3552         }
3553         else if (getRows() * getCols() < piMaxDim[0])
3554         {
3555             //free pArg content
3556             cleanIndexesArguments(_pArgs, &pArg);
3557             //out of range
3558             return NULL;
3559         }
3560     }
3561     else
3562     {
3563         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
3564         {
3565             bNeedToResize = true;
3566             iNewRows = std::max(getRows(), piMaxDim[0]);
3567             iNewCols = std::max(getCols(), piMaxDim[1]);
3568         }
3569     }
3570
3571     //check number of insertion
3572     if (pSource->isScalar() == false && pSource->getSize() != iSeqCount)
3573     {
3574         //free pArg content
3575         cleanIndexesArguments(_pArgs, &pArg);
3576         return NULL;
3577     }
3578
3579     //now you are sure to be able to insert values
3580     if (bNeedToResize)
3581     {
3582         if (resize(iNewRows, iNewCols) == NULL)
3583         {
3584             //free pArg content
3585             cleanIndexesArguments(_pArgs, &pArg);
3586             return NULL;
3587         }
3588     }
3589
3590     if (iDims == 1)
3591     {
3592         double* pIdx = pArg[0]->getAs<Double>()->get();
3593         for (int i = 0; i < iSeqCount; i++)
3594         {
3595             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
3596             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
3597             if (pSource->isScalar())
3598             {
3599                 set(iRow, iCol, pSource->get(0) != 0, false);
3600             }
3601             else
3602             {
3603                 set(iRow, iCol, pSource->get(i) != 0, false);
3604             }
3605         }
3606     }
3607     else
3608     {
3609         double* pIdxRow = pArg[0]->getAs<Double>()->get();
3610         int iRowSize = pArg[0]->getAs<Double>()->getSize();
3611         double* pIdxCol = pArg[1]->getAs<Double>()->get();
3612
3613         for (int i = 0; i < iSeqCount; i++)
3614         {
3615             if (pSource->isScalar())
3616             {
3617                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(0) != 0, false);
3618             }
3619             else
3620             {
3621                 int iRowOrig = i % pSource->getRows();
3622                 int iColOrig = i / pSource->getRows();
3623
3624                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(iRowOrig, iColOrig) != 0, false);
3625             }
3626         }
3627     }
3628
3629     finalize();
3630
3631     //free pArg content
3632     cleanIndexesArguments(_pArgs, &pArg);
3633     return this;
3634 }
3635
3636 GenericType* SparseBool::remove(typed_list* _pArgs)
3637 {
3638     SparseBool* pOut = NULL;
3639     int iDims = (int)_pArgs->size();
3640     if (iDims > 2)
3641     {
3642         //sparse are only in 2 dims
3643         return NULL;
3644     }
3645
3646     typed_list pArg;
3647
3648     int piMaxDim[2];
3649     int piCountDim[2];
3650
3651     //evaluate each argument and replace by appropriate value and compute the count of combinations
3652     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3653     if (iSeqCount == 0)
3654     {
3655         //free pArg content
3656         cleanIndexesArguments(_pArgs, &pArg);
3657         return this;
3658     }
3659
3660     bool* pbFull = new bool[iDims];
3661     //coord must represent all values on a dimension
3662     for (int i = 0; i < iDims; i++)
3663     {
3664         pbFull[i] = false;
3665         int iDimToCheck = getVarMaxDim(i, iDims);
3666         int iIndexSize = pArg[i]->getAs<GenericType>()->getSize();
3667
3668         //we can have index more than once
3669         if (iIndexSize >= iDimToCheck)
3670         {
3671             //size is good, now check datas
3672             double* pIndexes = getDoubleArrayFromDouble(pArg[i]);
3673             for (int j = 0; j < iDimToCheck; j++)
3674             {
3675                 bool bFind = false;
3676                 for (int k = 0; k < iIndexSize; k++)
3677                 {
3678                     if ((int)pIndexes[k] == j + 1)
3679                     {
3680                         bFind = true;
3681                         break;
3682                     }
3683                 }
3684                 pbFull[i] = bFind;
3685             }
3686         }
3687     }
3688
3689     //only one dims can be not full/entire
3690     bool bNotEntire = false;
3691     int iNotEntire = 0;
3692     bool bTooMuchNotEntire = false;
3693     for (int i = 0; i < iDims; i++)
3694     {
3695         if (pbFull[i] == false)
3696         {
3697             if (bNotEntire == false)
3698             {
3699                 bNotEntire = true;
3700                 iNotEntire = i;
3701             }
3702             else
3703             {
3704                 bTooMuchNotEntire = true;
3705                 break;
3706             }
3707         }
3708     }
3709
3710     delete[] pbFull;
3711
3712     if (bTooMuchNotEntire == true)
3713     {
3714         //free pArg content
3715         cleanIndexesArguments(_pArgs, &pArg);
3716         return NULL;
3717     }
3718
3719     //find index to keep
3720     int iNotEntireSize = pArg[iNotEntire]->getAs<GenericType>()->getSize();
3721     double* piNotEntireIndex = getDoubleArrayFromDouble(pArg[iNotEntire]);
3722     int iKeepSize = getVarMaxDim(iNotEntire, iDims);
3723     bool* pbKeep = new bool[iKeepSize];
3724
3725     //fill pbKeep with true value
3726     for (int i = 0; i < iKeepSize; i++)
3727     {
3728         pbKeep[i] = true;
3729     }
3730
3731     for (int i = 0; i < iNotEntireSize; i++)
3732     {
3733         int idx = (int)piNotEntireIndex[i] - 1;
3734
3735         //don't care of value out of bounds
3736         if (idx < iKeepSize)
3737         {
3738             pbKeep[idx] = false;
3739         }
3740     }
3741
3742     int iNewDimSize = 0;
3743     for (int i = 0; i < iKeepSize; i++)
3744     {
3745         if (pbKeep[i] == true)
3746         {
3747             iNewDimSize++;
3748         }
3749     }
3750     delete[] pbKeep;
3751
3752     int* piNewDims = new int[iDims];
3753     for (int i = 0; i < iDims; i++)
3754     {
3755         if (i == iNotEntire)
3756         {
3757             piNewDims[i] = iNewDimSize;
3758         }
3759         else
3760         {
3761             piNewDims[i] = getVarMaxDim(i, iDims);
3762         }
3763     }
3764
3765     //remove last dimension if are == 1
3766     int iOrigDims = iDims;
3767     for (int i = (iDims - 1); i >= 2; i--)
3768     {
3769         if (piNewDims[i] == 1)
3770         {
3771             iDims--;
3772         }
3773         else
3774         {
3775             break;
3776         }
3777     }
3778
3779     if (iDims == 1)
3780     {
3781         if (iNewDimSize == 0)
3782         {
3783             //free pArg content
3784             cleanIndexesArguments(_pArgs, &pArg);
3785             return new SparseBool(0, 0);
3786         }
3787         else
3788         {
3789             //two cases, depends of original matrix/vector
3790             if ((*_pArgs)[0]->isColon() == false && m_iDims == 2 && m_piDims[0] == 1 && m_piDims[1] != 1)
3791             {
3792                 //special case for row vector
3793                 pOut = new SparseBool(1, iNewDimSize);
3794                 //in this case we have to care of 2nd dimension
3795                 //iNotEntire = 1;
3796             }
3797             else
3798             {
3799                 pOut = new SparseBool(iNewDimSize, 1);
3800             }
3801         }
3802     }
3803     else
3804     {
3805         pOut = new SparseBool(piNewDims[0], piNewDims[0]);
3806     }
3807
3808     delete[] piNewDims;
3809     //find a way to copy existing data to new variable ...
3810     int iNewPos = 0;
3811     int* piIndexes = new int[iOrigDims];
3812     int* piViewDims = new int[iOrigDims];
3813     for (int i = 0; i < iOrigDims; i++)
3814     {
3815         piViewDims[i] = getVarMaxDim(i, iOrigDims);
3816     }
3817
3818     for (int i = 0; i < getSize(); i++)
3819     {
3820         bool bByPass = false;
3821         getIndexesWithDims(i, piIndexes, piViewDims, iOrigDims);
3822
3823         //check if piIndexes use removed indexes
3824         for (int j = 0; j < iNotEntireSize; j++)
3825         {
3826             if ((piNotEntireIndex[j] - 1) == piIndexes[iNotEntire])
3827             {
3828                 //by pass this value
3829                 bByPass = true;
3830                 break;
3831             }
3832         }
3833
3834         if (bByPass == false)
3835         {
3836             //compute new index
3837             pOut->set(iNewPos, get(i));
3838             iNewPos++;
3839         }
3840     }
3841
3842     delete[] piIndexes;
3843     delete[] piViewDims;
3844
3845     //free pArg content
3846     cleanIndexesArguments(_pArgs, &pArg);
3847
3848     return pOut;
3849 }
3850
3851 SparseBool* SparseBool::append(int r, int c, SparseBool SPARSE_CONST* src)
3852 {
3853     SparseBool* pIT = checkRef(this, &SparseBool::append, r, c, src);
3854     if (pIT != this)
3855     {
3856         return pIT;
3857     }
3858
3859     doAppend(*src, r, c, *matrixBool);
3860     finalize();
3861     return this;
3862 }
3863
3864 GenericType* SparseBool::insertNew(typed_list* _pArgs)
3865 {
3866     typed_list pArg;
3867     SparseBool *pOut  = NULL;
3868
3869     int iDims = (int)_pArgs->size();
3870     int* piMaxDim = new int[iDims];
3871     int* piCountDim = new int[iDims];
3872     bool bUndefine = false;
3873
3874     //evaluate each argument and replace by appropriate value and compute the count of combinations
3875     int iSeqCount = checkIndexesArguments(NULL, _pArgs, &pArg, piMaxDim, piCountDim);
3876     if (iSeqCount == 0)
3877     {
3878         //free pArg content
3879         cleanIndexesArguments(_pArgs, &pArg);
3880         return createEmptyDouble();
3881     }
3882
3883     if (iSeqCount < 0)
3884     {
3885         iSeqCount = -iSeqCount;
3886         bUndefine = true;
3887     }
3888
3889     if (bUndefine)
3890     {
3891         //manage : and $ in creation by insertion
3892         int iSource = 0;
3893         int *piSourceDims = getDimsArray();
3894
3895         for (int i = 0; i < iDims; i++)
3896         {
3897             if (pArg[i] == NULL)
3898             {
3899                 //undefine value
3900                 if (isScalar())
3901                 {
3902                     piMaxDim[i] = 1;
3903                     piCountDim[i] = 1;
3904                 }
3905                 else
3906                 {
3907                     piMaxDim[i] = piSourceDims[iSource];
3908                     piCountDim[i] = piSourceDims[iSource];
3909                 }
3910                 iSource++;
3911                 //replace pArg value by the new one
3912                 pArg[i] = createDoubleVector(piMaxDim[i]);
3913             }
3914             //else
3915             //{
3916             //    piMaxDim[i] = piCountDim[i];
3917             //}
3918         }
3919     }
3920
3921     //remove last dimension at size 1
3922     //remove last dimension if are == 1
3923     for (int i = (iDims - 1); i >= 2; i--)
3924     {
3925         if (piMaxDim[i] == 1)
3926         {
3927             iDims--;
3928             pArg.pop_back();
3929         }
3930         else
3931         {
3932             break;
3933         }
3934     }
3935
3936     if (checkArgValidity(pArg) == false)
3937     {
3938         //free pArg content
3939         cleanIndexesArguments(_pArgs, &pArg);
3940         //contain bad index, like <= 0, ...
3941         return NULL;
3942     }
3943
3944     if (iDims == 1)
3945     {
3946         if (getCols() == 1)
3947         {
3948             pOut = new SparseBool(piCountDim[0], 1);
3949         }
3950         else
3951         {
3952             //rows == 1
3953             pOut = new SparseBool(1, piCountDim[0]);
3954         }
3955     }
3956     else
3957     {
3958         pOut = new SparseBool(piMaxDim[0], piMaxDim[1]);
3959     }
3960
3961     //insert values in new matrix
3962     SparseBool* pOut2 = pOut->insert(&pArg, this);
3963     if (pOut != pOut2)
3964     {
3965         delete pOut;
3966     }
3967
3968     //free pArg content
3969     cleanIndexesArguments(_pArgs, &pArg);
3970
3971     return pOut2;
3972 }
3973
3974 SparseBool* SparseBool::extract(int nbCoords, int SPARSE_CONST* coords, int SPARSE_CONST* maxCoords, int SPARSE_CONST* resSize, bool asVector) SPARSE_CONST
3975 {
3976     if ((asVector && maxCoords[0] > getSize()) ||
3977             (asVector == false && maxCoords[0] > getRows()) ||
3978             (asVector == false && maxCoords[1] > getCols()))
3979     {
3980         return 0;
3981     }
3982
3983     SparseBool * pSp(0);
3984     if (asVector)
3985     {
3986         pSp = (getRows() == 1) ? new SparseBool(1, resSize[0]) : new SparseBool(resSize[0], 1);
3987         mycopy_n(makeMatrixIterator<bool>(*this, Coords<true>(coords, getRows())), nbCoords
3988                  , makeMatrixIterator<bool>(*(pSp->matrixBool), RowWiseFullIterator(pSp->getRows(), pSp->getCols())));
3989     }
3990     else
3991     {
3992         pSp = new SparseBool(resSize[0], resSize[1]);
3993         mycopy_n(makeMatrixIterator<bool>(*this, Coords<false>(coords, getRows())), nbCoords
3994                  , makeMatrixIterator<bool>(*(pSp->matrixBool), RowWiseFullIterator(pSp->getRows(), pSp->getCols())));
3995
3996     }
3997     return pSp;
3998 }
3999
4000 /*
4001 * create a new SparseBool of dims according to resSize and fill it from currentSparseBool (along coords)
4002 */
4003 GenericType* SparseBool::extract(typed_list* _pArgs)
4004 {
4005     SparseBool* pOut = NULL;
4006     int iDims = (int)_pArgs->size();
4007     typed_list pArg;
4008
4009     int* piMaxDim = new int[iDims];
4010     int* piCountDim = new int[iDims];
4011
4012     //evaluate each argument and replace by appropriate value and compute the count of combinations
4013     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
4014     if (iSeqCount == 0)
4015     {
4016         //free pArg content
4017         cleanIndexesArguments(_pArgs, &pArg);
4018         if (_pArgs->size() == 0)
4019         {
4020             delete[] piMaxDim;
4021             delete[] piCountDim;
4022             //a()
4023             return this;
4024         }
4025         else
4026         {
4027             delete[] piMaxDim;
4028             delete[] piCountDim;
4029             //a([])
4030             return Double::Empty();
4031         }
4032     }
4033
4034     if (iDims < 2)
4035     {
4036         // Check that we stay inside the input size.
4037         if (piMaxDim[0] <= getSize())
4038         {
4039             int iNewRows = 0;
4040             int iNewCols = 0;
4041
4042             if (getRows() == 1 && getCols() != 1 && (*_pArgs)[0]->isColon() == false)
4043             {
4044                 //special case for row vector
4045                 iNewRows = 1;
4046                 iNewCols = piCountDim[0];
4047             }
4048             else
4049             {
4050                 iNewRows = piCountDim[0];
4051                 iNewCols = 1;
4052             }
4053
4054             pOut = new SparseBool(iNewRows, iNewCols);
4055             double* pIdx = pArg[0]->getAs<Double>()->get();
4056             // Write in output all elements extract from input.
4057             for (int i = 0; i < iSeqCount; i++)
4058             {
4059                 if (pIdx[i] < 1)
4060                 {
4061                     delete pOut;
4062                     pOut = NULL;
4063                     break;
4064                 }
4065                 int iRowRead = static_cast<int>(pIdx[i] - 1) % getRows();
4066                 int iColRead = static_cast<int>(pIdx[i] - 1) / getRows();
4067
4068                 int iRowWrite = static_cast<int>(i) % iNewRows;
4069                 int iColWrite = static_cast<int>(i) / iNewRows;
4070
4071                 bool bValue = get(iRowRead, iColRead);
4072                 if (bValue)
4073                 {
4074                     //only non zero values
4075                     pOut->set(iRowWrite, iColWrite, true, false);
4076                 }
4077             }
4078         }
4079         else
4080         {
4081             delete[] piMaxDim;
4082             delete[] piCountDim;
4083             //free pArg content
4084             cleanIndexesArguments(_pArgs, &pArg);
4085             return NULL;
4086         }
4087     }
4088     else
4089     {
4090         // Check that we stay inside the input size.
4091         if (piMaxDim[0] <= getRows() && piMaxDim[1] <= getCols())
4092         {
4093             double* pIdxRow = pArg[0]->getAs<Double>()->get();
4094             double* pIdxCol = pArg[1]->getAs<Double>()->get();
4095
4096             int iNewRows = pArg[0]->getAs<Double>()->getSize();
4097             int iNewCols = pArg[1]->getAs<Double>()->getSize();
4098
4099             pOut = new SparseBool(iNewRows, iNewCols);
4100
4101             int iPos = 0;
4102             // Write in output all elements extract from input.
4103             for (int iRow = 0; iRow < iNewRows; iRow++)
4104             {
4105                 for (int iCol = 0; iCol < iNewCols; iCol++)
4106                 {
4107                     if ((pIdxRow[iRow] < 1) || (pIdxCol[iCol] < 1))
4108                     {
4109                         delete pOut;
4110                         pOut = NULL;
4111                         delete[] piMaxDim;
4112                         delete[] piCountDim;
4113                         //free pArg content
4114                         cleanIndexesArguments(_pArgs, &pArg);
4115                         return NULL;
4116                     }
4117                     bool bValue = get((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
4118                     if (bValue)
4119                     {
4120                         //only non zero values
4121                         pOut->set(iRow, iCol, true, false);
4122                     }
4123                     iPos++;
4124                 }
4125             }
4126         }
4127         else
4128         {
4129             delete[] piMaxDim;
4130             delete[] piCountDim;
4131             //free pArg content
4132             cleanIndexesArguments(_pArgs, &pArg);
4133             return NULL;
4134         }
4135     }
4136
4137     if (pOut)
4138     {
4139         pOut->finalize();
4140     }
4141
4142     delete[] piMaxDim;
4143     delete[] piCountDim;
4144     //free pArg content
4145     cleanIndexesArguments(_pArgs, &pArg);
4146
4147     return pOut;
4148 }
4149
4150 bool SparseBool::invoke(typed_list & in, optional_list &/*opt*/, int /*_iRetCount*/, typed_list & out, const ast::Exp & e)
4151 {
4152     if (in.size() == 0)
4153     {
4154         out.push_back(this);
4155     }
4156     else
4157     {
4158         InternalType * _out = extract(&in);
4159         if (!_out)
4160         {
4161             std::wostringstream os;
4162             os << _W("Invalid index.\n");
4163             throw ast::InternalError(os.str(), 999, e.getLocation());
4164         }
4165         out.push_back(_out);
4166     }
4167
4168     return true;
4169 }
4170
4171 bool SparseBool::isInvokable() const
4172 {
4173     return true;
4174 }
4175
4176 bool SparseBool::hasInvokeOption() const
4177 {
4178     return false;
4179 }
4180
4181 int SparseBool::getInvokeNbIn()
4182 {
4183     return -1;
4184 }
4185
4186 int SparseBool::getInvokeNbOut()
4187 {
4188     return 1;
4189 }
4190
4191 std::size_t SparseBool::nbTrue() const
4192 {
4193     return  matrixBool->nonZeros();
4194 }
4195 std::size_t SparseBool::nbTrue(std::size_t r) const
4196 {
4197     int* piIndex = matrixBool->outerIndexPtr();
4198     return piIndex[r + 1] - piIndex[r];
4199 }
4200
4201
4202 void SparseBool::setTrue(bool finalize)
4203 {
4204     int rows = getRows();
4205     int cols = getCols();
4206
4207     std::vector<BoolTriplet_t> tripletList;
4208
4209     for (int i = 0; i < rows; ++i)
4210     {
4211         for (int j = 0; j < cols; ++j)
4212         {
4213             tripletList.emplace_back(i, j, true);
4214         }
4215     }
4216
4217     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<bool>());
4218
4219     if (finalize)
4220     {
4221         matrixBool->finalize();
4222     }
4223 }
4224
4225 void SparseBool::setFalse(bool finalize)
4226 {
4227     int rows = getRows();
4228     int cols = getCols();
4229
4230     std::vector<BoolTriplet_t> tripletList;
4231
4232     for (int i = 0; i < rows; ++i)
4233     {
4234         for (int j = 0; j < cols; ++j)
4235         {
4236             tripletList.emplace_back(i, j, false);
4237         }
4238     }
4239
4240     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<bool>());
4241
4242     if (finalize)
4243     {
4244         matrixBool->finalize();
4245     }
4246 }
4247
4248 int* SparseBool::getNbItemByRow(int* _piNbItemByRows)
4249 {
4250     int* piNbItemByRows = new int[getRows() + 1];
4251     mycopy_n(matrixBool->outerIndexPtr(), getRows() + 1, piNbItemByRows);
4252
4253     for (int i = 0; i < getRows(); i++)
4254     {
4255         _piNbItemByRows[i] = piNbItemByRows[i + 1] - piNbItemByRows[i];
4256     }
4257
4258     delete[] piNbItemByRows;
4259     return _piNbItemByRows;
4260 }
4261
4262 int* SparseBool::getColPos(int* _piColPos)
4263 {
4264     mycopy_n(matrixBool->innerIndexPtr(), nbTrue(), _piColPos);
4265     for (size_t i = 0; i < nbTrue(); i++)
4266     {
4267         _piColPos[i]++;
4268     }
4269
4270     return _piColPos;
4271 }
4272
4273 int* SparseBool::outputRowCol(int* out)const
4274 {
4275     return sparseTransform(*matrixBool, sparseTransform(*matrixBool, out, GetRow<BoolSparse_t>()), GetCol<BoolSparse_t>());
4276 }
4277
4278 int* SparseBool::getInnerPtr(int* count)
4279 {
4280     *count = static_cast<int>(matrixBool->innerSize());
4281     return matrixBool->innerIndexPtr();
4282 }
4283
4284 int* SparseBool::getOuterPtr(int* count)
4285 {
4286     *count = static_cast<int>(matrixBool->outerSize());
4287     return matrixBool->outerIndexPtr();
4288 }
4289
4290
4291 bool SparseBool::operator==(const InternalType& it) SPARSE_CONST
4292 {
4293     SparseBool* otherSparse = const_cast<SparseBool*>(dynamic_cast<SparseBool const*>(&it));/* types::GenericType is not const-correct :( */
4294     return (otherSparse
4295             && (otherSparse->getRows() == getRows())
4296             && (otherSparse->getCols() == getCols())
4297             && equal(*matrixBool, *otherSparse->matrixBool));
4298 }
4299
4300 bool SparseBool::operator!=(const InternalType& it) SPARSE_CONST
4301 {
4302     return !(*this == it);
4303 }
4304
4305 void SparseBool::finalize()
4306 {
4307     matrixBool->prune(&keepForSparse<bool>);
4308     matrixBool->finalize();
4309 }
4310
4311 bool SparseBool::get(int r, int c) SPARSE_CONST
4312 {
4313     return matrixBool->coeff(r, c);
4314 }
4315
4316 SparseBool* SparseBool::set(int _iRows, int _iCols, bool _bVal, bool _bFinalize) SPARSE_CONST
4317 {
4318     typedef SparseBool* (SparseBool::*set_t)(int, int, bool, bool);
4319     SparseBool* pIT = checkRef(this, (set_t)&SparseBool::set, _iRows, _iCols, _bVal, _bFinalize);
4320     if (pIT != this)
4321     {
4322         return pIT;
4323     }
4324
4325     if (matrixBool->isCompressed() && matrixBool->coeff(_iRows, _iCols) == false)
4326     {
4327         matrixBool->reserve(1);
4328     }
4329
4330     matrixBool->coeffRef(_iRows, _iCols) = _bVal;
4331
4332     if (_bFinalize)
4333     {
4334         finalize();
4335     }
4336
4337     return this;
4338 }
4339
4340 void SparseBool::fill(Bool& dest, int r, int c) SPARSE_CONST
4341 {
4342     mycopy_n(makeMatrixIterator<bool >(*matrixBool, RowWiseFullIterator(getRows(), getCols())), getSize()
4343              , makeMatrixIterator<bool >(dest, RowWiseFullIterator(dest.getRows(), dest.getCols(), r, c)));
4344 }
4345
4346 Sparse* SparseBool::newOnes() const
4347 {
4348     return new Sparse(new types::Sparse::RealSparse_t(matrixBool->cast<double>()), 0);
4349 }
4350
4351 SparseBool* SparseBool::newNotEqualTo(SparseBool const&o) const
4352 {
4353     return cwiseOp<std::not_equal_to>(*this, o);
4354 }
4355
4356 SparseBool* SparseBool::newEqualTo(SparseBool& o)
4357 {
4358     int rowL = getRows();
4359     int colL = getCols();
4360
4361     int rowR = o.getRows();
4362     int colR = o.getCols();
4363     int row = std::max(rowL, rowR);
4364     int col = std::max(colL, colR);
4365
4366     //create a boolean sparse matrix with dims of sparses
4367     types::SparseBool* ret = new types::SparseBool(row, col);
4368
4369     if (isScalar() && o.isScalar())
4370     {
4371         bool l = get(0, 0);
4372         bool r = o.get(0, 0);
4373         ret->set(0, 0, l == r, false);
4374     }
4375     else if (isScalar())
4376     {
4377         int nnzR = static_cast<int>(o.nbTrue());
4378         std::vector<int> rowcolR(nnzR * 2, 0);
4379         o.outputRowCol(rowcolR.data());
4380
4381         //compare all items of R with R[0]
4382         bool l = get(0, 0);
4383         for (int i = 0; i < nnzR; ++i)
4384         {
4385             bool r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
4386             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
4387         }
4388     }
4389     else if (o.isScalar())
4390     {
4391         int nnzL = static_cast<int>(nbTrue());
4392         std::vector<int> rowcolL(nnzL * 2, 0);
4393         outputRowCol(rowcolL.data());
4394
4395         bool r = get(0, 0);
4396         for (int i = 0; i < nnzL; ++i)
4397         {
4398             bool l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
4399             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
4400         }
4401     }
4402     else
4403     {
4404         int nnzR = static_cast<int>(o.nbTrue());
4405         std::vector<int> rowcolR(nnzR * 2, 0);
4406         o.outputRowCol(rowcolR.data());
4407         int nnzL = static_cast<int>(nbTrue());
4408         std::vector<int> rowcolL(nnzL * 2, 0);
4409         outputRowCol(rowcolL.data());
4410         //set all values to %t
4411         ret->setTrue(false);
4412         //set %f in each pL values
4413         for (int i = 0; i < nnzL; ++i)
4414         {
4415             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, false, false);
4416         }
4417         ret->finalize();
4418
4419         //set _pR[i] == _pL[i] for each _pR values
4420         for (int i = 0; i < nnzR; ++i)
4421         {
4422             //get l and r following non zeros value of R
4423             bool l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
4424             bool r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
4425             //set value following non zeros value of R
4426             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
4427         }
4428     }
4429
4430     ret->finalize();
4431     return ret;
4432 }
4433
4434 SparseBool* SparseBool::newLogicalOr(SparseBool const&o) const
4435 {
4436     return cwiseOp<std::logical_or>(*this, o);
4437 }
4438
4439 SparseBool* SparseBool::newLogicalAnd(SparseBool const&o) const
4440 {
4441     return cwiseOp<std::logical_and>(*this, o);
4442 }
4443
4444 SparseBool* SparseBool::reshape(int* _piDims, int _iDims)
4445 {
4446     SparseBool* pSpBool = NULL;
4447     int iCols = 1;
4448
4449     if (_iDims == 2)
4450     {
4451         iCols = _piDims[1];
4452     }
4453
4454     if (_iDims <= 2)
4455     {
4456         pSpBool = reshape(_piDims[0], iCols);
4457     }
4458
4459     return pSpBool;
4460 }
4461
4462 SparseBool* SparseBool::reshape(int _iNewRows, int _iNewCols)
4463 {
4464     typedef SparseBool* (SparseBool::*reshape_t)(int, int);
4465     SparseBool* pIT = checkRef(this, (reshape_t)&SparseBool::reshape, _iNewRows, _iNewCols);
4466     if (pIT != this)
4467     {
4468         return pIT;
4469     }
4470
4471     if (_iNewRows * _iNewCols != getRows() * getCols())
4472     {
4473         return NULL;
4474     }
4475
4476     SparseBool* res = NULL;
4477     try
4478     {
4479         //item count
4480         size_t iNonZeros = matrixBool->nonZeros();
4481         BoolSparse_t *newBool = new BoolSparse_t(_iNewRows, _iNewCols);
4482         newBool->reserve((int)iNonZeros);
4483
4484         //coords
4485         int* pRows = new int[iNonZeros * 2];
4486         outputRowCol(pRows);
4487         int* pCols = pRows + iNonZeros;
4488
4489         std::vector<BoolTriplet_t> tripletList;
4490
4491         for (size_t i = 0; i < iNonZeros; i++)
4492         {
4493             int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
4494             tripletList.emplace_back((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), true);
4495         }
4496
4497         newBool->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<bool>());
4498
4499         delete matrixBool;
4500         matrixBool = newBool;
4501         delete[] pRows;
4502
4503         m_iRows = _iNewRows;
4504         m_iCols = _iNewCols;
4505         m_iSize = _iNewRows * _iNewCols;
4506
4507         m_iDims = 2;
4508         m_piDims[0] = m_iRows;
4509         m_piDims[1] = m_iCols;
4510
4511         finalize();
4512
4513         res = this;
4514     }
4515     catch (...)
4516     {
4517         res = NULL;
4518     }
4519     return res;
4520 }
4521
4522 bool SparseBool::transpose(InternalType *& out)
4523 {
4524     out = new SparseBool(new BoolSparse_t(matrixBool->transpose()));
4525     return true;
4526 }
4527
4528 template<typename T>
4529 void neg(const int r, const int c, const T * const in, Eigen::SparseMatrix<bool, 1> * const out)
4530 {
4531     for (int i = 0; i < r; i++)
4532     {
4533         for (int j = 0; j < c; j++)
4534         {
4535             out->coeffRef(i, j) = !in->coeff(i, j);
4536         }
4537     }
4538
4539     out->prune(&keepForSparse<bool>);
4540     out->finalize();
4541 }
4542 }