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