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