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     Sparse* res = NULL;
 876     try
 877     {
 878         if (matrixReal)
 879         {
 880             //item count
 881             size_t iNonZeros = nonZeros();
 882             RealSparse_t *newReal = new RealSparse_t(_iNewRows, _iNewCols);
 883             newReal->reserve((int)iNonZeros);
 884
 885
 886             //coords
 887             int* pRows = new int[iNonZeros * 2];
 888             outputRowCol(pRows);
 889             int* pCols = pRows + iNonZeros;
 890
 891             //values
 892             double* pNonZeroR = new double[iNonZeros];
 893             double* pNonZeroI = new double[iNonZeros];
 894             outputValues(pNonZeroR, pNonZeroI);
 895
 896             std::vector<RealTriplet_t> tripletList;
 897             for (size_t i = 0; i < iNonZeros; i++)
 898             {
 899                 tripletList.emplace_back((int)pRows[i] - 1, (int)pCols[i] - 1, pNonZeroR[i]);
 900             }
 901
 902             newReal->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<double>());
 903
 904             delete matrixReal;
 905             matrixReal = newReal;
 906             delete[] pRows;
 907             delete[] pNonZeroR;
 908             delete[] pNonZeroI;
 909         }
 910         else
 911         {
 912             //item count
 913             size_t iNonZeros = nonZeros();
 914             CplxSparse_t *newCplx = new CplxSparse_t(_iNewRows, _iNewCols);
 915             newCplx->reserve((int)iNonZeros);
 916
 917             //coords
 918             int* pRows = new int[iNonZeros * 2];
 919             outputRowCol(pRows);
 920             int* pCols = pRows + iNonZeros;
 921
 922             //values
 923             double* pNonZeroR = new double[iNonZeros];
 924             double* pNonZeroI = new double[iNonZeros];
 925             outputValues(pNonZeroR, pNonZeroI);
 926
 927             std::vector<CplxTriplet_t> tripletList;
 928             for (size_t i = 0; i < iNonZeros; i++)
 929             {
 930                 tripletList.emplace_back((int)pRows[i] - 1, (int)pCols[i] - 1, std::complex<double>(pNonZeroR[i], pNonZeroI[i]));
 931             }
 932
 933             newCplx->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<std::complex<double>>());
 934
 935
 936             delete matrixCplx;
 937             matrixCplx = newCplx;
 938             delete[] pRows;
 939             delete[] pNonZeroR;
 940             delete[] pNonZeroI;
 941         }
 942
 943         m_iRows = _iNewRows;
 944         m_iCols = _iNewCols;
 945         m_iSize = _iNewRows * _iNewCols;
 946         m_piDims[0] = m_iRows;
 947         m_piDims[1] = m_iCols;
 948
 949         res = this;
 950     }
 951     catch (...)
 952     {
 953         res = NULL;
 954     }
 955     return res;
 956 }
 957 // TODO decide if a complex matrix with 0 imag can be == to a real matrix
 958 // not true for dense (cf double.cpp)
 959 bool Sparse::operator==(const InternalType& it) SPARSE_CONST
 960 {
 961     Sparse* otherSparse = const_cast<Sparse*>(dynamic_cast<Sparse const*>(&it));/* types::GenericType is not const-correct :( */
 962     Sparse & cthis(const_cast<Sparse&>(*this));
 963
 964     if (otherSparse == NULL)
 965     {
 966         return false;
 967     }
 968
 969     if (otherSparse->getRows() != cthis.getRows())
 970     {
 971         return false;
 972     }
 973
 974     if (otherSparse->getCols() != cthis.getCols())
 975     {
 976         return false;
 977     }
 978
 979     if (otherSparse->isComplex() != isComplex())
 980     {
 981         return false;
 982     }
 983
 984     if (isComplex())
 985     {
 986         return equal(*matrixCplx, *otherSparse->matrixCplx);
 987     }
 988     else
 989     {
 990         return equal(*matrixReal, *otherSparse->matrixReal);
 991     }
 992 }
 993
 994 bool Sparse::one_set()
 995 {
 996     if (isComplex())
 997     {
 998         return setNonZero(*matrixCplx);
 999     }
1000     else
1001     {
1002         return setNonZero(*matrixReal);
1003     }
1004 }
1005
1006 void Sparse::toComplex()
1007 {
1008     if (!isComplex())
1009     {
1010         try
1011         {
1012             matrixCplx = new CplxSparse_t(matrixReal->cast<std::complex<double> >());
1013             delete matrixReal;
1014             matrixReal = NULL;
1015         }
1016         catch (...)
1017         {
1018             delete matrixCplx;
1019             matrixCplx = NULL;
1020             throw;
1021         }
1022     }
1023 }
1024
1025 GenericType* Sparse::insertNew(typed_list* _pArgs)
1026 {
1027     typed_list pArg;
1028     Sparse *pOut        = NULL;
1029
1030     int iDims           = (int)_pArgs->size();
1031     int* piMaxDim       = new int[iDims];
1032     int* piCountDim     = new int[iDims];
1033     bool bComplex       = isComplex();
1034     bool bUndefine      = false;
1035
1036     //evaluate each argument and replace by appropriate value and compute the count of combinations
1037     int iSeqCount = checkIndexesArguments(NULL, _pArgs, &pArg, piMaxDim, piCountDim);
1038     if (iSeqCount == 0)
1039     {
1040         //free pArg content
1041         cleanIndexesArguments(_pArgs, &pArg);
1042         return createEmptyDouble();
1043     }
1044
1045     if (iSeqCount < 0)
1046     {
1047         iSeqCount = -iSeqCount;
1048         bUndefine = true;
1049     }
1050
1051     if (bUndefine)
1052     {
1053         //manage : and $ in creation by insertion
1054         int iSource = 0;
1055         int *piSourceDims = getDimsArray();
1056
1057         for (int i = 0; i < iDims; i++)
1058         {
1059             if (pArg[i] == NULL)
1060             {
1061                 //undefine value
1062                 if (isScalar())
1063                 {
1064                     piMaxDim[i] = 1;
1065                     piCountDim[i] = 1;
1066                 }
1067                 else
1068                 {
1069                     piMaxDim[i] = piSourceDims[iSource];
1070                     piCountDim[i] = piSourceDims[iSource];
1071                 }
1072                 iSource++;
1073                 //replace pArg value by the new one
1074                 pArg[i] = createDoubleVector(piMaxDim[i]);
1075             }
1076             //else
1077             //{
1078             //    piMaxDim[i] = piCountDim[i];
1079             //}
1080         }
1081     }
1082
1083     //remove last dimension at size 1
1084     //remove last dimension if are == 1
1085     for (int i = (iDims - 1); i >= 2; i--)
1086     {
1087         if (piMaxDim[i] == 1)
1088         {
1089             iDims--;
1090             pArg.pop_back();
1091         }
1092         else
1093         {
1094             break;
1095         }
1096     }
1097
1098     if (checkArgValidity(pArg) == false)
1099     {
1100         //free pArg content
1101         cleanIndexesArguments(_pArgs, &pArg);
1102         //contain bad index, like <= 0, ...
1103         return NULL;
1104     }
1105
1106     if (iDims == 1)
1107     {
1108         if (getCols() == 1)
1109         {
1110             pOut = new Sparse(piCountDim[0], 1, bComplex);
1111         }
1112         else
1113         {
1114             //rows == 1
1115             pOut = new Sparse(1, piCountDim[0], bComplex);
1116         }
1117     }
1118     else
1119     {
1120         pOut = new Sparse(piMaxDim[0], piMaxDim[1], bComplex);
1121         //pOut = createEmpty(iDims, piMaxDim, bComplex);
1122     }
1123
1124     //insert values in new matrix
1125     Sparse* pOut2 = pOut->insert(&pArg, this);
1126     if (pOut != pOut2)
1127     {
1128         delete pOut;
1129     }
1130
1131     //free pArg content
1132     cleanIndexesArguments(_pArgs, &pArg);
1133
1134     return pOut2;
1135 }
1136
1137 Sparse* Sparse::insert(typed_list* _pArgs, InternalType* _pSource)
1138 {
1139     typedef Sparse* (Sparse::*insert_t)(typed_list*, InternalType*);
1140     Sparse* pIT = checkRef(this, (insert_t)&Sparse::insert, _pArgs, _pSource);
1141     if (pIT != this)
1142     {
1143         return pIT;
1144     }
1145
1146     if (_pSource->isSparse())
1147     {
1148         return insert(_pArgs, _pSource->getAs<Sparse>());
1149     }
1150
1151     bool bNeedToResize  = false;
1152     int iDims           = (int)_pArgs->size();
1153     if (iDims > 2)
1154     {
1155         //sparse are only in 2 dims
1156         return NULL;
1157     }
1158
1159     typed_list pArg;
1160
1161     int piMaxDim[2];
1162     int piCountDim[2];
1163
1164     //on case of resize
1165     int iNewRows = 0;
1166     int iNewCols = 0;
1167     Double* pSource = _pSource->getAs<Double>();
1168
1169     //evaluate each argument and replace by appropriate value and compute the count of combinations
1170     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1171     if (iSeqCount == 0)
1172     {
1173         //free pArg content
1174         cleanIndexesArguments(_pArgs, &pArg);
1175         return this;
1176     }
1177
1178     if (iDims < 2)
1179     {
1180         //see as vector
1181         if (getRows() == 1 || getCols() == 1)
1182         {
1183             //vector or scalar
1184             if (getRows() * getCols() < piMaxDim[0])
1185             {
1186                 bNeedToResize = true;
1187
1188                 //need to enlarge sparse dimensions
1189                 if (getCols() == 1 || getRows() * getCols() == 0)
1190                 {
1191                     //column vector
1192                     iNewRows = piMaxDim[0];
1193                     iNewCols = 1;
1194                 }
1195                 else if (getRows() == 1)
1196                 {
1197                     //row vector
1198                     iNewRows = 1;
1199                     iNewCols = piMaxDim[0];
1200                 }
1201             }
1202         }
1203         else if ((size_t)getRows() * (size_t)getCols() < (size_t)piMaxDim[0])
1204         {
1205             //free pArg content
1206             cleanIndexesArguments(_pArgs, &pArg);
1207             //out of range
1208             return NULL;
1209         }
1210     }
1211     else
1212     {
1213         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
1214         {
1215             bNeedToResize = true;
1216             iNewRows = std::max(getRows(), piMaxDim[0]);
1217             iNewCols = std::max(getCols(), piMaxDim[1]);
1218         }
1219     }
1220
1221     //check number of insertion
1222     if (pSource->isScalar() == false && pSource->getSize() != iSeqCount)
1223     {
1224         //free pArg content
1225         cleanIndexesArguments(_pArgs, &pArg);
1226         return NULL;
1227     }
1228
1229     //now you are sure to be able to insert values
1230     if (bNeedToResize)
1231     {
1232         if (resize(iNewRows, iNewCols) == NULL)
1233         {
1234             //free pArg content
1235             cleanIndexesArguments(_pArgs, &pArg);
1236             return NULL;
1237         }
1238     }
1239
1240     //update complexity
1241     if (pSource->isComplex() && isComplex() == false)
1242     {
1243         toComplex();
1244     }
1245
1246     int rows = getRows();
1247     int cols = getCols();
1248
1249     int nnz = static_cast<int>(nonZeros());
1250
1251     double ratio = 1;
1252     int inserted = 0;
1253     int updated = 0;
1254
1255     if (nnz != 0)
1256     {
1257         if (iDims != 1)
1258         {
1259             if (isComplex())
1260             {
1261                 getinsertedupdated(matrixCplx, pArg[0]->getAs<Double>(), pArg[1]->getAs<Double>(), updated, inserted);
1262             }
1263             else
1264             {
1265                 getinsertedupdated(matrixReal, pArg[0]->getAs<Double>(), pArg[1]->getAs<Double>(), updated, inserted);
1266             }
1267         }
1268         else
1269         {
1270             if (isComplex())
1271             {
1272                 getinsertedupdated(matrixCplx, pArg[0]->getAs<Double>(), updated, inserted);
1273             }
1274             else
1275             {
1276                 getinsertedupdated(matrixReal, pArg[0]->getAs<Double>(), updated, inserted);
1277             }
1278         }
1279
1280         ratio = (double)inserted / (double)nnz;
1281     }
1282
1283     if (ratio < 0.05) // less 5%
1284     {
1285         int nnzFinal = nnz + inserted;
1286         if (isComplex())
1287         {
1288             matrixCplx->reserve(nnzFinal);
1289         }
1290         else
1291         {
1292             matrixReal->reserve(nnzFinal);
1293         }
1294
1295         if (iDims == 1)
1296         {
1297             double* pIdx = pArg[0]->getAs<Double>()->get();
1298             int rows = getRows();
1299             double* pR = pSource->get();
1300             double* pI = pSource->getImg();
1301
1302             for (int i = 0; i < iSeqCount; i++)
1303             {
1304                 int iRow = static_cast<int>(pIdx[i] - 1) % rows;
1305                 int iCol = static_cast<int>(pIdx[i] - 1) / rows;
1306                 if (pSource->isScalar())
1307                 {
1308                     if (pSource->isComplex())
1309                     {
1310                         set(iRow, iCol, std::complex<double>(pR[0], pI[0]), false);
1311                     }
1312                     else
1313                     {
1314                         set(iRow, iCol, pR[0], false);
1315                     }
1316                 }
1317                 else
1318                 {
1319                     if (pSource->isComplex())
1320                     {
1321                         set(iRow, iCol, std::complex<double>(pR[i], pI[i]), false);
1322                     }
1323                     else
1324                     {
1325                         set(iRow, iCol, pR[i], false);
1326                     }
1327                 }
1328             }
1329         }
1330         else
1331         {
1332             double* pIdxRow = pArg[0]->getAs<Double>()->get();
1333             int iRowSize = pArg[0]->getAs<Double>()->getSize();
1334             double* pIdxCol = pArg[1]->getAs<Double>()->get();
1335             double* pR = pSource->get();
1336             double* pI = pSource->getImg();
1337             if (pSource->isScalar())
1338             {
1339                 if (isComplex())
1340                 {
1341                     //scalar complex
1342                     std::complex<double> val(pR[0], pI[0]);
1343                     for (int i = 0; i < iSeqCount; i++)
1344                     {
1345                         set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, val, false);
1346                     }
1347                 }
1348                 else
1349                 {
1350                     //scalar real
1351                     double val = pR[0];
1352                     for (int i = 0; i < iSeqCount; i++)
1353                     {
1354                         set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, val, false);
1355                     }
1356                 }
1357             }
1358             else
1359             {
1360                 if (isComplex())
1361                 {
1362                     //matrix complex
1363                     for (int i = 0; i < iSeqCount; i++)
1364                     {
1365                         set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, std::complex<double>(pR[i], pI[i]), false);
1366                     }
1367                 }
1368                 else
1369                 {
1370                     //matrix real
1371                     for (int i = 0; i < iSeqCount; i++)
1372                     {
1373                         set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pR[i], false);
1374                     }
1375                 }
1376             }
1377         }
1378     }
1379     else
1380     {
1381         if (iDims == 1)
1382         {
1383             if (isComplex())
1384             {
1385                 std::vector<CplxTriplet_t> tripletList;
1386
1387                 double* pIdx = pArg[0]->getAs<Double>()->get();
1388                 double* srcR = pSource->get();
1389                 double* srcI = NULL;
1390                 double zero = 0;
1391                 int incR = pSource->isScalar() ? 0 : 1;
1392
1393                 int incI = 0;
1394                 if (pSource->isComplex())
1395                 {
1396                     srcI = pSource->getImg();
1397                     incI = pSource->isScalar() ? 0 : 1;
1398                 }
1399                 else
1400                 {
1401                     srcI = &zero;
1402                     incI = 0;
1403                 }
1404
1405                 //save old values
1406                 if (nnz != 0)
1407                 {
1408                     std::complex<double>* val = matrixCplx->valuePtr();
1409
1410                     //save old values
1411                     for (int k = 0; k < matrixCplx->outerSize(); ++k)
1412                     {
1413                         for (CplxSparse_t::InnerIterator it(*matrixCplx, k); it; ++it)
1414                         {
1415                             //m[static_cast<size_t>(it.row()) + static_cast<size_t>(it.col()) * rows] = it.value();
1416                             tripletList.emplace_back(it.row(), it.col(), it.value());
1417                         }
1418                     }
1419
1420                     matrixCplx->setZero();
1421                 }
1422
1423                 for (int i = 0; i < iSeqCount; i++)
1424                 {
1425                     size_t idx = static_cast<size_t>(pIdx[i] - 1);
1426                     int iRow = static_cast<int>(idx % rows);
1427                     int iCol = static_cast<int>(idx / rows);
1428                     //m[static_cast<size_t>(pIdx[i]) - 1] = std::complex<double>(*srcR, *srcI);
1429                     tripletList.emplace_back(iRow, iCol, std::complex<double>(*srcR, *srcI));
1430                     srcR += incR;
1431                     srcI += incI;
1432                 }
1433
1434                 matrixCplx->reserve(static_cast<int>(tripletList.size()));
1435                 matrixCplx->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<std::complex<double>>());
1436
1437             }
1438             else
1439             {
1440                 std::vector<RealTriplet_t> tripletList;
1441
1442                 double* pIdx = pArg[0]->getAs<Double>()->get();
1443                 double* src = pSource->get();
1444                 int inc = pSource->isScalar() ? 0 : 1;
1445
1446                 if (nnz != 0)
1447                 {
1448                     //save old values
1449                     for (int k = 0; k < matrixReal->outerSize(); ++k)
1450                     {
1451                         for (RealSparse_t::InnerIterator it(*matrixReal, k); it; ++it)
1452                         {
1453                             //m[static_cast<size_t>(it.row()) + static_cast<size_t>(it.col()) * rows] = it.value();
1454                             tripletList.emplace_back(it.row(), it.col(), it.value());
1455                         }
1456                     }
1457
1458                     matrixReal->setZero();
1459                 }
1460
1461                 for (int i = 0; i < iSeqCount; i++)
1462                 {
1463                     size_t idx = static_cast<size_t>(pIdx[i] - 1);
1464                     int iRow = static_cast<int>(idx % rows);
1465                     int iCol = static_cast<int>(idx / rows);
1466                     //m[static_cast<size_t>(pIdx[i]) - 1] = *src;
1467                     tripletList.emplace_back(iRow, iCol, *src);
1468                     src += inc;
1469                 }
1470
1471                 matrixReal->reserve(static_cast<int>(tripletList.size()));
1472                 matrixReal->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<double>());
1473             }
1474
1475         }
1476         else
1477         {
1478             int iRowSize = pArg[0]->getAs<Double>()->getSize();
1479             double* pI = pArg[0]->getAs<Double>()->get();
1480             double* pJ = pArg[1]->getAs<Double>()->get();
1481
1482             if (isComplex())
1483             {
1484                 std::vector<CplxTriplet_t> tripletList;
1485                 double* srcR = pSource->get();
1486                 double* srcI = NULL;
1487                 double zero = 0;
1488                 int incR = pSource->isScalar() ? 0 : 1;
1489
1490                 int incI = 0;
1491                 if (pSource->isComplex())
1492                 {
1493                     srcI = pSource->getImg();
1494                     incI = pSource->isScalar() ? 0 : 1;
1495                 }
1496                 else
1497                 {
1498                     srcI = &zero;
1499                     incI = 0;
1500                 }
1501
1502                 if (nnz != 0)
1503                 {
1504                     //save old values
1505                     for (int k = 0; k < matrixCplx->outerSize(); ++k)
1506                     {
1507                         for (CplxSparse_t::InnerIterator it(*matrixCplx, k); it; ++it)
1508                         {
1509                             //m[static_cast<size_t>(it.row()) + static_cast<size_t>(it.col()) * rows] = it.value();
1510                             tripletList.emplace_back(it.row(), it.col(), it.value());
1511                         }
1512                     }
1513
1514                     matrixCplx->setZero();
1515                 }
1516
1517                 //add new values
1518                 for (int i = 0; i < iSeqCount; i++)
1519                 {
1520                     int iRow = static_cast<int>(i % iRowSize);
1521                     int iCol = static_cast<int>(i / iRowSize);
1522                     tripletList.emplace_back(static_cast<int>(pI[iRow] - 1), static_cast<int>(pJ[iCol] - 1), std::complex<double>(*srcR, *srcI));
1523                     srcR += incR;
1524                     srcI += incI;
1525                 }
1526
1527                 matrixCplx->reserve(static_cast<int>(tripletList.size()));
1528                 matrixCplx->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<std::complex<double>>());
1529             }
1530             else
1531             {
1532                 std::vector<RealTriplet_t> tripletList;
1533                 double* src = pSource->get();
1534                 int inc = pSource->isScalar() ? 0 : 1;
1535
1536                 if (nnz != 0)
1537                 {
1538                     double* val = matrixReal->valuePtr();
1539
1540                     //save old values
1541                     for (int k = 0; k < matrixReal->outerSize(); ++k)
1542                     {
1543                         for (RealSparse_t::InnerIterator it(*matrixReal, k); it; ++it)
1544                         {
1545                             //m[static_cast<size_t>(it.row()) + static_cast<size_t>(it.col()) * rows] = it.value();
1546                             tripletList.emplace_back(it.row(), it.col(), it.value());
1547                         }
1548                     }
1549                 }
1550
1551                 //add new values
1552                 for (int i = 0; i < iSeqCount; ++i)
1553                 {
1554                     int iRow = static_cast<int>(i % iRowSize);
1555                     int iCol = static_cast<int>(i / iRowSize);
1556                     tripletList.emplace_back(static_cast<int>(pI[iRow]) - 1, static_cast<int>(pJ[iCol]) - 1, *src);
1557                     src += inc;
1558                 }
1559
1560                 matrixReal->setZero();
1561                 matrixReal->reserve(static_cast<int>(tripletList.size()));
1562                 matrixReal->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<double>());
1563             }
1564         }
1565     }
1566
1567     finalize();
1568
1569     //free pArg content
1570     cleanIndexesArguments(_pArgs, &pArg);
1571     return this;
1572 }
1573
1574 Sparse* Sparse::insert(typed_list* _pArgs, Sparse* _pSource)
1575 {
1576     bool bNeedToResize = false;
1577     int iDims = (int)_pArgs->size();
1578     if (iDims > 2)
1579     {
1580         //sparse are only in 2 dims
1581         return NULL;
1582     }
1583
1584     typed_list pArg;
1585
1586     int piMaxDim[2];
1587     int piCountDim[2];
1588
1589     //on case of resize
1590     int iNewRows = 0;
1591     int iNewCols = 0;
1592
1593     //evaluate each argument and replace by appropriate value and compute the count of combinations
1594     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1595     if (iSeqCount == 0)
1596     {
1597         //free pArg content
1598         cleanIndexesArguments(_pArgs, &pArg);
1599         return this;
1600     }
1601
1602     if (iDims < 2)
1603     {
1604         //see as vector
1605         if (getRows() == 1 || getCols() == 1)
1606         {
1607             //vector or scalar
1608             bNeedToResize = true;
1609             if (getSize() < piMaxDim[0])
1610             {
1611                 //need to enlarge sparse dimensions
1612                 if (getCols() == 1 || getSize() == 0)
1613                 {
1614                     //column vector
1615                     iNewRows = piMaxDim[0];
1616                     iNewCols = 1;
1617                 }
1618                 else if (getRows() == 1)
1619                 {
1620                     //row vector
1621                     iNewRows = 1;
1622                     iNewCols = piMaxDim[0];
1623                 }
1624             }
1625         }
1626         else if (getSize() < piMaxDim[0])
1627         {
1628             //free pArg content
1629             cleanIndexesArguments(_pArgs, &pArg);
1630             //out of range
1631             return NULL;
1632         }
1633     }
1634     else
1635     {
1636         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
1637         {
1638             bNeedToResize = true;
1639             iNewRows = std::max(getRows(), piMaxDim[0]);
1640             iNewCols = std::max(getCols(), piMaxDim[1]);
1641         }
1642     }
1643
1644     //check number of insertion
1645     if (_pSource->isScalar() == false && _pSource->getSize() != iSeqCount)
1646     {
1647         //free pArg content
1648         cleanIndexesArguments(_pArgs, &pArg);
1649         return NULL;
1650     }
1651
1652     //now you are sure to be able to insert values
1653     if (bNeedToResize)
1654     {
1655         if (resize(iNewRows, iNewCols) == NULL)
1656         {
1657             //free pArg content
1658             cleanIndexesArguments(_pArgs, &pArg);
1659             return NULL;
1660         }
1661     }
1662
1663     //update complexity
1664     if (_pSource->isComplex() && isComplex() == false)
1665     {
1666         toComplex();
1667     }
1668
1669     if (iDims == 1)
1670     {
1671         double* pIdx = pArg[0]->getAs<Double>()->get();
1672         for (int i = 0; i < iSeqCount; i++)
1673         {
1674             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
1675             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
1676
1677             if (_pSource->isScalar())
1678             {
1679                 if (_pSource->isComplex())
1680                 {
1681                     set(iRow, iCol, _pSource->getImg(0, 0), false);
1682                 }
1683                 else
1684                 {
1685                     set(iRow, iCol, _pSource->get(0, 0), false);
1686                 }
1687             }
1688             else
1689             {
1690                 int iRowOrig = i % _pSource->getRows();
1691                 int iColOrig = i / _pSource->getRows();
1692                 if (_pSource->isComplex())
1693                 {
1694                     set(iRow, iCol, _pSource->getImg(iRowOrig, iColOrig), false);
1695                 }
1696                 else
1697                 {
1698                     set(iRow, iCol, _pSource->get(iRowOrig, iColOrig), false);
1699                 }
1700             }
1701         }
1702     }
1703     else
1704     {
1705         double* pIdxRow = pArg[0]->getAs<Double>()->get();
1706         int iRowSize = pArg[0]->getAs<Double>()->getSize();
1707         double* pIdxCol = pArg[1]->getAs<Double>()->get();
1708
1709         for (int i = 0; i < iSeqCount; i++)
1710         {
1711             if (_pSource->isScalar())
1712             {
1713                 if (_pSource->isComplex())
1714                 {
1715                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->getImg(0, 0), false);
1716                 }
1717                 else
1718                 {
1719                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(0, 0), false);
1720                 }
1721             }
1722             else
1723             {
1724                 int iRowOrig = i % _pSource->getRows();
1725                 int iColOrig = i / _pSource->getRows();
1726                 if (_pSource->isComplex())
1727                 {
1728                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->getImg(iRowOrig, iColOrig), false);
1729                 }
1730                 else
1731                 {
1732                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(iRowOrig, iColOrig), false);
1733                 }
1734             }
1735         }
1736     }
1737
1738     finalize();
1739
1740     //free pArg content
1741     cleanIndexesArguments(_pArgs, &pArg);
1742
1743     return this;
1744 }
1745
1746 GenericType* Sparse::remove(typed_list* _pArgs)
1747 {
1748     Sparse* pOut = NULL;
1749     int iDims = (int)_pArgs->size();
1750     if (iDims > 2)
1751     {
1752         //sparse are only in 2 dims
1753         return NULL;
1754     }
1755
1756     typed_list pArg;
1757
1758     int piMaxDim[2];
1759     int piCountDim[2];
1760
1761     //evaluate each argument and replace by appropriate value and compute the count of combinations
1762     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1763     if (iSeqCount == 0)
1764     {
1765         //free pArg content
1766         cleanIndexesArguments(_pArgs, &pArg);
1767         return this;
1768     }
1769
1770     bool* pbFull = new bool[iDims];
1771     //coord must represent all values on a dimension
1772     for (int i = 0; i < iDims; i++)
1773     {
1774         pbFull[i] = false;
1775         int iDimToCheck = getVarMaxDim(i, iDims);
1776         int iIndexSize = pArg[i]->getAs<GenericType>()->getSize();
1777
1778         //we can have index more than once
1779         if (iIndexSize >= iDimToCheck)
1780         {
1781             //size is good, now check datas
1782             double* pIndexes = getDoubleArrayFromDouble(pArg[i]);
1783             for (int j = 0; j < iDimToCheck; j++)
1784             {
1785                 bool bFind = false;
1786                 for (int k = 0; k < iIndexSize; k++)
1787                 {
1788                     if ((int)pIndexes[k] == j + 1)
1789                     {
1790                         bFind = true;
1791                         break;
1792                     }
1793                 }
1794                 pbFull[i] = bFind;
1795             }
1796         }
1797     }
1798
1799     //only one dims can be not full/entire
1800     bool bNotEntire = false;
1801     int iNotEntire = 0;
1802     bool bTooMuchNotEntire = false;
1803     for (int i = 0; i < iDims; i++)
1804     {
1805         if (pbFull[i] == false)
1806         {
1807             if (bNotEntire == false)
1808             {
1809                 bNotEntire = true;
1810                 iNotEntire = i;
1811             }
1812             else
1813             {
1814                 bTooMuchNotEntire = true;
1815                 break;
1816             }
1817         }
1818     }
1819
1820     delete[] pbFull;
1821
1822     if (bTooMuchNotEntire == true)
1823     {
1824         //free pArg content
1825         cleanIndexesArguments(_pArgs, &pArg);
1826         return NULL;
1827     }
1828
1829     //find index to keep
1830     int iNotEntireSize = pArg[iNotEntire]->getAs<GenericType>()->getSize();
1831     double* piNotEntireIndex = getDoubleArrayFromDouble(pArg[iNotEntire]);
1832     int iKeepSize = getVarMaxDim(iNotEntire, iDims);
1833     bool* pbKeep = new bool[iKeepSize];
1834
1835     //fill pbKeep with true value
1836     for (int i = 0; i < iKeepSize; i++)
1837     {
1838         pbKeep[i] = true;
1839     }
1840
1841     for (int i = 0; i < iNotEntireSize; i++)
1842     {
1843         int idx = (int)piNotEntireIndex[i] - 1;
1844
1845         //don't care of value out of bounds
1846         if (idx < iKeepSize)
1847         {
1848             pbKeep[idx] = false;
1849         }
1850     }
1851
1852     int iNewDimSize = 0;
1853     for (int i = 0; i < iKeepSize; i++)
1854     {
1855         if (pbKeep[i] == true)
1856         {
1857             iNewDimSize++;
1858         }
1859     }
1860     delete[] pbKeep;
1861
1862     int* piNewDims = new int[iDims];
1863     for (int i = 0; i < iDims; i++)
1864     {
1865         if (i == iNotEntire)
1866         {
1867             piNewDims[i] = iNewDimSize;
1868         }
1869         else
1870         {
1871             piNewDims[i] = getVarMaxDim(i, iDims);
1872         }
1873     }
1874
1875     //remove last dimension if are == 1
1876     int iOrigDims = iDims;
1877     for (int i = (iDims - 1); i >= 2; i--)
1878     {
1879         if (piNewDims[i] == 1)
1880         {
1881             iDims--;
1882         }
1883         else
1884         {
1885             break;
1886         }
1887     }
1888
1889     if (iNewDimSize == 0)
1890     {
1891         delete[] piNewDims;
1892         //free pArg content
1893         cleanIndexesArguments(_pArgs, &pArg);
1894         return new Sparse(0, 0);
1895     }
1896
1897     if (iDims == 1)
1898     {
1899         //two cases, depends of original matrix/vector
1900         if ((*_pArgs)[0]->isColon() == false && m_iDims == 2 && m_piDims[0] == 1 && m_piDims[1] != 1)
1901         {
1902             //special case for row vector
1903             pOut = new Sparse(1, iNewDimSize, isComplex());
1904             //in this case we have to care of 2nd dimension
1905             //iNotEntire = 1;
1906         }
1907         else
1908         {
1909             pOut = new Sparse(iNewDimSize, 1, isComplex());
1910         }
1911     }
1912     else
1913     {
1914         pOut = new Sparse(piNewDims[0], piNewDims[1], isComplex());
1915     }
1916
1917     delete[] piNewDims;
1918     //find a way to copy existing data to new variable ...
1919     int iNewPos = 0;
1920     int* piIndexes = new int[iOrigDims];
1921     int* piViewDims = new int[iOrigDims];
1922     for (int i = 0; i < iOrigDims; i++)
1923     {
1924         piViewDims[i] = getVarMaxDim(i, iOrigDims);
1925     }
1926
1927     for (int i = 0; i < getSize(); i++)
1928     {
1929         bool bByPass = false;
1930         getIndexesWithDims(i, piIndexes, piViewDims, iOrigDims);
1931
1932         //check if piIndexes use removed indexes
1933         for (int j = 0; j < iNotEntireSize; j++)
1934         {
1935             if ((piNotEntireIndex[j] - 1) == piIndexes[iNotEntire])
1936             {
1937                 //by pass this value
1938                 bByPass = true;
1939                 break;
1940             }
1941         }
1942
1943         if (bByPass == false)
1944         {
1945             //compute new index
1946             if (isComplex())
1947             {
1948                 pOut->set(iNewPos, getImg(i));
1949             }
1950             else
1951             {
1952                 pOut->set(iNewPos, get(i));
1953             }
1954             iNewPos++;
1955         }
1956     }
1957
1958     delete[] piIndexes;
1959     delete[] piViewDims;
1960
1961     //free pArg content
1962     cleanIndexesArguments(_pArgs, &pArg);
1963
1964     return pOut;
1965 }
1966
1967 Sparse* Sparse::append(int r, int c, types::Sparse SPARSE_CONST* src)
1968 {
1969     Sparse* pIT = checkRef(this, &Sparse::append, r, c, src);
1970     if (pIT != this)
1971     {
1972         return pIT;
1973     }
1974
1975     //        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;
1976     if (src->isComplex())
1977     {
1978         toComplex();
1979     }
1980     if (isComplex())
1981     {
1982         if (src->isComplex())
1983         {
1984             doAppend(*(src->matrixCplx), r, c, *matrixCplx);
1985         }
1986         else
1987         {
1988             doAppend(*(src->matrixReal), r, c, *matrixCplx);
1989         }
1990     }
1991     else
1992     {
1993         doAppend(*(src->matrixReal), r, c, *matrixReal);
1994     }
1995
1996     finalize();
1997
1998     return this; // realloc is meaningless for sparse matrices
1999 }
2000
2001 /*
2002 * create a new Sparse of dims according to resSize and fill it from currentSparse (along coords)
2003 */
2004 GenericType* Sparse::extract(typed_list* _pArgs)
2005 {
2006     Sparse* pOut = NULL;
2007     int iDims = (int)_pArgs->size();
2008     typed_list pArg;
2009
2010     int* piMaxDim = new int[iDims];
2011     int* piCountDim = new int[iDims];
2012
2013     //evaluate each argument and replace by appropriate value and compute the count of combinations
2014     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
2015     if (iSeqCount == 0)
2016     {
2017         //free pArg content
2018         cleanIndexesArguments(_pArgs, &pArg);
2019         if (_pArgs->size() == 0)
2020         {
2021             //a()
2022             delete[] piMaxDim;
2023             delete[] piCountDim;
2024             //free pArg content
2025             cleanIndexesArguments(_pArgs, &pArg);
2026             return this;
2027         }
2028         else
2029         {
2030             //a([])
2031             delete[] piMaxDim;
2032             delete[] piCountDim;
2033             //free pArg content
2034             cleanIndexesArguments(_pArgs, &pArg);
2035             return new types::Sparse(0,0,false);
2036         }
2037     }
2038
2039     if (iDims < 2)
2040     {
2041         if (piMaxDim[0] <= getSize())
2042         {
2043             int iNewRows = 0;
2044             int iNewCols = 0;
2045
2046             if (getRows() == 1 && getCols() != 1 && (*_pArgs)[0]->isColon() == false)
2047             {
2048                 //special case for row vector
2049                 iNewRows = 1;
2050                 iNewCols = piCountDim[0];
2051             }
2052             else
2053             {
2054                 iNewRows = piCountDim[0];
2055                 iNewCols = 1;
2056             }
2057
2058             double* pIdx = pArg[0]->getAs<Double>()->get();
2059             if (isComplex())
2060             {
2061                 bool bIsReal = true;
2062                 std::vector<CplxTriplet_t> tripletList;
2063                 std::vector<RealTriplet_t> realTripletList;
2064                 int row = getRows();
2065                 for (int i = 0; i < iSeqCount; i++)
2066                 {
2067                     int iRowRead = static_cast<int>(pIdx[i] - 1) % row;
2068                     int iColRead = static_cast<int>(pIdx[i] - 1) / row;
2069                     int iRowWrite = i % iNewRows;
2070                     int iColWrite = i / iNewRows;
2071
2072                     std::complex<double> dbl = getImg(iRowRead, iColRead);
2073                     if (dbl != 0.)
2074                     {
2075                         //only non zero values
2076                         tripletList.emplace_back(iRowWrite, iColWrite, dbl);
2077                         bIsReal &= (dbl.imag() == 0);
2078                         if (bIsReal)
2079                         {
2080                             realTripletList.emplace_back(iRowWrite, iColWrite, dbl.real());
2081                         }
2082                     }
2083                 }
2084                 if (bIsReal)
2085                 {
2086                     RealSparse_t* real = new RealSparse_t(iNewRows, iNewCols);
2087                     real->setFromTriplets(realTripletList.begin(), realTripletList.end(), DupFunctor<double>());
2088                     pOut = new Sparse(real, nullptr);
2089                 }
2090                 else
2091                 {
2092                     CplxSparse_t* cplx = new CplxSparse_t(iNewRows, iNewCols);
2093                     cplx->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<std::complex<double>>());
2094                     pOut = new Sparse(nullptr, cplx);
2095                 }
2096             }
2097             else
2098             {
2099                 RealSparse_t* real = new RealSparse_t(iNewRows, iNewCols);
2100                 std::vector<RealTriplet_t> tripletList;
2101                 int row = getRows();
2102                 for (int i = 0; i < iSeqCount; i++)
2103                 {
2104                     int iRowRead = static_cast<int>(pIdx[i] - 1) % row;
2105                     int iColRead = static_cast<int>(pIdx[i] - 1) / row;
2106                     int iRowWrite = i % iNewRows;
2107                     int iColWrite = i / iNewRows;
2108
2109                     double dbl = get(iRowRead, iColRead);
2110                     if (dbl != 0)
2111                     {
2112                         //only non zero values
2113                         tripletList.emplace_back(iRowWrite, iColWrite, dbl);
2114                     }
2115                 }
2116
2117                 real->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<double>());
2118                 pOut = new Sparse(real, nullptr);
2119             }
2120         }
2121         else
2122         {
2123             delete[] piMaxDim;
2124             delete[] piCountDim;
2125             //free pArg content
2126             cleanIndexesArguments(_pArgs, &pArg);
2127             return NULL;
2128         }
2129     }
2130     else
2131     {
2132         if (piMaxDim[0] <= getRows() && piMaxDim[1] <= getCols())
2133         {
2134             double* pIdxRow = pArg[0]->getAs<Double>()->get();
2135             double* pIdxCol = pArg[1]->getAs<Double>()->get();
2136
2137             int iNewRows = pArg[0]->getAs<Double>()->getSize();
2138             int iNewCols = pArg[1]->getAs<Double>()->getSize();
2139
2140             if (isComplex())
2141             {
2142                 bool bIsReal = true;
2143                 std::vector<CplxTriplet_t> tripletList;
2144                 std::vector<RealTriplet_t> realTripletList;
2145
2146                 for (int iRow = 0; iRow < iNewRows; iRow++)
2147                 {
2148                     for (int iCol = 0; iCol < iNewCols; iCol++)
2149                     {
2150                         std::complex<double> dbl = getImg((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
2151                         if (dbl != 0.)
2152                         {
2153                             //only non zero values
2154                             tripletList.emplace_back(iRow, iCol, dbl);
2155                             bIsReal &= (dbl.imag() == 0.);
2156                             if (bIsReal)
2157                             {
2158                                 realTripletList.emplace_back(iRow, iCol, dbl.real());
2159                             }
2160                         }
2161                     }
2162                 }
2163                 if (bIsReal)
2164                 {
2165                     RealSparse_t* real = new RealSparse_t(iNewRows, iNewCols);
2166                     real->setFromTriplets(realTripletList.begin(), realTripletList.end(), DupFunctor<double>());
2167                     pOut = new Sparse(real, nullptr);
2168                 }
2169                 else
2170                 {
2171                     CplxSparse_t* cplx = new CplxSparse_t(iNewRows, iNewCols);
2172                     cplx->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<std::complex<double>>());
2173                     pOut = new Sparse(nullptr, cplx);
2174                 }
2175             }
2176             else
2177             {
2178                 RealSparse_t* real = new RealSparse_t(iNewRows, iNewCols);
2179                 std::vector<RealTriplet_t> tripletList;
2180                 for (int iRow = 0; iRow < iNewRows; iRow++)
2181                 {
2182                     for (int iCol = 0; iCol < iNewCols; iCol++)
2183                     {
2184                         double dbl = get((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
2185                         if (dbl != 0.)
2186                         {
2187                             //only non zero values
2188                             tripletList.emplace_back(iRow, iCol, dbl);
2189                         }
2190                     }
2191                 }
2192
2193                 real->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<double>());
2194                 pOut = new Sparse(real, nullptr);
2195             }
2196         }
2197         else
2198         {
2199             delete[] piMaxDim;
2200             delete[] piCountDim;
2201             //free pArg content
2202             cleanIndexesArguments(_pArgs, &pArg);
2203             return NULL;
2204         }
2205     }
2206
2207     pOut->finalize();
2208
2209     delete[] piMaxDim;
2210     delete[] piCountDim;
2211     //free pArg content
2212     cleanIndexesArguments(_pArgs, &pArg);
2213
2214     return pOut;
2215 }
2216
2217 Sparse* Sparse::extract(int nbCoords, int SPARSE_CONST* coords, int SPARSE_CONST* maxCoords, int SPARSE_CONST* resSize, bool asVector) SPARSE_CONST
2218 {
2219     if ((asVector && maxCoords[0] > getSize()) ||
2220             (asVector == false && maxCoords[0] > getRows()) ||
2221             (asVector == false && maxCoords[1] > getCols()))
2222     {
2223         return 0;
2224     }
2225
2226     bool const cplx(isComplex());
2227     Sparse * pSp(0);
2228     if (asVector)
2229     {
2230         pSp = (getRows() == 1) ? new Sparse(1, resSize[0], cplx) : new Sparse(resSize[0], 1, cplx);
2231     }
2232     else
2233     {
2234         pSp = new Sparse(resSize[0], resSize[1], cplx);
2235     }
2236     //        std::cerr<<"extracted sparse:"<<pSp->getRows()<<", "<<pSp->getCols()<<"seqCount="<<nbCoords<<"maxDim="<<maxCoords[0] <<","<< maxCoords[1]<<std::endl;
2237     if (!(asVector
2238             ? copyToSparse(*this, Coords<true>(coords, getRows()), nbCoords
2239                            , *pSp, RowWiseFullIterator(pSp->getRows(), pSp->getCols()))
2240             : copyToSparse(*this, Coords<false>(coords), nbCoords
2241                            , *pSp, RowWiseFullIterator(pSp->getRows(), pSp->getCols()))))
2242     {
2243         delete pSp;
2244         pSp = NULL;
2245     }
2246     return pSp;
2247 }
2248
2249 bool Sparse::invoke(typed_list & in, optional_list & /*opt*/, int /*_iRetCount*/, typed_list & out, const ast::Exp & e)
2250 {
2251     if (in.size() == 0)
2252     {
2253         out.push_back(this);
2254     }
2255     else
2256     {
2257         InternalType * _out = extract(&in);
2258         if (!_out)
2259         {
2260             std::wostringstream os;
2261             os << _W("Invalid index.\n");
2262             throw ast::InternalError(os.str(), 999, e.getLocation());
2263         }
2264         out.push_back(_out);
2265     }
2266
2267     return true;
2268 }
2269
2270
2271 bool Sparse::isInvokable() const
2272 {
2273     return true;
2274 }
2275
2276 bool Sparse::hasInvokeOption() const
2277 {
2278     return false;
2279 }
2280
2281 int Sparse::getInvokeNbIn()
2282 {
2283     return -1;
2284 }
2285
2286 int Sparse::getInvokeNbOut()
2287 {
2288     return 1;
2289 }
2290
2291 /*
2292 coords are Scilab 1-based
2293 extract std::make_pair(coords, asVector), rowIter
2294 */
2295 template<typename Src, typename SrcTraversal, typename Sz, typename DestTraversal>
2296 bool Sparse::copyToSparse(Src SPARSE_CONST& src, SrcTraversal srcTrav, Sz n, Sparse& sp, DestTraversal destTrav)
2297 {
2298     if (!(src.isComplex() || sp.isComplex()))
2299     {
2300         mycopy_n(makeMatrixIterator<double>(src, srcTrav), n
2301                  , makeMatrixIterator<double>(*sp.matrixReal, destTrav));
2302     }
2303     else
2304     {
2305         sp.toComplex();
2306         mycopy_n(makeMatrixIterator<std::complex<double> >(src, srcTrav), n
2307                  , makeMatrixIterator<std::complex<double> >(*sp.matrixCplx, destTrav));
2308     }
2309
2310     sp.finalize();
2311     return true;
2312 }
2313
2314 // GenericType because we might return a Double* for scalar operand
2315 Sparse* Sparse::add(Sparse const& o) const
2316 {
2317     RealSparse_t* realSp(0);
2318     CplxSparse_t* cplxSp(0);
2319     if (isComplex() == false && o.isComplex() == false)
2320     {
2321         //R + R -> R
2322         realSp = new RealSparse_t(*matrixReal + * (o.matrixReal));
2323     }
2324     else if (isComplex() == false && o.isComplex() == true)
2325     {
2326         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >() + * (o.matrixCplx));
2327     }
2328     else if (isComplex() == true && o.isComplex() == false)
2329     {
2330         cplxSp = new CplxSparse_t(*matrixCplx + o.matrixReal->cast<std::complex<double> >());
2331     }
2332     else if (isComplex() == true && o.isComplex() == true)
2333     {
2334         cplxSp = new CplxSparse_t(*matrixCplx + * (o.matrixCplx));
2335     }
2336
2337     return new Sparse(realSp, cplxSp);
2338 }
2339
2340 Sparse* Sparse::substract(Sparse const& o) const
2341 {
2342     RealSparse_t* realSp(0);
2343     CplxSparse_t* cplxSp(0);
2344     if (isComplex() == false && o.isComplex() == false)
2345     {
2346         //R - R -> R
2347         realSp = new RealSparse_t(*matrixReal - * (o.matrixReal));
2348     }
2349     else if (isComplex() == false && o.isComplex() == true)
2350     {
2351         //R - C -> C
2352         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >() - * (o.matrixCplx));
2353     }
2354     else if (isComplex() == true && o.isComplex() == false)
2355     {
2356         //C - R -> C
2357         cplxSp = new CplxSparse_t(*matrixCplx - o.matrixReal->cast<std::complex<double> >());
2358     }
2359     else if (isComplex() == true && o.isComplex() == true)
2360     {
2361         //C - C -> C
2362         cplxSp = new CplxSparse_t(*matrixCplx - * (o.matrixCplx));
2363     }
2364
2365     return new Sparse(realSp, cplxSp);
2366 }
2367
2368 Sparse* Sparse::multiply(double s) const
2369 {
2370     return new Sparse(isComplex() ? 0 : new RealSparse_t((*matrixReal)*s)
2371                       , isComplex() ? new CplxSparse_t((*matrixCplx)*s) : 0);
2372 }
2373
2374 Sparse* Sparse::multiply(std::complex<double> s) const
2375 {
2376     return new Sparse(0
2377                       , isComplex() ? new CplxSparse_t((*matrixCplx) * s) : new CplxSparse_t((*matrixReal) * s));
2378 }
2379
2380 Sparse* Sparse::multiply(Sparse const& o) const
2381 {
2382     RealSparse_t* realSp(0);
2383     CplxSparse_t* cplxSp(0);
2384
2385     if (isComplex() == false && o.isComplex() == false)
2386     {
2387         realSp = new RealSparse_t(*matrixReal **(o.matrixReal));
2388     }
2389     else if (isComplex() == false && o.isComplex() == true)
2390     {
2391         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >() **(o.matrixCplx));
2392     }
2393     else if (isComplex() == true && o.isComplex() == false)
2394     {
2395         cplxSp = new CplxSparse_t(*matrixCplx * o.matrixReal->cast<std::complex<double> >());
2396     }
2397     else if (isComplex() == true && o.isComplex() == true)
2398     {
2399         cplxSp = new CplxSparse_t(*matrixCplx **(o.matrixCplx));
2400     }
2401
2402     return new Sparse(realSp, cplxSp);
2403 }
2404
2405 Sparse* Sparse::dotMultiply(Sparse SPARSE_CONST& o) const
2406 {
2407     RealSparse_t* realSp(0);
2408     CplxSparse_t* cplxSp(0);
2409     if (isComplex() == false && o.isComplex() == false)
2410     {
2411         realSp = new RealSparse_t(matrixReal->cwiseProduct(*(o.matrixReal)));
2412     }
2413     else if (isComplex() == false && o.isComplex() == true)
2414     {
2415         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >().cwiseProduct(*(o.matrixCplx)));
2416     }
2417     else if (isComplex() == true && o.isComplex() == false)
2418     {
2419         cplxSp = new CplxSparse_t(matrixCplx->cwiseProduct(o.matrixReal->cast<std::complex<double> >()));
2420     }
2421     else if (isComplex() == true && o.isComplex() == true)
2422     {
2423         cplxSp = new CplxSparse_t(matrixCplx->cwiseProduct(*(o.matrixCplx)));
2424     }
2425
2426     return new Sparse(realSp, cplxSp);
2427 }
2428
2429 Sparse* Sparse::dotDivide(Sparse SPARSE_CONST& o) const
2430 {
2431     RealSparse_t* realSp(0);
2432     CplxSparse_t* cplxSp(0);
2433     if (isComplex() == false && o.isComplex() == false)
2434     {
2435         realSp = new RealSparse_t(matrixReal->cwiseQuotient(*(o.matrixReal)));
2436     }
2437     else if (isComplex() == false && o.isComplex() == true)
2438     {
2439         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >().cwiseQuotient(*(o.matrixCplx)));
2440     }
2441     else if (isComplex() == true && o.isComplex() == false)
2442     {
2443         cplxSp = new CplxSparse_t(matrixCplx->cwiseQuotient(o.matrixReal->cast<std::complex<double> >()));
2444     }
2445     else if (isComplex() == true && o.isComplex() == true)
2446     {
2447         cplxSp = new CplxSparse_t(matrixCplx->cwiseQuotient(*(o.matrixCplx)));
2448     }
2449
2450     return new Sparse(realSp, cplxSp);
2451 }
2452
2453 int Sparse::newCholLLT(Sparse** _SpPermut, Sparse** _SpFactor) const
2454 {
2455     typedef Eigen::SparseMatrix<double, Eigen::ColMajor> RealSparseCol_t;
2456     RealSparseCol_t spColMajor = RealSparseCol_t((const RealSparse_t&) * matrixReal);
2457
2458     // Constructs and performs the LLT factorization of sparse
2459     Eigen::SimplicialLLT<RealSparseCol_t> pLLT(spColMajor);
2460     int iInfo = pLLT.info();
2461     if (iInfo != Eigen::Success)
2462     {
2463         *_SpFactor = NULL;
2464         *_SpPermut = NULL;
2465         return iInfo;
2466     }
2467
2468     // Get the lower matrix of factorization.
2469     // The new RealSparse_t will be set in Sparse without copy.
2470     *_SpFactor = new Sparse(new RealSparse_t(pLLT.matrixL()), NULL);
2471
2472     // Get the permutation matrix.
2473     Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic, int> p = pLLT.permutationP();
2474     *_SpPermut = new Sparse(static_cast<int>(p.rows()), static_cast<int>(p.cols()));
2475     for (int i = 0; i < p.rows(); i++)
2476     {
2477         (*_SpPermut)->set(i, p.indices()[i], 1, false);
2478     }
2479
2480     (*_SpPermut)->finalize();
2481
2482     return iInfo;
2483 }
2484
2485 bool Sparse::transpose(InternalType *& out)
2486 {
2487     out = new Sparse(matrixReal ? new RealSparse_t(matrixReal->transpose()) : 0, matrixCplx ? new CplxSparse_t(matrixCplx->transpose()) : 0);
2488     return true;
2489 }
2490
2491 bool Sparse::adjoint(InternalType *& out)
2492 {
2493     out = new Sparse(matrixReal ? new RealSparse_t(matrixReal->adjoint()) : 0, matrixCplx ? new CplxSparse_t(matrixCplx->adjoint()) : 0);
2494     return true;
2495 }
2496
2497 struct BoolCast
2498 {
2499     BoolCast(std::complex<double> const& c) : b(c.real() || c.imag()) {}
2500     operator bool() const
2501     {
2502         return b;
2503     }
2504     operator double() const
2505     {
2506         return b ? 1. : 0.;
2507     }
2508     bool b;
2509 };
2510 Sparse* Sparse::newOnes() const
2511 {
2512     // result is never cplx
2513     return new Sparse(matrixReal
2514                       ? new RealSparse_t(matrixReal->cast<bool>().cast<double>())
2515                       : new RealSparse_t(matrixCplx->cast<BoolCast>().cast<double>())
2516                       , 0);
2517 }
2518
2519 struct RealCast
2520 {
2521     RealCast(std::complex<double> const& c) : b(c.real()) {}
2522     operator bool() const
2523     {
2524         return b != 0;
2525     }
2526     operator double() const
2527     {
2528         return b;
2529     }
2530     double b;
2531 };
2532 Sparse* Sparse::newReal() const
2533 {
2534     return new Sparse(matrixReal
2535                       ? matrixReal
2536                       : new RealSparse_t(matrixCplx->cast<RealCast>().cast<double>())
2537                       , 0);
2538 }
2539
2540 std::size_t Sparse::nonZeros() const
2541 {
2542     if (isComplex())
2543     {
2544         return matrixCplx->nonZeros();
2545     }
2546     else
2547     {
2548         return matrixReal->nonZeros();
2549     }
2550 }
2551 std::size_t Sparse::nonZeros(std::size_t r) const
2552 {
2553     std::size_t res;
2554     if (matrixReal)
2555     {
2556         int* piIndex = matrixReal->outerIndexPtr();
2557         res = piIndex[r + 1] - piIndex[r];
2558     }
2559     else
2560     {
2561         int* piIndex = matrixCplx->outerIndexPtr();
2562         res = piIndex[r + 1] - piIndex[r];
2563     }
2564
2565     return res;
2566 }
2567
2568 int* Sparse::getNbItemByRow(int* _piNbItemByRows)
2569 {
2570     int* piNbItemByCols = new int[getRows() + 1];
2571     if (isComplex())
2572     {
2573         mycopy_n(matrixCplx->outerIndexPtr(), getRows() + 1, piNbItemByCols);
2574     }
2575     else
2576     {
2577         mycopy_n(matrixReal->outerIndexPtr(), getRows() + 1, piNbItemByCols);
2578     }
2579
2580     for (int i = 0; i < getRows(); i++)
2581     {
2582         _piNbItemByRows[i] = piNbItemByCols[i + 1] - piNbItemByCols[i];
2583     }
2584
2585     delete[] piNbItemByCols;
2586     return _piNbItemByRows;
2587 }
2588
2589 int* Sparse::getColPos(int* _piColPos)
2590 {
2591     if (isComplex())
2592     {
2593         mycopy_n(matrixCplx->innerIndexPtr(), nonZeros(), _piColPos);
2594     }
2595     else
2596     {
2597         mycopy_n(matrixReal->innerIndexPtr(), nonZeros(), _piColPos);
2598     }
2599
2600     for (size_t i = 0; i < nonZeros(); i++)
2601     {
2602         _piColPos[i]++;
2603     }
2604
2605     return _piColPos;
2606 }
2607
2608 int* Sparse::getInnerPtr(int* count)
2609 {
2610     int* ret = nullptr;
2611     if (isComplex())
2612     {
2613         ret = matrixCplx->innerIndexPtr();
2614         *count = static_cast<int>(matrixCplx->innerSize());
2615     }
2616     else
2617     {
2618         ret = matrixReal->innerIndexPtr();
2619         *count = static_cast<int>(matrixReal->innerSize());
2620     }
2621
2622     return ret;
2623 }
2624
2625 int* Sparse::getOuterPtr(int* count)
2626 {
2627     int* ret = nullptr;
2628     if (isComplex())
2629     {
2630         ret = matrixCplx->outerIndexPtr();
2631         *count = static_cast<int>(matrixCplx->outerSize());
2632     }
2633     else
2634     {
2635         ret = matrixReal->outerIndexPtr();
2636         *count = static_cast<int>(matrixReal->outerSize());
2637     }
2638
2639     return ret;
2640 }
2641
2642 namespace
2643 {
2644 template<typename S> struct GetReal : std::unary_function<typename S::InnerIterator, double>
2645 {
2646     double operator()(typename S::InnerIterator it) const
2647     {
2648         return it.value();
2649     }
2650 };
2651 template<> struct GetReal< Eigen::SparseMatrix<std::complex<double >, Eigen::RowMajor > >
2652     : std::unary_function<Sparse::CplxSparse_t::InnerIterator, double>
2653 {
2654     double operator()(Sparse::CplxSparse_t::InnerIterator it) const
2655     {
2656         return it.value().real();
2657     }
2658 };
2659 template<typename S> struct GetImag : std::unary_function<typename S::InnerIterator, double>
2660 {
2661     double operator()(typename S::InnerIterator it) const
2662     {
2663         return it.value().imag();
2664     }
2665 };
2666 template<typename S> struct GetRow : std::unary_function<typename S::InnerIterator, int>
2667 {
2668     int operator()(typename S::InnerIterator it) const
2669     {
2670         return static_cast<int>(it.row() + 1);
2671     }
2672 };
2673 template<typename S> struct GetCol : std::unary_function<typename S::InnerIterator, int>
2674 {
2675     int operator()(typename S::InnerIterator it) const
2676     {
2677         return static_cast<int>(it.col() + 1);
2678     }
2679 };
2680
2681 template<typename S, typename Out, typename F> Out sparseTransform(S& s, Out o, F f)
2682 {
2683     int size = static_cast<int>(s.outerSize());
2684     for (std::size_t k(0); k < size; ++k)
2685     {
2686         for (typename S::InnerIterator it(s, (int)k); it; ++it, ++o)
2687         {
2688             *o = f(it);
2689         }
2690     }
2691     return o;
2692 }
2693 }
2694
2695 std::pair<double*, double*> Sparse::outputValues(double* outReal, double* outImag)const
2696 {
2697     return matrixReal
2698            ? std::make_pair(sparseTransform(*matrixReal, outReal, GetReal<RealSparse_t>()), outImag)
2699            : std::make_pair(sparseTransform(*matrixCplx, outReal, GetReal<CplxSparse_t>())
2700                             , sparseTransform(*matrixCplx, outImag, GetImag<CplxSparse_t>()));
2701 }
2702
2703 int* Sparse::outputRowCol(int* out)const
2704 {
2705     return matrixReal
2706            ? sparseTransform(*matrixReal, sparseTransform(*matrixReal, out, GetRow<RealSparse_t>()), GetCol<RealSparse_t>())
2707            : sparseTransform(*matrixCplx, sparseTransform(*matrixCplx, out, GetRow<CplxSparse_t>()), GetCol<CplxSparse_t>());
2708 }
2709 double* Sparse::outputCols(double* out) const
2710 {
2711     if (isComplex())
2712     {
2713         mycopy_n(matrixCplx->innerIndexPtr(), nonZeros(), out);
2714     }
2715     else
2716     {
2717         mycopy_n(matrixReal->innerIndexPtr(), nonZeros(), out);
2718     }
2719
2720     return std::transform(out, out, out, std::bind(std::plus<double>(), std::placeholders::_1, 1));
2721
2722 }
2723
2724 void Sparse::opposite(void)
2725 {
2726     if (isComplex())
2727     {
2728         std::complex<double>* data = matrixCplx->valuePtr();
2729         std::transform(data, data + matrixCplx->nonZeros(), data, std::negate<std::complex<double> >());
2730     }
2731     else
2732     {
2733         double* data = matrixReal->valuePtr();
2734         std::transform(data, data + matrixReal->nonZeros(), data, std::negate<double>());
2735     }
2736 }
2737
2738 SparseBool* Sparse::newLessThan(Sparse &o)
2739 {
2740     //only real values !
2741
2742     //return cwiseOp<std::less>(*this, o);
2743     int rowL = getRows();
2744     int colL = getCols();
2745
2746     int rowR = o.getRows();
2747     int colR = o.getCols();
2748     int row = std::max(rowL, rowR);
2749     int col = std::max(colL, colR);
2750
2751     //create a boolean sparse matrix with dims of sparses
2752     types::SparseBool* ret = new types::SparseBool(row, col);
2753     if (isScalar() && o.isScalar())
2754     {
2755         double l = get(0, 0);
2756         double r = o.get(0, 0);
2757         ret->set(0, 0, l < r, false);
2758     }
2759     else if (isScalar())
2760     {
2761         int nnzR = static_cast<int>(o.nonZeros());
2762         std::vector<int> rowcolR(nnzR * 2, 0);
2763         o.outputRowCol(rowcolR.data());
2764
2765         //compare all items of R with R[0]
2766         double l = get(0, 0);
2767         if (l < 0)
2768         {
2769             //set true
2770             ret->setTrue(false);
2771         }
2772
2773         for (int i = 0; i < nnzR; ++i)
2774         {
2775             double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2776             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l < r, false);
2777         }
2778     }
2779     else if (o.isScalar())
2780     {
2781         int nnzL = static_cast<int>(nonZeros());
2782         std::vector<int> rowcolL(nnzL * 2, 0);
2783         outputRowCol(rowcolL.data());
2784
2785         double r = o.get(0, 0);
2786         if (r > 0)
2787         {
2788             ret->setTrue(true);
2789         }
2790
2791         for (int i = 0; i < nnzL; ++i)
2792         {
2793             double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2794             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l < r, false);
2795         }
2796     }
2797     else
2798     {
2799         int nnzR = static_cast<int>(o.nonZeros());
2800         std::vector<int> rowcolR(nnzR * 2, 0);
2801         o.outputRowCol(rowcolR.data());
2802         int nnzL = static_cast<int>(nonZeros());
2803         std::vector<int> rowcolL(nnzL * 2, 0);
2804         outputRowCol(rowcolL.data());
2805         //set all values to %t
2806         ret->setFalse(false);
2807
2808         for (int i = 0; i < nnzL; ++i)
2809         {
2810             double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2811             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l < 0, false);
2812         }
2813         ret->finalize();
2814
2815         //set _pR[i] == _pL[i] for each _pR values
2816         for (int i = 0; i < nnzR; ++i)
2817         {
2818             //get l and r following non zeros value of R
2819             double l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2820             double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2821             //set value following non zeros value of R
2822             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l < r, false);
2823         }
2824     }
2825
2826     ret->finalize();
2827     return ret;
2828 }
2829
2830 SparseBool* Sparse::newNotEqualTo(Sparse const&o) const
2831 {
2832     return cwiseOp<std::not_equal_to>(*this, o);
2833 }
2834
2835 SparseBool* Sparse::newLessOrEqual(Sparse &o)
2836 {
2837     //only real values !
2838
2839     //return cwiseOp<std::less>(*this, o);
2840     int rowL = getRows();
2841     int colL = getCols();
2842
2843     int rowR = o.getRows();
2844     int colR = o.getCols();
2845     int row = std::max(rowL, rowR);
2846     int col = std::max(colL, colR);
2847
2848     //create a boolean sparse matrix with dims of sparses
2849     types::SparseBool* ret = new types::SparseBool(row, col);
2850     if (isScalar() && o.isScalar())
2851     {
2852         double l = get(0, 0);
2853         double r = o.get(0, 0);
2854         ret->set(0, 0, l <= r, false);
2855     }
2856     else if (isScalar())
2857     {
2858         int nnzR = static_cast<int>(o.nonZeros());
2859         std::vector<int> rowcolR(nnzR * 2, 0);
2860         o.outputRowCol(rowcolR.data());
2861
2862         //compare all items of R with R[0]
2863         double l = get(0, 0);
2864         if (l <= 0)
2865         {
2866             //set true
2867             ret->setTrue(false);
2868         }
2869
2870         for (int i = 0; i < nnzR; ++i)
2871         {
2872             double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2873             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l <= r, false);
2874         }
2875     }
2876     else if (o.isScalar())
2877     {
2878         int nnzL = static_cast<int>(nonZeros());
2879         std::vector<int> rowcolL(nnzL * 2, 0);
2880         outputRowCol(rowcolL.data());
2881
2882         double r = o.get(0, 0);
2883         if (r >= 0)
2884         {
2885             ret->setTrue(true);
2886         }
2887
2888         for (int i = 0; i < nnzL; ++i)
2889         {
2890             double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2891             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l <= r, false);
2892         }
2893     }
2894     else
2895     {
2896         int nnzR = static_cast<int>(o.nonZeros());
2897         std::vector<int> rowcolR(nnzR * 2, 0);
2898         o.outputRowCol(rowcolR.data());
2899         int nnzL = static_cast<int>(nonZeros());
2900         std::vector<int> rowcolL(nnzL * 2, 0);
2901         outputRowCol(rowcolL.data());
2902         //set all values to %t
2903         ret->setTrue(false);
2904
2905         for (int i = 0; i < nnzL; ++i)
2906         {
2907             double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2908             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l <= 0, false);
2909         }
2910         ret->finalize();
2911
2912         //set _pR[i] == _pL[i] for each _pR values
2913         for (int i = 0; i < nnzR; ++i)
2914         {
2915             //get l and r following non zeros value of R
2916             double l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2917             double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2918             //set value following non zeros value of R
2919             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l <= r, false);
2920         }
2921     }
2922
2923     ret->finalize();
2924     return ret;
2925 }
2926
2927 SparseBool* Sparse::newEqualTo(Sparse &o)
2928 {
2929     int rowL = getRows();
2930     int colL = getCols();
2931
2932     int rowR = o.getRows();
2933     int colR = o.getCols();
2934     int row = std::max(rowL, rowR);
2935     int col = std::max(colL, colR);
2936
2937     //create a boolean sparse matrix with dims of sparses
2938     types::SparseBool* ret = new types::SparseBool(row, col);
2939     if (isScalar() && o.isScalar())
2940     {
2941         if (isComplex() || o.isComplex())
2942         {
2943             std::complex<double> l = getImg(0, 0);
2944             std::complex<double> r = o.getImg(0, 0);
2945             ret->set(0, 0, l == r, false);
2946         }
2947         else
2948         {
2949             double l = get(0, 0);
2950             double r = o.get(0, 0);
2951             ret->set(0, 0, l == r, false);
2952         }
2953     }
2954     else if (isScalar())
2955     {
2956         int nnzR = static_cast<int>(o.nonZeros());
2957         std::vector<int> rowcolR(nnzR * 2, 0);
2958         o.outputRowCol(rowcolR.data());
2959
2960         //compare all items of R with R[0]
2961         if (isComplex() || o.isComplex())
2962         {
2963             std::complex<double> l = getImg(0, 0);
2964             for (int i = 0; i < nnzR; ++i)
2965             {
2966                 std::complex<double> r = o.getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2967                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2968             }
2969         }
2970         else
2971         {
2972             double l = get(0, 0);
2973             for (int i = 0; i < nnzR; ++i)
2974             {
2975                 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2976                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2977             }
2978         }
2979     }
2980     else if (o.isScalar())
2981     {
2982         int nnzL = static_cast<int>(nonZeros());
2983         std::vector<int> rowcolL(nnzL * 2, 0);
2984         outputRowCol(rowcolL.data());
2985
2986         if (isComplex() || o.isComplex())
2987         {
2988             std::complex<double> r = o.getImg(0, 0);
2989             for (int i = 0; i < nnzL; ++i)
2990             {
2991                 std::complex<double> l = getImg(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2992                 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
2993             }
2994         }
2995         else
2996         {
2997             double r = get(0, 0);
2998             for (int i = 0; i < nnzL; ++i)
2999             {
3000                 double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
3001                 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
3002             }
3003         }
3004     }
3005     else
3006     {
3007         int nnzR = static_cast<int>(o.nonZeros());
3008         std::vector<int> rowcolR(nnzR * 2, 0);
3009         o.outputRowCol(rowcolR.data());
3010         int nnzL = static_cast<int>(nonZeros());
3011         std::vector<int> rowcolL(nnzL * 2, 0);
3012         outputRowCol(rowcolL.data());
3013         //set all values to %t
3014         ret->setTrue(false);
3015         //set %f in each pL values
3016         for (int i = 0; i < nnzL; ++i)
3017         {
3018             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, false, false);
3019         }
3020         ret->finalize();
3021
3022         //set _pR[i] == _pL[i] for each _pR values
3023         if (isComplex() || o.isComplex())
3024         {
3025             for (int i = 0; i < nnzR; ++i)
3026             {
3027                 //get l and r following non zeros value of R
3028                 std::complex<double> l = getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
3029                 std::complex<double> r = o.getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
3030                 //set value following non zeros value of R
3031                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
3032             }
3033         }
3034         else
3035         {
3036             for (int i = 0; i < nnzR; ++i)
3037             {
3038                 //get l and r following non zeros value of R
3039                 double l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
3040                 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
3041                 //set value following non zeros value of R
3042                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
3043             }
3044         }
3045     }
3046
3047     ret->finalize();
3048     return ret;
3049 }
3050
3051 Sparse* Sparse::reshape(int* _piDims, int _iDims)
3052 {
3053     Sparse* pSp = NULL;
3054     int iCols = 1;
3055
3056     if (_iDims == 2)
3057     {
3058         iCols = _piDims[1];
3059     }
3060
3061     if (_iDims <= 2)
3062     {
3063         pSp = reshape(_piDims[0], iCols);
3064     }
3065
3066     return pSp;
3067 }
3068
3069 Sparse* Sparse::reshape(int _iNewRows, int _iNewCols)
3070 {
3071     typedef Sparse* (Sparse::*reshape_t)(int, int);
3072     Sparse* pIT = checkRef(this, (reshape_t)&Sparse::reshape, _iNewRows, _iNewCols);
3073     if (pIT != this)
3074     {
3075         return pIT;
3076     }
3077
3078     if (_iNewRows * _iNewCols != getRows() * getCols())
3079     {
3080         return NULL;
3081     }
3082
3083     Sparse* res = NULL;
3084     try
3085     {
3086         if (matrixReal)
3087         {
3088             //item count
3089             size_t iNonZeros = nonZeros();
3090             RealSparse_t *newReal = new RealSparse_t(_iNewRows, _iNewCols);
3091             newReal->reserve((int)iNonZeros);
3092
3093             //coords
3094             int* pRows = new int[iNonZeros * 2];
3095             outputRowCol(pRows);
3096             int* pCols = pRows + iNonZeros;
3097
3098             //values
3099             double* pNonZeroR = new double[iNonZeros];
3100             double* pNonZeroI = new double[iNonZeros];
3101             outputValues(pNonZeroR, pNonZeroI);
3102
3103             std::vector<RealTriplet_t> tripletList;
3104             for (size_t i = 0; i < iNonZeros; i++)
3105             {
3106                 int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
3107                 tripletList.emplace_back((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), pNonZeroR[i]);
3108             }
3109
3110             newReal->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<double>());
3111
3112             delete matrixReal;
3113             matrixReal = newReal;
3114             delete[] pRows;
3115             delete[] pNonZeroR;
3116             delete[] pNonZeroI;
3117         }
3118         else
3119         {
3120             //item count
3121             size_t iNonZeros = nonZeros();
3122             CplxSparse_t *newCplx = new CplxSparse_t(_iNewRows, _iNewCols);
3123             newCplx->reserve((int)iNonZeros);
3124
3125             //coords
3126             int* pRows = new int[iNonZeros * 2];
3127             outputRowCol(pRows);
3128             int* pCols = pRows + iNonZeros;
3129
3130             //values
3131             double* pNonZeroR = new double[iNonZeros];
3132             double* pNonZeroI = new double[iNonZeros];
3133             outputValues(pNonZeroR, pNonZeroI);
3134
3135             std::vector<CplxTriplet_t> tripletList;
3136
3137             for (size_t i = 0; i < iNonZeros; i++)
3138             {
3139                 int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
3140                 tripletList.emplace_back((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), std::complex<double>(pNonZeroR[i], pNonZeroI[i]));
3141             }
3142
3143             newCplx->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<std::complex<double>>());
3144
3145             delete matrixCplx;
3146             matrixCplx = newCplx;
3147             delete[] pRows;
3148             delete[] pNonZeroR;
3149             delete[] pNonZeroI;
3150         }
3151
3152         m_iRows = _iNewRows;
3153         m_iCols = _iNewCols;
3154         m_iSize = _iNewRows * _iNewCols;
3155
3156         m_iDims = 2;
3157         m_piDims[0] = m_iRows;
3158         m_piDims[1] = m_iCols;
3159
3160         finalize();
3161
3162         res = this;
3163     }
3164     catch (...)
3165     {
3166         res = NULL;
3167     }
3168     return res;
3169 }
3170
3171 //    SparseBool* SparseBool::new
3172
3173 SparseBool::SparseBool(Bool SPARSE_CONST& src)
3174 {
3175     //compute idx
3176     int size = src.getSize();
3177     int row = src.getRows();
3178     Double* idx = new Double(src.getSize(), 2);
3179     double* p = idx->get();
3180     for (int i = 0; i < size; ++i)
3181     {
3182         p[i] = (double)(i % row) + 1;
3183         p[i + size] = (double)(i / row) + 1;
3184     }
3185     create2(src.getRows(), src.getCols(), src, *idx);
3186     idx->killMe();
3187 #ifndef NDEBUG
3188     Inspector::addItem(this);
3189 #endif
3190 }
3191 /* @param src : Bool matrix to copy into a new sparse matrix
3192 @param idx : Double matrix to use as indexes to get values from the src
3193 **/
3194 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
3195 {
3196     int idxrow = idx.getRows();
3197     int rows = static_cast<int>(*std::max_element(idx.get(), idx.get() + idxrow));
3198     int cols = static_cast<int>(*std::max_element(idx.get() + idxrow, idx.get() + idxrow * 2));
3199     create2(rows, cols, src, idx);
3200 #ifndef NDEBUG
3201     Inspector::addItem(this);
3202 #endif
3203 }
3204
3205 /* @param src : Bool matrix to copy into a new sparse matrix
3206 @param idx : Double matrix to use as indexes to get values from the src
3207 @param dims : Double matrix containing the dimensions of the new matrix
3208 **/
3209 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx, Double SPARSE_CONST& dims)
3210 {
3211     create2(static_cast<int>(dims.get(0)), static_cast<int>(dims.get(1)), src, idx);
3212 #ifndef NDEBUG
3213     Inspector::addItem(this);
3214 #endif
3215 }
3216
3217 SparseBool::SparseBool(int _iRows, int _iCols) : matrixBool(new BoolSparse_t(_iRows, _iCols))
3218 {
3219     m_iRows = _iRows;
3220     m_iCols = _iCols;
3221     m_iSize = _iRows * _iCols;
3222     m_iDims = 2;
3223     m_piDims[0] = _iRows;
3224     m_piDims[1] = _iCols;
3225 #ifndef NDEBUG
3226     Inspector::addItem(this);
3227 #endif
3228 }
3229
3230 SparseBool::SparseBool(SparseBool const& src) : matrixBool(new BoolSparse_t(*src.matrixBool))
3231 {
3232     m_iDims = 2;
3233     m_iRows = const_cast<SparseBool*>(&src)->getRows();
3234     m_iCols = const_cast<SparseBool*>(&src)->getCols();
3235     m_iSize = m_iRows * m_iCols;
3236     m_piDims[0] = m_iRows;
3237     m_piDims[1] = m_iCols;
3238 #ifndef NDEBUG
3239     Inspector::addItem(this);
3240 #endif
3241 }
3242
3243 SparseBool::SparseBool(BoolSparse_t* src) : matrixBool(src)
3244 {
3245     m_iRows = static_cast<int>(src->rows());
3246     m_iCols = static_cast<int>(src->cols());
3247     m_iSize = m_iRows * m_iCols;
3248     m_iDims = 2;
3249     m_piDims[0] = m_iRows;
3250     m_piDims[1] = m_iCols;
3251 #ifndef NDEBUG
3252     Inspector::addItem(this);
3253 #endif
3254 }
3255
3256 SparseBool::SparseBool(int rows, int cols, int trues, int* inner, int* outer)
3257 {
3258     int* out = nullptr;
3259     int* in = nullptr;
3260
3261     matrixBool = new BoolSparse_t(rows, cols);
3262     matrixBool->reserve((int)trues);
3263     out = matrixBool->outerIndexPtr();
3264     in = matrixBool->innerIndexPtr();
3265
3266     //update outerIndexPtr
3267     memcpy(out, outer, sizeof(int) * (rows + 1));
3268     //update innerIndexPtr
3269     memcpy(in, inner, sizeof(int) * trues);
3270
3271     bool* data = matrixBool->valuePtr();
3272     for (int i = 0; i < trues; ++i)
3273     {
3274         data[i] = true;
3275     }
3276
3277     m_iCols = cols;
3278     m_iRows = rows;
3279     m_iSize = cols * rows;
3280     m_iDims = 2;
3281     m_piDims[0] = m_iRows;
3282     m_piDims[1] = m_iCols;
3283
3284     matrixBool->resizeNonZeros(trues);
3285 }
3286
3287 bool SparseBool::getMemory(long long *_piSize, long long* _piSizePlusType)
3288 {
3289     *_piSize = nbTrue() * sizeof(bool);
3290     *_piSizePlusType = *_piSize + sizeof(*this);
3291     return true;
3292 }
3293
3294 void SparseBool::create2(int rows, int cols, Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
3295 {
3296     int nnz = src.getSize();
3297     double* i = idx.get();
3298     double* j = i + idx.getRows();
3299     int* val = src.get();
3300
3301     std::vector<BoolTriplet_t> tripletList;
3302     tripletList.reserve((int)nnz);
3303
3304     for (int k = 0; k < nnz; ++k)
3305     {
3306         tripletList.emplace_back(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, val[k] == 1);
3307     }
3308
3309     matrixBool = new BoolSparse_t(rows, cols);
3310     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
3311
3312     m_iRows = static_cast<int>(matrixBool->rows());
3313     m_iCols = static_cast<int>(matrixBool->cols());
3314     m_iSize = cols * rows;
3315     m_iDims = 2;
3316     m_piDims[0] = m_iRows;
3317     m_piDims[1] = m_iCols;
3318     finalize();
3319 }
3320
3321 SparseBool::~SparseBool()
3322 {
3323     delete matrixBool;
3324 #ifndef NDEBUG
3325     Inspector::removeItem(this);
3326 #endif
3327 }
3328
3329 bool SparseBool::toString(std::wostringstream& ostr)
3330 {
3331     ostr << ::toString(*matrixBool, 0);
3332     return true;
3333 }
3334
3335 void SparseBool::whoAmI() SPARSE_CONST
3336 {
3337     std::cout << "types::SparseBool";
3338 }
3339
3340 SparseBool* SparseBool::clone(void)
3341 {
3342     return new SparseBool(*this);
3343 }
3344
3345 SparseBool* SparseBool::resize(int _iNewRows, int _iNewCols)
3346 {
3347     typedef SparseBool* (SparseBool::*resize_t)(int, int);
3348     SparseBool* pIT = checkRef(this, (resize_t)&SparseBool::resize, _iNewRows, _iNewCols);
3349     if (pIT != this)
3350     {
3351         return pIT;
3352     }
3353
3354     if (_iNewRows <= getRows() && _iNewCols <= getCols())
3355     {
3356         //nothing to do: hence we do NOT fail
3357         return this;
3358     }
3359
3360     SparseBool* res = NULL;
3361     try
3362     {
3363         //item count
3364         size_t iNonZeros = nbTrue();
3365
3366         BoolSparse_t *newBool = new BoolSparse_t(_iNewRows, _iNewCols);
3367         newBool->reserve((int)iNonZeros);
3368
3369         //coords
3370         int* pRows = new int[iNonZeros * 2];
3371         outputRowCol(pRows);
3372         int* pCols = pRows + iNonZeros;
3373
3374         std::vector<BoolTriplet_t> tripletList;
3375
3376         for (size_t i = 0; i < iNonZeros; i++)
3377         {
3378             tripletList.emplace_back((int)pRows[i] - 1, (int)pCols[i] - 1, true);
3379         }
3380
3381         newBool->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<bool>());
3382
3383         delete matrixBool;
3384         matrixBool = newBool;
3385         delete[] pRows;
3386
3387         m_iRows = _iNewRows;
3388         m_iCols = _iNewCols;
3389         m_iSize = _iNewRows * _iNewCols;
3390         m_piDims[0] = m_iRows;
3391         m_piDims[1] = m_iCols;
3392
3393         res = this;
3394     }
3395     catch (...)
3396     {
3397         res = NULL;
3398     }
3399     return res;
3400 }
3401
3402 SparseBool* SparseBool::insert(typed_list* _pArgs, SparseBool* _pSource)
3403 {
3404     bool bNeedToResize = false;
3405     int iDims = (int)_pArgs->size();
3406     if (iDims > 2)
3407     {
3408         //sparse are only in 2 dims
3409         return NULL;
3410     }
3411
3412     typed_list pArg;
3413
3414     int piMaxDim[2];
3415     int piCountDim[2];
3416
3417     //on case of resize
3418     int iNewRows = 0;
3419     int iNewCols = 0;
3420
3421     //evaluate each argument and replace by appropriate value and compute the count of combinations
3422     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3423     if (iSeqCount == 0)
3424     {
3425         //free pArg content
3426         cleanIndexesArguments(_pArgs, &pArg);
3427         return this;
3428     }
3429
3430     if (iDims < 2)
3431     {
3432         //see as vector
3433         if (getRows() == 1 || getCols() == 1)
3434         {
3435             //vector or scalar
3436             if (getRows() * getCols() < piMaxDim[0])
3437             {
3438                 bNeedToResize = true;
3439
3440                 //need to enlarge sparse dimensions
3441                 if (getCols() == 1 || getSize() == 0)
3442                 {
3443                     //column vector
3444                     iNewRows = piMaxDim[0];
3445                     iNewCols = 1;
3446                 }
3447                 else if (getRows() == 1)
3448                 {
3449                     //row vector
3450                     iNewRows = 1;
3451                     iNewCols = piMaxDim[0];
3452                 }
3453             }
3454         }
3455         else if (getRows() * getCols() < piMaxDim[0])
3456         {
3457             //free pArg content
3458             cleanIndexesArguments(_pArgs, &pArg);
3459             //out of range
3460             return NULL;
3461         }
3462     }
3463     else
3464     {
3465         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
3466         {
3467             bNeedToResize = true;
3468             iNewRows = std::max(getRows(), piMaxDim[0]);
3469             iNewCols = std::max(getCols(), piMaxDim[1]);
3470         }
3471     }
3472
3473     //check number of insertion
3474     if (_pSource->isScalar() == false && _pSource->getSize() != iSeqCount)
3475     {
3476         //free pArg content
3477         cleanIndexesArguments(_pArgs, &pArg);
3478         return NULL;
3479     }
3480
3481     //now you are sure to be able to insert values
3482     if (bNeedToResize)
3483     {
3484         if (resize(iNewRows, iNewCols) == NULL)
3485         {
3486             //free pArg content
3487             cleanIndexesArguments(_pArgs, &pArg);
3488             return NULL;
3489         }
3490     }
3491
3492     if (iDims == 1)
3493     {
3494         double* pIdx = pArg[0]->getAs<Double>()->get();
3495         for (int i = 0; i < iSeqCount; i++)
3496         {
3497             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
3498             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
3499
3500             if (_pSource->isScalar())
3501             {
3502                 set(iRow, iCol, _pSource->get(0, 0), false);
3503             }
3504             else
3505             {
3506                 int iRowOrig = i % _pSource->getRows();
3507                 int iColOrig = i / _pSource->getRows();
3508                 set(iRow, iCol, _pSource->get(iRowOrig, iColOrig), false);
3509             }
3510         }
3511     }
3512     else
3513     {
3514         double* pIdxRow = pArg[0]->getAs<Double>()->get();
3515         int iRowSize = pArg[0]->getAs<Double>()->getSize();
3516         double* pIdxCol = pArg[1]->getAs<Double>()->get();
3517
3518         for (int i = 0; i < iSeqCount; i++)
3519         {
3520             if (_pSource->isScalar())
3521             {
3522                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(0, 0), false);
3523             }
3524             else
3525             {
3526                 int iRowOrig = i % _pSource->getRows();
3527                 int iColOrig = i / _pSource->getRows();
3528                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(iRowOrig, iColOrig), false);
3529             }
3530         }
3531     }
3532
3533     finalize();
3534
3535     //free pArg content
3536     cleanIndexesArguments(_pArgs, &pArg);
3537
3538     return this;
3539 }
3540
3541 SparseBool* SparseBool::insert(typed_list* _pArgs, InternalType* _pSource)
3542 {
3543     typedef SparseBool* (SparseBool::*insert_t)(typed_list*, InternalType*);
3544     SparseBool* pIT = checkRef(this, (insert_t)&SparseBool::insert, _pArgs, _pSource);
3545
3546
3547     if (pIT != this)
3548     {
3549         return pIT;
3550     }
3551
3552     if (_pSource->isSparseBool())
3553     {
3554         return insert(_pArgs, _pSource->getAs<SparseBool>());
3555     }
3556
3557     bool bNeedToResize  = false;
3558     int iDims           = (int)_pArgs->size();
3559     if (iDims > 2)
3560     {
3561         //sparse are only in 2 dims
3562         return NULL;
3563     }
3564
3565     typed_list pArg;
3566
3567     int piMaxDim[2];
3568     int piCountDim[2];
3569
3570     //on case of resize
3571     int iNewRows = 0;
3572     int iNewCols = 0;
3573     Bool* pSource = _pSource->getAs<Bool>();
3574
3575     //evaluate each argument and replace by appropriate value and compute the count of combinations
3576     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3577     if (iSeqCount == 0)
3578     {
3579         //free pArg content
3580         cleanIndexesArguments(_pArgs, &pArg);
3581         return this;
3582     }
3583
3584     if (iDims < 2)
3585     {
3586         //see as vector
3587         if (getRows() == 1 || getCols() == 1)
3588         {
3589             //vector or scalar
3590             bNeedToResize = true;
3591             if (getRows() * getCols() < piMaxDim[0])
3592             {
3593                 //need to enlarge sparse dimensions
3594                 if (getCols() == 1 || getRows() * getCols() == 0)
3595                 {
3596                     //column vector
3597                     iNewRows = piMaxDim[0];
3598                     iNewCols = 1;
3599                 }
3600                 else if (getRows() == 1)
3601                 {
3602                     //row vector
3603                     iNewRows = 1;
3604                     iNewCols = piMaxDim[0];
3605                 }
3606             }
3607         }
3608         else if (getRows() * getCols() < piMaxDim[0])
3609         {
3610             //free pArg content
3611             cleanIndexesArguments(_pArgs, &pArg);
3612             //out of range
3613             return NULL;
3614         }
3615     }
3616     else
3617     {
3618         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
3619         {
3620             bNeedToResize = true;
3621             iNewRows = std::max(getRows(), piMaxDim[0]);
3622             iNewCols = std::max(getCols(), piMaxDim[1]);
3623         }
3624     }
3625
3626     //check number of insertion
3627     if (pSource->isScalar() == false && pSource->getSize() != iSeqCount)
3628     {
3629         //free pArg content
3630         cleanIndexesArguments(_pArgs, &pArg);
3631         return NULL;
3632     }
3633
3634     //now you are sure to be able to insert values
3635     if (bNeedToResize)
3636     {
3637         if (resize(iNewRows, iNewCols) == NULL)
3638         {
3639             //free pArg content
3640             cleanIndexesArguments(_pArgs, &pArg);
3641             return NULL;
3642         }
3643     }
3644
3645     if (iDims == 1)
3646     {
3647         double* pIdx = pArg[0]->getAs<Double>()->get();
3648         for (int i = 0; i < iSeqCount; i++)
3649         {
3650             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
3651             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
3652             if (pSource->isScalar())
3653             {
3654                 set(iRow, iCol, pSource->get(0) != 0, false);
3655             }
3656             else
3657             {
3658                 set(iRow, iCol, pSource->get(i) != 0, false);
3659             }
3660         }
3661     }
3662     else
3663     {
3664         double* pIdxRow = pArg[0]->getAs<Double>()->get();
3665         int iRowSize = pArg[0]->getAs<Double>()->getSize();
3666         double* pIdxCol = pArg[1]->getAs<Double>()->get();
3667
3668         for (int i = 0; i < iSeqCount; i++)
3669         {
3670             if (pSource->isScalar())
3671             {
3672                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(0) != 0, false);
3673             }
3674             else
3675             {
3676                 int iRowOrig = i % pSource->getRows();
3677                 int iColOrig = i / pSource->getRows();
3678
3679                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(iRowOrig, iColOrig) != 0, false);
3680             }
3681         }
3682     }
3683
3684     finalize();
3685
3686     //free pArg content
3687     cleanIndexesArguments(_pArgs, &pArg);
3688     return this;
3689 }
3690
3691 GenericType* SparseBool::remove(typed_list* _pArgs)
3692 {
3693     SparseBool* pOut = NULL;
3694     int iDims = (int)_pArgs->size();
3695     if (iDims > 2)
3696     {
3697         //sparse are only in 2 dims
3698         return NULL;
3699     }
3700
3701     typed_list pArg;
3702
3703     int piMaxDim[2];
3704     int piCountDim[2];
3705
3706     //evaluate each argument and replace by appropriate value and compute the count of combinations
3707     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3708     if (iSeqCount == 0)
3709     {
3710         //free pArg content
3711         cleanIndexesArguments(_pArgs, &pArg);
3712         return this;
3713     }
3714
3715     bool* pbFull = new bool[iDims];
3716     //coord must represent all values on a dimension
3717     for (int i = 0; i < iDims; i++)
3718     {
3719         pbFull[i] = false;
3720         int iDimToCheck = getVarMaxDim(i, iDims);
3721         int iIndexSize = pArg[i]->getAs<GenericType>()->getSize();
3722
3723         //we can have index more than once
3724         if (iIndexSize >= iDimToCheck)
3725         {
3726             //size is good, now check datas
3727             double* pIndexes = getDoubleArrayFromDouble(pArg[i]);
3728             for (int j = 0; j < iDimToCheck; j++)
3729             {
3730                 bool bFind = false;
3731                 for (int k = 0; k < iIndexSize; k++)
3732                 {
3733                     if ((int)pIndexes[k] == j + 1)
3734                     {
3735                         bFind = true;
3736                         break;
3737                     }
3738                 }
3739                 pbFull[i] = bFind;
3740             }
3741         }
3742     }
3743
3744     //only one dims can be not full/entire
3745     bool bNotEntire = false;
3746     int iNotEntire = 0;
3747     bool bTooMuchNotEntire = false;
3748     for (int i = 0; i < iDims; i++)
3749     {
3750         if (pbFull[i] == false)
3751         {
3752             if (bNotEntire == false)
3753             {
3754                 bNotEntire = true;
3755                 iNotEntire = i;
3756             }
3757             else
3758             {
3759                 bTooMuchNotEntire = true;
3760                 break;
3761             }
3762         }
3763     }
3764
3765     delete[] pbFull;
3766
3767     if (bTooMuchNotEntire == true)
3768     {
3769         //free pArg content
3770         cleanIndexesArguments(_pArgs, &pArg);
3771         return NULL;
3772     }
3773
3774     //find index to keep
3775     int iNotEntireSize = pArg[iNotEntire]->getAs<GenericType>()->getSize();
3776     double* piNotEntireIndex = getDoubleArrayFromDouble(pArg[iNotEntire]);
3777     int iKeepSize = getVarMaxDim(iNotEntire, iDims);
3778     bool* pbKeep = new bool[iKeepSize];
3779
3780     //fill pbKeep with true value
3781     for (int i = 0; i < iKeepSize; i++)
3782     {
3783         pbKeep[i] = true;
3784     }
3785
3786     for (int i = 0; i < iNotEntireSize; i++)
3787     {
3788         int idx = (int)piNotEntireIndex[i] - 1;
3789
3790         //don't care of value out of bounds
3791         if (idx < iKeepSize)
3792         {
3793             pbKeep[idx] = false;
3794         }
3795     }
3796
3797     int iNewDimSize = 0;
3798     for (int i = 0; i < iKeepSize; i++)
3799     {
3800         if (pbKeep[i] == true)
3801         {
3802             iNewDimSize++;
3803         }
3804     }
3805     delete[] pbKeep;
3806
3807     int* piNewDims = new int[iDims];
3808     for (int i = 0; i < iDims; i++)
3809     {
3810         if (i == iNotEntire)
3811         {
3812             piNewDims[i] = iNewDimSize;
3813         }
3814         else
3815         {
3816             piNewDims[i] = getVarMaxDim(i, iDims);
3817         }
3818     }
3819
3820     //remove last dimension if are == 1
3821     int iOrigDims = iDims;
3822     for (int i = (iDims - 1); i >= 2; i--)
3823     {
3824         if (piNewDims[i] == 1)
3825         {
3826             iDims--;
3827         }
3828         else
3829         {
3830             break;
3831         }
3832     }
3833
3834     if (iDims == 1)
3835     {
3836         if (iNewDimSize == 0)
3837         {
3838             //free pArg content
3839             cleanIndexesArguments(_pArgs, &pArg);
3840             delete[] piNewDims;
3841             return new SparseBool(0, 0);
3842         }
3843         else
3844         {
3845             //two cases, depends of original matrix/vector
3846             if ((*_pArgs)[0]->isColon() == false && m_iDims == 2 && m_piDims[0] == 1 && m_piDims[1] != 1)
3847             {
3848                 //special case for row vector
3849                 pOut = new SparseBool(1, iNewDimSize);
3850                 //in this case we have to care of 2nd dimension
3851                 //iNotEntire = 1;
3852             }
3853             else
3854             {
3855                 pOut = new SparseBool(iNewDimSize, 1);
3856             }
3857         }
3858     }
3859     else
3860     {
3861         pOut = new SparseBool(piNewDims[0], piNewDims[0]);
3862     }
3863
3864     delete[] piNewDims;
3865     //find a way to copy existing data to new variable ...
3866     int iNewPos = 0;
3867     int* piIndexes = new int[iOrigDims];
3868     int* piViewDims = new int[iOrigDims];
3869     for (int i = 0; i < iOrigDims; i++)
3870     {
3871         piViewDims[i] = getVarMaxDim(i, iOrigDims);
3872     }
3873
3874     for (int i = 0; i < getSize(); i++)
3875     {
3876         bool bByPass = false;
3877         getIndexesWithDims(i, piIndexes, piViewDims, iOrigDims);
3878
3879         //check if piIndexes use removed indexes
3880         for (int j = 0; j < iNotEntireSize; j++)
3881         {
3882             if ((piNotEntireIndex[j] - 1) == piIndexes[iNotEntire])
3883             {
3884                 //by pass this value
3885                 bByPass = true;
3886                 break;
3887             }
3888         }
3889
3890         if (bByPass == false)
3891         {
3892             //compute new index
3893             pOut->set(iNewPos, get(i));
3894             iNewPos++;
3895         }
3896     }
3897
3898     delete[] piIndexes;
3899     delete[] piViewDims;
3900
3901     //free pArg content
3902     cleanIndexesArguments(_pArgs, &pArg);
3903
3904     return pOut;
3905 }
3906
3907 SparseBool* SparseBool::append(int r, int c, SparseBool SPARSE_CONST* src)
3908 {
3909     SparseBool* pIT = checkRef(this, &SparseBool::append, r, c, src);
3910     if (pIT != this)
3911     {
3912         return pIT;
3913     }
3914
3915     doAppend(*src, r, c, *matrixBool);
3916     finalize();
3917     return this;
3918 }
3919
3920 GenericType* SparseBool::insertNew(typed_list* _pArgs)
3921 {
3922     typed_list pArg;
3923     SparseBool *pOut  = NULL;
3924
3925     int iDims = (int)_pArgs->size();
3926     int* piMaxDim = new int[iDims];
3927     int* piCountDim = new int[iDims];
3928     bool bUndefine = false;
3929
3930     //evaluate each argument and replace by appropriate value and compute the count of combinations
3931     int iSeqCount = checkIndexesArguments(NULL, _pArgs, &pArg, piMaxDim, piCountDim);
3932     if (iSeqCount == 0)
3933     {
3934         //free pArg content
3935         cleanIndexesArguments(_pArgs, &pArg);
3936         return createEmptyDouble();
3937     }
3938
3939     if (iSeqCount < 0)
3940     {
3941         iSeqCount = -iSeqCount;
3942         bUndefine = true;
3943     }
3944
3945     if (bUndefine)
3946     {
3947         //manage : and $ in creation by insertion
3948         int iSource = 0;
3949         int *piSourceDims = getDimsArray();
3950
3951         for (int i = 0; i < iDims; i++)
3952         {
3953             if (pArg[i] == NULL)
3954             {
3955                 //undefine value
3956                 if (isScalar())
3957                 {
3958                     piMaxDim[i] = 1;
3959                     piCountDim[i] = 1;
3960                 }
3961                 else
3962                 {
3963                     piMaxDim[i] = piSourceDims[iSource];
3964                     piCountDim[i] = piSourceDims[iSource];
3965                 }
3966                 iSource++;
3967                 //replace pArg value by the new one
3968                 pArg[i] = createDoubleVector(piMaxDim[i]);
3969             }
3970             //else
3971             //{
3972             //    piMaxDim[i] = piCountDim[i];
3973             //}
3974         }
3975     }
3976
3977     //remove last dimension at size 1
3978     //remove last dimension if are == 1
3979     for (int i = (iDims - 1); i >= 2; i--)
3980     {
3981         if (piMaxDim[i] == 1)
3982         {
3983             iDims--;
3984             pArg.pop_back();
3985         }
3986         else
3987         {
3988             break;
3989         }
3990     }
3991
3992     if (checkArgValidity(pArg) == false)
3993     {
3994         //free pArg content
3995         cleanIndexesArguments(_pArgs, &pArg);
3996         //contain bad index, like <= 0, ...
3997         return NULL;
3998     }
3999
4000     if (iDims == 1)
4001     {
4002         if (getCols() == 1)
4003         {
4004             pOut = new SparseBool(piCountDim[0], 1);
4005         }
4006         else
4007         {
4008             //rows == 1
4009             pOut = new SparseBool(1, piCountDim[0]);
4010         }
4011     }
4012     else
4013     {
4014         pOut = new SparseBool(piMaxDim[0], piMaxDim[1]);
4015     }
4016
4017     //insert values in new matrix
4018     SparseBool* pOut2 = pOut->insert(&pArg, this);
4019     if (pOut != pOut2)
4020     {
4021         delete pOut;
4022     }
4023
4024     //free pArg content
4025     cleanIndexesArguments(_pArgs, &pArg);
4026
4027     return pOut2;
4028 }
4029
4030 SparseBool* SparseBool::extract(int nbCoords, int SPARSE_CONST* coords, int SPARSE_CONST* maxCoords, int SPARSE_CONST* resSize, bool asVector) SPARSE_CONST
4031 {
4032     if ((asVector && maxCoords[0] > getSize()) ||
4033             (asVector == false && maxCoords[0] > getRows()) ||
4034             (asVector == false && maxCoords[1] > getCols()))
4035     {
4036         return 0;
4037     }
4038
4039     SparseBool * pSp(0);
4040     if (asVector)
4041     {
4042         pSp = (getRows() == 1) ? new SparseBool(1, resSize[0]) : new SparseBool(resSize[0], 1);
4043         mycopy_n(makeMatrixIterator<bool>(*this, Coords<true>(coords, getRows())), nbCoords
4044                  , makeMatrixIterator<bool>(*(pSp->matrixBool), RowWiseFullIterator(pSp->getRows(), pSp->getCols())));
4045     }
4046     else
4047     {
4048         pSp = new SparseBool(resSize[0], resSize[1]);
4049         mycopy_n(makeMatrixIterator<bool>(*this, Coords<false>(coords, getRows())), nbCoords
4050                  , makeMatrixIterator<bool>(*(pSp->matrixBool), RowWiseFullIterator(pSp->getRows(), pSp->getCols())));
4051
4052     }
4053     return pSp;
4054 }
4055
4056 /*
4057 * create a new SparseBool of dims according to resSize and fill it from currentSparseBool (along coords)
4058 */
4059 GenericType* SparseBool::extract(typed_list* _pArgs)
4060 {
4061     SparseBool* pOut = NULL;
4062     int iDims = (int)_pArgs->size();
4063     typed_list pArg;
4064
4065     int* piMaxDim = new int[iDims];
4066     int* piCountDim = new int[iDims];
4067
4068     //evaluate each argument and replace by appropriate value and compute the count of combinations
4069     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
4070     if (iSeqCount == 0)
4071     {
4072         //free pArg content
4073         cleanIndexesArguments(_pArgs, &pArg);
4074         if (_pArgs->size() == 0)
4075         {
4076             delete[] piMaxDim;
4077             delete[] piCountDim;
4078             //a()
4079             return this;
4080         }
4081         else
4082         {
4083             delete[] piMaxDim;
4084             delete[] piCountDim;
4085             //a([])
4086             return new types::SparseBool(0,0);
4087         }
4088     }
4089
4090     if (iDims < 2)
4091     {
4092         // Check that we stay inside the input size.
4093         if (piMaxDim[0] <= getSize())
4094         {
4095             int iNewRows = 0;
4096             int iNewCols = 0;
4097
4098             if (getRows() == 1 && getCols() != 1 && (*_pArgs)[0]->isColon() == false)
4099             {
4100                 //special case for row vector
4101                 iNewRows = 1;
4102                 iNewCols = piCountDim[0];
4103             }
4104             else
4105             {
4106                 iNewRows = piCountDim[0];
4107                 iNewCols = 1;
4108             }
4109
4110             pOut = new SparseBool(iNewRows, iNewCols);
4111             double* pIdx = pArg[0]->getAs<Double>()->get();
4112             // Write in output all elements extract from input.
4113             for (int i = 0; i < iSeqCount; i++)
4114             {
4115                 if (pIdx[i] < 1)
4116                 {
4117                     delete pOut;
4118                     pOut = NULL;
4119                     break;
4120                 }
4121                 int iRowRead = static_cast<int>(pIdx[i] - 1) % getRows();
4122                 int iColRead = static_cast<int>(pIdx[i] - 1) / getRows();
4123
4124                 int iRowWrite = static_cast<int>(i) % iNewRows;
4125                 int iColWrite = static_cast<int>(i) / iNewRows;
4126
4127                 bool bValue = get(iRowRead, iColRead);
4128                 if (bValue)
4129                 {
4130                     //only non zero values
4131                     pOut->set(iRowWrite, iColWrite, true, false);
4132                 }
4133             }
4134         }
4135         else
4136         {
4137             delete[] piMaxDim;
4138             delete[] piCountDim;
4139             //free pArg content
4140             cleanIndexesArguments(_pArgs, &pArg);
4141             return NULL;
4142         }
4143     }
4144     else
4145     {
4146         // Check that we stay inside the input size.
4147         if (piMaxDim[0] <= getRows() && piMaxDim[1] <= getCols())
4148         {
4149             double* pIdxRow = pArg[0]->getAs<Double>()->get();
4150             double* pIdxCol = pArg[1]->getAs<Double>()->get();
4151
4152             int iNewRows = pArg[0]->getAs<Double>()->getSize();
4153             int iNewCols = pArg[1]->getAs<Double>()->getSize();
4154
4155             pOut = new SparseBool(iNewRows, iNewCols);
4156
4157             int iPos = 0;
4158             // Write in output all elements extract from input.
4159             for (int iRow = 0; iRow < iNewRows; iRow++)
4160             {
4161                 for (int iCol = 0; iCol < iNewCols; iCol++)
4162                 {
4163                     if ((pIdxRow[iRow] < 1) || (pIdxCol[iCol] < 1))
4164                     {
4165                         delete pOut;
4166                         pOut = NULL;
4167                         delete[] piMaxDim;
4168                         delete[] piCountDim;
4169                         //free pArg content
4170                         cleanIndexesArguments(_pArgs, &pArg);
4171                         return NULL;
4172                     }
4173                     bool bValue = get((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
4174                     if (bValue)
4175                     {
4176                         //only non zero values
4177                         pOut->set(iRow, iCol, true, false);
4178                     }
4179                     iPos++;
4180                 }
4181             }
4182         }
4183         else
4184         {
4185             delete[] piMaxDim;
4186             delete[] piCountDim;
4187             //free pArg content
4188             cleanIndexesArguments(_pArgs, &pArg);
4189             return NULL;
4190         }
4191     }
4192
4193     if (pOut)
4194     {
4195         pOut->finalize();
4196     }
4197
4198     delete[] piMaxDim;
4199     delete[] piCountDim;
4200     //free pArg content
4201     cleanIndexesArguments(_pArgs, &pArg);
4202
4203     return pOut;
4204 }
4205
4206 bool SparseBool::invoke(typed_list & in, optional_list &/*opt*/, int /*_iRetCount*/, typed_list & out, const ast::Exp & e)
4207 {
4208     if (in.size() == 0)
4209     {
4210         out.push_back(this);
4211     }
4212     else
4213     {
4214         InternalType * _out = extract(&in);
4215         if (!_out)
4216         {
4217             std::wostringstream os;
4218             os << _W("Invalid index.\n");
4219             throw ast::InternalError(os.str(), 999, e.getLocation());
4220         }
4221         out.push_back(_out);
4222     }
4223
4224     return true;
4225 }
4226
4227 bool SparseBool::isInvokable() const
4228 {
4229     return true;
4230 }
4231
4232 bool SparseBool::hasInvokeOption() const
4233 {
4234     return false;
4235 }
4236
4237 int SparseBool::getInvokeNbIn()
4238 {
4239     return -1;
4240 }
4241
4242 int SparseBool::getInvokeNbOut()
4243 {
4244     return 1;
4245 }
4246
4247 std::size_t SparseBool::nbTrue() const
4248 {
4249     return  matrixBool->nonZeros();
4250 }
4251 std::size_t SparseBool::nbTrue(std::size_t r) const
4252 {
4253     int* piIndex = matrixBool->outerIndexPtr();
4254     return piIndex[r + 1] - piIndex[r];
4255 }
4256
4257
4258 void SparseBool::setTrue(bool finalize)
4259 {
4260     int rows = getRows();
4261     int cols = getCols();
4262
4263     std::vector<BoolTriplet_t> tripletList;
4264
4265     for (int i = 0; i < rows; ++i)
4266     {
4267         for (int j = 0; j < cols; ++j)
4268         {
4269             tripletList.emplace_back(i, j, true);
4270         }
4271     }
4272
4273     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<bool>());
4274
4275     if (finalize)
4276     {
4277         matrixBool->finalize();
4278     }
4279 }
4280
4281 void SparseBool::setFalse(bool finalize)
4282 {
4283     int rows = getRows();
4284     int cols = getCols();
4285
4286     std::vector<BoolTriplet_t> tripletList;
4287
4288     for (int i = 0; i < rows; ++i)
4289     {
4290         for (int j = 0; j < cols; ++j)
4291         {
4292             tripletList.emplace_back(i, j, false);
4293         }
4294     }
4295
4296     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<bool>());
4297
4298     if (finalize)
4299     {
4300         matrixBool->finalize();
4301     }
4302 }
4303
4304 int* SparseBool::getNbItemByRow(int* _piNbItemByRows)
4305 {
4306     int* piNbItemByRows = new int[getRows() + 1];
4307     mycopy_n(matrixBool->outerIndexPtr(), getRows() + 1, piNbItemByRows);
4308
4309     for (int i = 0; i < getRows(); i++)
4310     {
4311         _piNbItemByRows[i] = piNbItemByRows[i + 1] - piNbItemByRows[i];
4312     }
4313
4314     delete[] piNbItemByRows;
4315     return _piNbItemByRows;
4316 }
4317
4318 int* SparseBool::getColPos(int* _piColPos)
4319 {
4320     mycopy_n(matrixBool->innerIndexPtr(), nbTrue(), _piColPos);
4321     for (size_t i = 0; i < nbTrue(); i++)
4322     {
4323         _piColPos[i]++;
4324     }
4325
4326     return _piColPos;
4327 }
4328
4329 int* SparseBool::outputRowCol(int* out)const
4330 {
4331     return sparseTransform(*matrixBool, sparseTransform(*matrixBool, out, GetRow<BoolSparse_t>()), GetCol<BoolSparse_t>());
4332 }
4333
4334 int* SparseBool::getInnerPtr(int* count)
4335 {
4336     *count = static_cast<int>(matrixBool->innerSize());
4337     return matrixBool->innerIndexPtr();
4338 }
4339
4340 int* SparseBool::getOuterPtr(int* count)
4341 {
4342     *count = static_cast<int>(matrixBool->outerSize());
4343     return matrixBool->outerIndexPtr();
4344 }
4345
4346
4347 bool SparseBool::operator==(const InternalType& it) SPARSE_CONST
4348 {
4349     SparseBool* otherSparse = const_cast<SparseBool*>(dynamic_cast<SparseBool const*>(&it));/* types::GenericType is not const-correct :( */
4350     return (otherSparse
4351             && (otherSparse->getRows() == getRows())
4352             && (otherSparse->getCols() == getCols())
4353             && equal(*matrixBool, *otherSparse->matrixBool));
4354 }
4355
4356 bool SparseBool::operator!=(const InternalType& it) SPARSE_CONST
4357 {
4358     return !(*this == it);
4359 }
4360
4361 void SparseBool::finalize()
4362 {
4363     matrixBool->prune(&keepForSparse<bool>);
4364     matrixBool->finalize();
4365 }
4366
4367 bool SparseBool::get(int r, int c) SPARSE_CONST
4368 {
4369     return matrixBool->coeff(r, c);
4370 }
4371
4372 SparseBool* SparseBool::set(int _iRows, int _iCols, bool _bVal, bool _bFinalize) SPARSE_CONST
4373 {
4374     typedef SparseBool* (SparseBool::*set_t)(int, int, bool, bool);
4375     SparseBool* pIT = checkRef(this, (set_t)&SparseBool::set, _iRows, _iCols, _bVal, _bFinalize);
4376     if (pIT != this)
4377     {
4378         return pIT;
4379     }
4380
4381     if (matrixBool->isCompressed() && matrixBool->coeff(_iRows, _iCols) == false)
4382     {
4383         matrixBool->reserve(1);
4384     }
4385
4386     matrixBool->coeffRef(_iRows, _iCols) = _bVal;
4387
4388     if (_bFinalize)
4389     {
4390         finalize();
4391     }
4392
4393     return this;
4394 }
4395
4396 void SparseBool::fill(Bool& dest, int r, int c) SPARSE_CONST
4397 {
4398     mycopy_n(makeMatrixIterator<bool >(*matrixBool, RowWiseFullIterator(getRows(), getCols())), getSize()
4399              , makeMatrixIterator<bool >(dest, RowWiseFullIterator(dest.getRows(), dest.getCols(), r, c)));
4400 }
4401
4402 Sparse* SparseBool::newOnes() const
4403 {
4404     return new Sparse(new types::Sparse::RealSparse_t(matrixBool->cast<double>()), 0);
4405 }
4406
4407 SparseBool* SparseBool::newNotEqualTo(SparseBool const&o) const
4408 {
4409     return cwiseOp<std::not_equal_to>(*this, o);
4410 }
4411
4412 SparseBool* SparseBool::newEqualTo(SparseBool& o)
4413 {
4414     int rowL = getRows();
4415     int colL = getCols();
4416
4417     int rowR = o.getRows();
4418     int colR = o.getCols();
4419     int row = std::max(rowL, rowR);
4420     int col = std::max(colL, colR);
4421
4422     //create a boolean sparse matrix with dims of sparses
4423     types::SparseBool* ret = new types::SparseBool(row, col);
4424
4425     if (isScalar() && o.isScalar())
4426     {
4427         bool l = get(0, 0);
4428         bool r = o.get(0, 0);
4429         ret->set(0, 0, l == r, false);
4430     }
4431     else if (isScalar())
4432     {
4433         int nnzR = static_cast<int>(o.nbTrue());
4434         std::vector<int> rowcolR(nnzR * 2, 0);
4435         o.outputRowCol(rowcolR.data());
4436
4437         //compare all items of R with R[0]
4438         bool l = get(0, 0);
4439         for (int i = 0; i < nnzR; ++i)
4440         {
4441             bool r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
4442             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
4443         }
4444     }
4445     else if (o.isScalar())
4446     {
4447         int nnzL = static_cast<int>(nbTrue());
4448         std::vector<int> rowcolL(nnzL * 2, 0);
4449         outputRowCol(rowcolL.data());
4450
4451         bool r = get(0, 0);
4452         for (int i = 0; i < nnzL; ++i)
4453         {
4454             bool l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
4455             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
4456         }
4457     }
4458     else
4459     {
4460         int nnzR = static_cast<int>(o.nbTrue());
4461         std::vector<int> rowcolR(nnzR * 2, 0);
4462         o.outputRowCol(rowcolR.data());
4463         int nnzL = static_cast<int>(nbTrue());
4464         std::vector<int> rowcolL(nnzL * 2, 0);
4465         outputRowCol(rowcolL.data());
4466         //set all values to %t
4467         ret->setTrue(false);
4468         //set %f in each pL values
4469         for (int i = 0; i < nnzL; ++i)
4470         {
4471             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, false, false);
4472         }
4473         ret->finalize();
4474
4475         //set _pR[i] == _pL[i] for each _pR values
4476         for (int i = 0; i < nnzR; ++i)
4477         {
4478             //get l and r following non zeros value of R
4479             bool l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
4480             bool r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
4481             //set value following non zeros value of R
4482             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
4483         }
4484     }
4485
4486     ret->finalize();
4487     return ret;
4488 }
4489
4490 SparseBool* SparseBool::newLogicalOr(SparseBool const&o) const
4491 {
4492     return cwiseOp<std::logical_or>(*this, o);
4493 }
4494
4495 SparseBool* SparseBool::newLogicalAnd(SparseBool const&o) const
4496 {
4497     return cwiseOp<std::logical_and>(*this, o);
4498 }
4499
4500 SparseBool* SparseBool::reshape(int* _piDims, int _iDims)
4501 {
4502     SparseBool* pSpBool = NULL;
4503     int iCols = 1;
4504
4505     if (_iDims == 2)
4506     {
4507         iCols = _piDims[1];
4508     }
4509
4510     if (_iDims <= 2)
4511     {
4512         pSpBool = reshape(_piDims[0], iCols);
4513     }
4514
4515     return pSpBool;
4516 }
4517
4518 SparseBool* SparseBool::reshape(int _iNewRows, int _iNewCols)
4519 {
4520     typedef SparseBool* (SparseBool::*reshape_t)(int, int);
4521     SparseBool* pIT = checkRef(this, (reshape_t)&SparseBool::reshape, _iNewRows, _iNewCols);
4522     if (pIT != this)
4523     {
4524         return pIT;
4525     }
4526
4527     if (_iNewRows * _iNewCols != getRows() * getCols())
4528     {
4529         return NULL;
4530     }
4531
4532     SparseBool* res = NULL;
4533     try
4534     {
4535         //item count
4536         size_t iNonZeros = matrixBool->nonZeros();
4537         BoolSparse_t *newBool = new BoolSparse_t(_iNewRows, _iNewCols);
4538         newBool->reserve((int)iNonZeros);
4539
4540         //coords
4541         int* pRows = new int[iNonZeros * 2];
4542         outputRowCol(pRows);
4543         int* pCols = pRows + iNonZeros;
4544
4545         std::vector<BoolTriplet_t> tripletList;
4546
4547         for (size_t i = 0; i < iNonZeros; i++)
4548         {
4549             int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
4550             tripletList.emplace_back((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), true);
4551         }
4552
4553         newBool->setFromTriplets(tripletList.begin(), tripletList.end(), DupFunctor<bool>());
4554
4555         delete matrixBool;
4556         matrixBool = newBool;
4557         delete[] pRows;
4558
4559         m_iRows = _iNewRows;
4560         m_iCols = _iNewCols;
4561         m_iSize = _iNewRows * _iNewCols;
4562
4563         m_iDims = 2;
4564         m_piDims[0] = m_iRows;
4565         m_piDims[1] = m_iCols;
4566
4567         finalize();
4568
4569         res = this;
4570     }
4571     catch (...)
4572     {
4573         res = NULL;
4574     }
4575     return res;
4576 }
4577
4578 bool SparseBool::transpose(InternalType *& out)
4579 {
4580     out = new SparseBool(new BoolSparse_t(matrixBool->transpose()));
4581     return true;
4582 }
4583
4584 template<typename T>
4585 void neg(const int r, const int c, const T * const in, Eigen::SparseMatrix<bool, 1> * const out)
4586 {
4587     for (int i = 0; i < r; i++)
4588     {
4589         for (int j = 0; j < c; j++)
4590         {
4591             out->coeffRef(i, j) = !in->coeff(i, j);
4592         }
4593     }
4594
4595     out->prune(&keepForSparse<bool>);
4596     out->finalize();
4597 }
4598 }