2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2010 - DIGITEO - Bernard HUGUENEY
5 * Copyright (C) 2012 - 2016 - Scilab Enterprises
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.
18 #include <Eigen/Sparse>
24 #include <Eigen/IterativeLinearSolvers>
25 #include <Eigen/SparseCholesky>
29 #include "tostring_common.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"
38 #include "types_tools.hxx"
40 #include "sparseOp.hxx"
44 #include "elem_common.h"
49 /* used for debuging output
51 template<typename Os, typename In, typename Sz> Os& writeData(wchar_t const* title, In beg, Sz n, Os& os)
54 /* TODO: use tostring_common (with a kind of std::boolalpha for boolean output)
56 mycopy_n(beg, n, std::ostream_iterator<typename std::iterator_traits<In>::value_type, char>(os, L" "));
63 Printer (int precision) : p(precision)
67 std::wstring emptyName( /* */) const
73 std::wstring operator()(T const& t) const
76 std::wostringstream ostr;
85 std::wstring Printer::operator()(bool const& b) const
98 std::wstring Printer::operator()(double const& d) const
100 std::wostringstream ostr;
102 getDoubleFormat(d, &df);
103 addDoubleValue(&ostr, d, &df);
108 std::wstring Printer::operator()(std::complex<double > const& c) const
110 std::wostringstream ostr;
112 DoubleFormat dfR, dfI;
113 getComplexFormat(c.real(), c.imag(), &iLen, &dfR, &dfI);
114 addDoubleComplexValue(&ostr, c.real(), c.imag(), iLen, &dfR, &dfI);
119 std::wstring Printer::emptyName<bool>() const
125 template<typename T> std::wstring toString(T const& m, int precision)
127 std::wostringstream ostr;
131 getSignedIntFormat(m.rows(), &iWidthRows);
132 getSignedIntFormat(m.cols(), &iWidthCols);
135 addUnsignedIntValue<unsigned long long>(&ostr, m.rows(), iWidthRows);
137 addUnsignedIntValue<unsigned long long>(&ostr, m.cols(), iWidthCols);
140 Printer p(precision);
143 ostr << ( p.emptyName<typename Eigen::internal::traits<T>::Scalar>());
145 ostr << L" sparse matrix\n\n";
147 auto * pIColPos = m.innerIndexPtr();
148 auto * pINbItemByRow = m.outerIndexPtr();
152 for (size_t j = 1 ; j < m.rows() + 1 ; j++)
154 for (size_t i = pINbItemByRow[j - 1] ; i < pINbItemByRow[j] ; i++)
157 addUnsignedIntValue<unsigned long long>(&ostr, (int)j, iWidthRows);
159 addUnsignedIntValue<unsigned long long>(&ostr, pIColPos[iPos] + 1, iWidthCols);
160 ostr << L")\t" << p(m.valuePtr()[iPos]) << std::endl;
169 /** utility function to compare two Eigen::Sparse matrices to equality
171 template<typename T> bool equal(T const& s1, T const& s2)
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);
178 for (int k = 0; res && k != s1.outerSize(); ++k)
180 for (typename T::InnerIterator it1(s1, k), it2(s2, k); res && it1 && it2 ; ++it1, ++it2, ++nbElts)
182 res = (it1.value() == it2.value()
183 && it1.row() == it2.row()
184 && it1.col() == it2.col());
187 return res && (nbElts == s1.nonZeros()) && (nbElts == s2.nonZeros());
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.)
194 template<typename T> bool setNonZero(T& s, typename Eigen::internal::traits<T>::Scalar v = 1.)
196 for (auto j = 0; j < s.outerSize(); ++j)
198 for (typename T::InnerIterator it(s, j); it; ++it)
208 template<typename Src, typename Sp>
209 void doAppend(Src SPARSE_CONST& src, int r, int c, Sp& dest)
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))));
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)
219 typedef typename Eigen::SparseMatrix<Scalar1, Eigen::RowMajor>::InnerIterator srcIt_t;
220 for (std::size_t k = 0; k != src.outerSize(); ++k)
222 for (srcIt_t it(src, (int)k); it; ++it)
224 dest.insert( it.row() + r, it.col() + c) = it.value();
229 Sp is an Eigen::SparseMatrix
231 template<typename Sp, typename M>
232 void cwiseInPlaceProduct(Sp& sp, M SPARSE_CONST& m)
234 // should be a transform_n() over makeNonZerosIterator(src)
235 for (std::size_t k = 0; k != sp.outerSize(); ++k)
237 for (typename Sp::InnerIterator it(sp, k); it; ++it)
239 it.valueRef() *= get<typename Eigen::internal::traits<Sp>::Scalar >(m, it.row(), it.col());
248 template<typename T, typename Arg>
249 T* create_new(Arg const& a)
255 Double* create_new(double const& d)
257 Double* res(new Double(1, 1, false));
263 Double* create_new(std::complex<double>const& c)
265 Double* res(new Double(1, 1, true));
266 res->set(0, 0, c.real());
267 res->setImg(0, 0, c.imag());
272 Double* create_new(Sparse const& s)
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);
286 Inspector::removeItem(this);
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)
295 m_iRows = const_cast<Sparse*>(&src)->getRows();
296 m_iCols = const_cast<Sparse*>(&src)->getCols();
297 m_iSize = m_iRows * m_iCols;
299 m_piDims[0] = m_iRows;
300 m_piDims[1] = m_iCols;
302 Inspector::addItem(this);
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)
312 m_iSize = _iRows * _iCols;
314 m_piDims[0] = _iRows;
315 m_piDims[1] = _iCols;
317 Inspector::addItem(this);
321 Sparse::Sparse(Double SPARSE_CONST& src)
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)
330 p[i] = (double)(i % row) + 1;
331 p[i + size] = (double)(i / row) + 1;
333 create2(src.getRows(), src.getCols(), src, *idx);
336 Inspector::addItem(this);
340 Sparse::Sparse(Double SPARSE_CONST& src, Double SPARSE_CONST& idx)
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));
346 create2(rows, cols, src, idx);
348 Inspector::removeItem(this);
352 Sparse::Sparse(Double SPARSE_CONST& src, Double SPARSE_CONST& idx, Double SPARSE_CONST& dims)
354 create2(static_cast<int>(dims.get(0)), static_cast<int>(dims.get(1)), src, idx);
356 Inspector::addItem(this);
360 Sparse::Sparse(RealSparse_t* realSp, CplxSparse_t* cplxSp): matrixReal(realSp), matrixCplx(cplxSp)
364 m_iCols = realSp->cols();
365 m_iRows = realSp->rows();
369 m_iCols = cplxSp->cols();
370 m_iRows = cplxSp->rows();
372 m_iSize = m_iCols * m_iRows;
374 m_piDims[0] = m_iRows;
375 m_piDims[1] = m_iCols;
379 Inspector::addItem(this);
383 Sparse::Sparse(Double SPARSE_CONST& xadj, Double SPARSE_CONST& adjncy, Double SPARSE_CONST& src, std::size_t r, std::size_t c)
385 Adjacency a(xadj.get(), adjncy.get());
386 create(static_cast<int>(r), static_cast<int>(c), src, makeIteratorFromVar(a), src.getSize());
388 Inspector::addItem(this);
392 Sparse::Sparse(int rows, int cols, int nonzeros, int* inner, int* outer, double* real, double* img)
399 matrixCplx = new CplxSparse_t(rows, cols);
400 matrixCplx->reserve((int)nonzeros);
401 out = matrixCplx->outerIndexPtr();
402 in = matrixCplx->innerIndexPtr();
403 matrixReal = nullptr;
407 matrixReal = new RealSparse_t(rows, cols);
408 matrixReal->reserve((int)nonzeros);
409 out = matrixReal->outerIndexPtr();
410 in = matrixReal->innerIndexPtr();
411 matrixCplx = nullptr;
414 //update outerIndexPtr
415 memcpy(out, outer, sizeof(int) * (rows + 1));
416 //update innerIndexPtr
417 memcpy(in, inner, sizeof(int) * nonzeros);
421 std::complex<double>* data = matrixCplx->valuePtr();
422 for (int i = 0; i < nonzeros; ++i)
424 data[i] = std::complex<double>(real[i], img[i]);
429 double* data = matrixReal->valuePtr();
430 for (int i = 0; i < nonzeros; ++i)
439 m_iSize = cols * rows;
441 m_piDims[0] = m_iRows;
442 m_piDims[1] = m_iCols;
444 matrixCplx ? matrixCplx->resizeNonZeros(nonzeros) : matrixReal->resizeNonZeros(nonzeros);
448 template<typename DestIter>
449 void Sparse::create(int rows, int cols, Double SPARSE_CONST& src, DestIter o, std::size_t n)
453 m_iSize = cols * rows;
455 m_piDims[0] = m_iRows;
456 m_piDims[1] = m_iCols;
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));
467 matrixReal = new RealSparse_t(rows, cols);
468 matrixReal->reserve((int)n);
470 mycopy_n(makeMatrixIterator<double >(src, RowWiseFullIterator(src.getRows(), src.getCols())), n
471 , makeMatrixIterator<double>(*matrixReal, o));
476 void Sparse::create2(int rows, int cols, Double SPARSE_CONST& src, Double SPARSE_CONST& idx)
478 int nnz = src.getSize();
479 double* i = idx.get();
480 double* j = i + idx.getRows();
481 double* valR = src.get();
487 typedef Eigen::Triplet<std::complex<double> > T;
488 std::vector<T> tripletList;
489 tripletList.reserve((int)nnz);
491 double* valI = src.getImg();
493 for (int k = 0; k < nnz; ++k)
495 tripletList.push_back(T(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, std::complex<double>(valR[k], valI[k])));
498 matrixCplx = new CplxSparse_t(rows, cols);
499 matrixCplx->setFromTriplets(tripletList.begin(), tripletList.end());
500 m_iRows = matrixCplx->rows();
501 m_iCols = matrixCplx->cols();
507 typedef Eigen::Triplet<double> T;
508 std::vector<T> tripletList;
509 tripletList.reserve((int)nnz);
511 for (int k = 0; k < nnz; ++k)
513 tripletList.push_back(T(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, valR[k]));
516 matrixReal = new RealSparse_t(rows, cols);
517 matrixReal->setFromTriplets(tripletList.begin(), tripletList.end());
519 m_iRows = matrixReal->rows();
520 m_iCols = matrixReal->cols();
523 m_iSize = m_iCols * m_iRows;
525 m_piDims[0] = m_iRows;
526 m_piDims[1] = m_iCols;
530 void Sparse::fill(Double& dest, int r, int c) SPARSE_CONST
532 Sparse & cthis(const_cast<Sparse&>(*this));
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)));
540 mycopy_n( makeMatrixIterator<double>(*matrixReal, RowWiseFullIterator(cthis.getRows(), cthis.getCols())), cthis.getSize()
541 , makeMatrixIterator<double >(dest, RowWiseFullIterator(dest.getRows(), dest.getCols(), r, c)));
545 Sparse* Sparse::set(int _iRows, int _iCols, std::complex<double> v, bool _bFinalize)
547 if (_iRows >= getRows() || _iCols >= getCols())
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);
561 matrixReal->coeffRef(_iRows, _iCols) = v.real();
565 matrixCplx->coeffRef(_iRows, _iCols) = v;
575 Sparse* Sparse::set(int _iRows, int _iCols, double _dblReal, bool _bFinalize)
577 if (_iRows >= getRows() || _iCols >= getCols())
582 typedef Sparse* (Sparse::*set_t)(int, int, double, bool);
583 Sparse* pIT = checkRef(this, (set_t)&Sparse::set, _iRows, _iCols, _dblReal, _bFinalize);
591 matrixReal->coeffRef(_iRows, _iCols) = _dblReal;
595 matrixCplx->coeffRef(_iRows, _iCols) = std::complex<double>(_dblReal, 0);
607 void Sparse::finalize()
611 matrixCplx->prune(&keepForSparse<std::complex<double> >);
612 matrixCplx->finalize();
616 matrixReal->prune(&keepForSparse<double>);
617 matrixReal->finalize();
622 bool Sparse::neg(InternalType *& out)
624 SparseBool * _out = new SparseBool(getRows(), getCols());
625 types::neg(getRows(), getCols(), matrixReal, _out->matrixBool);
632 bool Sparse::isComplex() const
634 return static_cast<bool>(matrixCplx != NULL);
637 // TODO: should have both a bounds checking and a non-checking interface to elt access
638 double* Sparse::get()
640 if (isComplex() == false)
642 return matrixReal->valuePtr();
648 double Sparse::get(int _iRows, int _iCols) const
650 return getReal(_iRows, _iCols);
653 double Sparse::getReal(int _iRows, int _iCols) const
658 res = matrixReal->coeff(_iRows, _iCols);
662 res = matrixCplx->coeff(_iRows, _iCols).real();
667 std::complex<double>* Sparse::getImg()
671 return matrixCplx->valuePtr();
677 std::complex<double> Sparse::getImg(int _iRows, int _iCols) const
679 std::complex<double> res;
682 res = matrixCplx->coeff(_iRows, _iCols);
686 res = std::complex<double>(matrixReal->coeff(_iRows, _iCols), 0.);
692 void Sparse::whoAmI() SPARSE_CONST
694 std::cout << "types::Sparse";
697 Sparse* Sparse::clone(void)
699 return new Sparse(*this);
702 bool Sparse::zero_set()
706 matrixReal->setZero();
710 matrixCplx->setZero();
716 // TODO: handle precision and line length
717 bool Sparse::toString(std::wostringstream& ostr)
719 int iPrecision = ConfigVariable::getFormatSize();
723 res = ::toString(*matrixReal, iPrecision);
727 res = ::toString(*matrixCplx, iPrecision);
734 Sparse* Sparse::resize(int _iNewRows, int _iNewCols)
736 typedef Sparse* (Sparse::*resize_t)(int, int);
737 Sparse* pIT = checkRef(this, (resize_t)&Sparse::resize, _iNewRows, _iNewCols);
743 if (_iNewRows <= getRows() && _iNewCols <= getCols())
745 //nothing to do: hence we do NOT fail
755 size_t iNonZeros = nonZeros();
756 RealSparse_t *newReal = new RealSparse_t(_iNewRows, _iNewCols);
757 newReal->reserve((int)iNonZeros);
761 int* pRows = new int[iNonZeros * 2];
763 int* pCols = pRows + iNonZeros;
766 double* pNonZeroR = new double[iNonZeros];
767 double* pNonZeroI = new double[iNonZeros];
768 outputValues(pNonZeroR, pNonZeroI);
770 typedef Eigen::Triplet<double> triplet;
771 std::vector<triplet> tripletList;
773 for (size_t i = 0 ; i < iNonZeros ; i++)
775 tripletList.push_back(triplet((int)pRows[i] - 1, (int)pCols[i] - 1, pNonZeroR[i]));
778 newReal->setFromTriplets(tripletList.begin(), tripletList.end());
781 matrixReal = newReal;
789 size_t iNonZeros = nonZeros();
790 CplxSparse_t *newCplx = new CplxSparse_t(_iNewRows, _iNewCols);
791 newCplx->reserve((int)iNonZeros);
794 int* pRows = new int[iNonZeros * 2];
796 int* pCols = pRows + iNonZeros;
799 double* pNonZeroR = new double[iNonZeros];
800 double* pNonZeroI = new double[iNonZeros];
801 outputValues(pNonZeroR, pNonZeroI);
803 typedef Eigen::Triplet<std::complex<double> > triplet;
804 std::vector<triplet> tripletList;
806 for (size_t i = 0 ; i < iNonZeros ; i++)
808 tripletList.push_back(triplet((int)pRows[i] - 1, (int)pCols[i] - 1, std::complex<double>(pNonZeroR[i], pNonZeroI[i])));
811 newCplx->setFromTriplets(tripletList.begin(), tripletList.end());
815 matrixCplx = newCplx;
823 m_iSize = _iNewRows * _iNewCols;
824 m_piDims[0] = m_iRows;
825 m_piDims[1] = m_iCols;
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
839 Sparse* otherSparse = const_cast<Sparse*>(dynamic_cast<Sparse const*>(&it));/* types::GenericType is not const-correct :( */
840 Sparse & cthis (const_cast<Sparse&>(*this));
842 if (otherSparse == NULL)
847 if (otherSparse->getRows() != cthis.getRows())
852 if (otherSparse->getCols() != cthis.getCols())
857 if (otherSparse->isComplex() != isComplex())
864 return equal(*matrixCplx, *otherSparse->matrixCplx);
868 return equal(*matrixReal, *otherSparse->matrixReal);
872 bool Sparse::one_set()
876 return setNonZero(*matrixCplx);
880 return setNonZero(*matrixReal);
884 void Sparse::toComplex()
890 matrixCplx = new CplxSparse_t(matrixReal->cast<std::complex<double> >());
903 GenericType* Sparse::insertNew(typed_list* _pArgs)
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;
914 //evaluate each argument and replace by appropriate value and compute the count of combinations
915 int iSeqCount = checkIndexesArguments(NULL, _pArgs, &pArg, piMaxDim, piCountDim);
919 cleanIndexesArguments(_pArgs, &pArg);
920 return createEmptyDouble();
925 iSeqCount = -iSeqCount;
931 //manage : and $ in creation by insertion
933 int *piSourceDims = getDimsArray();
935 for (int i = 0 ; i < iDims ; i++)
947 piMaxDim[i] = piSourceDims[iSource];
948 piCountDim[i] = piSourceDims[iSource];
951 //replace pArg value by the new one
952 pArg[i] = createDoubleVector(piMaxDim[i]);
956 // piMaxDim[i] = piCountDim[i];
961 //remove last dimension at size 1
962 //remove last dimension if are == 1
963 for (int i = (iDims - 1) ; i >= 2 ; i--)
965 if (piMaxDim[i] == 1)
976 if (checkArgValidity(pArg) == false)
979 cleanIndexesArguments(_pArgs, &pArg);
980 //contain bad index, like <= 0, ...
988 pOut = new Sparse(piCountDim[0], 1, bComplex);
993 pOut = new Sparse(1, piCountDim[0], bComplex);
998 pOut = new Sparse(piMaxDim[0], piMaxDim[1], bComplex);
999 //pOut = createEmpty(iDims, piMaxDim, bComplex);
1002 //insert values in new matrix
1003 Sparse* pOut2 = pOut->insert(&pArg, this);
1010 cleanIndexesArguments(_pArgs, &pArg);
1015 Sparse* Sparse::insert(typed_list* _pArgs, InternalType* _pSource)
1017 typedef Sparse* (Sparse::*insert_t)(typed_list*, InternalType*);
1018 Sparse* pIT = checkRef(this, (insert_t)&Sparse::insert, _pArgs, _pSource);
1024 if (_pSource->isSparse())
1026 return insert(_pArgs, _pSource->getAs<Sparse>());
1029 bool bNeedToResize = false;
1030 int iDims = (int)_pArgs->size();
1033 //sparse are only in 2 dims
1045 Double* pSource = _pSource->getAs<Double>();
1047 //evaluate each argument and replace by appropriate value and compute the count of combinations
1048 int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1052 cleanIndexesArguments(_pArgs, &pArg);
1059 if (getRows() == 1 || getCols() == 1)
1062 if (getSize() < piMaxDim[0])
1064 bNeedToResize = true;
1066 //need to enlarge sparse dimensions
1067 if (getCols() == 1 || getSize() == 0)
1070 iNewRows = piMaxDim[0];
1073 else if (getRows() == 1)
1077 iNewCols = piMaxDim[0];
1081 else if (getSize() < piMaxDim[0])
1084 cleanIndexesArguments(_pArgs, &pArg);
1091 if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
1093 bNeedToResize = true;
1094 iNewRows = std::max(getRows(), piMaxDim[0]);
1095 iNewCols = std::max(getCols(), piMaxDim[1]);
1099 //check number of insertion
1100 if (pSource->isScalar() == false && pSource->getSize() != iSeqCount)
1103 cleanIndexesArguments(_pArgs, &pArg);
1107 //now you are sure to be able to insert values
1110 if (resize(iNewRows, iNewCols) == false)
1113 cleanIndexesArguments(_pArgs, &pArg);
1119 if (pSource->isComplex() && isComplex() == false)
1127 double* pIdx = pArg[0]->getAs<Double>()->get();
1128 for (int i = 0 ; i < iSeqCount ; i++)
1130 int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
1131 int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
1132 if (pSource->isScalar())
1134 if (pSource->isComplex())
1136 set(iRow, iCol, std::complex<double>(pSource->get(0), pSource->getImg(0)), false);
1140 set(iRow, iCol, pSource->get(0), false);
1145 if (pSource->isComplex())
1147 set(iRow, iCol, std::complex<double>(pSource->get(i), pSource->getImg(i)), false);
1151 set(iRow, iCol, pSource->get(i), false);
1158 double* pIdxRow = pArg[0]->getAs<Double>()->get();
1159 int iRowSize = pArg[0]->getAs<Double>()->getSize();
1160 double* pIdxCol = pArg[1]->getAs<Double>()->get();
1162 for (int i = 0 ; i < iSeqCount ; i++)
1164 if (pSource->isScalar())
1166 if (pSource->isComplex())
1168 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, std::complex<double>(pSource->get(0), pSource->getImg(0)), false);
1172 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(0), false);
1177 int iRowOrig = i % pSource->getRows();
1178 int iColOrig = i / pSource->getRows();
1180 if (pSource->isComplex())
1182 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, std::complex<double>(pSource->get(iRowOrig, iColOrig), pSource->getImg(iRowOrig, iColOrig)), false);
1186 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(iRowOrig, iColOrig), false);
1195 cleanIndexesArguments(_pArgs, &pArg);
1200 Sparse* Sparse::insert(typed_list* _pArgs, Sparse* _pSource)
1202 bool bNeedToResize = false;
1203 int iDims = (int)_pArgs->size();
1206 //sparse are only in 2 dims
1219 //evaluate each argument and replace by appropriate value and compute the count of combinations
1220 int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1224 cleanIndexesArguments(_pArgs, &pArg);
1231 if (getRows() == 1 || getCols() == 1)
1234 bNeedToResize = true;
1235 if (getSize() < piMaxDim[0])
1237 //need to enlarge sparse dimensions
1238 if (getCols() == 1 || getSize() == 0)
1241 iNewRows = piMaxDim[0];
1244 else if (getRows() == 1)
1248 iNewCols = piMaxDim[0];
1252 else if (getSize() < piMaxDim[0])
1255 cleanIndexesArguments(_pArgs, &pArg);
1262 if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
1264 bNeedToResize = true;
1265 iNewRows = std::max(getRows(), piMaxDim[0]);
1266 iNewCols = std::max(getCols(), piMaxDim[1]);
1270 //check number of insertion
1271 if (_pSource->isScalar() == false && _pSource->getSize() != iSeqCount)
1274 cleanIndexesArguments(_pArgs, &pArg);
1278 //now you are sure to be able to insert values
1281 if (resize(iNewRows, iNewCols) == false)
1284 cleanIndexesArguments(_pArgs, &pArg);
1290 if (_pSource->isComplex() && isComplex() == false)
1297 double* pIdx = pArg[0]->getAs<Double>()->get();
1298 for (int i = 0 ; i < iSeqCount ; i++)
1300 int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
1301 int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
1303 if (_pSource->isScalar())
1305 if (_pSource->isComplex())
1307 set(iRow, iCol, _pSource->getImg(0, 0), false);
1311 set(iRow, iCol, _pSource->get(0, 0), false);
1316 int iRowOrig = i % _pSource->getRows();
1317 int iColOrig = i / _pSource->getRows();
1318 if (_pSource->isComplex())
1320 set(iRow, iCol, _pSource->getImg(iRowOrig, iColOrig), false);
1324 set(iRow, iCol, _pSource->get(iRowOrig, iColOrig), false);
1331 double* pIdxRow = pArg[0]->getAs<Double>()->get();
1332 int iRowSize = pArg[0]->getAs<Double>()->getSize();
1333 double* pIdxCol = pArg[1]->getAs<Double>()->get();
1335 for (int i = 0 ; i < iSeqCount ; i++)
1337 if (_pSource->isScalar())
1339 if (_pSource->isComplex())
1341 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->getImg(0, 0), false);
1345 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(0, 0), false);
1350 int iRowOrig = i % _pSource->getRows();
1351 int iColOrig = i / _pSource->getRows();
1352 if (_pSource->isComplex())
1354 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->getImg(iRowOrig, iColOrig), false);
1358 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(iRowOrig, iColOrig), false);
1367 cleanIndexesArguments(_pArgs, &pArg);
1372 GenericType* Sparse::remove(typed_list* _pArgs)
1374 Sparse* pOut = NULL;
1375 int iDims = (int)_pArgs->size();
1378 //sparse are only in 2 dims
1387 //evaluate each argument and replace by appropriate value and compute the count of combinations
1388 int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1392 cleanIndexesArguments(_pArgs, &pArg);
1396 bool* pbFull = new bool[iDims];
1397 //coord must represent all values on a dimension
1398 for (int i = 0 ; i < iDims ; i++)
1401 int iDimToCheck = getVarMaxDim(i, iDims);
1402 int iIndexSize = pArg[i]->getAs<GenericType>()->getSize();
1404 //we can have index more than once
1405 if (iIndexSize >= iDimToCheck)
1407 //size is good, now check datas
1408 double* pIndexes = getDoubleArrayFromDouble(pArg[i]);
1409 for (int j = 0 ; j < iDimToCheck ; j++)
1412 for (int k = 0 ; k < iIndexSize ; k++)
1414 if ((int)pIndexes[k] == j + 1)
1425 //only one dims can be not full/entire
1426 bool bNotEntire = false;
1428 bool bTooMuchNotEntire = false;
1429 for (int i = 0 ; i < iDims ; i++)
1431 if (pbFull[i] == false)
1433 if (bNotEntire == false)
1440 bTooMuchNotEntire = true;
1448 if (bTooMuchNotEntire == true)
1451 cleanIndexesArguments(_pArgs, &pArg);
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];
1461 //fill pbKeep with true value
1462 for (int i = 0 ; i < iKeepSize ; i++)
1467 for (int i = 0 ; i < iNotEntireSize ; i++)
1469 int idx = (int)piNotEntireIndex[i] - 1;
1471 //don't care of value out of bounds
1472 if (idx < iKeepSize)
1474 pbKeep[idx] = false;
1478 int iNewDimSize = 0;
1479 for (int i = 0 ; i < iKeepSize ; i++)
1481 if (pbKeep[i] == true)
1488 int* piNewDims = new int[iDims];
1489 for (int i = 0 ; i < iDims ; i++)
1491 if (i == iNotEntire)
1493 piNewDims[i] = iNewDimSize;
1497 piNewDims[i] = getVarMaxDim(i, iDims);
1501 //remove last dimension if are == 1
1502 int iOrigDims = iDims;
1503 for (int i = (iDims - 1) ; i >= 2 ; i--)
1505 if (piNewDims[i] == 1)
1517 if (iNewDimSize == 0)
1521 cleanIndexesArguments(_pArgs, &pArg);
1522 return new Sparse(0, 0);
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)
1529 //special case for row vector
1530 pOut = new Sparse(1, iNewDimSize, isComplex());
1531 //in this case we have to care of 2nd dimension
1536 pOut = new Sparse(iNewDimSize, 1, isComplex());
1542 pOut = new Sparse(piNewDims[0], piNewDims[0], isComplex());
1546 //find a way to copy existing data to new variable ...
1548 int* piIndexes = new int[iOrigDims];
1549 int* piViewDims = new int[iOrigDims];
1550 for (int i = 0 ; i < iOrigDims ; i++)
1552 piViewDims[i] = getVarMaxDim(i, iOrigDims);
1555 for (int i = 0 ; i < getSize() ; i++)
1557 bool bByPass = false;
1558 getIndexesWithDims(i, piIndexes, piViewDims, iOrigDims);
1560 //check if piIndexes use removed indexes
1561 for (int j = 0 ; j < iNotEntireSize ; j++)
1563 if ((piNotEntireIndex[j] - 1) == piIndexes[iNotEntire])
1565 //by pass this value
1571 if (bByPass == false)
1576 pOut->set(iNewPos, getImg(i));
1580 pOut->set(iNewPos, get(i));
1586 //free allocated data
1587 for (int i = 0 ; i < iDims ; i++)
1589 if (pArg[i] != (*_pArgs)[i])
1596 delete[] piViewDims;
1599 cleanIndexesArguments(_pArgs, &pArg);
1604 Sparse* Sparse::append(int r, int c, types::Sparse SPARSE_CONST* src)
1606 Sparse* pIT = checkRef(this, &Sparse::append, r, c, src);
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())
1619 if (src->isComplex())
1621 doAppend(*(src->matrixCplx), r, c, *matrixCplx);
1625 doAppend(*(src->matrixReal), r, c, *matrixCplx);
1630 doAppend(*(src->matrixReal), r, c, *matrixReal);
1635 return this; // realloc is meaningless for sparse matrices
1639 * create a new Sparse of dims according to resSize and fill it from currentSparse (along coords)
1641 GenericType* Sparse::extract(typed_list* _pArgs)
1643 Sparse* pOut = NULL;
1644 int iDims = (int)_pArgs->size();
1647 int* piMaxDim = new int[iDims];
1648 int* piCountDim = new int[iDims];
1650 //evaluate each argument and replace by appropriate value and compute the count of combinations
1651 int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1655 cleanIndexesArguments(_pArgs, &pArg);
1656 if (_pArgs->size() == 0)
1660 delete[] piCountDim;
1662 cleanIndexesArguments(_pArgs, &pArg);
1669 delete[] piCountDim;
1671 cleanIndexesArguments(_pArgs, &pArg);
1672 return Double::Empty();
1678 if (piMaxDim[0] <= getSize())
1683 if (getRows() == 1 && getCols() != 1 && (*_pArgs)[0]->isColon() == false)
1685 //special case for row vector
1687 iNewCols = piCountDim[0];
1691 iNewRows = piCountDim[0];
1695 pOut = new Sparse(iNewRows, iNewCols, isComplex());
1696 double* pIdx = pArg[0]->getAs<Double>()->get();
1697 for (int i = 0 ; i < iSeqCount ; i++)
1704 delete[] piCountDim;
1705 cleanIndexesArguments(_pArgs, &pArg);
1708 int iRowRead = static_cast<int>(pIdx[i] - 1) % getRows();
1709 int iColRead = static_cast<int>(pIdx[i] - 1) / getRows();
1711 int iRowWrite = static_cast<int>(i) % iNewRows;
1712 int iColWrite = static_cast<int>(i) / iNewRows;
1715 std::complex<double> dbl = getImg(iRowRead, iColRead);
1716 if (dbl.real() != 0 || dbl.imag() != 0)
1718 //only non zero values
1719 pOut->set(iRowWrite, iColWrite, dbl, false);
1724 double dbl = get(iRowRead, iColRead);
1727 //only non zero values
1728 pOut->set(iRowWrite, iColWrite, dbl, false);
1736 delete[] piCountDim;
1738 cleanIndexesArguments(_pArgs, &pArg);
1744 if (piMaxDim[0] <= getRows() && piMaxDim[1] <= getCols())
1746 double* pIdxRow = pArg[0]->getAs<Double>()->get();
1747 double* pIdxCol = pArg[1]->getAs<Double>()->get();
1749 int iNewRows = pArg[0]->getAs<Double>()->getSize();
1750 int iNewCols = pArg[1]->getAs<Double>()->getSize();
1752 pOut = new Sparse(iNewRows, iNewCols, isComplex());
1755 for (int iRow = 0 ; iRow < iNewRows ; iRow++)
1757 for (int iCol = 0 ; iCol < iNewCols ; iCol++)
1759 if ((pIdxRow[iRow] < 1) || (pIdxCol[iCol] < 1))
1764 delete[] piCountDim;
1766 cleanIndexesArguments(_pArgs, &pArg);
1771 std::complex<double> dbl = getImg((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
1772 if (dbl.real() != 0 || dbl.imag() != 0)
1774 //only non zero values
1775 pOut->set(iRow, iCol, dbl, false);
1780 double dbl = get((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
1783 //only non zero values
1784 pOut->set(iRow, iCol, dbl, false);
1794 delete[] piCountDim;
1796 cleanIndexesArguments(_pArgs, &pArg);
1804 delete[] piCountDim;
1806 cleanIndexesArguments(_pArgs, &pArg);
1811 Sparse* Sparse::extract(int nbCoords, int SPARSE_CONST* coords, int SPARSE_CONST* maxCoords, int SPARSE_CONST* resSize, bool asVector) SPARSE_CONST
1813 if ( (asVector && maxCoords[0] > getSize()) ||
1814 (asVector == false && maxCoords[0] > getRows()) ||
1815 (asVector == false && maxCoords[1] > getCols()))
1820 bool const cplx(isComplex());
1824 pSp = (getRows() == 1) ? new Sparse(1, resSize[0], cplx) : new Sparse(resSize[0], 1, cplx);
1828 pSp = new Sparse(resSize[0], resSize[1], cplx);
1830 // std::cerr<<"extracted sparse:"<<pSp->getRows()<<", "<<pSp->getCols()<<"seqCount="<<nbCoords<<"maxDim="<<maxCoords[0] <<","<< maxCoords[1]<<std::endl;
1832 ? copyToSparse(*this, Coords<true>(coords, getRows()), nbCoords
1833 , *pSp, RowWiseFullIterator(pSp->getRows(), pSp->getCols()))
1834 : copyToSparse(*this, Coords<false>(coords), nbCoords
1835 , *pSp, RowWiseFullIterator(pSp->getRows(), pSp->getCols()))))
1843 bool Sparse::invoke(typed_list & in, optional_list & /*opt*/, int /*_iRetCount*/, typed_list & out, const ast::Exp & e)
1847 out.push_back(this);
1851 InternalType * _out = extract(&in);
1854 std::wostringstream os;
1855 os << _W("Invalid index.\n");
1856 throw ast::InternalError(os.str(), 999, e.getLocation());
1858 out.push_back(_out);
1865 bool Sparse::isInvokable() const
1870 bool Sparse::hasInvokeOption() const
1875 int Sparse::getInvokeNbIn()
1880 int Sparse::getInvokeNbOut()
1886 coords are Scilab 1-based
1887 extract std::make_pair(coords, asVector), rowIter
1889 template<typename Src, typename SrcTraversal, typename Sz, typename DestTraversal>
1890 bool Sparse::copyToSparse(Src SPARSE_CONST& src, SrcTraversal srcTrav, Sz n, Sparse& sp, DestTraversal destTrav)
1892 if (!(src.isComplex() || sp.isComplex()))
1894 mycopy_n(makeMatrixIterator<double>(src, srcTrav), n
1895 , makeMatrixIterator<double>(*sp.matrixReal, destTrav));
1900 mycopy_n(makeMatrixIterator<std::complex<double> >(src, srcTrav), n
1901 , makeMatrixIterator<std::complex<double> >(*sp.matrixCplx, destTrav));
1908 // GenericType because we might return a Double* for scalar operand
1909 Sparse* Sparse::add(Sparse const& o) const
1911 RealSparse_t* realSp(0);
1912 CplxSparse_t* cplxSp(0);
1913 if (isComplex() == false && o.isComplex() == false)
1916 realSp = new RealSparse_t(*matrixReal + * (o.matrixReal));
1918 else if (isComplex() == false && o.isComplex() == true)
1920 cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >() + * (o.matrixCplx));
1922 else if (isComplex() == true && o.isComplex() == false)
1924 cplxSp = new CplxSparse_t(*matrixCplx + o.matrixReal->cast<std::complex<double> >());
1926 else if (isComplex() == true && o.isComplex() == true)
1928 cplxSp = new CplxSparse_t(*matrixCplx + * (o.matrixCplx));
1931 return new Sparse(realSp, cplxSp);
1934 Sparse* Sparse::substract(Sparse const& o) const
1936 RealSparse_t* realSp(0);
1937 CplxSparse_t* cplxSp(0);
1938 if (isComplex() == false && o.isComplex() == false)
1941 realSp = new RealSparse_t(*matrixReal - * (o.matrixReal));
1943 else if (isComplex() == false && o.isComplex() == true)
1946 cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >() - * (o.matrixCplx));
1948 else if (isComplex() == true && o.isComplex() == false)
1951 cplxSp = new CplxSparse_t(*matrixCplx - o.matrixReal->cast<std::complex<double> >());
1953 else if (isComplex() == true && o.isComplex() == true)
1956 cplxSp = new CplxSparse_t(*matrixCplx - * (o.matrixCplx));
1959 return new Sparse(realSp, cplxSp);
1962 Sparse* Sparse::multiply(double s) const
1964 return new Sparse( isComplex() ? 0 : new RealSparse_t((*matrixReal)*s)
1965 , isComplex() ? new CplxSparse_t((*matrixCplx)*s) : 0);
1968 Sparse* Sparse::multiply(std::complex<double> s) const
1970 return new Sparse( 0
1971 , isComplex() ? new CplxSparse_t((*matrixCplx) * s) : new CplxSparse_t((*matrixReal) * s));
1974 Sparse* Sparse::multiply(Sparse const& o) const
1976 RealSparse_t* realSp(0);
1977 CplxSparse_t* cplxSp(0);
1979 if (isComplex() == false && o.isComplex() == false)
1981 realSp = new RealSparse_t(*matrixReal **(o.matrixReal));
1983 else if (isComplex() == false && o.isComplex() == true)
1985 cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >() **(o.matrixCplx));
1987 else if (isComplex() == true && o.isComplex() == false)
1989 cplxSp = new CplxSparse_t(*matrixCplx * o.matrixReal->cast<std::complex<double> >());
1991 else if (isComplex() == true && o.isComplex() == true)
1993 cplxSp = new CplxSparse_t(*matrixCplx **(o.matrixCplx));
1996 return new Sparse(realSp, cplxSp);
1999 Sparse* Sparse::dotMultiply(Sparse SPARSE_CONST& o) const
2001 RealSparse_t* realSp(0);
2002 CplxSparse_t* cplxSp(0);
2003 if (isComplex() == false && o.isComplex() == false)
2005 realSp = new RealSparse_t(matrixReal->cwiseProduct(*(o.matrixReal)));
2007 else if (isComplex() == false && o.isComplex() == true)
2009 cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >().cwiseProduct( *(o.matrixCplx)));
2011 else if (isComplex() == true && o.isComplex() == false)
2013 cplxSp = new CplxSparse_t(matrixCplx->cwiseProduct(o.matrixReal->cast<std::complex<double> >()));
2015 else if (isComplex() == true && o.isComplex() == true)
2017 cplxSp = new CplxSparse_t(matrixCplx->cwiseProduct(*(o.matrixCplx)));
2020 return new Sparse(realSp, cplxSp);
2023 Sparse* Sparse::dotDivide(Sparse SPARSE_CONST& o) const
2025 RealSparse_t* realSp(0);
2026 CplxSparse_t* cplxSp(0);
2027 if (isComplex() == false && o.isComplex() == false)
2029 realSp = new RealSparse_t(matrixReal->cwiseQuotient(*(o.matrixReal)));
2031 else if (isComplex() == false && o.isComplex() == true)
2033 cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >().cwiseQuotient( *(o.matrixCplx)));
2035 else if (isComplex() == true && o.isComplex() == false)
2037 cplxSp = new CplxSparse_t(matrixCplx->cwiseQuotient(o.matrixReal->cast<std::complex<double> >()));
2039 else if (isComplex() == true && o.isComplex() == true)
2041 cplxSp = new CplxSparse_t(matrixCplx->cwiseQuotient(*(o.matrixCplx)));
2044 return new Sparse(realSp, cplxSp);
2047 int Sparse::newCholLLT(Sparse** _SpPermut, Sparse** _SpFactor) const
2049 typedef Eigen::SparseMatrix<double, Eigen::ColMajor> RealSparseCol_t;
2050 RealSparseCol_t spColMajor = RealSparseCol_t((const RealSparse_t&) * matrixReal);
2052 // Constructs and performs the LLT factorization of sparse
2053 Eigen::SimplicialLLT<RealSparseCol_t> pLLT(spColMajor);
2054 int iInfo = pLLT.info();
2055 if (iInfo != Eigen::Success)
2062 // Get the lower matrix of factorization.
2063 // The new RealSparse_t will be setted in Sparse without copy.
2064 *_SpFactor = new Sparse(new RealSparse_t(pLLT.matrixL()), NULL);
2066 // Get the permutation matrix.
2067 Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic, int> p = pLLT.permutationP();
2068 *_SpPermut = new Sparse(p.rows(), p.cols());
2069 for (int i = 0; i < p.rows(); i++)
2071 (*_SpPermut)->set(i, p.indices()[i], 1, false);
2074 (*_SpPermut)->finalize();
2079 bool Sparse::transpose(InternalType *& out)
2081 out = new Sparse(matrixReal ? new RealSparse_t(matrixReal->transpose()) : 0, matrixCplx ? new CplxSparse_t(matrixCplx->transpose()) : 0);
2085 bool Sparse::adjoint(InternalType *& out)
2087 out = new Sparse(matrixReal ? new RealSparse_t(matrixReal->adjoint()) : 0, matrixCplx ? new CplxSparse_t(matrixCplx->adjoint()) : 0);
2093 BoolCast(std::complex<double> const& c): b(c.real() || c.imag()) {}
2094 operator bool () const
2098 operator double() const
2104 Sparse* Sparse::newOnes() const
2106 // result is never cplx
2107 return new Sparse( matrixReal
2108 ? new RealSparse_t(matrixReal->cast<bool>().cast<double>())
2109 : new RealSparse_t(matrixCplx->cast<BoolCast>().cast<double>())
2115 RealCast(std::complex<double> const& c): b(c.real()) {}
2116 operator bool () const
2120 operator double() const
2126 Sparse* Sparse::newReal() const
2128 return new Sparse( matrixReal
2130 : new RealSparse_t(matrixCplx->cast<RealCast>().cast<double>())
2134 std::size_t Sparse::nonZeros() const
2138 return matrixCplx->nonZeros();
2142 return matrixReal->nonZeros();
2145 std::size_t Sparse::nonZeros(std::size_t r) const
2150 int* piIndex = matrixReal->outerIndexPtr();
2151 res = piIndex[r + 1] - piIndex[r];
2155 int* piIndex = matrixCplx->outerIndexPtr();
2156 res = piIndex[r + 1] - piIndex[r];
2162 int* Sparse::getNbItemByRow(int* _piNbItemByRows)
2164 int* piNbItemByCols = new int[getRows() + 1];
2167 mycopy_n(matrixCplx->outerIndexPtr(), getRows() + 1, piNbItemByCols);
2171 mycopy_n(matrixReal->outerIndexPtr(), getRows() + 1, piNbItemByCols);
2174 for (int i = 0 ; i < getRows() ; i++)
2176 _piNbItemByRows[i] = piNbItemByCols[i + 1] - piNbItemByCols[i];
2179 delete[] piNbItemByCols;
2180 return _piNbItemByRows;
2183 int* Sparse::getColPos(int* _piColPos)
2187 mycopy_n(matrixCplx->innerIndexPtr(), nonZeros(), _piColPos);
2191 mycopy_n(matrixReal->innerIndexPtr(), nonZeros(), _piColPos);
2194 for (int i = 0; i < nonZeros(); i++)
2202 int* Sparse::getInnerPtr(int* count)
2207 ret = matrixCplx->innerIndexPtr();
2208 *count = matrixCplx->innerSize();
2212 ret = matrixReal->innerIndexPtr();
2213 *count = matrixReal->innerSize();
2219 int* Sparse::getOuterPtr(int* count)
2224 ret = matrixCplx->outerIndexPtr();
2225 *count = matrixCplx->outerSize();
2229 ret = matrixReal->outerIndexPtr();
2230 *count = matrixReal->outerSize();
2238 template<typename S> struct GetReal: std::unary_function<typename S::InnerIterator, double>
2240 double operator()(typename S::InnerIterator it) const
2245 template<> struct GetReal< Eigen::SparseMatrix<std::complex<double >, Eigen::RowMajor > >
2246 : std::unary_function<Sparse::CplxSparse_t::InnerIterator, double>
2248 double operator()( Sparse::CplxSparse_t::InnerIterator it) const
2250 return it.value().real();
2253 template<typename S> struct GetImag: std::unary_function<typename S::InnerIterator, double>
2255 double operator()(typename S::InnerIterator it) const
2257 return it.value().imag();
2260 template<typename S> struct GetRow: std::unary_function<typename S::InnerIterator, int>
2262 int operator()(typename S::InnerIterator it) const
2264 return it.row() + 1;
2267 template<typename S> struct GetCol: std::unary_function<typename S::InnerIterator, int>
2269 int operator()(typename S::InnerIterator it) const
2271 return it.col() + 1;
2275 template<typename S, typename Out, typename F> Out sparseTransform(S& s, Out o, F f)
2277 for (std::size_t k(0); k < s.outerSize(); ++k)
2279 for (typename S::InnerIterator it(s, (int)k); it; ++it, ++o)
2288 std::pair<double*, double*> Sparse::outputValues(double* outReal, double* outImag)const
2291 ? std::make_pair(sparseTransform(*matrixReal, outReal, GetReal<RealSparse_t>()), outImag)
2292 : std::make_pair(sparseTransform(*matrixCplx, outReal, GetReal<CplxSparse_t>())
2293 , sparseTransform(*matrixCplx, outImag, GetImag<CplxSparse_t>()));
2296 int* Sparse::outputRowCol(int* out)const
2299 ? sparseTransform(*matrixReal, sparseTransform(*matrixReal, out, GetRow<RealSparse_t>()), GetCol<RealSparse_t>())
2300 : sparseTransform(*matrixCplx, sparseTransform(*matrixCplx, out, GetRow<CplxSparse_t>()), GetCol<CplxSparse_t>());
2302 double* Sparse::outputCols(double* out) const
2306 mycopy_n(matrixCplx->innerIndexPtr(), nonZeros(), out);
2310 mycopy_n(matrixReal->innerIndexPtr(), nonZeros(), out);
2313 return std::transform(out, out, out, std::bind2nd(std::plus<double>(), 1));
2317 void Sparse::opposite(void)
2321 std::complex<double>* data = matrixCplx->valuePtr();
2322 std::transform(data, data + matrixCplx->nonZeros(), data, std::negate<std::complex<double> >());
2326 double* data = matrixReal->valuePtr();
2327 std::transform(data, data + matrixReal->nonZeros(), data, std::negate<double>());
2331 SparseBool* Sparse::newLessThan(Sparse &o)
2333 //only real values !
2335 //return cwiseOp<std::less>(*this, o);
2336 int rowL = getRows();
2337 int colL = getCols();
2339 int rowR = o.getRows();
2340 int colR = o.getCols();
2341 int row = std::max(rowL, rowR);
2342 int col = std::max(colL, colR);
2344 //create a boolean sparse matrix with dims of sparses
2345 types::SparseBool* ret = new types::SparseBool(row, col);
2346 if (isScalar() && o.isScalar())
2348 double l = get(0, 0);
2349 double r = o.get(0, 0);
2350 ret->set(0, 0, l < r, false);
2352 else if (isScalar())
2354 int nnzR = static_cast<int>(o.nonZeros());
2355 std::vector<int> rowcolR(nnzR * 2, 0);
2356 o.outputRowCol(rowcolR.data());
2358 //compare all items of R with R[0]
2359 double l = get(0, 0);
2363 ret->setTrue(false);
2366 for (int i = 0; i < nnzR; ++i)
2368 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2369 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l < r, false);
2372 else if (o.isScalar())
2374 int nnzL = static_cast<int>(nonZeros());
2375 std::vector<int> rowcolL(nnzL * 2, 0);
2376 outputRowCol(rowcolL.data());
2378 double r = o.get(0, 0);
2384 for (int i = 0; i < nnzL; ++i)
2386 double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2387 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l < r, false);
2392 int nnzR = static_cast<int>(o.nonZeros());
2393 std::vector<int> rowcolR(nnzR * 2, 0);
2394 o.outputRowCol(rowcolR.data());
2395 int nnzL = static_cast<int>(nonZeros());
2396 std::vector<int> rowcolL(nnzL * 2, 0);
2397 outputRowCol(rowcolL.data());
2398 //set all values to %t
2399 ret->setFalse(false);
2401 for (int i = 0; i < nnzL; ++i)
2403 double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2404 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l < 0, false);
2408 //set _pR[i] == _pL[i] for each _pR values
2409 for (int i = 0; i < nnzR; ++i)
2411 //get l and r following non zeros value of R
2412 double l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2413 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2414 //set value following non zeros value of R
2415 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l < r, false);
2423 SparseBool* Sparse::newNotEqualTo(Sparse const&o) const
2425 return cwiseOp<std::not_equal_to>(*this, o);
2428 SparseBool* Sparse::newLessOrEqual(Sparse &o)
2430 //only real values !
2432 //return cwiseOp<std::less>(*this, o);
2433 int rowL = getRows();
2434 int colL = getCols();
2436 int rowR = o.getRows();
2437 int colR = o.getCols();
2438 int row = std::max(rowL, rowR);
2439 int col = std::max(colL, colR);
2441 //create a boolean sparse matrix with dims of sparses
2442 types::SparseBool* ret = new types::SparseBool(row, col);
2443 if (isScalar() && o.isScalar())
2445 double l = get(0, 0);
2446 double r = o.get(0, 0);
2447 ret->set(0, 0, l <= r, false);
2449 else if (isScalar())
2451 int nnzR = static_cast<int>(o.nonZeros());
2452 std::vector<int> rowcolR(nnzR * 2, 0);
2453 o.outputRowCol(rowcolR.data());
2455 //compare all items of R with R[0]
2456 double l = get(0, 0);
2460 ret->setTrue(false);
2463 for (int i = 0; i < nnzR; ++i)
2465 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2466 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l <= r, false);
2469 else if (o.isScalar())
2471 int nnzL = static_cast<int>(nonZeros());
2472 std::vector<int> rowcolL(nnzL * 2, 0);
2473 outputRowCol(rowcolL.data());
2475 double r = o.get(0, 0);
2481 for (int i = 0; i < nnzL; ++i)
2483 double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2484 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l <= r, false);
2489 int nnzR = static_cast<int>(o.nonZeros());
2490 std::vector<int> rowcolR(nnzR * 2, 0);
2491 o.outputRowCol(rowcolR.data());
2492 int nnzL = static_cast<int>(nonZeros());
2493 std::vector<int> rowcolL(nnzL * 2, 0);
2494 outputRowCol(rowcolL.data());
2495 //set all values to %t
2496 ret->setTrue(false);
2498 for (int i = 0; i < nnzL; ++i)
2500 double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2501 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l <= 0, false);
2505 //set _pR[i] == _pL[i] for each _pR values
2506 for (int i = 0; i < nnzR; ++i)
2508 //get l and r following non zeros value of R
2509 double l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2510 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2511 //set value following non zeros value of R
2512 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l <= r, false);
2520 SparseBool* Sparse::newEqualTo(Sparse &o)
2522 int rowL = getRows();
2523 int colL = getCols();
2525 int rowR = o.getRows();
2526 int colR = o.getCols();
2527 int row = std::max(rowL, rowR);
2528 int col = std::max(colL, colR);
2530 //create a boolean sparse matrix with dims of sparses
2531 types::SparseBool* ret = new types::SparseBool(row, col);
2532 if (isScalar() && o.isScalar())
2534 if (isComplex() || o.isComplex())
2536 std::complex<double> l = getImg(0, 0);
2537 std::complex<double> r = o.getImg(0, 0);
2538 ret->set(0, 0, l == r, false);
2542 double l = get(0, 0);
2543 double r = o.get(0, 0);
2544 ret->set(0, 0, l == r, false);
2547 else if (isScalar())
2549 int nnzR = static_cast<int>(o.nonZeros());
2550 std::vector<int> rowcolR(nnzR * 2, 0);
2551 o.outputRowCol(rowcolR.data());
2553 //compare all items of R with R[0]
2554 if (isComplex() || o.isComplex())
2556 std::complex<double> l = getImg(0, 0);
2557 for (int i = 0; i < nnzR; ++i)
2559 std::complex<double> r = o.getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2560 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2565 double l = get(0, 0);
2566 for (int i = 0; i < nnzR; ++i)
2568 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2569 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2573 else if (o.isScalar())
2575 int nnzL = static_cast<int>(nonZeros());
2576 std::vector<int> rowcolL(nnzL * 2, 0);
2577 outputRowCol(rowcolL.data());
2579 if (isComplex() || o.isComplex())
2581 std::complex<double> r = o.getImg(0, 0);
2582 for (int i = 0; i < nnzL; ++i)
2584 std::complex<double> l = getImg(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2585 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
2590 double r = get(0, 0);
2591 for (int i = 0; i < nnzL; ++i)
2593 double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2594 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
2600 int nnzR = static_cast<int>(o.nonZeros());
2601 std::vector<int> rowcolR(nnzR * 2, 0);
2602 o.outputRowCol(rowcolR.data());
2603 int nnzL = static_cast<int>(nonZeros());
2604 std::vector<int> rowcolL(nnzL * 2, 0);
2605 outputRowCol(rowcolL.data());
2606 //set all values to %t
2607 ret->setTrue(false);
2608 //set %f in each pL values
2609 for (int i = 0; i < nnzL; ++i)
2611 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, false, false);
2615 //set _pR[i] == _pL[i] for each _pR values
2616 if (isComplex() || o.isComplex())
2618 for (int i = 0; i < nnzR; ++i)
2620 //get l and r following non zeros value of R
2621 std::complex<double> l = getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2622 std::complex<double> r = o.getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2623 //set value following non zeros value of R
2624 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2629 for (int i = 0; i < nnzR; ++i)
2631 //get l and r following non zeros value of R
2632 double l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2633 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2634 //set value following non zeros value of R
2635 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2644 Sparse* Sparse::reshape(int* _piDims, int _iDims)
2656 pSp = reshape(_piDims[0], iCols);
2662 Sparse* Sparse::reshape(int _iNewRows, int _iNewCols)
2664 typedef Sparse* (Sparse::*reshape_t)(int, int);
2665 Sparse* pIT = checkRef(this, (reshape_t)&Sparse::reshape, _iNewRows, _iNewCols);
2671 if (_iNewRows * _iNewCols != getRows() * getCols())
2682 size_t iNonZeros = nonZeros();
2683 RealSparse_t *newReal = new RealSparse_t(_iNewRows, _iNewCols);
2684 newReal->reserve((int)iNonZeros);
2687 int* pRows = new int[iNonZeros * 2];
2688 outputRowCol(pRows);
2689 int* pCols = pRows + iNonZeros;
2692 double* pNonZeroR = new double[iNonZeros];
2693 double* pNonZeroI = new double[iNonZeros];
2694 outputValues(pNonZeroR, pNonZeroI);
2696 typedef Eigen::Triplet<double> triplet;
2697 std::vector<triplet> tripletList;
2699 for (size_t i = 0 ; i < iNonZeros ; i++)
2701 int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
2702 tripletList.push_back(triplet((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), pNonZeroR[i]));
2705 newReal->setFromTriplets(tripletList.begin(), tripletList.end());
2708 matrixReal = newReal;
2716 size_t iNonZeros = nonZeros();
2717 CplxSparse_t *newCplx = new CplxSparse_t(_iNewRows, _iNewCols);
2718 newCplx->reserve((int)iNonZeros);
2721 int* pRows = new int[iNonZeros * 2];
2722 outputRowCol(pRows);
2723 int* pCols = pRows + iNonZeros;
2726 double* pNonZeroR = new double[iNonZeros];
2727 double* pNonZeroI = new double[iNonZeros];
2728 outputValues(pNonZeroR, pNonZeroI);
2730 typedef Eigen::Triplet<std::complex<double> > triplet;
2731 std::vector<triplet> tripletList;
2733 for (size_t i = 0 ; i < iNonZeros ; i++)
2735 int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
2736 tripletList.push_back(triplet((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), std::complex<double>(pNonZeroR[i], pNonZeroI[i])));
2739 newCplx->setFromTriplets(tripletList.begin(), tripletList.end());
2742 matrixCplx = newCplx;
2748 m_iRows = _iNewRows;
2749 m_iCols = _iNewCols;
2750 m_iSize = _iNewRows * _iNewCols;
2753 m_piDims[0] = m_iRows;
2754 m_piDims[1] = m_iCols;
2767 // SparseBool* SparseBool::new
2769 SparseBool::SparseBool(Bool SPARSE_CONST& src)
2772 int size = src.getSize();
2773 int row = src.getRows();
2774 Double* idx = new Double(src.getSize(), 2);
2775 double* p = idx->get();
2776 for (int i = 0; i < size; ++i)
2778 p[i] = (double)(i % row) + 1;
2779 p[i + size] = (double)(i / row) + 1;
2781 create2(src.getRows(), src.getCols(), src, *idx);
2784 Inspector::addItem(this);
2787 /* @param src : Bool matrix to copy into a new sparse matrix
2788 @param idx : Double matrix to use as indexes to get values from the src
2790 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
2792 int idxrow = idx.getRows();
2793 int rows = static_cast<int>(*std::max_element(idx.get(), idx.get() + idxrow));
2794 int cols = static_cast<int>(*std::max_element(idx.get() + idxrow, idx.get() + idxrow * 2));
2795 create2(rows, cols, src, idx);
2797 Inspector::addItem(this);
2801 /* @param src : Bool matrix to copy into a new sparse matrix
2802 @param idx : Double matrix to use as indexes to get values from the src
2803 @param dims : Double matrix containing the dimensions of the new matrix
2805 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx, Double SPARSE_CONST& dims)
2807 create2(static_cast<int>(dims.get(0)), static_cast<int>(dims.get(1)), src, idx);
2809 Inspector::addItem(this);
2813 SparseBool::SparseBool(int _iRows, int _iCols) : matrixBool(new BoolSparse_t(_iRows, _iCols))
2817 m_iSize = _iRows * _iCols;
2819 m_piDims[0] = _iRows;
2820 m_piDims[1] = _iCols;
2822 Inspector::addItem(this);
2826 SparseBool::SparseBool(SparseBool const& src) : matrixBool(new BoolSparse_t(*src.matrixBool))
2829 m_iRows = const_cast<SparseBool*>(&src)->getRows();
2830 m_iCols = const_cast<SparseBool*>(&src)->getCols();
2831 m_iSize = m_iRows * m_iCols;
2832 m_piDims[0] = m_iRows;
2833 m_piDims[1] = m_iCols;
2835 Inspector::addItem(this);
2839 SparseBool::SparseBool(BoolSparse_t* src) : matrixBool(src)
2841 m_iRows = src->rows();
2842 m_iCols = src->cols();
2843 m_iSize = m_iRows * m_iCols;
2845 m_piDims[0] = m_iRows;
2846 m_piDims[1] = m_iCols;
2848 Inspector::addItem(this);
2852 SparseBool::SparseBool(int rows, int cols, int trues, int* inner, int* outer)
2857 matrixBool = new BoolSparse_t(rows, cols);
2858 matrixBool->reserve((int)trues);
2859 out = matrixBool->outerIndexPtr();
2860 in = matrixBool->innerIndexPtr();
2862 //update outerIndexPtr
2863 memcpy(out, outer, sizeof(int) * (rows + 1));
2864 //update innerIndexPtr
2865 memcpy(in, inner, sizeof(int) * trues);
2867 bool* data = matrixBool->valuePtr();
2868 for (int i = 0; i < trues; ++i)
2875 m_iSize = cols * rows;
2877 m_piDims[0] = m_iRows;
2878 m_piDims[1] = m_iCols;
2880 matrixBool->resizeNonZeros(trues);
2883 void SparseBool::create2(int rows, int cols, Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
2885 int nnz = src.getSize();
2886 double* i = idx.get();
2887 double* j = i + idx.getRows();
2888 int* val = src.get();
2890 typedef Eigen::Triplet<bool> T;
2891 std::vector<T> tripletList;
2892 tripletList.reserve((int)nnz);
2894 for (int k = 0; k < nnz; ++k)
2896 tripletList.push_back(T(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, val[k] == 1));
2899 matrixBool = new BoolSparse_t(rows, cols);
2900 matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
2902 m_iRows = matrixBool->rows();
2903 m_iCols = matrixBool->cols();
2904 m_iSize = cols * rows;
2906 m_piDims[0] = m_iRows;
2907 m_piDims[1] = m_iCols;
2911 SparseBool::~SparseBool()
2915 Inspector::removeItem(this);
2919 bool SparseBool::toString(std::wostringstream& ostr)
2921 ostr << ::toString(*matrixBool, 0);
2925 void SparseBool::whoAmI() SPARSE_CONST
2927 std::cout << "types::SparseBool";
2930 SparseBool* SparseBool::clone(void)
2932 return new SparseBool(*this);
2935 SparseBool* SparseBool::resize(int _iNewRows, int _iNewCols)
2937 typedef SparseBool* (SparseBool::*resize_t)(int, int);
2938 SparseBool* pIT = checkRef(this, (resize_t)&SparseBool::resize, _iNewRows, _iNewCols);
2944 if (_iNewRows <= getRows() && _iNewCols <= getCols())
2946 //nothing to do: hence we do NOT fail
2950 SparseBool* res = NULL;
2954 size_t iNonZeros = nbTrue();
2956 BoolSparse_t *newBool = new BoolSparse_t(_iNewRows, _iNewCols);
2957 newBool->reserve((int)iNonZeros);
2960 int* pRows = new int[iNonZeros * 2];
2961 outputRowCol(pRows);
2962 int* pCols = pRows + iNonZeros;
2964 typedef Eigen::Triplet<bool> triplet;
2965 std::vector<triplet> tripletList;
2967 for (size_t i = 0 ; i < iNonZeros ; i++)
2969 tripletList.push_back(triplet((int)pRows[i] - 1, (int)pCols[i] - 1, true));
2972 newBool->setFromTriplets(tripletList.begin(), tripletList.end());
2975 matrixBool = newBool;
2978 m_iRows = _iNewRows;
2979 m_iCols = _iNewCols;
2980 m_iSize = _iNewRows * _iNewCols;
2981 m_piDims[0] = m_iRows;
2982 m_piDims[1] = m_iCols;
2993 SparseBool* SparseBool::insert(typed_list* _pArgs, SparseBool* _pSource)
2995 bool bNeedToResize = false;
2996 int iDims = (int)_pArgs->size();
2999 //sparse are only in 2 dims
3012 //evaluate each argument and replace by appropriate value and compute the count of combinations
3013 int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3017 cleanIndexesArguments(_pArgs, &pArg);
3024 if (getRows() == 1 || getCols() == 1)
3027 if (getSize() < piMaxDim[0])
3029 bNeedToResize = true;
3031 //need to enlarge sparse dimensions
3032 if (getCols() == 1 || getSize() == 0)
3035 iNewRows = piMaxDim[0];
3038 else if (getRows() == 1)
3042 iNewCols = piMaxDim[0];
3046 else if (getSize() < piMaxDim[0])
3049 cleanIndexesArguments(_pArgs, &pArg);
3056 if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
3058 bNeedToResize = true;
3059 iNewRows = std::max(getRows(), piMaxDim[0]);
3060 iNewCols = std::max(getCols(), piMaxDim[1]);
3064 //check number of insertion
3065 if (_pSource->isScalar() == false && _pSource->getSize() != iSeqCount)
3068 cleanIndexesArguments(_pArgs, &pArg);
3072 //now you are sure to be able to insert values
3075 if (resize(iNewRows, iNewCols) == false)
3078 cleanIndexesArguments(_pArgs, &pArg);
3085 double* pIdx = pArg[0]->getAs<Double>()->get();
3086 for (int i = 0 ; i < iSeqCount ; i++)
3088 int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
3089 int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
3091 if (_pSource->isScalar())
3093 set(iRow, iCol, _pSource->get(0, 0), false);
3097 int iRowOrig = i % _pSource->getRows();
3098 int iColOrig = i / _pSource->getRows();
3099 set(iRow, iCol, _pSource->get(iRowOrig, iColOrig), false);
3105 double* pIdxRow = pArg[0]->getAs<Double>()->get();
3106 int iRowSize = pArg[0]->getAs<Double>()->getSize();
3107 double* pIdxCol = pArg[1]->getAs<Double>()->get();
3109 for (int i = 0 ; i < iSeqCount ; i++)
3111 if (_pSource->isScalar())
3113 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(0, 0), false);
3117 int iRowOrig = i % _pSource->getRows();
3118 int iColOrig = i / _pSource->getRows();
3119 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(iRowOrig, iColOrig), false);
3127 cleanIndexesArguments(_pArgs, &pArg);
3132 SparseBool* SparseBool::insert(typed_list* _pArgs, InternalType* _pSource)
3134 typedef SparseBool* (SparseBool::*insert_t)(typed_list*, InternalType*);
3135 SparseBool* pIT = checkRef(this, (insert_t)&SparseBool::insert, _pArgs, _pSource);
3141 if (_pSource->isSparseBool())
3143 return insert(_pArgs, _pSource->getAs<SparseBool>());
3146 bool bNeedToResize = false;
3147 int iDims = (int)_pArgs->size();
3150 //sparse are only in 2 dims
3162 Bool* pSource = _pSource->getAs<Bool>();
3164 //evaluate each argument and replace by appropriate value and compute the count of combinations
3165 int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3169 cleanIndexesArguments(_pArgs, &pArg);
3176 if (getRows() == 1 || getCols() == 1)
3179 bNeedToResize = true;
3180 if (getSize() < piMaxDim[0])
3182 //need to enlarge sparse dimensions
3183 if (getCols() == 1 || getSize() == 0)
3186 iNewRows = piMaxDim[0];
3189 else if (getRows() == 1)
3193 iNewCols = piMaxDim[0];
3197 else if (getSize() < piMaxDim[0])
3200 cleanIndexesArguments(_pArgs, &pArg);
3207 if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
3209 bNeedToResize = true;
3210 iNewRows = std::max(getRows(), piMaxDim[0]);
3211 iNewCols = std::max(getCols(), piMaxDim[1]);
3215 //check number of insertion
3216 if (pSource->isScalar() == false && pSource->getSize() != iSeqCount)
3219 cleanIndexesArguments(_pArgs, &pArg);
3223 //now you are sure to be able to insert values
3226 if (resize(iNewRows, iNewCols) == false)
3229 cleanIndexesArguments(_pArgs, &pArg);
3236 double* pIdx = pArg[0]->getAs<Double>()->get();
3237 for (int i = 0 ; i < iSeqCount ; i++)
3239 int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
3240 int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
3241 if (pSource->isScalar())
3243 set(iRow, iCol, pSource->get(0) != 0, false);
3247 set(iRow, iCol, pSource->get(i) != 0, false);
3253 double* pIdxRow = pArg[0]->getAs<Double>()->get();
3254 int iRowSize = pArg[0]->getAs<Double>()->getSize();
3255 double* pIdxCol = pArg[1]->getAs<Double>()->get();
3257 for (int i = 0 ; i < iSeqCount ; i++)
3259 if (pSource->isScalar())
3261 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(0) != 0, false);
3265 int iRowOrig = i % pSource->getRows();
3266 int iColOrig = i / pSource->getRows();
3268 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(iRowOrig, iColOrig) != 0, false);
3276 cleanIndexesArguments(_pArgs, &pArg);
3280 GenericType* SparseBool::remove(typed_list* _pArgs)
3282 SparseBool* pOut = NULL;
3283 int iDims = (int)_pArgs->size();
3286 //sparse are only in 2 dims
3295 //evaluate each argument and replace by appropriate value and compute the count of combinations
3296 int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3300 cleanIndexesArguments(_pArgs, &pArg);
3304 bool* pbFull = new bool[iDims];
3305 //coord must represent all values on a dimension
3306 for (int i = 0 ; i < iDims ; i++)
3309 int iDimToCheck = getVarMaxDim(i, iDims);
3310 int iIndexSize = pArg[i]->getAs<GenericType>()->getSize();
3312 //we can have index more than once
3313 if (iIndexSize >= iDimToCheck)
3315 //size is good, now check datas
3316 double* pIndexes = getDoubleArrayFromDouble(pArg[i]);
3317 for (int j = 0 ; j < iDimToCheck ; j++)
3320 for (int k = 0 ; k < iIndexSize ; k++)
3322 if ((int)pIndexes[k] == j + 1)
3333 //only one dims can be not full/entire
3334 bool bNotEntire = false;
3336 bool bTooMuchNotEntire = false;
3337 for (int i = 0 ; i < iDims ; i++)
3339 if (pbFull[i] == false)
3341 if (bNotEntire == false)
3348 bTooMuchNotEntire = true;
3356 if (bTooMuchNotEntire == true)
3359 cleanIndexesArguments(_pArgs, &pArg);
3363 //find index to keep
3364 int iNotEntireSize = pArg[iNotEntire]->getAs<GenericType>()->getSize();
3365 double* piNotEntireIndex = getDoubleArrayFromDouble(pArg[iNotEntire]);
3366 int iKeepSize = getVarMaxDim(iNotEntire, iDims);
3367 bool* pbKeep = new bool[iKeepSize];
3369 //fill pbKeep with true value
3370 for (int i = 0 ; i < iKeepSize ; i++)
3375 for (int i = 0 ; i < iNotEntireSize ; i++)
3377 int idx = (int)piNotEntireIndex[i] - 1;
3379 //don't care of value out of bounds
3380 if (idx < iKeepSize)
3382 pbKeep[idx] = false;
3386 int iNewDimSize = 0;
3387 for (int i = 0 ; i < iKeepSize ; i++)
3389 if (pbKeep[i] == true)
3396 int* piNewDims = new int[iDims];
3397 for (int i = 0 ; i < iDims ; i++)
3399 if (i == iNotEntire)
3401 piNewDims[i] = iNewDimSize;
3405 piNewDims[i] = getVarMaxDim(i, iDims);
3409 //remove last dimension if are == 1
3410 int iOrigDims = iDims;
3411 for (int i = (iDims - 1) ; i >= 2 ; i--)
3413 if (piNewDims[i] == 1)
3425 if (iNewDimSize == 0)
3428 cleanIndexesArguments(_pArgs, &pArg);
3429 return new SparseBool(0, 0);
3433 //two cases, depends of original matrix/vector
3434 if ((*_pArgs)[0]->isColon() == false && m_iDims == 2 && m_piDims[0] == 1 && m_piDims[1] != 1)
3436 //special case for row vector
3437 pOut = new SparseBool(1, iNewDimSize);
3438 //in this case we have to care of 2nd dimension
3443 pOut = new SparseBool(iNewDimSize, 1);
3449 pOut = new SparseBool(piNewDims[0], piNewDims[0]);
3453 //find a way to copy existing data to new variable ...
3455 int* piIndexes = new int[iOrigDims];
3456 int* piViewDims = new int[iOrigDims];
3457 for (int i = 0 ; i < iOrigDims ; i++)
3459 piViewDims[i] = getVarMaxDim(i, iOrigDims);
3462 for (int i = 0 ; i < getSize() ; i++)
3464 bool bByPass = false;
3465 getIndexesWithDims(i, piIndexes, piViewDims, iOrigDims);
3467 //check if piIndexes use removed indexes
3468 for (int j = 0 ; j < iNotEntireSize ; j++)
3470 if ((piNotEntireIndex[j] - 1) == piIndexes[iNotEntire])
3472 //by pass this value
3478 if (bByPass == false)
3481 pOut->set(iNewPos, get(i));
3486 //free allocated data
3487 for (int i = 0 ; i < iDims ; i++)
3489 if (pArg[i] != (*_pArgs)[i])
3496 delete[] piViewDims;
3499 cleanIndexesArguments(_pArgs, &pArg);
3504 SparseBool* SparseBool::append(int r, int c, SparseBool SPARSE_CONST* src)
3506 SparseBool* pIT = checkRef(this, &SparseBool::append, r, c, src);
3512 doAppend(*src, r, c, *matrixBool);
3517 GenericType* SparseBool::insertNew(typed_list* _pArgs)
3520 SparseBool *pOut = NULL;
3522 int iDims = (int)_pArgs->size();
3523 int* piMaxDim = new int[iDims];
3524 int* piCountDim = new int[iDims];
3525 bool bUndefine = false;
3527 //evaluate each argument and replace by appropriate value and compute the count of combinations
3528 int iSeqCount = checkIndexesArguments(NULL, _pArgs, &pArg, piMaxDim, piCountDim);
3532 cleanIndexesArguments(_pArgs, &pArg);
3533 return createEmptyDouble();
3538 iSeqCount = -iSeqCount;
3544 //manage : and $ in creation by insertion
3546 int *piSourceDims = getDimsArray();
3548 for (int i = 0 ; i < iDims ; i++)
3550 if (pArg[i] == NULL)
3560 piMaxDim[i] = piSourceDims[iSource];
3561 piCountDim[i] = piSourceDims[iSource];
3564 //replace pArg value by the new one
3565 pArg[i] = createDoubleVector(piMaxDim[i]);
3569 // piMaxDim[i] = piCountDim[i];
3574 //remove last dimension at size 1
3575 //remove last dimension if are == 1
3576 for (int i = (iDims - 1) ; i >= 2 ; i--)
3578 if (piMaxDim[i] == 1)
3589 if (checkArgValidity(pArg) == false)
3592 cleanIndexesArguments(_pArgs, &pArg);
3593 //contain bad index, like <= 0, ...
3601 pOut = new SparseBool(piCountDim[0], 1);
3606 pOut = new SparseBool(1, piCountDim[0]);
3611 pOut = new SparseBool(piMaxDim[0], piMaxDim[1]);
3614 //insert values in new matrix
3615 SparseBool* pOut2 = pOut->insert(&pArg, this);
3622 cleanIndexesArguments(_pArgs, &pArg);
3627 SparseBool* SparseBool::extract(int nbCoords, int SPARSE_CONST* coords, int SPARSE_CONST* maxCoords, int SPARSE_CONST* resSize, bool asVector) SPARSE_CONST
3629 if ( (asVector && maxCoords[0] > getSize()) ||
3630 (asVector == false && maxCoords[0] > getRows()) ||
3631 (asVector == false && maxCoords[1] > getCols()))
3636 SparseBool * pSp (0);
3639 pSp = (getRows() == 1) ? new SparseBool(1, resSize[0]) : new SparseBool(resSize[0], 1);
3640 mycopy_n(makeMatrixIterator<bool>(*this, Coords<true>(coords, getRows())), nbCoords
3641 , makeMatrixIterator<bool>(*(pSp->matrixBool), RowWiseFullIterator(pSp->getRows(), pSp->getCols())));
3645 pSp = new SparseBool(resSize[0], resSize[1]);
3646 mycopy_n(makeMatrixIterator<bool>(*this, Coords<false>(coords, getRows())), nbCoords
3647 , makeMatrixIterator<bool>(*(pSp->matrixBool), RowWiseFullIterator(pSp->getRows(), pSp->getCols())));
3654 * create a new SparseBool of dims according to resSize and fill it from currentSparseBool (along coords)
3656 GenericType* SparseBool::extract(typed_list* _pArgs)
3658 SparseBool* pOut = NULL;
3659 int iDims = (int)_pArgs->size();
3662 int* piMaxDim = new int[iDims];
3663 int* piCountDim = new int[iDims];
3665 //evaluate each argument and replace by appropriate value and compute the count of combinations
3666 int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3670 cleanIndexesArguments(_pArgs, &pArg);
3671 if (_pArgs->size() == 0)
3674 delete[] piCountDim;
3681 delete[] piCountDim;
3683 return Double::Empty();
3689 // Check that we stay inside the input size.
3690 if (piMaxDim[0] <= getSize())
3695 if (getRows() == 1 && getCols() != 1 && (*_pArgs)[0]->isColon() == false)
3697 //special case for row vector
3699 iNewCols = piCountDim[0];
3703 iNewRows = piCountDim[0];
3707 pOut = new SparseBool(iNewRows, iNewCols);
3708 double* pIdx = pArg[0]->getAs<Double>()->get();
3709 // Write in output all elements extract from input.
3710 for (int i = 0 ; i < iSeqCount ; i++)
3718 int iRowRead = static_cast<int>(pIdx[i] - 1) % getRows();
3719 int iColRead = static_cast<int>(pIdx[i] - 1) / getRows();
3721 int iRowWrite = static_cast<int>(i) % iNewRows;
3722 int iColWrite = static_cast<int>(i) / iNewRows;
3724 bool bValue = get(iRowRead, iColRead);
3727 //only non zero values
3728 pOut->set(iRowWrite, iColWrite, true, false);
3735 delete[] piCountDim;
3737 cleanIndexesArguments(_pArgs, &pArg);
3743 // Check that we stay inside the input size.
3744 if (piMaxDim[0] <= getRows() && piMaxDim[1] <= getCols())
3746 double* pIdxRow = pArg[0]->getAs<Double>()->get();
3747 double* pIdxCol = pArg[1]->getAs<Double>()->get();
3749 int iNewRows = pArg[0]->getAs<Double>()->getSize();
3750 int iNewCols = pArg[1]->getAs<Double>()->getSize();
3752 pOut = new SparseBool(iNewRows, iNewCols);
3755 // Write in output all elements extract from input.
3756 for (int iRow = 0 ; iRow < iNewRows ; iRow++)
3758 for (int iCol = 0 ; iCol < iNewCols ; iCol++)
3760 if ((pIdxRow[iRow] < 1) || (pIdxCol[iCol] < 1))
3765 delete[] piCountDim;
3767 cleanIndexesArguments(_pArgs, &pArg);
3770 bool bValue = get((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
3773 //only non zero values
3774 pOut->set(iRow, iCol, true, false);
3783 delete[] piCountDim;
3785 cleanIndexesArguments(_pArgs, &pArg);
3793 delete[] piCountDim;
3795 cleanIndexesArguments(_pArgs, &pArg);
3800 bool SparseBool::invoke(typed_list & in, optional_list &/*opt*/, int /*_iRetCount*/, typed_list & out, const ast::Exp & e)
3804 out.push_back(this);
3808 InternalType * _out = extract(&in);
3811 std::wostringstream os;
3812 os << _W("Invalid index.\n");
3813 throw ast::InternalError(os.str(), 999, e.getLocation());
3815 out.push_back(_out);
3821 bool SparseBool::isInvokable() const
3826 bool SparseBool::hasInvokeOption() const
3831 int SparseBool::getInvokeNbIn()
3836 int SparseBool::getInvokeNbOut()
3841 std::size_t SparseBool::nbTrue() const
3843 return matrixBool->nonZeros() ;
3845 std::size_t SparseBool::nbTrue(std::size_t r) const
3847 int* piIndex = matrixBool->outerIndexPtr();
3848 return piIndex[r + 1] - piIndex[r];
3852 void SparseBool::setTrue(bool finalize)
3854 int rows = getRows();
3855 int cols = getCols();
3857 typedef Eigen::Triplet<bool> triplet;
3858 std::vector<triplet> tripletList;
3860 for (int i = 0; i < rows; ++i)
3862 for (int j = 0; j < cols; ++j)
3864 tripletList.push_back(triplet(i, j, true));
3868 matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
3872 matrixBool->finalize();
3876 void SparseBool::setFalse(bool finalize)
3878 int rows = getRows();
3879 int cols = getCols();
3881 typedef Eigen::Triplet<bool> triplet;
3882 std::vector<triplet> tripletList;
3884 for (int i = 0; i < rows; ++i)
3886 for (int j = 0; j < cols; ++j)
3888 tripletList.push_back(triplet(i, j, false));
3892 matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
3896 matrixBool->finalize();
3900 int* SparseBool::getNbItemByRow(int* _piNbItemByRows)
3902 int* piNbItemByRows = new int[getRows() + 1];
3903 mycopy_n(matrixBool->outerIndexPtr(), getRows() + 1, piNbItemByRows);
3905 for (int i = 0 ; i < getRows() ; i++)
3907 _piNbItemByRows[i] = piNbItemByRows[i + 1] - piNbItemByRows[i];
3910 delete[] piNbItemByRows;
3911 return _piNbItemByRows;
3914 int* SparseBool::getColPos(int* _piColPos)
3916 mycopy_n(matrixBool->innerIndexPtr(), nbTrue(), _piColPos);
3917 for (int i = 0; i < nbTrue(); i++)
3925 int* SparseBool::outputRowCol(int* out)const
3927 return sparseTransform(*matrixBool, sparseTransform(*matrixBool, out, GetRow<BoolSparse_t>()), GetCol<BoolSparse_t>());
3930 int* SparseBool::getInnerPtr(int* count)
3932 *count = matrixBool->innerSize();
3933 return matrixBool->innerIndexPtr();
3936 int* SparseBool::getOuterPtr(int* count)
3938 *count = matrixBool->outerSize();
3939 return matrixBool->outerIndexPtr();
3943 bool SparseBool::operator==(const InternalType& it) SPARSE_CONST
3945 SparseBool* otherSparse = const_cast<SparseBool*>(dynamic_cast<SparseBool const*>(&it));/* types::GenericType is not const-correct :( */
3947 && (otherSparse->getRows() == getRows())
3948 && (otherSparse->getCols() == getCols())
3949 && equal(*matrixBool, *otherSparse->matrixBool));
3952 bool SparseBool::operator!=(const InternalType& it) SPARSE_CONST
3954 return !(*this == it);
3957 void SparseBool::finalize()
3959 matrixBool->prune(&keepForSparse<bool>);
3960 matrixBool->finalize();
3963 bool SparseBool::get(int r, int c) SPARSE_CONST
3965 return matrixBool->coeff(r, c);
3968 SparseBool* SparseBool::set(int _iRows, int _iCols, bool _bVal, bool _bFinalize) SPARSE_CONST
3970 typedef SparseBool* (SparseBool::*set_t)(int, int, bool, bool);
3971 SparseBool* pIT = checkRef(this, (set_t)&SparseBool::set, _iRows, _iCols, _bVal, _bFinalize);
3977 matrixBool->coeffRef(_iRows, _iCols) = _bVal;
3987 void SparseBool::fill(Bool& dest, int r, int c) SPARSE_CONST
3989 mycopy_n(makeMatrixIterator<bool >(*matrixBool, RowWiseFullIterator(getRows(), getCols())), getSize()
3990 , makeMatrixIterator<bool >(dest, RowWiseFullIterator(dest.getRows(), dest.getCols(), r, c)));
3993 Sparse* SparseBool::newOnes() const
3995 return new Sparse(new types::Sparse::RealSparse_t(matrixBool->cast<double>()), 0);
3998 SparseBool* SparseBool::newNotEqualTo(SparseBool const&o) const
4000 return cwiseOp<std::not_equal_to>(*this, o);
4003 SparseBool* SparseBool::newEqualTo(SparseBool& o)
4005 int rowL = getRows();
4006 int colL = getCols();
4008 int rowR = o.getRows();
4009 int colR = o.getCols();
4010 int row = std::max(rowL, rowR);
4011 int col = std::max(colL, colR);
4013 //create a boolean sparse matrix with dims of sparses
4014 types::SparseBool* ret = new types::SparseBool(row, col);
4016 if (isScalar() && o.isScalar())
4019 bool r = o.get(0, 0);
4020 ret->set(0, 0, l == r, false);
4022 else if (isScalar())
4024 int nnzR = static_cast<int>(o.nbTrue());
4025 std::vector<int> rowcolR(nnzR * 2, 0);
4026 o.outputRowCol(rowcolR.data());
4028 //compare all items of R with R[0]
4030 for (int i = 0; i < nnzR; ++i)
4032 bool r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
4033 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
4036 else if (o.isScalar())
4038 int nnzL = static_cast<int>(nbTrue());
4039 std::vector<int> rowcolL(nnzL * 2, 0);
4040 outputRowCol(rowcolL.data());
4043 for (int i = 0; i < nnzL; ++i)
4045 bool l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
4046 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
4051 int nnzR = static_cast<int>(o.nbTrue());
4052 std::vector<int> rowcolR(nnzR * 2, 0);
4053 o.outputRowCol(rowcolR.data());
4054 int nnzL = static_cast<int>(nbTrue());
4055 std::vector<int> rowcolL(nnzL * 2, 0);
4056 outputRowCol(rowcolL.data());
4057 //set all values to %t
4058 ret->setTrue(false);
4059 //set %f in each pL values
4060 for (int i = 0; i < nnzL; ++i)
4062 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, false, false);
4066 //set _pR[i] == _pL[i] for each _pR values
4067 for (int i = 0; i < nnzR; ++i)
4069 //get l and r following non zeros value of R
4070 bool l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
4071 bool r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
4072 //set value following non zeros value of R
4073 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
4081 SparseBool* SparseBool::newLogicalOr(SparseBool const&o) const
4083 return cwiseOp<std::logical_or>(*this, o);
4086 SparseBool* SparseBool::newLogicalAnd(SparseBool const&o) const
4088 return cwiseOp<std::logical_and>(*this, o);
4091 SparseBool* SparseBool::reshape(int* _piDims, int _iDims)
4093 SparseBool* pSpBool = NULL;
4103 pSpBool = reshape(_piDims[0], iCols);
4109 SparseBool* SparseBool::reshape(int _iNewRows, int _iNewCols)
4111 typedef SparseBool* (SparseBool::*reshape_t)(int, int);
4112 SparseBool* pIT = checkRef(this, (reshape_t)&SparseBool::reshape, _iNewRows, _iNewCols);
4118 if (_iNewRows * _iNewCols != getRows() * getCols())
4123 SparseBool* res = NULL;
4127 size_t iNonZeros = matrixBool->nonZeros();
4128 BoolSparse_t *newBool = new BoolSparse_t(_iNewRows, _iNewCols);
4129 newBool->reserve((int)iNonZeros);
4132 int* pRows = new int[iNonZeros * 2];
4133 outputRowCol(pRows);
4134 int* pCols = pRows + iNonZeros;
4136 typedef Eigen::Triplet<bool> triplet;
4137 std::vector<triplet> tripletList;
4139 for (size_t i = 0 ; i < iNonZeros ; i++)
4141 int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
4142 tripletList.push_back(triplet((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), true));
4145 newBool->setFromTriplets(tripletList.begin(), tripletList.end());
4148 matrixBool = newBool;
4151 m_iRows = _iNewRows;
4152 m_iCols = _iNewCols;
4153 m_iSize = _iNewRows * _iNewCols;
4156 m_piDims[0] = m_iRows;
4157 m_piDims[1] = m_iCols;
4170 bool SparseBool::transpose(InternalType *& out)
4172 out = new SparseBool(new BoolSparse_t(matrixBool->transpose()));
4176 template<typename T>
4177 void neg(const int r, const int c, const T * const in, Eigen::SparseMatrix<bool, 1> * const out)
4179 for (int i = 0; i < r; i++)
4181 for (int j = 0; j < c; j++)
4183 out->coeffRef(i, j) = !in->coeff(i, j);
4187 out->prune(&keepForSparse<bool>);