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

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