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

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