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