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