MFEM  v3.4
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
blockmatrix.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "../general/array.hpp"
13 #include "../general/globals.hpp"
14 #include "matrix.hpp"
15 #include "sparsemat.hpp"
16 #include "blockvector.hpp"
17 #include "blockmatrix.hpp"
18 
19 namespace mfem
20 {
21 
23  AbstractSparseMatrix(offsets.Last()),
24  owns_blocks(false),
25  nRowBlocks(offsets.Size()-1),
26  nColBlocks(offsets.Size()-1),
27  row_offsets(const_cast< Array<int>& >(offsets).GetData(), offsets.Size()),
28  col_offsets(const_cast< Array<int>& >(offsets).GetData(), offsets.Size()),
29  Aij(nRowBlocks, nColBlocks)
30 {
31  Aij = (SparseMatrix *)NULL;
32 }
33 
34 BlockMatrix::BlockMatrix(const Array<int> & row_offsets_,
35  const Array<int> & col_offsets_):
36  AbstractSparseMatrix(row_offsets_.Last(), col_offsets_.Last()),
37  owns_blocks(false),
38  nRowBlocks(row_offsets_.Size()-1),
39  nColBlocks(col_offsets_.Size()-1),
40  row_offsets(const_cast< Array<int>& >(row_offsets_).GetData(),
41  row_offsets_.Size()),
42  col_offsets(const_cast< Array<int>& >(col_offsets_).GetData(),
43  col_offsets_.Size()),
44  Aij(nRowBlocks, nColBlocks)
45 {
46  Aij = (SparseMatrix *)NULL;
47 }
48 
49 
51 {
52  if (owns_blocks)
53  for (SparseMatrix ** it = Aij.GetRow(0);
54  it != Aij.GetRow(0)+(Aij.NumRows()*Aij.NumCols()); ++it)
55  {
56  delete *it;
57  }
58 }
59 
60 void BlockMatrix::SetBlock(int i, int j, SparseMatrix * mat)
61 {
62 #ifdef MFEM_DEBUG
63  if (nRowBlocks <= i || nColBlocks <= j)
64  {
65  mfem_error("BlockMatrix::SetBlock #0");
66  }
67 
68  if (mat->Height() != row_offsets[i+1] - row_offsets[i])
69  {
70  mfem_error("BlockMatrix::SetBlock #1");
71  }
72 
73  if (mat->Width() != col_offsets[j+1] - col_offsets[j])
74  {
75  mfem_error("BlockMatrix::SetBlock #2");
76  }
77 #endif
78  Aij(i,j) = mat;
79 }
80 
82 {
83 #ifdef MFEM_DEBUG
84  if (nRowBlocks <= i || nColBlocks <= j)
85  {
86  mfem_error("BlockMatrix::Block #0");
87  }
88 
89  if (IsZeroBlock(i,j))
90  {
91  mfem_error("BlockMatrix::Block #1");
92  }
93 #endif
94  return *Aij(i,j);
95 }
96 
97 const SparseMatrix & BlockMatrix::GetBlock(int i, int j) const
98 {
99 #ifdef MFEM_DEBUG
100  if (nRowBlocks <= i || nColBlocks <= j)
101  {
102  mfem_error("BlockMatrix::Block const #0");
103  }
104 
105  if (IsZeroBlock(i,j))
106  {
107  mfem_error("BlockMatrix::Block const #1");
108  }
109 #endif
110 
111  return *Aij(i,j);
112 }
113 
114 
116 {
117  int nnz_elem = 0;
118  for (int jcol = 0; jcol != nColBlocks; ++jcol)
119  for (int irow = 0; irow != nRowBlocks; ++irow)
120  {
121  if (Aij(irow,jcol))
122  {
123  nnz_elem+= Aij(irow,jcol)->NumNonZeroElems();
124  }
125  }
126  return nnz_elem;
127 }
128 
129 
130 double& BlockMatrix::Elem (int i, int j)
131 {
132  int iloc, iblock;
133  int jloc, jblock;
134 
135  findGlobalRow(i, iblock, iloc);
136  findGlobalCol(j, jblock, jloc);
137 
138  if (IsZeroBlock(i, j))
139  {
140  mfem_error("BlockMatrix::Elem");
141  }
142 
143  return Aij(iblock, jblock)->Elem(iloc, jloc);
144 }
145 
146 const double& BlockMatrix::Elem (int i, int j) const
147 {
148  int iloc, iblock;
149  int jloc, jblock;
150 
151  findGlobalRow(i, iblock, iloc);
152  findGlobalCol(j, jblock, jloc);
153 
154  if (IsZeroBlock(i, j))
155  {
156  mfem_error("BlockMatrix::Elem");
157  }
158 
159  return Aij(iblock, jblock)->Elem(iloc, jloc);
160 }
161 
162 int BlockMatrix::RowSize(const int i) const
163 {
164  int rowsize = 0;
165 
166  int iblock, iloc;
167  findGlobalRow(i, iblock, iloc);
168 
169  for (int jblock = 0; jblock < nColBlocks; ++jblock)
170  if (Aij(iblock,jblock) != NULL)
171  {
172  rowsize += Aij(iblock,jblock)->RowSize(iloc);
173  }
174 
175  return rowsize;
176 }
177 
178 int BlockMatrix::GetRow(const int row, Array<int> &cols, Vector &srow) const
179 {
180  int iblock, iloc, rowsize;
181  findGlobalRow(row, iblock, iloc);
182  rowsize = RowSize(row);
183  cols.SetSize(rowsize);
184  srow.SetSize(rowsize);
185 
186  Array<int> bcols;
187  Vector bsrow;
188 
189  int * it_cols = cols.GetData();
190  double *it_srow = srow.GetData();
191 
192  for (int jblock = 0; jblock < nColBlocks; ++jblock)
193  if (Aij(iblock,jblock) != NULL)
194  {
195  Aij(iblock,jblock)->GetRow(iloc, bcols, bsrow);
196  for (int i = 0; i < bcols.Size(); ++i)
197  {
198  *(it_cols++) = bcols[i] + col_offsets[jblock];
199  *(it_srow++) = bsrow(i);
200  }
201  }
202 
203  return 0;
204 }
205 
207 {
208  // Find the block to which the dof belongs and its local number
209  int idx, iiblock;
210  for (iiblock = 0; iiblock < nRowBlocks; ++iiblock)
211  {
212  idx = rc - row_offsets[iiblock];
213  if (idx < 0 ) { break; }
214  }
215  iiblock--;
216  idx = rc - row_offsets[iiblock];
217 
218  // Asserts
219  MFEM_ASSERT(nRowBlocks == nColBlocks,
220  "BlockMatrix::EliminateRowCol: nRowBlocks != nColBlocks");
221 
222  MFEM_ASSERT(row_offsets[iiblock] == col_offsets[iiblock],
223  "BlockMatrix::EliminateRowCol: row_offests["
224  << iiblock << "] != col_offsets["<<iiblock<<"]");
225 
226  MFEM_ASSERT(Aij(iiblock, iiblock),
227  "BlockMatrix::EliminateRowCol: Null diagonal block");
228 
229  // Apply the constraint idx to the iiblock
230  for (int jjblock = 0; jjblock < nRowBlocks; ++jjblock)
231  {
232  if (iiblock == jjblock) { continue; }
233  if (Aij(iiblock,jjblock)) { Aij(iiblock,jjblock)->EliminateRow(idx); }
234  }
235  for (int jjblock = 0; jjblock < nRowBlocks; ++jjblock)
236  {
237  if (iiblock == jjblock) { continue; }
238  if (Aij(jjblock,iiblock)) { Aij(jjblock,iiblock)->EliminateCol(idx); }
239  }
240  Aij(iiblock, iiblock)->EliminateRowCol(idx,dpolicy);
241 }
242 
244  Vector & rhs)
245 {
246  if (nRowBlocks != nColBlocks)
247  {
248  mfem_error("BlockMatrix::EliminateRowCol: nRowBlocks != nColBlocks");
249  }
250 
251  for (int iiblock = 0; iiblock < nRowBlocks; ++iiblock)
252  if (row_offsets[iiblock] != col_offsets[iiblock])
253  {
254  mfem::out << "BlockMatrix::EliminateRowCol: row_offests["
255  << iiblock << "] != col_offsets["<<iiblock<<"]\n";
256  mfem_error();
257  }
258 
259  // We also have to do the same for each Aij
260  Array<int> block_dofs;
261  Vector block_sol, block_rhs;
262 
263  for (int iiblock = 0; iiblock < nRowBlocks; ++iiblock)
264  {
265  int dsize = row_offsets[iiblock+1] - row_offsets[iiblock];
266  block_dofs.MakeRef(ess_bc_dofs.GetData()+row_offsets[iiblock], dsize);
267  block_sol.SetDataAndSize(sol.GetData()+row_offsets[iiblock], dsize);
268  block_rhs.SetDataAndSize(rhs.GetData()+row_offsets[iiblock], dsize);
269 
270  if (Aij(iiblock, iiblock))
271  {
272  for (int i = 0; i < block_dofs.Size(); ++i)
273  if (block_dofs[i])
274  {
275  Aij(iiblock, iiblock)->EliminateRowCol(i,block_sol(i), block_rhs);
276  }
277  }
278  else
279  {
280  for (int i = 0; i < block_dofs.Size(); ++i)
281  if (block_dofs[i])
282  {
283  mfem_error("BlockMatrix::EliminateRowCol: Null diagonal block \n");
284  }
285  }
286 
287  for (int jjblock = 0; jjblock < nRowBlocks; ++jjblock)
288  {
289  if (jjblock != iiblock && Aij(iiblock, jjblock))
290  {
291  for (int i = 0; i < block_dofs.Size(); ++i)
292  if (block_dofs[i])
293  {
294  Aij(iiblock, jjblock)->EliminateRow(i);
295  }
296  }
297  if (jjblock != iiblock && Aij(jjblock, iiblock))
298  {
299  block_rhs.SetDataAndSize(rhs.GetData()+row_offsets[jjblock],
300  row_offsets[jjblock+1] - row_offsets[jjblock]);
301  Aij(jjblock, iiblock)->EliminateCols(block_dofs, &block_sol, &block_rhs);
302  }
303  }
304  }
305 }
306 
307 void BlockMatrix::EliminateZeroRows(const double threshold)
308 {
309  MFEM_VERIFY(nRowBlocks == nColBlocks, "not a square matrix");
310 
311  for (int iblock = 0; iblock < nRowBlocks; ++iblock)
312  {
313  if (Aij(iblock,iblock))
314  {
315  double norm;
316  for (int i = 0; i < Aij(iblock, iblock)->NumRows(); ++i)
317  {
318  norm = 0.;
319  for (int jblock = 0; jblock < nColBlocks; ++jblock)
320  if (Aij(iblock,jblock))
321  {
322  norm += Aij(iblock,jblock)->GetRowNorml1(i);
323  }
324 
325  if (norm <= threshold)
326  {
327  for (int jblock = 0; jblock < nColBlocks; ++jblock)
328  if (Aij(iblock,jblock))
329  {
330  Aij(iblock,jblock)->EliminateRow(
331  i, (iblock==jblock) ? DIAG_ONE : DIAG_ZERO);
332  }
333  }
334  }
335  }
336  else
337  {
338  double norm;
339  for (int i = 0; i < row_offsets[iblock+1] - row_offsets[iblock]; ++i)
340  {
341  norm = 0.;
342  for (int jblock = 0; jblock < nColBlocks; ++jblock)
343  if (Aij(iblock,jblock))
344  {
345  norm += Aij(iblock,jblock)->GetRowNorml1(i);
346  }
347 
348  MFEM_VERIFY(!(norm <= threshold), "diagonal block is NULL:"
349  " iblock = " << iblock << ", i = " << i << ", norm = "
350  << norm);
351  }
352  }
353  }
354 }
355 
356 void BlockMatrix::Finalize(int skip_zeros, bool fix_empty_rows)
357 {
358  for (int iblock = 0; iblock < nRowBlocks; ++iblock)
359  {
360  for (int jblock = 0; jblock < nColBlocks; ++jblock)
361  {
362  if (!Aij(iblock,jblock)) { continue; }
363  if (!Aij(iblock,jblock)->Finalized())
364  {
365  Aij(iblock,jblock)->Finalize(skip_zeros, fix_empty_rows);
366  }
367  }
368  }
369 }
370 
371 void BlockMatrix::Mult(const Vector & x, Vector & y) const
372 {
373  if (x.GetData() == y.GetData())
374  {
375  mfem_error("Error: x and y can't point to the same datas \n");
376  }
377 
378  MFEM_ASSERT(width == x.Size(), "Input vector size (" << x.Size()
379  << ") must match matrix width (" << width << ")");
380  MFEM_ASSERT(height == y.Size(), "Output vector size (" << y.Size()
381  << ") must match matrix height (" << height << ")");
382  y = 0.;
383  AddMult(x, y, 1.0);
384 }
385 
386 void BlockMatrix::AddMult(const Vector & x, Vector & y, const double val) const
387 {
388  if (x.GetData() == y.GetData())
389  {
390  mfem_error("Error: x and y can't point to the same datas \n");
391  }
392 
393  Vector xblockview, yblockview;
394 
395  for (int iblock = 0; iblock != nRowBlocks; ++iblock)
396  {
397  yblockview.SetDataAndSize(y.GetData() + row_offsets[iblock],
398  row_offsets[iblock+1] - row_offsets[iblock]);
399 
400  for (int jblock = 0; jblock != nColBlocks; ++jblock)
401  {
402  if (Aij(iblock, jblock) != NULL)
403  {
404  xblockview.SetDataAndSize(
405  x.GetData() + col_offsets[jblock],
406  col_offsets[jblock+1] - col_offsets[jblock]);
407 
408  Aij(iblock, jblock)->AddMult(xblockview, yblockview, val);
409  }
410  }
411  }
412 }
413 
414 void BlockMatrix::MultTranspose(const Vector & x, Vector & y) const
415 {
416  if (x.GetData() == y.GetData())
417  {
418  mfem_error("Error: x and y can't point to the same datas \n");
419  }
420 
421  y = 0.;
422  AddMultTranspose(x, y, 1.0);
423 }
424 
426  const double val) const
427 {
428  if (x.GetData() == y.GetData())
429  {
430  mfem_error("Error: x and y can't point to the same datas \n");
431  }
432 
433  Vector xblockview, yblockview;
434 
435  for (int iblock = 0; iblock != nColBlocks; ++iblock)
436  {
437  yblockview.SetDataAndSize(y.GetData() + col_offsets[iblock],
438  col_offsets[iblock+1] - col_offsets[iblock]);
439 
440  for (int jblock = 0; jblock != nRowBlocks; ++jblock)
441  {
442  if (Aij(jblock, iblock) != NULL)
443  {
444  xblockview.SetDataAndSize(
445  x.GetData() + row_offsets[jblock],
446  row_offsets[jblock+1] - row_offsets[jblock]);
447 
448  Aij(jblock, iblock)->AddMultTranspose(xblockview, yblockview, val);
449  }
450  }
451  }
452 }
453 
455 {
456  int nnz = NumNonZeroElems();
457 
458  int * i_amono = new int[ row_offsets[nRowBlocks]+2 ];
459  int * j_amono = new int[ nnz ];
460  double * data = new double[ nnz ];
461 
462  for (int i = 0; i < row_offsets[nRowBlocks]+2; i++)
463  {
464  i_amono[i] = 0;
465  }
466 
467  int * i_amono_construction = i_amono+1;
468 
469  for (int iblock = 0; iblock != nRowBlocks; ++iblock)
470  {
471  for (int irow(row_offsets[iblock]); irow < row_offsets[iblock+1]; ++irow)
472  {
473  int local_row = irow - row_offsets[iblock];
474  int ind = i_amono_construction[irow];
475  for (int jblock = 0; jblock < nColBlocks; ++jblock)
476  {
477  if (Aij(iblock,jblock) != NULL)
478  ind += Aij(iblock, jblock)->GetI()[local_row+1]
479  - Aij(iblock, jblock)->GetI()[local_row];
480  }
481  i_amono_construction[irow+1] = ind;
482  }
483  }
484 
485  // Fill in the jarray and copy the data
486  for (int iblock = 0; iblock != nRowBlocks; ++iblock)
487  {
488  for (int jblock = 0; jblock != nColBlocks; ++jblock)
489  {
490  if (Aij(iblock,jblock) != NULL)
491  {
492  int nrow = row_offsets[iblock+1]-row_offsets[iblock];
493  int * i_aij = Aij(iblock, jblock)->GetI();
494  int * j_aij = Aij(iblock, jblock)->GetJ();
495  double * data_aij = Aij(iblock, jblock)->GetData();
496  int *i_it = i_amono_construction+row_offsets[iblock];
497 
498  int loc_start_index = 0;
499  int loc_end_index = 0;
500  int glob_start_index = 0;
501 
502  int shift(col_offsets[jblock]);
503  for (int * i_it_aij(i_aij+1); i_it_aij != i_aij+nrow+1; ++i_it_aij)
504  {
505  glob_start_index = *i_it;
506 
507 #ifdef MFEM_DEBUG
508  if (glob_start_index > nnz)
509  {
510  mfem::out<<"glob_start_index = " << glob_start_index << "\n";
511  mfem::out<<"Block:" << iblock << " " << jblock << "\n";
512  mfem::out<<std::endl;
513  }
514 #endif
515  loc_end_index = *(i_it_aij);
516  for (int cnt = 0; cnt < loc_end_index-loc_start_index; cnt++)
517  {
518  data[glob_start_index+cnt] = data_aij[loc_start_index+cnt];
519  j_amono[glob_start_index+cnt] = j_aij[loc_start_index+cnt] + shift;
520  }
521 
522  *i_it += loc_end_index-loc_start_index;
523  ++i_it;
524  loc_start_index = loc_end_index;
525  }
526  }
527  }
528  }
529 
530  return new SparseMatrix(i_amono, j_amono, data, row_offsets[nRowBlocks],
531  col_offsets[nColBlocks]);
532 }
533 
534 void BlockMatrix::PrintMatlab(std::ostream & os) const
535 {
536 
537  Vector row_data;
538  Array<int> row_ind;
539  int nnz_elem = NumNonZeroElems();
540  os<<"% size " << row_offsets.Last() << " " << col_offsets.Last() << "\n";
541  os<<"% Non Zeros " << nnz_elem << "\n";
542  int i, j;
543  std::ios::fmtflags old_fmt = os.flags();
544  os.setf(std::ios::scientific);
545  std::streamsize old_prec = os.precision(14);
546  for (i = 0; i < row_offsets.Last(); i++)
547  {
548  GetRow(i, row_ind, row_data);
549  for (j = 0; j < row_ind.Size(); j++)
550  {
551  os << i+1 << " " << row_ind[j]+1 << " " << row_data[j] << std::endl;
552  }
553  }
554 
555  os.precision(old_prec);
556  os.flags(old_fmt);
557 }
558 
560 {
561  BlockMatrix * At = new BlockMatrix(A.ColOffsets(), A.RowOffsets());
562  At->owns_blocks = 1;
563 
564  for (int irowAt = 0; irowAt < At->NumRowBlocks(); ++irowAt)
565  for (int jcolAt = 0; jcolAt < At->NumColBlocks(); ++jcolAt)
566  if (!A.IsZeroBlock(jcolAt, irowAt))
567  {
568  At->SetBlock(irowAt, jcolAt, Transpose(A.GetBlock(jcolAt, irowAt)));
569  }
570  return At;
571 }
572 
573 BlockMatrix * Mult(const BlockMatrix & A, const BlockMatrix & B)
574 {
575  BlockMatrix * C= new BlockMatrix(A.RowOffsets(), B.ColOffsets());
576  C->owns_blocks = 1;
577  Array<SparseMatrix *> CijPieces(A.NumColBlocks());
578 
579  for (int irowC = 0; irowC < A.NumRowBlocks(); ++irowC)
580  for (int jcolC = 0; jcolC < B.NumColBlocks(); ++jcolC)
581  {
582  CijPieces.SetSize(0, static_cast<SparseMatrix *>(NULL));
583  for (int k = 0; k < A.NumColBlocks(); ++k)
584  if (!A.IsZeroBlock(irowC, k) && !B.IsZeroBlock(k, jcolC))
585  {
586  CijPieces.Append(Mult(A.GetBlock(irowC, k), B.GetBlock(k, jcolC)));
587  }
588 
589  if (CijPieces.Size() > 1)
590  {
591  C->SetBlock(irowC, jcolC, Add(CijPieces));
592  for (SparseMatrix ** it = CijPieces.GetData();
593  it != CijPieces.GetData()+CijPieces.Size(); ++it)
594  {
595  delete *it;
596  }
597  }
598  else if (CijPieces.Size() == 1)
599  {
600  C->SetBlock(irowC, jcolC, CijPieces[0]);
601  }
602  }
603 
604  return C;
605 }
606 
607 }
int Size() const
Logical size of the array.
Definition: array.hpp:133
SparseMatrix * CreateMonolithic() const
Returns a monolithic CSR matrix that represents this operator.
int owns_blocks
If owns_blocks the SparseMatrix objects Aij will be deallocated.
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:328
Set the diagonal value to zero.
Definition: matrix.hpp:34
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:478
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
Abstract data type for sparse matrices.
Definition: matrix.hpp:77
T * GetData()
Returns the data.
Definition: array.hpp:115
int Size() const
Returns the size of the vector.
Definition: vector.hpp:120
BlockMatrix(const Array< int > &offsets)
Constructor for square block matrices.
Definition: blockmatrix.cpp:22
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:129
virtual void AddMultTranspose(const Vector &x, Vector &y, const double val=1.) const
MatrixTranspose-Vector Multiplication y = y + val*A&#39;*x.
virtual int GetRow(const int row, Array< int > &cols, Vector &srow) const
Gets the columns indexes and values for row row.
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2942
void SetBlock(int i, int j, SparseMatrix *mat)
Set A(i,j) = mat.
Definition: blockmatrix.cpp:60
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
virtual void Finalize(int skip_zeros=1)
Finalize all the submatrices.
Definition: blockmatrix.hpp:85
Data type sparse matrix.
Definition: sparsemat.hpp:38
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
int IsZeroBlock(int i, int j) const
Check if block (i,j) is a zero block.
Definition: blockmatrix.hpp:55
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:146
int NumRowBlocks() const
Return the number of row blocks.
Definition: blockmatrix.hpp:46
Array< int > & ColOffsets()
Return the columns offsets for block starts.
Definition: blockmatrix.hpp:59
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:418
void PrintMatlab(std::ostream &os=mfem::out) const
Export the monolithic matrix to file.
Array< int > & RowOffsets()
Return the row offsets for block starts.
Definition: blockmatrix.hpp:57
virtual int NumNonZeroElems() const
Returns the total number of non zeros in the matrix.
virtual ~BlockMatrix()
Destructor.
Definition: blockmatrix.cpp:50
double * GetData() const
Return element data, i.e. array A.
Definition: sparsemat.hpp:138
Set the diagonal value to one.
Definition: matrix.hpp:35
virtual void MultTranspose(const Vector &x, Vector &y) const
MatrixTranspose-Vector Multiplication y = A&#39;*x.
virtual void Mult(const Vector &x, Vector &y) const
Matrix-Vector Multiplication y = A*x.
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:569
int NumColBlocks() const
Return the number of column blocks.
Definition: blockmatrix.hpp:48
void EliminateRowCol(int rc, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate the row and column rc from the matrix.
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:101
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:24
virtual void AddMult(const Vector &x, Vector &y, const double val=1.) const
Matrix-Vector Multiplication y = y + val*A*x.
SparseMatrix & GetBlock(int i, int j)
Return a reference to block (i,j). Reference may be invalid if Aij(i,j) == NULL.
Definition: blockmatrix.cpp:81
T & Last()
Return the last element in the array.
Definition: array.hpp:647
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:720
Vector data type.
Definition: vector.hpp:48
int RowSize(const int i) const
Return the number of non zeros in row i.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:64
virtual void EliminateZeroRows(const double threshold=1e-12)
If the matrix is square, this method will place 1 on the diagonal (i,i) if row i has &quot;almost&quot; zero l1...
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:25