escript Revision_
SparseMatrix.h
Go to the documentation of this file.
1
2/*****************************************************************************
3*
4* Copyright (c) 2003-2020 by The University of Queensland
5* http://www.uq.edu.au
6*
7* Primary Business: Queensland, Australia
8* Licensed under the Apache License, version 2.0
9* http://www.apache.org/licenses/LICENSE-2.0
10*
11* Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12* Development 2012-2013 by School of Earth Sciences
13* Development from 2014-2017 by Centre for Geoscience Computing (GeoComp)
14* Development from 2019 by School of Earth and Environmental Sciences
15**
16*****************************************************************************/
17
18
19/****************************************************************************/
20
21/* Paso: SparseMatrix */
22
23/****************************************************************************/
24
25/* Author: lgross@uq.edu.au */
26
27/****************************************************************************/
28
29#ifndef __PASO_SPARSEMATRIX_H__
30#define __PASO_SPARSEMATRIX_H__
31
32#include "Pattern.h"
33
34namespace paso {
35
36template <typename T> struct SparseMatrix;
37template <typename T> using SparseMatrix_ptr = boost::shared_ptr<SparseMatrix<T> >;
38template <typename T> using const_SparseMatrix_ptr = boost::shared_ptr<const SparseMatrix<T> >;
39
40typedef int SparseMatrixType;
41
42// this struct holds a sparse matrix
43template <typename T>
44struct SparseMatrix : boost::enable_shared_from_this<SparseMatrix<T> >
45{
47 dim_t rowBlockSize, dim_t colBlockSize,
48 bool patternIsUnrolled);
49
51
52 void setValues(T value);
53
54 void copyFromMainDiagonal(double* out) const;
55
56 void copyToMainDiagonal(const double* in);
57
58 void copyBlockFromMainDiagonal(double* out) const;
59
60 void copyBlockToMainDiagonal(const double* in);
61
62 void applyBlockMatrix(double* block_diag, index_t* pivot, double* x,
63 const double* b) const;
64
65 void invMain(double* inv_diag, index_t* pivot) const;
66
69 dim_t n_col_sub,
70 const index_t* row_list,
71 const index_t* new_col_index) const;
72
74
76
77 void saveHB_CSC(const char* filename) const;
78
79 void saveMM(const char* filename) const;
80
82 {
83 return pattern->borrowMainDiagonalPointer();
84 }
85
87 {
88 return pattern->borrowColoringPointer();
89 }
90
91 inline dim_t getNumColors() const
92 {
93 return pattern->getNumColors();
94 }
95
96 inline dim_t maxDeg() const
97 {
98 return pattern->maxDeg();
99 }
100
101 inline dim_t getTotalNumRows() const
102 {
103 return numRows * row_block_size;
104 }
105
106 inline dim_t getTotalNumCols() const
107 {
108 return numCols * col_block_size;
109 }
110
111 inline dim_t getNumRows() const
112 {
113 return numRows;
114 }
115
116 inline dim_t getNumCols() const
117 {
118 return numCols;
119 }
120
121 inline double getSize() const
122 {
123 return (double)len;
124 }
125
126 inline double getSparsity() const
127 {
128 return getSize() / ((double)getTotalNumRows()*getTotalNumCols());
129 }
130
131 static SparseMatrix_ptr<double> loadMM_toCSR(const char* filename);
132
133
134 void nullifyRowsAndCols_CSC_BLK1(const double* mask_row,
135 const double* mask_col,
136 double main_diagonal_value);
137
138 void nullifyRowsAndCols_CSR_BLK1(const double* mask_row,
139 const double* mask_col,
140 double main_diagonal_value);
141
142 void nullifyRowsAndCols_CSC(const double* mask_row, const double* mask_col,
143 double main_diagonal_value);
144
145 void nullifyRowsAndCols_CSR(const double* mask_row, const double* mask_col,
146 double main_diagonal_value);
147
148 void nullifyRows_CSR_BLK1(const double* mask_row,
149 double main_diagonal_value);
150
151 void nullifyRows_CSR(const double* mask_row, double main_diagonal_value);
152
153 void maxAbsRow_CSR_OFFSET0(double* array) const;
154
155 void addAbsRow_CSR_OFFSET0(double* array) const;
156
157 void addRow_CSR_OFFSET0(double* array) const;
158
159 void applyDiagonal_CSR_OFFSET0(const double* left, const double* right);
160
169
171 T* val;
172
175
177 void* solver_p;
178};
179
180// interfaces:
181
182void SparseMatrix_MatrixVector_CSC_OFFSET0(const double alpha,
184 const double* in,
185 const double beta, double* out);
186
187void SparseMatrix_MatrixVector_CSC_OFFSET1(const double alpha,
189 const double* in,
190 const double beta, double* out);
191
192void SparseMatrix_MatrixVector_CSR_OFFSET0(const double alpha,
194 const double* in,
195 const double beta, double* out);
196
197template <typename T>
198void SparseMatrix_MatrixVector_CSR_OFFSET1(const double alpha,
200 const T* in,
201 const double beta, T* out);
202
203void SparseMatrix_MatrixVector_CSR_OFFSET0_DIAG(const double alpha,
205 const double* in,
206 const double beta, double* out);
207
210
214
215} // namespace paso
216
217#include "Options.h"
218#include "MKL.h"
219#include "MUMPS.h"
220#include "UMFPACK.h"
221
222namespace paso {
223
224struct Preconditioner_LocalSmoother;
225void PASO_DLL_API Preconditioner_LocalSmoother_free(Preconditioner_LocalSmoother * in);
226
227template <>
228void PASO_DLL_API SparseMatrix<double>::saveHB_CSC(const char* filename) const;
229template <>
230void PASO_DLL_API SparseMatrix<double>::saveMM(const char* filename) const;
231
232template <>
233void PASO_DLL_API SparseMatrix<cplx_t>::saveHB_CSC(const char* filename) const;
234template <>
235void PASO_DLL_API SparseMatrix<cplx_t>::saveMM(const char* filename) const;
236
237/* Allocates a SparseMatrix of given type using the given matrix pattern.
238 Values are initialized with zero.
239 If patternIsUnrolled and type & MATRIX_FORMAT_BLK1, it is assumed that the
240 pattern is already unrolled to match the requested block size
241 and offsets. Otherwise unrolling and offset adjustment will be performed.
242*/
243template <typename T>
245 dim_t rowBlockSize, dim_t colBlockSize,
246 bool patternIsUnrolled) :
247 type(ntype),
248 val(NULL),
249 solver_package(PASO_PASO),
250 solver_p(NULL)
251{
252 if (patternIsUnrolled) {
253 if ((ntype & MATRIX_FORMAT_OFFSET1) != (npattern->type & MATRIX_FORMAT_OFFSET1)) {
254 throw PasoException("SparseMatrix: requested offset and pattern offset do not match.");
255 }
256 }
257 // do we need to apply unrolling?
258 bool unroll
259 // we don't like non-square blocks
260 = (rowBlockSize != colBlockSize)
261#ifndef ESYS_HAVE_LAPACK
262 // or any block size bigger than 3
263 || (colBlockSize > 3)
264#endif
265 // or if block size one requested and the block size is not 1
266 || ((ntype & MATRIX_FORMAT_BLK1) && (colBlockSize > 1))
267 // or if offsets don't match
268 || ((ntype & MATRIX_FORMAT_OFFSET1) != (npattern->type & MATRIX_FORMAT_OFFSET1));
269
270 SparseMatrixType pattern_format_out = (ntype & MATRIX_FORMAT_OFFSET1)
272
273 // === compressed sparse columns ===
274 if (ntype & MATRIX_FORMAT_CSC) {
275 if (unroll) {
276 if (patternIsUnrolled) {
277 pattern = npattern;
278 } else {
279 pattern = npattern->unrollBlocks(pattern_format_out,
280 colBlockSize, rowBlockSize);
281 }
282 row_block_size = 1;
283 col_block_size = 1;
284 } else {
285 pattern = npattern->unrollBlocks(pattern_format_out, 1, 1);
286 row_block_size = rowBlockSize;
287 col_block_size = colBlockSize;
288 }
289 numRows = pattern->numInput;
290 numCols = pattern->numOutput;
291 } else {
292 // === compressed sparse row ===
293 if (unroll) {
294 if (patternIsUnrolled) {
295 pattern = npattern;
296 } else {
297 pattern = npattern->unrollBlocks(pattern_format_out,
298 rowBlockSize, colBlockSize);
299 }
300 row_block_size = 1;
301 col_block_size = 1;
302 } else {
303 pattern = npattern->unrollBlocks(pattern_format_out, 1, 1);
304 row_block_size = rowBlockSize;
305 col_block_size = colBlockSize;
306 }
307 numRows = pattern->numOutput;
308 numCols = pattern->numInput;
309 }
310 if (ntype & MATRIX_FORMAT_DIAGONAL_BLOCK) {
312 } else {
314 }
315 len = (size_t)(pattern->len)*(size_t)(block_size);
316
317 val = new T[len];
318 setValues(0.);
319}
320
321template <typename T>
323{
324 switch (solver_package) {
325 case PASO_SMOOTHER:
327 break;
328
329 case PASO_MKL:
330 MKL_free(this);
331 break;
332
333 case PASO_UMFPACK:
334 UMFPACK_free(this);
335 break;
336
337 case PASO_MUMPS:
338 MUMPS_free(this);
339 break;
340 }
341 delete[] val;
342}
343
344template <typename T>
346{
347 const index_t index_offset=(type & MATRIX_FORMAT_OFFSET1 ? 1:0);
348 if (!pattern->isEmpty()) {
349 const dim_t nOut = pattern->numOutput;
350#pragma omp parallel for
351 for (dim_t i=0; i < nOut; ++i) {
352 for (index_t iptr=pattern->ptr[i]-index_offset; iptr < pattern->ptr[i+1]-index_offset; ++iptr) {
353 for (dim_t j=0; j<block_size; ++j)
354 val[iptr*block_size+j] = value;
355 }
356 }
357 }
358}
359
360/* CSR format with offset 1 */
361template <typename T>
364 const T* in,
365 double beta, T* out)
366{
367 const int totalRowSize = A->numRows * A->row_block_size;
368 if (std::abs(beta) > 0) {
369 if (beta != 1.) {
370#pragma omp parallel for schedule(static)
371 for (index_t irow=0; irow < totalRowSize; irow++) {
372 out[irow] *= beta;
373 }
374 }
375 } else {
376#pragma omp parallel for schedule(static)
377 for (index_t irow=0; irow < totalRowSize; irow++) {
378 out[irow] = 0.;
379 }
380 }
381
382 if (std::abs(alpha) > 0) {
383 const int nRows = A->pattern->numOutput;
384 if (A->col_block_size==1 && A->row_block_size==1) {
385#pragma omp parallel for schedule(static)
386 for (index_t irow=0; irow < nRows; irow++) {
387 T reg = 0.;
388 #pragma ivdep
389 for (index_t iptr=A->pattern->ptr[irow]-1;
390 iptr < A->pattern->ptr[irow+1]-1; ++iptr) {
391 reg += A->val[iptr] * in[A->pattern->index[iptr]-1];
392 }
393 out[irow] += alpha * reg;
394 }
395 } else if (A->col_block_size==2 && A->row_block_size==2) {
396#pragma omp parallel for schedule(static)
397 for (index_t ir=0; ir < nRows; ir++) {
398 T reg1 = 0.;
399 T reg2 = 0.;
400 #pragma ivdep
401 for (index_t iptr=A->pattern->ptr[ir]-1;
402 iptr < A->pattern->ptr[ir+1]-1; iptr++) {
403 const index_t ic=2*(A->pattern->index[iptr]-1);
404 reg1 += A->val[iptr*4 ]*in[ic] + A->val[iptr*4+2]*in[1+ic];
405 reg2 += A->val[iptr*4+1]*in[ic] + A->val[iptr*4+3]*in[1+ic];
406 }
407 out[ 2*ir] += alpha * reg1;
408 out[1+2*ir] += alpha * reg2;
409 }
410 } else if (A->col_block_size==3 && A->row_block_size==3) {
411#pragma omp parallel for schedule(static)
412 for (index_t ir=0; ir < nRows; ir++) {
413 T reg1 = 0.;
414 T reg2 = 0.;
415 T reg3 = 0.;
416 #pragma ivdep
417 for (index_t iptr=A->pattern->ptr[ir]-1;
418 iptr < A->pattern->ptr[ir+1]-1; iptr++) {
419 const index_t ic=3*(A->pattern->index[iptr]-1);
420 reg1 += A->val[iptr*9 ]*in[ic] + A->val[iptr*9+3]*in[1+ic] + A->val[iptr*9+6]*in[2+ic];
421 reg2 += A->val[iptr*9+1]*in[ic] + A->val[iptr*9+4]*in[1+ic] + A->val[iptr*9+7]*in[2+ic];
422 reg3 += A->val[iptr*9+2]*in[ic] + A->val[iptr*9+5]*in[1+ic] + A->val[iptr*9+8]*in[2+ic];
423 }
424 out[ 3*ir] += alpha * reg1;
425 out[1+3*ir] += alpha * reg2;
426 out[2+3*ir] += alpha * reg3;
427 }
428 } else {
429#pragma omp parallel for schedule(static)
430 for (index_t ir=0; ir < nRows; ir++) {
431 for (index_t iptr=A->pattern->ptr[ir]-1;
432 iptr < A->pattern->ptr[ir+1]-1; iptr++) {
433 for (index_t irb=0; irb < A->row_block_size; irb++) {
434 T reg = 0.;
435 #pragma ivdep
436 for (index_t icb=0; icb < A->col_block_size; icb++) {
437 const index_t icol=icb+A->col_block_size*(A->pattern->index[iptr]-1);
438 reg += A->val[iptr*A->block_size+irb+A->row_block_size*icb] * in[icol];
439 }
440 const index_t irow=irb+A->row_block_size*ir;
441 out[irow] += alpha * reg;
442 }
443 }
444 }
445 } // blocksizes
446 } // alpha > 0
447}
448
449template <typename T>
451 const double* mask_col,
452 double main_diagonal_value)
453{
454 const index_t index_offset=(type & MATRIX_FORMAT_OFFSET1 ? 1:0);
455 const int nOut = pattern->numOutput;
456#pragma omp parallel for
457 for (index_t icol=0; icol < nOut; icol++) {
458 #pragma ivdep
459 for (index_t iptr=pattern->ptr[icol]-index_offset; iptr < pattern->ptr[icol+1]-index_offset; iptr++) {
460 const index_t irow = pattern->index[iptr]-index_offset;
461 if (mask_col[icol]>0. || mask_row[irow]>0.) {
462 val[iptr] = (irow==icol ? main_diagonal_value : 0);
463 }
464 }
465 }
466}
467
468template <typename T>
469void SparseMatrix<T>::nullifyRowsAndCols_CSC(const double* mask_row,
470 const double* mask_col,
471 double main_diagonal_value)
472{
473 const index_t index_offset=(type & MATRIX_FORMAT_OFFSET1 ? 1:0);
474 const int nOut = pattern->numOutput;
475#pragma omp parallel for
476 for (index_t ic=0; ic < nOut; ic++) {
477 for (index_t iptr=pattern->ptr[ic]-index_offset; iptr < pattern->ptr[ic+1]-index_offset; iptr++) {
478 for (index_t irb=0; irb < row_block_size; irb++) {
479 const index_t irow=irb+row_block_size*(pattern->index[iptr]-index_offset);
480 #pragma ivdep
481 for (index_t icb=0; icb < col_block_size; icb++) {
482 const index_t icol=icb+col_block_size*ic;
483 if (mask_col[icol]>0. || mask_row[irow]>0.) {
484 const index_t l=iptr*block_size+irb+row_block_size*icb;
485 val[l] = (irow==icol ? main_diagonal_value : 0);
486 }
487 }
488 }
489 }
490 }
491}
492
493template <typename T>
495 const double* mask_col,
496 double main_diagonal_value)
497{
498 const index_t index_offset=(type & MATRIX_FORMAT_OFFSET1 ? 1:0);
499 const int nOut = pattern->numOutput;
500#pragma omp parallel for
501 for (index_t irow=0; irow < nOut; irow++) {
502 #pragma ivdep
503 for (index_t iptr=pattern->ptr[irow]-index_offset; iptr < pattern->ptr[irow+1]-index_offset; iptr++) {
504 const index_t icol = pattern->index[iptr]-index_offset;
505 if (mask_col[icol]>0. || mask_row[irow]>0.) {
506 val[iptr] = (irow==icol ? main_diagonal_value : 0);
507 }
508 }
509 }
510}
511
512template <typename T>
513void SparseMatrix<T>::nullifyRowsAndCols_CSR(const double* mask_row,
514 const double* mask_col,
515 double main_diagonal_value)
516{
517 const index_t index_offset=(type & MATRIX_FORMAT_OFFSET1 ? 1:0);
518 const int nOut = pattern->numOutput;
519#pragma omp parallel for
520 for (index_t ir=0; ir < nOut; ir++) {
521 for (index_t iptr=pattern->ptr[ir]-index_offset; iptr < pattern->ptr[ir+1]-index_offset; iptr++) {
522 for (index_t irb=0; irb < row_block_size; irb++) {
523 const index_t irow=irb+row_block_size*ir;
524 #pragma ivdep
525 for (index_t icb=0; icb < col_block_size; icb++) {
526 const index_t icol=icb+col_block_size*(pattern->index[iptr]-index_offset);
527 if (mask_col[icol]>0. || mask_row[irow]>0.) {
528 const index_t l=iptr*block_size+irb+row_block_size*icb;
529 val[l] = (irow==icol ? main_diagonal_value : 0);
530 }
531 }
532 }
533 }
534 }
535}
536
537} // namespace paso
538
539#endif // __PASO_SPARSEMATRIX_H__
540
#define PASO_SMOOTHER
Definition: Options.h:75
#define PASO_MUMPS
Definition: Options.h:57
#define PASO_UMFPACK
Definition: Options.h:51
#define PASO_PASO
Definition: Options.h:56
#define PASO_MKL
Definition: Options.h:50
#define MATRIX_FORMAT_BLK1
Definition: Paso.h:63
#define MATRIX_FORMAT_DEFAULT
Definition: Paso.h:61
#define MATRIX_FORMAT_OFFSET1
Definition: Paso.h:64
#define MATRIX_FORMAT_CSC
Definition: Paso.h:62
#define MATRIX_FORMAT_DIAGONAL_BLOCK
Definition: Paso.h:65
PasoException exception class.
Definition: PasoException.h:34
index_t dim_t
Definition: DataTypes.h:66
int index_t
type for array/matrix indices used both globally and on each rank
Definition: DataTypes.h:61
Definition: BiCGStab.cpp:25
void SparseMatrix_MatrixVector_CSR_OFFSET0(const double alpha, const_SparseMatrix_ptr< double > A, const double *in, const double beta, double *out)
Definition: SparseMatrix_MatrixVector.cpp:191
boost::shared_ptr< SparseMatrix< T > > SparseMatrix_ptr
Definition: SparseMatrix.h:37
boost::shared_ptr< Pattern > Pattern_ptr
Definition: Pattern.h:40
void SparseMatrix_MatrixVector_CSR_OFFSET0_DIAG(const double alpha, const_SparseMatrix_ptr< double > A, const double *in, const double beta, double *out)
Definition: SparseMatrix_MatrixVector.cpp:338
void SparseMatrix_MatrixVector_CSR_OFFSET1(const double alpha, const_SparseMatrix_ptr< T > A, const T *in, const double beta, T *out)
Definition: SparseMatrix.h:362
void SparseMatrix_MatrixVector_CSC_OFFSET0(const double alpha, const_SparseMatrix_ptr< double > A, const double *in, const double beta, double *out)
Definition: SparseMatrix_MatrixVector.cpp:43
void UMFPACK_free(SparseMatrix< double > *A)
frees any UMFPACK related data from the matrix
Definition: UMFPACK.cpp:35
SparseMatrix_ptr< double > SparseMatrix_MatrixMatrixTranspose(const_SparseMatrix_ptr< double > A, const_SparseMatrix_ptr< double > B, const_SparseMatrix_ptr< double > T)
Definition: SparseMatrix_MatrixMatrixTranspose.cpp:52
void SparseMatrix_MatrixVector_CSC_OFFSET1(const double alpha, const_SparseMatrix_ptr< double > A, const double *in, const double beta, double *out)
Definition: SparseMatrix_MatrixVector.cpp:119
boost::shared_ptr< const SparseMatrix< T > > const_SparseMatrix_ptr
Definition: SparseMatrix.h:38
void Preconditioner_LocalSmoother_free(Preconditioner_LocalSmoother *in)
Definition: Smoother.cpp:42
void MUMPS_free(SparseMatrix< T > *A)
frees any MUMPS related data from the matrix
Definition: MUMPS.h:118
void MKL_free(SparseMatrix< double > *A)
Definition: MKL.cpp:35
SparseMatrix_ptr< double > SparseMatrix_MatrixMatrix(const_SparseMatrix_ptr< double > A, const_SparseMatrix_ptr< double > B)
Definition: SparseMatrix_MatrixMatrix.cpp:44
int SparseMatrixType
Definition: SparseMatrix.h:40
#define PASO_DLL_API
Definition: paso/src/system_dep.h:29
Definition: Preconditioner.h:57
Definition: SparseMatrix.h:45
index_t solver_package
package controlling the solver pointer
Definition: SparseMatrix.h:174
void invMain(double *inv_diag, index_t *pivot) const
SparseMatrixType type
Definition: SparseMatrix.h:161
void copyBlockToMainDiagonal(const double *in)
void copyFromMainDiagonal(double *out) const
dim_t len
Definition: SparseMatrix.h:168
double getSparsity() const
Definition: SparseMatrix.h:126
dim_t block_size
Definition: SparseMatrix.h:164
dim_t getTotalNumCols() const
Definition: SparseMatrix.h:106
SparseMatrix_ptr< double > unroll(SparseMatrixType type) const
void nullifyRowsAndCols_CSC_BLK1(const double *mask_row, const double *mask_col, double main_diagonal_value)
Definition: SparseMatrix.h:450
Pattern_ptr pattern
Definition: SparseMatrix.h:167
dim_t numRows
Definition: SparseMatrix.h:165
T * val
this is used for classical CSR or CSC
Definition: SparseMatrix.h:171
dim_t getTotalNumRows() const
Definition: SparseMatrix.h:101
void applyBlockMatrix(double *block_diag, index_t *pivot, double *x, const double *b) const
void saveMM(const char *filename) const
void applyDiagonal_CSR_OFFSET0(const double *left, const double *right)
dim_t numCols
Definition: SparseMatrix.h:166
void addRow_CSR_OFFSET0(double *array) const
void saveHB_CSC(const char *filename) const
SparseMatrix_ptr< double > getSubmatrix(dim_t n_row_sub, dim_t n_col_sub, const index_t *row_list, const index_t *new_col_index) const
static SparseMatrix_ptr< double > loadMM_toCSR(const char *filename)
void copyBlockFromMainDiagonal(double *out) const
dim_t maxDeg() const
Definition: SparseMatrix.h:96
void nullifyRows_CSR(const double *mask_row, double main_diagonal_value)
void * solver_p
pointer to data needed by a solver
Definition: SparseMatrix.h:177
SparseMatrix_ptr< double > getTranspose() const
void nullifyRowsAndCols_CSR_BLK1(const double *mask_row, const double *mask_col, double main_diagonal_value)
Definition: SparseMatrix.h:494
SparseMatrix(SparseMatrixType type, Pattern_ptr pattern, dim_t rowBlockSize, dim_t colBlockSize, bool patternIsUnrolled)
Definition: SparseMatrix.h:244
dim_t col_block_size
Definition: SparseMatrix.h:163
index_t * borrowMainDiagonalPointer() const
Definition: SparseMatrix.h:81
void maxAbsRow_CSR_OFFSET0(double *array) const
void nullifyRows_CSR_BLK1(const double *mask_row, double main_diagonal_value)
~SparseMatrix()
Definition: SparseMatrix.h:322
dim_t getNumRows() const
Definition: SparseMatrix.h:111
dim_t row_block_size
Definition: SparseMatrix.h:162
index_t * borrowColoringPointer() const
Definition: SparseMatrix.h:86
void nullifyRowsAndCols_CSC(const double *mask_row, const double *mask_col, double main_diagonal_value)
Definition: SparseMatrix.h:469
void nullifyRowsAndCols_CSR(const double *mask_row, const double *mask_col, double main_diagonal_value)
Definition: SparseMatrix.h:513
void setValues(T value)
Definition: SparseMatrix.h:345
double getSize() const
Definition: SparseMatrix.h:121
void copyToMainDiagonal(const double *in)
void addAbsRow_CSR_OFFSET0(double *array) const
dim_t getNumColors() const
Definition: SparseMatrix.h:91
dim_t getNumCols() const
Definition: SparseMatrix.h:116
SparseMatrix_ptr< double > getBlock(int blockid) const