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