1 #ifndef SPARSEMATRIX_INL_H
2 #define SPARSEMATRIX_INL_H
3
4 #include "mpi.h"
5 #include "DenseVector.h"
6
7 namespace ATC_matrix {
8
9 template <typename T>
TRI_COORD(INDEX row,INDEX col)10 TRI_COORD<T>::TRI_COORD(INDEX row, INDEX col) : i(row), j(col) {}
11 template <typename T>
TRI_COORD(INDEX row,INDEX col,T val,bool add_to)12 TRI_COORD<T>::TRI_COORD(INDEX row, INDEX col, T val, bool add_to)
13 : i(row), j(col), v(val), add(add_to) {}
14
15 //-----------------------------------------------------------------------------
16 // default constructor - creates an empty sparsematrix with specified size
17 //-----------------------------------------------------------------------------
18 template<typename T>
SparseMatrix(INDEX rows,INDEX cols)19 SparseMatrix<T>::SparseMatrix(INDEX rows, INDEX cols)
20 : _val(NULL), _ia(NULL), _ja(NULL), _size(0), _nRowsCRS(0), hasTemplate_(false),
21 _nRows(rows),_nCols(cols) {}
22 //-----------------------------------------------------------------------------
23 // copy constructor
24 //-----------------------------------------------------------------------------
25 template<typename T>
SparseMatrix(const SparseMatrix<T> & C)26 SparseMatrix<T>::SparseMatrix(const SparseMatrix<T>& C)
27 : _val(NULL), _ia(NULL), _ja(NULL), hasTemplate_(false)
28 {
29 _copy(C);
30 }
31 //-----------------------------------------------------------------------------
32 // copy constructor - converts from DenseMatrix
33 //-----------------------------------------------------------------------------
34 template<typename T>
SparseMatrix(const DenseMatrix<T> & C)35 SparseMatrix<T>::SparseMatrix(const DenseMatrix<T>& C)
36 : _val(NULL), _ia(NULL), _ja(NULL), hasTemplate_(false)
37 {
38 reset(C);
39 }
40
41 //-----------------------------------------------------------------------------
42 // constructor - creates a sparse matrix given an array of row indeces,
43 // an array of col indeces, and an array of nonzero values.
44 //-----------------------------------------------------------------------------
45 template<typename T>
SparseMatrix(INDEX * rows,INDEX * cols,T * vals,INDEX size,INDEX nRows,INDEX nCols,INDEX nRowsCRS)46 SparseMatrix<T>::SparseMatrix(INDEX* rows, INDEX* cols, T* vals,
47 INDEX size, INDEX nRows, INDEX nCols, INDEX nRowsCRS)
48 : hasTemplate_(true)
49 {
50 _val = vals;
51 _ia = rows;
52 _ja = cols;
53 _size = size;
54 _nRows = nRows;
55 _nCols = nCols;
56 _nRowsCRS = nRowsCRS;
57 }
58
59 //-----------------------------------------------------------------------------
60 // assigns internal storage for CRS
61 //-----------------------------------------------------------------------------
62 template<typename T>
_create(INDEX size,INDEX nrows)63 void SparseMatrix<T>::_create(INDEX size, INDEX nrows)
64 {
65 _size = size;
66 _nRowsCRS = nrows;
67 // assign memory to hold matrix
68 try
69 {
70 _val = (_size*nrows) ? new T [_size] : NULL;
71 _ia = (_size*nrows) ? new INDEX [_nRowsCRS+1] : NULL;
72 _ja = (_size*nrows) ? new INDEX [_size] : NULL;
73 }
74 catch (std::exception &e)
75 {
76 std::cout << "Could not allocate SparseMatrix of "<< _size << " nonzeros.\n";
77 ERROR_FOR_BACKTRACE
78 exit(EXIT_FAILURE);
79 }
80 if (!_ia) return;
81 // automatically handle the ends of rowpointer
82 *_ia = 0; // first non-zero is the zero index
83 _ia[_nRowsCRS] = _size; // last row pointer is the size
84 }
85 //-----------------------------------------------------------------------------
86 // cleans up internal storage, but retains nRows & nCols
87 //-----------------------------------------------------------------------------
88 template<typename T>
_delete()89 void SparseMatrix<T>::_delete()
90 {
91
92 std::vector<TRI_COORD<T> >().swap(_tri); // completely deletes _tri
93 if (_val) delete [] _val;
94 if (_ia) delete [] _ia;
95 if (_ja) delete [] _ja;
96 _size = _nRowsCRS = 0;
97 _val = NULL;
98 _ia = _ja = NULL;
99 }
100 //-----------------------------------------------------------------------------
101 // full memory copy of C into this
102 //-----------------------------------------------------------------------------
103 template<typename T>
_copy(const SparseMatrix<T> & C)104 void SparseMatrix<T>::_copy(const SparseMatrix<T> &C)
105 {
106 compress(C);
107 _delete();
108 _create(C.size(), C._nRowsCRS);
109 if (_size) {
110 std::copy(C._val, C._val+_size, _val);
111 std::copy(C._ja, C._ja+_size, _ja);
112 }
113 if (_nRowsCRS) {
114 std::copy(C._ia, C._ia+_nRowsCRS+1, _ia);
115 }
116 _nCols = C._nCols;
117 _nRows = C._nRows;
118 if (_nCols > 0 && _nRows > 0) hasTemplate_ = true; // needs if since map seems to call the copy instead of the default constructor
119 }
120 // this version is accessible to derived classes
121 template<typename T>
copy(const SparseMatrix<T> & C)122 void SparseMatrix<T>::copy(const SparseMatrix<T> &C)
123 {
124 _copy(C);
125 }
126 //----------------------------------------------------------------------------
127 // general sparse matrix assignment
128 //----------------------------------------------------------------------------
129 template<typename T>
_set_equal(const Matrix<T> & r)130 void SparseMatrix<T>::_set_equal(const Matrix<T> &r)
131 {
132 this->resize(r.nRows(), r.nCols());
133 const Matrix<T> *ptr_r = &r;
134
135 const SparseMatrix<T> *s_ptr = dynamic_cast<const SparseMatrix<T>*>(ptr_r);
136 if (s_ptr) this->reset(*s_ptr);
137 else if (dynamic_cast<const DiagonalMatrix<T>*>(ptr_r))
138 for (INDEX i=0; i<r.size(); i++) set(i,i,r[i]);
139 else if (dynamic_cast<const DenseMatrix<T>*>(ptr_r)) this->reset(r);
140 else
141 {
142 std::cout <<"Error in general sparse matrix assignment\n";
143 exit(1);
144 }
145 }
146 // General flat index by value operator (by nth nonzero)
147 template <typename T> inline T SparseMatrix<T>::operator[](INDEX i) const
148 {
149 VICK(i); return _val[i];
150 }
151
152 // General flat index by reference operator (by nth nonzero)
153 template <typename T> inline T& SparseMatrix<T>::operator[](INDEX i)
154 {
155 VICK(i); return _val[i];
156 }
157
158 template<typename T>
159 T SparseMatrix<T>::_zero = T(0);
160
161 //-----------------------------------------------------------------------------
162 // triplet comparison operator returns true if x < y
163 //-----------------------------------------------------------------------------
164 template <typename T>
triplet_comparision(const TRI_COORD<T> & x,const TRI_COORD<T> & y)165 bool triplet_comparision(const TRI_COORD<T> &x, const TRI_COORD<T> &y)
166 {
167 const bool row_less = (x.i) < (y.i);
168 const bool row_equal = (x.i) == (y.i);
169 const bool col_less = (x.j) < (y.j);
170 return (row_less || (row_equal && col_less));
171 }
172 //-----------------------------------------------------------------------------
173 // triplet comparison operator returns true if x == y
174 //-----------------------------------------------------------------------------
175 template <typename T>
triplets_equal(const TRI_COORD<T> & x,const TRI_COORD<T> & y)176 bool triplets_equal(const TRI_COORD<T> &x, const TRI_COORD<T> &y)
177 {
178 return x.i==y.i && x.j==y.j;
179 }
180 //-----------------------------------------------------------------------------
181 // multiply sparse matrix by a vector
182 //-----------------------------------------------------------------------------
183 template<typename T>
184 DenseVector<T> operator*(const SparseMatrix<T> &A, const Vector<T>& x)
185 {
186 DenseVector<T> y(A.nRows(), true);
187
188 A.MultMv(x, y);
189
190 return y;
191 }
192 //-----------------------------------------------------------------------------
193 // multiply a vector by a sparse matrix
194 //-----------------------------------------------------------------------------
195 template<typename T>
196 DenseVector<T> operator*(const Vector<T>& x, const SparseMatrix<T> &A)
197 {
198 return A.transMat(x);
199 }
200 //-----------------------------------------------------------------------------
201 // multiply sparse matrix by dense matrix
202 //-----------------------------------------------------------------------------
203 template<typename T>
204 DenseMatrix<T> operator*(const SparseMatrix<T> &A, const Matrix<T>& D)
205 {
206 DenseMatrix<T> C(A.nRows(), D.nCols(), true); // initialized to zero
207
208 A.MultAB(D, C);
209
210 return C;
211 }
212 //-----------------------------------------------------------------------------
213 // multiply sparse matrix by a diagonal matrix - scales each column
214 //-----------------------------------------------------------------------------
215 template<typename T>
216 SparseMatrix<T> operator*(const SparseMatrix<T> &A, const DiagonalMatrix<T>& D)
217 {
218 GCK(A, D, A.nCols()!=D.nRows(),"SparseMatrix * DiagonalMatrix")
219 SparseMatrix<T> C(A); // C has same sparcity as A
220
221 // C(i,j) = A(i,k) * D(k, j) * j==k
222 INDEX i, ij;
223 for (i=0; i<A._nRowsCRS; i++)
224 for (ij=A._ia[i]; ij<A._ia[i+1]; ij++)
225 C[ij] = A._val[ij]*D(A._ja[ij],A._ja[ij]);
226
227 return C;
228 }
229 //-----------------------------------------------------------------------------
230 // multiplies two sparse matrices - assumes their output is sparse
231 //-----------------------------------------------------------------------------
232 template<typename T>
233 SparseMatrix<T> operator*(const SparseMatrix<T> &A, const SparseMatrix<T> &B)
234 {
235 SparseMatrix<T> At(A.transpose());
236 SparseMatrix<T>::compress(B);
237
238 GCK(A, B, A.nCols()!=B.nRows(), "SparseMatrix * SparseMatrix");
239
240 SparseMatrix<T> C(A.nRows(), B.nCols());
241 if (At.empty() || B.empty()) return C;
242
243 INDEX k, ki, kj;
244 INDEX K = std::min(At._nRowsCRS, B._nRowsCRS);
245 for (k=0; k<K; k++) // loop over rows of A or B (smallest)
246 for (ki=At._ia[k]; ki<At._ia[k+1]; ki++) // loop over row nonzeros of A
247 for (kj=B._ia[k]; kj<B._ia[k+1]; kj++) // loop over row nonzeros of B
248 C.add(At._ja[ki], B._ja[kj], At[ki]*B[kj]); // C(i,j) = At(k,i)*B(k, j)
249
250 C.compress();
251 return C;
252 }
253 //-----------------------------------------------------------------------------
254 // returns the first row number with a nonzero entry or -1 if no rows
255 //-----------------------------------------------------------------------------
256 template<typename T>
_first_nonzero_row_crs()257 int SparseMatrix<T>::_first_nonzero_row_crs() const
258 {
259 if (!_nRowsCRS) return -1;
260 INDEX r;
261 for (r=0; r<_nRowsCRS; r++)
262 if (_ia[r+1]>0) return r;
263 return -1;
264 }
265 //-----------------------------------------------------------------------------
266 // converts T to CRS
267 //-----------------------------------------------------------------------------
268 template<typename T>
compress(const SparseMatrix<T> & C)269 void SparseMatrix<T>::compress(const SparseMatrix<T> &C)
270 {
271 const_cast<SparseMatrix<T>*>(&C)->compress();
272 }
273 //-----------------------------------------------------------------------------
274 // merges all the _tri triples with CRS storage
275 //-----------------------------------------------------------------------------
276 template<typename T>
compress()277 void SparseMatrix<T>::compress()
278 {
279 if (_tri.empty()) return;
280
281 // Sort and find the number of unique triplets.
282 // Triplet values will all be not present in existing CRS structure.
283 const INDEX nUnique = CountUniqueTriplets();
284
285 // Max number of rows in new CRS structure.
286 const INDEX nRows = std::max((INDEX)_tri.back().i+1, _nRowsCRS);
287
288 // make a new CRS structure
289 INDEX *ia = new INDEX [nRows+1];
290 INDEX *ja = new INDEX [nUnique];
291 T *val = new T [nUnique];
292
293 // Set first and last row ptr to 0 and nnz respectively.
294 // Set all else to a flagvalue MAX_UNSIGNED (~0).
295 ia[0] = 0;
296 INDEX i;
297 for (i=1; i<nRows; i++) ia[i]=~0; // ~0 is max(INDEX)
298 ia[nRows] = nUnique;
299
300 INDEX crs_pt, crs_row;
301 unsigned tri_ct; // must be unsigned to interface with std::vector without warnings
302
303 // Get the first CRS and triplet coordinates (if they exist).
304 TRI_COORD<T> nextCRS, nextTRI(_tri[0]), next;
305 int first_row = _first_nonzero_row_crs();
306 if (first_row != -1) nextCRS = TRI_COORD<T>(first_row, _ja[0], _val[0]);
307
308 // merge sorted triplets into a new CRS structure
309 crs_pt = crs_row = tri_ct = 0; // initialize counters
310 for (i=0; i<nUnique; i++)
311 {
312 // is the next non-zero in the new triplet vector
313 if (tri_ct < _tri.size()
314 && (triplet_comparision(nextTRI, nextCRS) || crs_pt>=_size)) {
315 next = nextTRI;
316 // advance the triplet counter, and skip voided TRIPLET entries
317 do tri_ct++;
318 while ( tri_ct<_tri.size() && _tri[tri_ct].j == ~0 );
319
320 // if not at the end of the vector, set the next triplet
321 if (tri_ct<_tri.size()) nextTRI = _tri[tri_ct];
322 }
323 // is the next nonzero in the old CRS data
324 else if (crs_pt < _size) {
325 next = nextCRS;
326 // Advance the CRS counter, don't set next if we are at the end.
327 if (++crs_pt < _size) {
328 // advance to the row corresponding to this value
329 while (crs_pt >= _ia[crs_row+1]) {
330 crs_row++;
331 }
332 nextCRS = TRI_COORD<T>(crs_row, _ja[crs_pt], _val[crs_pt]);
333 }
334 }
335 else std::cout << "SparseMatrix - Error in compressing CRS\n";
336
337 // Add next to the new CRS structure.
338 // Is this a new row (is j>0 and is ja[j] == 0)?
339 if (ia[next.i]==~0) ia[next.i] = i;
340 ja[i] = next.j;
341 val[i] = next.v;
342 }
343 // sweep backwards through row pointers and check for skipped rows
344 for (i=nRows-1; i>0; i--) ia[i] = (ia[i]==~0) ? ia[i+1] : ia[i];
345
346 _delete();
347 _val = val;
348 _ia = ia;
349 _ja = ja;
350 _size = nUnique;
351 _nRowsCRS = nRows;
352 hasTemplate_=true;
353 }
354 //-----------------------------------------------------------------------------
355 // Sorts the triplets, condenses duplicates, and returns the # of unique values
356 //-----------------------------------------------------------------------------
357 template<typename T>
CountUniqueTriplets()358 INDEX SparseMatrix<T>::CountUniqueTriplets()
359 {
360 if (_tri.empty()) return _size;
361 std::sort(_tri.begin(), _tri.end(), triplet_comparision<T>);
362 INDEX nUnique=1 + _size;
363
364 typename std::vector<TRI_COORD<T> >::reverse_iterator t;
365 // Loop backwards over all new triplets.
366 for (t = _tri.rbegin(); t+1!=_tri.rend(); ++t) {
367 // If this triplet is the same as the preceding one.
368 if (triplets_equal(*(t+1), *t)) {
369 if (t->add) (t+1)->v += t->v; // Add to previous
370 else (t+1)->v = t->v; // Replace previous -- DOES THIS WORK?
371 t->j = ~0; // Void this entry's column pointer
372 }
373 else nUnique++;
374 }
375 return nUnique;
376 }
377 //-----------------------------------------------------------------------------
378 // Checks if a value has been set
379 //-----------------------------------------------------------------------------
380 template<typename T>
has_entry(INDEX i,INDEX j)381 bool SparseMatrix<T>::has_entry(INDEX i, INDEX j) const
382 {
383 if (has_entry_compressed(i,j)) return true;
384 if (has_entry_uncompressed(i,j)) return true;
385 return false;
386 }
387
388 template<typename T>
has_entry_uncompressed(INDEX i,INDEX j)389 bool SparseMatrix<T>::has_entry_uncompressed(INDEX i, INDEX j) const
390 {
391 for (unsigned k=0; k<_tri.size() ; k++) {
392 if (_tri[k].i == i && _tri[k].j == j) return true;
393 }
394 return false;
395 }
396
397 template<typename T>
has_entry_compressed(INDEX i,INDEX j)398 bool SparseMatrix<T>::has_entry_compressed(INDEX i, INDEX j) const
399 {
400 if (_size == 0) return false;
401 if (i >= _nRowsCRS) return false;
402 if (_ia[i] < _ia[i+1]) {
403 return -1 < ATC_Utility::search_sorted(_ja, j, _ia[i], _ia[i+1]);
404 }
405 return false;
406 }
407 //-----------------------------------------------------------------------------
408 // check if the matrix has been compressed at least once
409 //-----------------------------------------------------------------------------
410 template<typename T>
has_template(void)411 bool SparseMatrix<T>::has_template(void) const
412 {
413 return hasTemplate_;
414 }
415 //-----------------------------------------------------------------------------
416 // Index by copy operator - return zero if not found
417 //-----------------------------------------------------------------------------
418 template<typename T>
operator()419 T SparseMatrix<T>::operator()(INDEX i, INDEX j) const
420 {
421 MICK(i,j); // Matrix Index ChecKing
422 compress(*this);
423 if (i>=_nRowsCRS || _ia[i+1]==_ia[i]) return 0.0;
424 INDEX f = std::lower_bound(_ja+_ia[i], _ja+_ia[i+1]-1, j) - _ja;
425 if (f>=_ia[i] && f<_ia[i+1] && _ja[f] == j) return _val[f];
426 return 0.0;
427 }
428 //-----------------------------------------------------------------------------
429 // Index by reference operator - add to _tri if not found
430 //-----------------------------------------------------------------------------
431 template<typename T>
operator()432 T& SparseMatrix<T>::operator()(INDEX i, INDEX j)
433 {
434 MICK(i,j); // Matrix Index ChecKing
435 compress(*this);
436 if (i < _nRowsCRS && _ia[i+1]>_ia[i]) {
437 INDEX f = std::lower_bound(_ja+_ia[i], _ja+_ia[i+1]-1, j) - _ja;
438 if (f>=_ia[i] && f<_ia[i+1] && _ja[f] == j) return _val[f];
439 }
440 // NEVER use index operator as LHS to modify values not already in the
441 // sparcity pattern - the crude check below will only catch this on the
442 // second infraction.
443 if (_zero != T(0)) std::cout << "Use add or set for SparseMatrix\n";
444 return _zero;
445 }
446 //-----------------------------------------------------------------------------
447 // Sets (i,j) to value
448 //-----------------------------------------------------------------------------
449 template<typename T>
set(INDEX i,INDEX j,T v)450 void SparseMatrix<T>::set(INDEX i, INDEX j, T v)
451 {
452 MICK(i,j); // Matrix Index ChecKing
453 if (i < _nRowsCRS)
454 {
455 const int loc = ATC_Utility::search_sorted(_ja, j, _ia[i], _ia[i+1]);
456 if (loc >=0 )
457 {
458 _val[loc] = v;
459 return;
460 }
461 }
462 _tri.push_back(TRI_COORD<T>(i,j,v,false));
463 }
464 //-----------------------------------------------------------------------------
465 // Adds (i,j) to value
466 //-----------------------------------------------------------------------------
467 template<typename T>
add(INDEX i,INDEX j,T v)468 void SparseMatrix<T>::add(INDEX i, INDEX j, T v)
469 {
470 MICK(i,j); // Matrix Index ChecKing
471 if (i < _nRowsCRS)
472 {
473 const int loc = ATC_Utility::search_sorted(_ja, j, _ia[i], _ia[i+1]);
474 if (loc >=0 )
475 {
476 _val[loc] += v;
477 return;
478 }
479 }
480 _tri.push_back(TRI_COORD<T>(i,j,v,true));
481 }
482 //-----------------------------------------------------------------------------
483 // returns a triplet value of the ith nonzero
484 //-----------------------------------------------------------------------------
485 template<typename T>
triplet(INDEX i)486 TRIPLET<T> SparseMatrix<T>::triplet(INDEX i) const
487 {
488 compress(*this);
489 if (i >= _ia[_nRowsCRS]) {
490 gerror("ERROR: tried indexing triplet of sparse matrix beyond range");
491 }
492
493 INDEX row(std::lower_bound(_ia, _ia+_nRowsCRS, i)-_ia);
494 row -= _ia[row] != i;
495 return TRIPLET<T>(row, _ja[i], _val[i]);
496 }
497 //-----------------------------------------------------------------------------
498 // full reset - completely wipes out all SparseMatrix data, zero is ignored
499 //-----------------------------------------------------------------------------
500 template<typename T>
reset(INDEX rows,INDEX cols,bool zero)501 void SparseMatrix<T>::reset(INDEX rows, INDEX cols, bool zero)
502 {
503 _delete();
504 _nRows = rows;
505 _nCols = cols;
506 }
507 //-----------------------------------------------------------------------------
508 // resize - changes the _nRows and _nCols without changing anything else if
509 // the matrix is being enlarged, other wise wipes it
510 //-----------------------------------------------------------------------------
511 template<typename T>
resize(INDEX rows,INDEX cols,bool copy)512 void SparseMatrix<T>::resize(INDEX rows, INDEX cols, bool copy)
513 {
514 //if (copy) throw;
515 if (_nRowsCRS>rows) {
516 _delete();
517 }
518 if (copy)
519 _nRows = rows;
520 _nCols = cols; // a check on this would be expensive
521 }
522 //-----------------------------------------------------------------------------
523 // get sparsity from DenseMatrix, if TOL < 0, then only zero values are added
524 //-----------------------------------------------------------------------------
525 template<typename T>
reset(const DenseMatrix<T> & D,double TOL)526 void SparseMatrix<T>::reset(const DenseMatrix<T>& D, double TOL)
527 {
528 _delete(); // clears all values
529 // if TOL is specified then TOL = TOL^2 * max(abs(D))^2
530 if (TOL > 0.0)
531 {
532 TOL *= D.maxabs();
533 TOL *= TOL;
534 }
535 _nRows = D.nRows();
536 _nCols = D.nCols();
537 for (INDEX i=0; i<D.nRows(); i++)
538 for (INDEX j=0; j<D.nCols(); j++)
539 if (D(i,j)*D(i,j) >= TOL) // if TOL wasn't specified then TOL < 0
540 set(i, j, D(i,j));
541
542 compress();
543 }
544 //-----------------------------------------------------------------------------
545 // copy - dangerous: ignores rows & columns
546 //-----------------------------------------------------------------------------
547 template<typename T>
copy(const T * ptr,INDEX rows,INDEX cols)548 void SparseMatrix<T>::copy(const T * ptr, INDEX rows, INDEX cols)
549 {
550 std::cout << "SparseMatrix<T>::copy() has no effect.\n";
551 throw;
552 }
553 //-----------------------------------------------------------------------------
554 // dense_copy - copy to dense matrix
555 //-----------------------------------------------------------------------------
556 template<typename T>
dense_copy(DenseMatrix<T> & D)557 void SparseMatrix<T>::dense_copy(DenseMatrix <T> & D ) const
558 {
559 SparseMatrix<T>::compress(*this);
560 D.reset(nRows(),nCols());
561 for (INDEX i=0; i<_nRowsCRS; i++)
562 for (INDEX j=_ia[i]; j<_ia[i+1]; j++)
563 D(i, _ja[j]) = _val[j];
564 }
565 template<typename T>
dense_copy(void)566 DenseMatrix <T> SparseMatrix<T>::dense_copy(void) const
567 {
568 DenseMatrix<T> D;
569 dense_copy(D);
570 return D;
571 }
572 //-----------------------------------------------------------------------------
573 // returns true if the matrix has no non-zero elements
574 //-----------------------------------------------------------------------------
575 template<typename T>
empty()576 bool SparseMatrix<T>::empty() const
577 {
578 return _size==0 && _tri.empty();
579 }
580 //-----------------------------------------------------------------------------
581 // returns the number of rows specified by the user
582 //-----------------------------------------------------------------------------
583 template<typename T>
nRows()584 inline INDEX SparseMatrix<T>::nRows() const
585 {
586 return _nRows;
587 }
588 //-----------------------------------------------------------------------------
589 // returns ??????????????????????
590 //-----------------------------------------------------------------------------
591 template<typename T>
nRowsCRS()592 inline INDEX SparseMatrix<T>::nRowsCRS() const
593 {
594 return _nRowsCRS;
595 }
596 //-----------------------------------------------------------------------------
597 // returns the number of columns specified by the user
598 //-----------------------------------------------------------------------------
599 template<typename T>
nCols()600 inline INDEX SparseMatrix<T>::nCols() const
601 {
602 return _nCols;
603 }
604 //-----------------------------------------------------------------------------
605 // returns the number of non-zeros in the matrix
606 //-----------------------------------------------------------------------------
607 template<typename T>
size()608 INDEX SparseMatrix<T>::size() const
609 {
610 compress(*this);
611 return _size;
612 }
613 //-----------------------------------------------------------------------------
614 // returns the number of nonzero elements in a row
615 //-----------------------------------------------------------------------------
616 template<typename T>
RowSize(INDEX r)617 INDEX SparseMatrix<T>::RowSize(INDEX r) const
618 {
619 compress(*this);
620 GCHK(r>=_nRows, "Rowsize: invalid row");
621 if (r >= _nRowsCRS) return 0;
622 return _ia[r+1]-_ia[r];
623 }
624 //-----------------------------------------------------------------------------
625 // returns a pointer to the data, causes a compress
626 //-----------------------------------------------------------------------------
627 template<typename T>
ptr()628 T* SparseMatrix<T>::ptr() const
629 {
630 compress(*this);
631 return _val;
632 }
633 template<typename T>
rows()634 INDEX* SparseMatrix<T>::rows() const
635 {
636 compress(*this);
637 return _ia;
638 }
639 template<typename T>
cols()640 INDEX* SparseMatrix<T>::cols() const
641 {
642 compress(*this);
643 return _ja;
644 }
645 //-----------------------------------------------------------------------------
646 // returns true if (i,j) falls in the user specified range
647 //-----------------------------------------------------------------------------
648 template<typename T>
in_range(INDEX i,INDEX j)649 bool SparseMatrix<T>::in_range(INDEX i, INDEX j) const
650 {
651 return i < nRows() && j < nCols();
652 }
653 //-----------------------------------------------------------------------------
654 // assigns this sparsematrix from another one - full memory copy
655 //-----------------------------------------------------------------------------
656 template<typename T>
657 SparseMatrix<T>& SparseMatrix<T>::operator=(const SparseMatrix<T> &C)
658 {
659 _delete();
660 _copy(C);
661 return *this;
662 }
663 //-----------------------------------------------------------------------------
664 // assigns existing sparsematrix to a value, preserving structure
665 //-----------------------------------------------------------------------------
666 template<typename T>
667 SparseMatrix<T>& SparseMatrix<T>::operator=(const T v)
668 {
669 this->set_all_elements_to(v);
670 return *this;
671 }
672 //-----------------------------------------------------------------------------
673 // scales this sparse matrix by a constant
674 //-----------------------------------------------------------------------------
675 template<typename T>
set_all_elements_to(const T & a)676 void SparseMatrix<T>::set_all_elements_to(const T &a)
677 {
678 compress(*this);
679 for (INDEX i=0; i<size(); i++) _val[i] = a;
680 }
681 //-----------------------------------------------------------------------------
682 // scales this sparse matrix by a constant
683 //-----------------------------------------------------------------------------
684 template<typename T>
685 SparseMatrix<T>& SparseMatrix<T>::operator*=(const T &a)
686 {
687 compress(*this);
688 for (INDEX i=0; i<size(); i++) _val[i] *= a;
689 return *this;
690 }
691
692 template<typename T>
693 SparseMatrix<T>& SparseMatrix<T>::operator*=(const SparseMatrix<T> &a)
694 {
695 compress(*this);
696 Matrix<T>::operator*=(a);
697 return *this;
698 }
699
700 //-----------------------------------------------------------------------------
701 // Adds two sparse matrices together.
702 //-----------------------------------------------------------------------------
703 template<typename T>
704 SparseMatrix<T>& SparseMatrix<T>::operator+=(const SparseMatrix & R)
705 {
706
707 compress(R);
708
709 int *Ria = R.rows();
710 int *Rja = R.cols();
711 T *Rval = R.ptr();
712
713 int nRowsCRS = R.nRowsCRS();
714
715 int rowR, colR;
716 T valR;
717 for (rowR = 0; rowR < nRowsCRS; ++rowR) {
718
719 for (int j = Ria[rowR]; j < Ria[rowR+1]; ++j) {
720 colR = Rja[j];
721 valR = Rval[j];
722
723 // Because we simply want to add the value, we call add and let compress
724 // take care of the rest--we don't have to worry about extant entries.
725 add(rowR, colR, valR);
726 }
727 }
728 return *this;
729 }
730
731 //-----------------------------------------------------------------------------
732 // Return matrix transpose
733 //-----------------------------------------------------------------------------
734 template<typename T>
transpose()735 SparseMatrix<T> SparseMatrix<T>::transpose() const
736 {
737 compress(*this);
738 SparseMatrix<T> At(nCols(), nRows());
739
740 for (INDEX i=0; i<_nRowsCRS; i++)
741 for (INDEX ij=_ia[i]; ij<_ia[i+1]; ij++)
742 At.set(_ja[ij], i, _val[ij]);
743 compress(At);
744 return At;
745 }
746
747 //-----------------------------------------------------------------------------
748 // multiplies each row by the corresponding element in Vector scale
749 //-----------------------------------------------------------------------------
750 template<typename T>
row_scale(const Vector<T> & v)751 SparseMatrix<T>& SparseMatrix<T>::row_scale(const Vector<T> &v)
752 {
753 compress(*this);
754 INDEX i,ij;
755 GCK(*this, v, v.size()!=nRows(), "Incompatible Vector length in row_scale.");
756 for(i=0; i<_nRowsCRS; i++)
757 for(ij=_ia[i]; ij<_ia[i+1]; ij++) _val[ij] *= v[i];
758 return *this;
759 }
760 //-----------------------------------------------------------------------------
761 // multiples each column by the corresponding element in Vector scale
762 //-----------------------------------------------------------------------------
763 template<typename T>
col_scale(const Vector<T> & v)764 SparseMatrix<T>& SparseMatrix<T>::col_scale(const Vector<T> &v)
765 {
766 compress(*this);
767 INDEX i,ij;
768 GCK(*this, v, v.size()!=nCols(), "Incompatible Vector length in col_scale.");
769 for(i=0; i<_nRowsCRS; i++)
770 for(ij=_ia[i]; ij<_ia[i+1]; ij++) _val[ij] *= v[_ja[ij]];
771 return *this;
772 }
773 //-----------------------------------------------------------------------------
774 // Returns a vector of the sums of each column
775 //-----------------------------------------------------------------------------
776 template<typename T>
col_sum()777 DenseVector<T> SparseMatrix<T>::col_sum() const
778 {
779 compress(*this);
780 INDEX i,ij;
781 GCHK(!nRows(), "SparseMatrix::Matrix not initialized in col_sum.")
782 DenseVector<T> csum(nCols());
783 for(i=0; i<_nRowsCRS; i++)
784 for(ij=_ia[i]; ij<_ia[i+1]; ij++) csum(_ja[ij]) += _val[ij];
785 return(csum);
786 }
787 //-----------------------------------------------------------------------------
788 // Returns a vector with the number of nonzeros in each column
789 //-----------------------------------------------------------------------------
790 template<typename T>
column_count()791 DenseVector<INDEX> SparseMatrix<T>::column_count() const
792 {
793 compress(*this);
794 INDEX i,j;
795 DenseVector<INDEX> counts(nCols());
796
797 for (i=0; i<_nRowsCRS; i++)
798 for(j=_ia[i]; j<_ia[i+1]; j++) counts(_ja[j])++;
799 return(counts);
800 }
801 //-----------------------------------------------------------------------------
802 // Writes a the nonzeros of a row to a vector
803 //-----------------------------------------------------------------------------
804 template<typename T>
row(INDEX i,DenseVector<T> & row,DenseVector<INDEX> & indx)805 void SparseMatrix<T>::row(INDEX i, DenseVector<T>& row, DenseVector<INDEX>& indx) const
806 {
807 compress(*this);
808 GCHK(i>=nRows(), "get_row() - invalid row number");
809 if (i >= _nRowsCRS) {
810 row.resize(0);
811 indx.resize(0);
812 return;
813 }
814 row.resize(RowSize(i));
815 indx.resize(row.size());
816 INDEX idx=0, ij;
817 for(ij=_ia[i]; ij<_ia[i+1]; ij++)
818 {
819 row(idx) = _val[ij];
820 indx(idx++) = _ja[ij];
821 }
822 }
823 //-----------------------------------------------------------------------------
824 // Computes the product of N'DN
825 //-----------------------------------------------------------------------------
826 template<typename T>
827 void SparseMatrix<T>::
weighted_least_squares(const SparseMatrix<T> & N,const DiagonalMatrix<T> & D)828 weighted_least_squares(const SparseMatrix<T> &N, const DiagonalMatrix<T> &D)
829 {
830 compress(N);
831 GCK(N,D,N.nRows()!=D.nRows(),"SparseMatrix::WeightedLeastSquares()");
832 INDEX k, ki, kj;
833
834 resize(N.nCols(), N.nCols()); // set size of this matrix
835 for (k=0; k<_size; k++) _val[k] = 0.0;
836 // compute R(i,j) = N(k,i) D(k,q) N(i,j) = N(k,i)*D(k,k)*N(k,j) (sum on k)
837 for (k=0; k<N._nRowsCRS; k++)
838 for (ki=N._ia[k]; ki<N._ia[k+1]; ki++)
839 for (kj=N._ia[k]; kj<N._ia[k+1]; kj++)
840 add(N._ja[ki],N._ja[kj], D[k]*N[kj]*N[ki]);
841 compress();
842 }
843 //-----------------------------------------------------------------------------
844 // Return a diagonal matrix containing the diagonal entries of this matrix
845 //-----------------------------------------------------------------------------
846 template<typename T>
diag()847 DiagonalMatrix<T> SparseMatrix<T>::diag() const
848 {
849 compress(*this);
850 DiagonalMatrix<T> D(nRows(), true); // initialized to zero
851 INDEX i, ij;
852 for (i=0; i<_nRowsCRS; i++)
853 {
854 for(ij=_ia[i]; ij<_ia[i+1]; ij++)
855 {
856 if (_ja[ij]>=i) // have we reached or passed the diagonal?
857 {
858 if (_ja[ij]==i) D[i]=_val[ij]; // this this the diagonal?
859 break; // D[i] is already zero if there is no diagonal
860 }
861 }
862 }
863 return D;
864 }
865 //-----------------------------------------------------------------------------
866 // Return a diagonal matrix containing row-sum lumped entries of the matrix
867 //-----------------------------------------------------------------------------
868 template<typename T>
row_sum_lump()869 DiagonalMatrix<T> SparseMatrix<T>::row_sum_lump() const
870 {
871 compress(*this);
872 DiagonalMatrix<T> D(nRows(), true); // initialized to zero
873 INDEX i, ij;
874 for (i=0; i<_nRowsCRS; i++)
875 {
876 for(ij=_ia[i]; ij<_ia[i+1]; ij++)
877 {
878 D(i,i) += _val[ij];
879 }
880 }
881 return D;
882 }
883 //-----------------------------------------------------------------------------
884 // output function - builds a string with each nonzero triplet value
885 //-----------------------------------------------------------------------------
886 template<typename T>
to_string()887 std::string SparseMatrix<T>::to_string() const
888 {
889 compress(*this);
890 std::string out;
891 INDEX i, ij;
892 for(i=0; i<_nRowsCRS; i++)
893 {
894 for(ij=_ia[i]; ij<_ia[i+1]; ij++)
895 {
896 if (ij) out += "\n"; // append newline if not first nonzero
897 out += "(" + ATC_Utility::to_string(i) + ", "; // append "(i,"
898 out += ATC_Utility::to_string(_ja[ij]) + ") = "; // append "j) = "
899 out += ATC_Utility::to_string(_val[ij]); // append "value"
900 }
901 }
902 return out; // return the completed string
903 }
904 //-----------------------------------------------------------------------------
905 // returns the maximum value in the row
906 //-----------------------------------------------------------------------------
907 template<typename T>
row_max(INDEX row)908 T SparseMatrix<T>::row_max(INDEX row) const
909 {
910 compress(*this);
911 if (!RowSize(row)) return (T)0; // if there are no nonzeros in the row
912 INDEX ij;
913 T max = _val[_ia[row]];
914 for(ij=_ia[row]+1; ij<_ia[row+1]; ij++) max = std::max(max,_val[ij]);
915 return max;
916 }
917 //-----------------------------------------------------------------------------
918 // returns the minimum value in the row
919 //-----------------------------------------------------------------------------
920 template<typename T>
row_min(INDEX row)921 T SparseMatrix<T>::row_min(INDEX row) const
922 {
923 compress(*this);
924 if (!RowSize(row)) return (T)0; // if there are no nonzeros in the row
925 INDEX ij;
926 T min = _val[_ia[row]];
927 for(ij=_ia[row]+1; ij<_ia[row+1]; ij++) min = std::min(min,_val[ij]);
928 return min;
929 }
930 //-----------------------------------------------------------------------------
931 // prints a histogram of the values of a row to the screen
932 //-----------------------------------------------------------------------------
933 template<typename T>
print_row_histogram(const std::string & name,INDEX nbins)934 void SparseMatrix<T>::print_row_histogram(const std::string &name, INDEX nbins) const
935 {
936 compress(*this);
937 std::cout << "Begin histogram " << name << "\n";
938 std::cout << "# rows: " << _nRows << " columns: " << _nCols
939 << " size: " << _size << "\n";
940 for(INDEX i=0; i<_nRows; i++)
941 {
942 print_row_histogram(i, nbins);
943 std::cout << "\n";
944 }
945 std::cout << "End histogram " << name << "\n";
946 }
947 //-----------------------------------------------------------------------------
948 // prints a histogram of the values of a row to the screen
949 //-----------------------------------------------------------------------------
950 template<typename T>
print_row_histogram(INDEX row,INDEX nbins)951 void SparseMatrix<T>::print_row_histogram(INDEX row, INDEX nbins) const
952 {
953 compress(*this);
954 if (!nbins) nbins++;
955 std::vector<INDEX> counts(nbins, 0);
956 const T min = row_min(row);
957 const T max = row_max(row);
958 const T range = max-min;
959 const double bin_size = range/double(nbins);
960 if (range<=0.0) counts[nbins-1]=RowSize(row);
961 else
962 {
963 for(INDEX ij=_ia[row]; ij<_ia[row+1]; ij++)
964 {
965 INDEX bin = INDEX((_val[ij]-min)/bin_size);
966 counts[bin-(bin==nbins)]++;
967 }
968 }
969 std::cout<<std::showbase<<std::scientific;
970 std::cout<<"# Histogram: row "<<row<<" min "<<min<<" max "<<max<<" cnt " <<RowSize(row)<<"\n";
971 T bin_start = min;
972 for(INDEX i=0; i<nbins; i++)
973 {
974 std::cout << "(" << bin_start << ",";
975 bin_start += bin_size;
976 std::cout << bin_start << ") " << counts[i] << "\n";
977 }
978 }
979 //-----------------------------------------------------------------------------
980 // prints the triplets the screen
981 //-----------------------------------------------------------------------------
982 template<typename T>
print_triplets()983 void SparseMatrix<T>::print_triplets() const
984 {
985 typename std::vector<TRI_COORD<T> >::const_iterator t;
986 std::string out;
987 out += "==================BEGIN TRIPLETS=======================\n";
988 // Loop backwards over all new triplets.
989 for (t = _tri.begin(); t!=_tri.end(); ++t) {
990 out += "(" + ATC_Utility::to_string(t->i) + ", "; // append "(i,"
991 out += ATC_Utility::to_string(t->j) + ") = "; // append "j) = "
992 out += ATC_Utility::to_string(t->v); // append "value"
993 out += "\n";
994 }
995 out += "===================END TRIPLETS========================\n";
996 std::cout << out;
997 }
998 //-----------------------------------------------------------------------------
999 // Outputs a string to a sparse Matlab type
1000 //-----------------------------------------------------------------------------
1001 template<typename T>
matlab(std::ostream & o,const std::string & s)1002 void SparseMatrix<T>::matlab(std::ostream &o, const std::string &s) const
1003 {
1004 compress(*this);
1005 INDEX i, ij;
1006 o << s <<" = sparse(" << nRows() << "," << nCols() << ");\n";
1007 o << std::showbase << std::scientific;
1008 for(i=0; i<_nRowsCRS; i++)
1009 for(ij=_ia[i]; ij<_ia[i+1]; ij++)
1010 o<<s<<"("<<i+1<<","<<_ja[ij]+1<<")="<<_val[ij]<<";\n";
1011 }
1012
1013 //-----------------------------------------------------------------------------
1014 // Writes the matrix to a binary file (after a compress).
1015 //-----------------------------------------------------------------------------
1016 template<typename T>
binary_write(std::fstream & f)1017 void SparseMatrix<T>::binary_write(std::fstream& f) const
1018 {
1019 compress(*this);
1020 f.write((char*)&_size, sizeof(INDEX)); // writes number of nonzeros
1021 f.write((char*)&_nRowsCRS, sizeof(INDEX)); // writes number of rows in crs
1022 f.write((char*)&_nRows, sizeof(INDEX)); // write matrix rows
1023 f.write((char*)&_nCols, sizeof(INDEX)); // write number of columns
1024 if (!_size) return;
1025 f.write((char*)_val, sizeof(T) *_size);
1026 f.write((char*)_ja, sizeof(INDEX)*_size);
1027 f.write((char*)_ia, sizeof(INDEX)*(_nRowsCRS+1));
1028 }
1029 //-----------------------------------------------------------------------------
1030 // Reads a SparseMatrix from a binary file. (wipes out any original data)
1031 //-----------------------------------------------------------------------------
1032 template<typename T>
binary_read(std::fstream & f)1033 void SparseMatrix<T>::binary_read(std::fstream& f)
1034 {
1035 _delete();
1036 f.read((char*)&_size, sizeof(INDEX));
1037 f.read((char*)&_nRowsCRS, sizeof(INDEX));
1038 f.read((char*)&_nRows, sizeof(INDEX));
1039 f.read((char*)&_nCols, sizeof(INDEX));
1040 if (!_size) return;
1041 _create(_size,_nRowsCRS);
1042 f.read((char*)_val, sizeof(T)*_size);
1043 f.read((char*)_ja, sizeof(INDEX)*_size);
1044 f.read((char*)_ia, sizeof(INDEX)*(_nRowsCRS+1));
1045 }
1046
1047 //-----------------------------------------------------------------------------
1048 // Writes the sparse matrix to a file in a binary format
1049 //-----------------------------------------------------------------------------
1050 template<typename T>
write_restart(FILE * f)1051 void SparseMatrix<T>::write_restart(FILE *f) const
1052 {
1053 compress(*this);
1054 fwrite(&_size, sizeof(INDEX), 1 ,f); // write number of nonzeros
1055 fwrite(&_nRowsCRS, sizeof(INDEX), 1 ,f); // write number of rows
1056 fwrite(&_nRows, sizeof(INDEX), 1 ,f); // write number of columns
1057 fwrite(&_nCols, sizeof(INDEX), 1 ,f); // write number of columns
1058 if (!_size) return;
1059 fwrite(_val, sizeof(T), _size ,f);
1060 fwrite(_ja, sizeof(T), _size ,f);
1061 fwrite(_ia, sizeof(INDEX), _nRowsCRS+1 ,f);
1062 }
1063
1064 } // end namespace
1065 #endif
1066
1067