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