COPASI API  4.16.103
CSparseMatrix.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2013 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 // Copyright (C) 2008 - 2009 by Pedro Mendes, Virginia Tech Intellectual
7 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
8 // and The University of Manchester.
9 // All rights reserved.
10 
11 // Copyright (C) 2005 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 ///////////////////////////////////////////////////////////
16 // CSparseMatrix.cpp
17 // Implementation of the Class CSparseMatrix
18 // Created on: 29-Jun-2004 16:52:05
19 ///////////////////////////////////////////////////////////
20 
21 #include <algorithm>
22 #include <cmath>
23 #include <limits>
24 
25 #include "copasi.h"
26 
27 #include "CSparseMatrix.h"
28 #include "CMatrix.h"
29 
30 #include "lapack/blaswrap.h"
31 #include "lapack/lapackwrap.h"
32 
33 // ---------- CSparseMatrixElement
34 
36  const CSparseMatrixElement * pRhs)
37 {return pLhs->mRow < pRhs->mRow;}
38 
40  const CSparseMatrixElement * pRhs)
41 {return pLhs->mCol < pRhs->mCol;}
42 
44  const size_t & row,
45  const size_t & col,
46  const C_FLOAT64 & value):
47  mMatrix(matrix),
48  mRow(row),
49  mCol(col),
50  mValue(value)
51 {}
52 
54 {}
55 
57 {
58  if (fabs(value) < mMatrix.getTreshold() &&
59  fabs(mValue) >= mMatrix.getTreshold())
61  else if (fabs(value) >= mMatrix.getTreshold() &&
62  fabs(mValue) < mMatrix.getTreshold())
63  mMatrix.insert(mRow, mCol, value);
64  else
65  mValue = value;
66 
67  return value;
68 }
69 
70 CSparseMatrixElement::operator const C_FLOAT64 & () const
71 {return mValue;}
72 
73 const size_t & CSparseMatrixElement::row() const {return mRow;}
74 const size_t & CSparseMatrixElement::col() const {return mCol;}
75 
76 // ---------- CSparseMatrix
77 
78 CSparseMatrix::CSparseMatrix(const size_t & rows,
79  const size_t & cols):
80  mThreshold(std::numeric_limits< C_FLOAT64 >::epsilon()),
81  mNumRows(0),
82  mNumCols(0),
83  mRowIndex(),
84  mColIndex(),
85  mRows(),
86  mCols(),
87  mSearchRow(0),
88  mSearchCol(0),
89  mElement(*this, this->mSearchRow, this->mSearchCol, 0.0)
90 {resize(rows, cols);}
91 
93  mThreshold(std::numeric_limits< C_FLOAT64 >::epsilon()),
94  mNumRows(0),
95  mNumCols(0),
96  mRowIndex(),
97  mColIndex(),
98  mRows(),
99  mCols(),
100  mSearchRow(0),
101  mSearchCol(0),
102  mElement(*this, this->mSearchRow, this->mSearchCol, 0.0)
103 {*this = ccf;}
104 
106 {cleanup();}
107 
109 {
110  std::vector< std::vector< CSparseMatrixElement * > >::iterator itRow;
111  std::vector< std::vector< CSparseMatrixElement * > >::iterator endRow;
112 
113  std::vector< CSparseMatrixElement * >::iterator itElement;
114  std::vector< CSparseMatrixElement * >::iterator endElement;
115 
116  for (itRow = mRows.begin(), endRow = mRows.end(); itRow != endRow; ++itRow)
117  for (itElement = itRow->begin(), endElement = itRow->end(); itElement != endElement; ++itElement)
118  delete *itElement;
119 }
120 
121 bool CSparseMatrix::resize(const size_t & rows, const size_t & cols)
122 {
123  size_t i;
124  mNumRows = rows;
125  mNumCols = cols;
126 
127  mRowIndex.resize(mNumRows);
128 
129  for (i = 0; i < mNumRows; i++)
130  mRowIndex[i] = i;
131 
132  mColIndex.resize(mNumCols);
133 
134  for (i = 0; i < mNumCols; i++)
135  mColIndex[i] = i;
136 
137  mRows.resize(mNumRows);
138  mCols.resize(mNumCols);
139 
140  return true;
141 }
142 
143 const size_t & CSparseMatrix::numRows() const
144 {return mNumRows;}
145 
146 const size_t & CSparseMatrix::numCols() const
147 {return mNumCols;}
148 
150 {
151  size_t NonZeros = 0;
152 
153  std::vector< std::vector< CSparseMatrixElement * > >::const_iterator itRow;
154  std::vector< std::vector< CSparseMatrixElement * > >::const_iterator endRow;
155 
156  for (itRow = mRows.begin(), endRow = mRows.end(); itRow != endRow; ++itRow)
157  NonZeros += itRow->size();
158 
159  return NonZeros;
160 }
161 
162 const std::vector< std::vector< CSparseMatrixElement * > > & CSparseMatrix::getColumns() const
163 {return mCols;}
164 
165 std::vector< std::vector< CSparseMatrixElement * > > & CSparseMatrix::getColumns()
166 {return mCols;}
167 
168 const std::vector< std::vector< CSparseMatrixElement * > > & CSparseMatrix::getRows() const
169 {return mRows;}
170 
171 std::vector< std::vector< CSparseMatrixElement * > > & CSparseMatrix::getRows()
172 {return mRows;}
173 
175 {
176  resize(ccf.numRows(), ccf.numCols());
177 
178  const C_FLOAT64 * pValue = ccf.getValues();
179  const size_t * pRowIndex = ccf.getRowIndex();
180  const size_t * pColumnStart = ccf.getColumnStart() + 1;
181 
182  size_t i, j, k;
183  CSparseMatrixElement * pElement;
184 
185  for (j = 0, i = 0; j < mNumCols; j++, pColumnStart++)
186  {
187  mCols[j].resize(*pColumnStart - * (pColumnStart - 1));
188 
189  for (k = 0; i < *pColumnStart; i++, k++, pRowIndex++, pValue++)
190  {
191  pElement = new CSparseMatrixElement(*this, mRowIndex[*pRowIndex], mColIndex[j], *pValue);
192  mCols[j][k] = pElement;
193  mRows[*pRowIndex].push_back(pElement); // Note, these are out of order and must be sorted.
194  }
195  }
196 
197  for (i = 0; i < mNumRows; i++)
198  std::sort(mRows[i].begin(), mRows[i].end(), &CSparseMatrixElement::compareCol);
199 
200  return *this;
201 }
202 
204 {
205  size_t i, j;
206  CSparseMatrixElement * pElement;
207 
208  resize(matrix.numRows(), matrix.numCols());
209  const C_FLOAT64 * pTmp = matrix.array();
210 
211  for (i = 0; i < mNumRows; i++)
212  for (j = 0; j < mNumCols; j++, pTmp++)
213  if (fabs(*pTmp) >= mThreshold)
214  {
215  pElement = new CSparseMatrixElement(*this, mRowIndex[i], mColIndex[j], *pTmp);
216  mRows[i].push_back(pElement);
217  mCols[j].push_back(pElement); // Note, these are out of order and must be sorted.
218  }
219 
220  for (j = 0; j < mNumCols; j++)
221  std::sort(mCols[j].begin(), mCols[j].end(), &CSparseMatrixElement::compareRow);
222 
223  return *this;
224 }
225 
226 CSparseMatrix::operator const CMatrix< C_FLOAT64 > () const
227 {
229  M.resize(mNumRows, mNumCols);
230 
231  C_FLOAT64 * pTmp = M.array();
232  size_t i;
233 
234  std::vector< std::vector< CSparseMatrixElement * > >::const_iterator itRow;
235  std::vector< std::vector< CSparseMatrixElement * > >::const_iterator endRow;
236 
237  std::vector< CSparseMatrixElement * >::const_iterator itElement;
238  std::vector< CSparseMatrixElement * >::const_iterator endElement;
239 
240  for (itRow = mRows.begin(), endRow = mRows.end(); itRow != endRow; ++itRow)
241  {
242  for (i = 0, itElement = itRow->begin(), endElement = itRow->end();
243  itElement != endElement; ++itElement, i++, pTmp++)
244  {
245  for (; i < (*itElement)->col(); i++, pTmp++)
246  *pTmp = 0.0;
247 
248  *pTmp = **itElement;
249  }
250 
251  for (; i < mNumCols; i++, pTmp++)
252  *pTmp = 0.0;
253  }
254 
255  return M;
256 }
257 
258 bool CSparseMatrix::setTreshold(const C_FLOAT64 & threshold)
259 {
260  if (threshold < mThreshold)
261  {
262  std::vector< std::vector< CSparseMatrixElement * > >::iterator itRow;
263  std::vector< std::vector< CSparseMatrixElement * > >::iterator endRow;
264 
265  std::vector< CSparseMatrixElement * >::reverse_iterator itElement;
266  std::vector< CSparseMatrixElement * >::reverse_iterator endElement;
267 
268  for (itRow = mRows.begin(), endRow = mRows.end(); itRow != endRow; ++itRow)
269  for (itElement = itRow->rbegin(), endElement = itRow->rend(); itElement != endElement; ++itElement)
270  if (fabs(**itElement) < threshold)
271  remove((*itElement)->row(), (*itElement)->col());
272  }
273 
274  mThreshold = threshold;
275 
276  return true;
277 }
278 
280 {return mThreshold;}
281 
283  const size_t & col)
284 {
285  mSearchRow = row;
286  mSearchCol = col;
287 
288  std::vector< CSparseMatrixElement * >::iterator found
289  = std::lower_bound(mCols[col].begin(), mCols[col].end(), &mElement, CSparseMatrixElement::compareRow);
290 
291  if (found != mCols[col].end() &&
292  (*found)->row() == row) return **found;
293 
294  return mElement;
295 }
296 
298  const size_t & col) const
299 {
300  const_cast<CSparseMatrix *>(this)->mSearchRow = row;
301  const_cast<CSparseMatrix *>(this)->mSearchCol = col;
302 
303  std::vector< CSparseMatrixElement * >::const_iterator found
304  = std::lower_bound(mCols[col].begin(), mCols[col].end(), &mElement, CSparseMatrixElement::compareRow);
305 
306  if (found != mCols[col].end() &&
307  (*found)->row() == row) return **found;
308 
309  return mElement;
310 }
311 
312 bool CSparseMatrix::insert(const size_t & row,
313  const size_t & col,
314  const C_FLOAT64 & value)
315 {
316  mSearchRow = row;
317  mSearchCol = col;
318 
319  std::vector< CSparseMatrixElement * >::iterator found
320  = std::lower_bound(mCols[col].begin(), mCols[col].end(), &mElement, CSparseMatrixElement::compareRow);
321 
322  if (found != mCols[col].end() &&
323  (*found)->row() == row) return false; // The element already exists.
324 
325  CSparseMatrixElement * pElement = new CSparseMatrixElement(*this, mRowIndex[row], mColIndex[col], value);
326  mCols[col].insert(found, pElement);
327 
328  found = std::lower_bound(mRows[row].begin(), mRows[row].end(), &mElement, CSparseMatrixElement::compareCol);
329  mRows[row].insert(found, pElement);
330 
331  return true;
332 }
333 
334 bool CSparseMatrix::remove(const size_t & row,
335  const size_t & col)
336 {
337  mSearchRow = row;
338  mSearchCol = col;
339 
340  std::vector< CSparseMatrixElement * >::iterator found
341  = std::lower_bound(mCols[col].begin(), mCols[col].end(), &mElement, CSparseMatrixElement::compareRow);
342 
343  if (found == mCols[col].end()) return false; // The element does not exist.
344 
345  mCols[col].erase(found);
346 
347  found = std::lower_bound(mRows[row].begin(), mRows[row].end(), &mElement, CSparseMatrixElement::compareCol);
348  mRows[row].erase(found);
349 
350  delete *found;
351 
352  return true;
353 }
354 
355 std::ostream &operator<<(std::ostream &os, const CSparseMatrix & A)
356 {
357  os << "Matrix(" << A.mNumRows << "x" << A.mNumCols << ")" << std::endl;
358 
359  size_t i;
360 
361  std::vector< std::vector< CSparseMatrixElement * > >::const_iterator itRow;
362  std::vector< std::vector< CSparseMatrixElement * > >::const_iterator endRow;
363 
364  std::vector< CSparseMatrixElement * >::const_iterator itElement;
365  std::vector< CSparseMatrixElement * >::const_iterator endElement;
366 
367  for (itRow = A.mRows.begin(), endRow = A.mRows.end(); itRow != endRow; ++itRow)
368  {
369  for (i = 0, itElement = itRow->begin(), endElement = itRow->end();
370  itElement != endElement; ++itElement, i++)
371  {
372  for (; i < (*itElement)->col(); i++)
373  os << "\t" << 0.0;
374 
375  os << "\t" << **itElement;
376  }
377 
378  for (; i < A.mNumCols; i++)
379  os << "\t" << 0.0;
380 
381  os << std::endl;
382  }
383 
384  return os;
385 }
386 
387 // ---------- CCompressedColumnFormat
388 
390  const size_t & rowIndex):
391  mpMatrix(pMatrix),
392  mRowIndex(rowIndex),
393  mpRowIndex(NULL),
394  mColumnIndex(0),
395  mpColumnIndex(NULL),
396  mpCurrent(NULL)
397 {
398  if (mpMatrix &&
400  mRowIndex < mpMatrix->numRows())
401  {
404 
405  ++(*this);
406  }
407 }
408 
410  mpMatrix(src.mpMatrix),
411  mRowIndex(src.mRowIndex),
412  mpRowIndex(src.mpRowIndex),
413  mColumnIndex(src.mColumnIndex),
414  mpColumnIndex(src.mpColumnIndex),
415  mpCurrent(src.mpCurrent)
416 {}
417 
419 {}
420 
422 {return *mpCurrent;}
423 
425 {return *mpCurrent;}
426 
428 {return mpCurrent != rhs.mpCurrent;}
429 
431 {
432  mpMatrix = rhs.mpMatrix;
433  mRowIndex = rhs.mRowIndex;
434  mpRowIndex = rhs.mpRowIndex;
435  mColumnIndex = rhs.mColumnIndex;
436  mpColumnIndex = rhs.mpColumnIndex;
437  mpCurrent = rhs.mpCurrent;
438 
439  return *this;
440 }
441 
443 {
444  const size_t * pRowIndexEnd = mpMatrix->getRowIndex() + mpMatrix->numNonZeros();
445 
446  mpRowIndex++; // We need to make at least one step forward.
447 
448  while (mpRowIndex != pRowIndexEnd && *mpRowIndex != mRowIndex) mpRowIndex++;
449 
450  if (mpRowIndex != pRowIndexEnd)
451  {
452  size_t index = mpRowIndex - mpMatrix->getRowIndex();
453  mpCurrent = mpMatrix->getValues() + index;
454 
455  while (*mpColumnIndex <= index) mpColumnIndex++;
456 
457  mColumnIndex = mpColumnIndex - mpMatrix->getColumnStart() - 1;
458  }
459  else
460  mpCurrent = NULL;
461 
462  return *this;
463 }
464 
466 {return mColumnIndex;}
467 
469  mNumRows(0),
470  mNumCols(0),
471  mpColumnStart(new size_t[1]),
472  mpRowIndex(NULL),
473  mpValue(NULL)
474 {mpColumnStart[mNumCols] = 0;}
475 
477  const size_t & columns,
478  const size_t & nonZeros):
479  mNumRows(rows),
480  mNumCols(columns),
481  mpColumnStart(new size_t[columns + 1]),
482  mpRowIndex(nonZeros ? new size_t[nonZeros] : NULL),
483  mpValue(nonZeros ? new C_FLOAT64[nonZeros] : NULL)
484 {mpColumnStart[mNumCols] = nonZeros;}
485 
487  mNumRows(0),
488  mNumCols(0),
489  mpColumnStart(NULL),
490  mpRowIndex(NULL),
491  mpValue(NULL)
492 {*this = matrix;}
493 
495 {
496  pdelete(mpValue);
499 }
500 
502 {return mNumRows;}
503 
505 {return mNumCols;}
506 
508 {return mpColumnStart[mNumCols];}
509 
511 {return mpValue;}
512 
514 {return mpValue;}
515 
517 {return mpRowIndex;}
518 
520 {return mpRowIndex;}
521 
523 {return mpColumnStart;}
524 
526 {return mpColumnStart;}
527 
529 {
530  pdelete(mpValue);
533 
534  mNumRows = matrix.numRows();
535  mNumCols = matrix.numCols();
536 
537  mpColumnStart = new size_t[mNumCols + 1];
538  mpColumnStart[mNumCols] = matrix.numNonZeros();
539 
540  if (mpColumnStart[mNumCols])
541  {
542  mpRowIndex = new size_t[mpColumnStart[mNumCols]];
543  mpValue = new C_FLOAT64[mpColumnStart[mNumCols]];
544  }
545  else
546  {
547  mpRowIndex = NULL;
548  mpValue = NULL;
549  }
550 
551  size_t k = 0;
552  size_t i;
553 
554  std::vector< std::vector< CSparseMatrixElement * > >::const_iterator itCol;
555  std::vector< std::vector< CSparseMatrixElement * > >::const_iterator endCol;
556 
557  std::vector< CSparseMatrixElement * >::const_iterator itElement;
558  std::vector< CSparseMatrixElement * >::const_iterator endElement;
559 
560  size_t * pRowIndex = mpRowIndex;
561  C_FLOAT64 * pValue = mpValue;
562 
563  for (i = 0, itCol = matrix.getColumns().begin(), endCol = matrix.getColumns().end();
564  itCol != endCol; ++itCol, i++)
565  {
566  mpColumnStart[i] = k;
567 
568  for (itElement = itCol->begin(), endElement = itCol->end();
569  itElement != endElement; ++itElement, k++, pRowIndex++, pValue++)
570  {
571  *pRowIndex = (*itElement)->row();
572  *pValue = **itElement;
573  }
574  }
575 
576  return *this;
577 }
578 
580 {return const_row_iterator(this, row);}
582 {return const_row_iterator(this);}
583 
584 // ---------- SparseMatrixTest
585 
586 #ifdef COPASI_DEBUG
587 #include "randomGenerator/CRandom.h"
588 #include "report/CCopasiTimer.h"
589 
590 bool SparseMatrixTest(const size_t & size,
591  const C_FLOAT64 & sparseness,
592  const unsigned C_INT32 & seed,
593  const bool & RMP,
594  const bool & dgemmFlag,
595  const bool & SMP,
596  const bool & CCMP)
597 {
598  size_t i, j, l, loop = 1;
599  CRandom * pRandom =
601 
602  // If the sparseness is not specified we expect 4 metabolites per reaction
603  C_FLOAT64 Sparseness = sparseness;
604 
605  if (Sparseness == 0.0) Sparseness = 4.0 / size;
606 
607  CMatrix< C_FLOAT64 > M(size - 3, size);
608  CSparseMatrix S(size - 3, size);
609  CMatrix< C_FLOAT64 > MM(size, size + 3);
610  CSparseMatrix Ss(size, size + 3);
611  C_FLOAT64 tmp;
612 
613  for (i = 0; i < size - 3; i++)
614  for (j = 0; j < size; j++)
615  {
616  if (pRandom->getRandomCC() < Sparseness)
617  S(i, j) = (pRandom->getRandomCC() - 0.5) * 100.0;
618  }
619 
620  for (i = 0; i < size; i++)
621  for (j = 0; j < size + 3; j++)
622  {
623  if (pRandom->getRandomCC() < Sparseness)
624  Ss(i, j) = (pRandom->getRandomCC() - 0.5) * 100.0;
625  }
626 
627  M = S;
628  MM = Ss;
629 
632 
633  std::cout << "Memory requirements for sparseness:\t" << Sparseness << std::endl;
634 
635  tmp = (C_FLOAT64) sizeof(CMatrix< C_FLOAT64 >) + size * size * sizeof(C_FLOAT64);
636  std::cout << "Matrix(" << size << "x" << size << "):\t" << tmp << std::endl;
637 
638  C_FLOAT64 tmp2 = (C_FLOAT64) sizeof(CSparseMatrix)
639  + 2 * size * sizeof(std::vector<CSparseMatrixElement *>)
640  + 2 * size * sizeof(C_FLOAT64)
641  + S.numNonZeros() * sizeof(CSparseMatrixElement);
642  std::cout << "Sparse(" << size << "x" << size << "):\t" << tmp2 << std::endl;
643  std::cout << "Sparse/Matrix:\t" << tmp2 / tmp << std::endl;
644 
645  tmp2 = (C_FLOAT64) sizeof(CCompressedColumnFormat)
646  + 2 * C.numNonZeros() * sizeof(C_FLOAT64)
647  + (size + 1) * sizeof(C_FLOAT64);
648  std::cout << "CompressedColumnFormat(" << size << "x" << size << "):\t" << tmp2 << std::endl;
649  std::cout << "CompressedColumnFormat/Matrix:\t" << tmp2 / tmp << std::endl << std::endl;
650 
653 
654  if (RMP)
655  {
656  // Regular Matrix Product
657  CPU.start();
658  WALL.start();
659 
660  for (l = 0; l < loop; l++)
661  {
662  CMatrix< C_FLOAT64 > MR(M.numRows(), MM.numCols());
663  const C_FLOAT64 *pTmp1, *pTmp2, *pTmp4, *pTmp5;
664  const C_FLOAT64 *pEnd1, *pEnd2, *pEnd4;
665  C_FLOAT64 *pTmp3;
666 
667  size_t LDA = M.numCols();
668  size_t LDB = MM.numCols();
669 
670  pTmp1 = M.array();
671  pEnd1 = pTmp1 + M.numRows() * LDA;
672 
673  pEnd2 = MM.array() + LDB;
674  pTmp3 = MR.array();
675 
676  for (; pTmp1 < pEnd1; pTmp1 += LDA)
677  for (pTmp2 = MM.array(); pTmp2 < pEnd2; pTmp2++, pTmp3++)
678  {
679  *pTmp3 = 0.0;
680 
681  for (pTmp4 = pTmp1, pTmp5 = pTmp2, pEnd4 = pTmp4 + LDA;
682  pTmp4 < pEnd4; pTmp4++, pTmp5 += LDB)
683  * pTmp3 += *pTmp4 * *pTmp5;
684  }
685  }
686 
687  CPU.refresh();
688  WALL.refresh();
689  std::cout << "Matrix * Matrix:\t";
690  CPU.print(&std::cout);
691  std::cout << "\t";
692  WALL.print(&std::cout);
693  std::cout << std::endl;
694  }
695 
696  if (dgemmFlag)
697  {
698  CPU.start();
699  WALL.start();
700 
701  for (l = 0; l < loop; l++)
702  {
703  CMatrix< C_FLOAT64 > dgemmR(M.numRows(), MM.numCols());
704  char T = 'N';
705 
706  C_INT m = (C_INT) MM.numCols(); /* LDA, LDC */
707  C_INT n = (C_INT) M.numRows();
708  C_INT k = (C_INT) M.numCols(); /* LDB */
709 
710  C_FLOAT64 Alpha = 1.0;
711  C_FLOAT64 Beta = 0.0;
712 
713  dgemm_(&T, &T, &m, &n, &k, &Alpha, MM.array(), &m,
714  M.array(), &k, &Beta, dgemmR.array(), &m);
715  }
716 
717  /*
718  for (i = 0; i < MR.numRows(); i++)
719  for (j = 0; j < MR.numCols(); j++)
720  assert(fabs(MR(i, j) - dgemmR(i, j)) <= 100.0 * std::numeric_limits< C_FLOAT64 >::epsilon() * fabs(MR(i, j)));
721  */
722 
723  CPU.refresh();
724  WALL.refresh();
725  std::cout << "dgemm(Matrix, Matrix):\t";
726  CPU.print(&std::cout);
727  std::cout << "\t";
728  WALL.print(&std::cout);
729  std::cout << std::endl;
730  }
731 
732  // Sparse Matrix Product
733  if (SMP)
734  {
735  CPU.start();
736  WALL.start();
737 
738  for (l = 0; l < loop; l++)
739  {
740  CSparseMatrix SR(S.numRows(), Ss.numCols());
741  C_FLOAT64 Tmp;
742  std::vector< std::vector< CSparseMatrixElement * > >::const_iterator itRow;
743  std::vector< std::vector< CSparseMatrixElement * > >::const_iterator endRow;
744  std::vector< CSparseMatrixElement * >::const_iterator itRowElement;
745  std::vector< CSparseMatrixElement * >::const_iterator endRowElement;
746  std::vector< std::vector< CSparseMatrixElement * > >::const_iterator itCol;
747  std::vector< std::vector< CSparseMatrixElement * > >::const_iterator endCol;
748  std::vector< CSparseMatrixElement * >::const_iterator itColElement;
749  std::vector< CSparseMatrixElement * >::const_iterator endColElement;
750 
751  for (itRow = S.getRows().begin(), endRow = S.getRows().end(); itRow != endRow; ++itRow)
752  {
753  endRowElement = itRow->end();
754 
755  for (itCol = Ss.getColumns().begin(), endCol = Ss.getColumns().end(); itCol != endCol; ++itCol)
756  {
757  Tmp = 0;
758  itRowElement = itRow->begin();
759  itColElement = itCol->begin();
760  endColElement = itCol->end();
761 
762  while (itRowElement != endRowElement &&
763  itColElement != endColElement)
764  {
765  while (itRowElement != endRowElement &&
766  (*itRowElement)->col() < (*itColElement)->row()) ++itRowElement;
767 
768  if (itRowElement == endRowElement) break;
769 
770  while (itColElement != endColElement &&
771  (*itColElement)->row() < (*itRowElement)->col()) ++itColElement;
772 
773  if (itColElement == endColElement) break;
774 
775  if ((*itRowElement)->col() != (*itColElement)->row()) continue;
776 
777  Tmp += **itRowElement * **itColElement;
778  ++itRowElement;
779  ++itColElement;
780  }
781 
782  if (fabs(Tmp) < SR.getTreshold()) continue;
783 
784  SR.insert((*itRow->begin())->row(), (*itCol->begin())->col(), Tmp);
785  }
786  }
787  }
788 
789  CPU.refresh();
790  WALL.refresh();
791  std::cout << "Sparse * Sparse:\t";
792  CPU.print(&std::cout);
793  std::cout << "\t";
794  WALL.print(&std::cout);
795  std::cout << std::endl;
796 
797  /*
798  for (i = 0; i < MR.numRows(); i++)
799  for (j = 0; j < MR.numCols(); j++)
800  assert(fabs(MR(i, j) - SR(i, j)) < SR.getTreshold());
801  */
802  }
803 
804  // Compressed Column Format Product
805  if (CCMP)
806  {
807  CPU.start();
808  WALL.start();
809 
810  for (l = 0; l < loop; l++)
811  {
812  CSparseMatrix TmpR(C.numRows(), CC.numCols());
813  CCompressedColumnFormat CR(C.numRows(), CC.numCols(), 0);
814  C_FLOAT64 Tmp;
815  size_t imax = CR.numRows();
816  size_t jmax = CR.numCols();
817  C_FLOAT64 * pColElement, * pEndColElement;
818  size_t * pColElementRow, * pEndColElementRow;
819  size_t * pColStart;
821  CCompressedColumnFormat::const_row_iterator endRowElement = C.endRow(0);
822 
823  for (j = 0, pColStart = CC.getColumnStart(); j < jmax; j++, pColStart++)
824  {
825  for (i = 0; i < imax; i++)
826  {
827  Tmp = 0;
828 
829  itRowElement = C.beginRow(i);
830  pColElement = CC.getValues() + *pColStart;
831  pEndColElement = CC.getValues() + *(pColStart + 1);
832  pColElementRow = CC.getRowIndex() + *pColStart;
833  pEndColElementRow = CC.getRowIndex() + *(pColStart + 1);
834 
835  while (itRowElement != endRowElement &&
836  pColElement != pEndColElement)
837  {
838  while (itRowElement != endRowElement &&
839  itRowElement.getColumnIndex() < *pColElementRow) ++itRowElement;
840 
841  if (!(itRowElement != endRowElement)) break;
842 
843  while (pColElement != pEndColElement &&
844  *pColElementRow < itRowElement.getColumnIndex())
845  {
846  ++pColElement;
847  ++pColElementRow;
848  }
849 
850  if (pColElement == pEndColElement) break;
851 
852  if (itRowElement.getColumnIndex() != *pColElementRow) continue;
853 
854  Tmp += *itRowElement * *pColElement;
855  ++itRowElement;
856  ++pColElement;
857  ++pColElementRow;
858  }
859 
860  if (fabs(Tmp) < TmpR.getTreshold()) continue;
861 
862  TmpR.insert(i, j, Tmp);
863  }
864  }
865 
866  CR = TmpR;
867  }
868 
869  CPU.refresh();
870  WALL.refresh();
871  std::cout << "Compressed * Compressed:\t";
872  CPU.print(&std::cout);
873  std::cout << "\t";
874  WALL.print(&std::cout);
875  std::cout << std::endl;
876 
877  /*
878  for (i = 0; i < MR.numRows(); i++)
879  for (j = 0; j < MR.numCols(); j++)
880  assert(fabs(MR(i, j) - TmpR(i, j)) < SR.getTreshold());
881  */
882  }
883 
884  std::cout << std::endl;
885  std::cout << std::endl;
886 
887  return true;
888 }
889 #endif
const size_t & col() const
#define C_INT
Definition: copasi.h:115
#define pdelete(p)
Definition: copasi.h:215
const C_FLOAT64 & operator=(const C_FLOAT64 &value)
const_row_iterator(const CCompressedColumnFormat *pMatrix=NULL, const size_t &rowIndex=C_INVALID_INDEX)
bool resize(const size_t &rows, const size_t &columns)
const CCompressedColumnFormat * mpMatrix
CSparseMatrix & operator=(const CCompressedColumnFormat &ccf)
const size_t * getColumnStart() const
CCompressedColumnFormat & operator=(const CSparseMatrix &ccf)
const C_FLOAT64 & getTreshold() const
const size_t & numCols() const
static bool compareRow(const CSparseMatrixElement *pLhs, const CSparseMatrixElement *pRhs)
CSparseMatrix & mMatrix
Definition: CSparseMatrix.h:61
const size_t & mRow
Definition: CSparseMatrix.h:62
virtual size_t numRows() const
Definition: CMatrix.h:138
const size_t & numRows() const
std::vector< std::vector< CSparseMatrixElement * > > mRows
#define C_INVALID_INDEX
Definition: copasi.h:222
bool operator!=(const const_row_iterator &rhs)
const_row_iterator beginRow(const size_t &row) const
#define C_INT32
Definition: copasi.h:90
bool remove(const size_t &row, const size_t &col)
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
virtual C_FLOAT64 getRandomCC()
Definition: CRandom.cpp:235
static bool compareCol(const CSparseMatrixElement *pLhs, const CSparseMatrixElement *pRhs)
const std::vector< std::vector< CSparseMatrixElement * > > & getColumns() const
std::ostream & operator<<(std::ostream &os, const CSparseMatrix &A)
const std::vector< std::vector< CSparseMatrixElement * > > & getRows() const
const size_t * getRowIndex() const
const_row_iterator & operator=(const const_row_iterator &rhs)
int dgemm_(char *transa, char *transb, integer *m, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c, integer *ldc)
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
CSparseMatrix(const size_t &rows=0, const size_t &cols=0)
const size_t & row() const
C_FLOAT64 mThreshold
std::vector< size_t > mColIndex
std::vector< size_t > mRowIndex
bool insert(const size_t &row, const size_t &col, const C_FLOAT64 &value)
size_t numNonZeros() const
#define C_FLOAT64
Definition: copasi.h:92
const size_t & mCol
Definition: CSparseMatrix.h:63
std::vector< std::vector< CSparseMatrixElement * > > mCols
bool setTreshold(const C_FLOAT64 &threshold)
CSparseMatrixElement & operator()(const size_t &row, const size_t &col)
const C_FLOAT64 * getValues() const
virtual size_t numCols() const
Definition: CMatrix.h:144
CSparseMatrixElement mElement
const_row_iterator endRow(const size_t &row) const
virtual CType * array()
Definition: CMatrix.h:337
size_t numNonZeros() const