scilab/modules/ast/src/cpp/types/sparse.cpp

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