MFEM  v4.6.0
Finite element discretization library
blockmatrix.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 #include <vector>
19 
20 namespace mfem
21 {
22 
24  AbstractSparseMatrix(offsets.Last()),
25  owns_blocks(false),
26  nRowBlocks(offsets.Size()-1),
27  nColBlocks(offsets.Size()-1),
28  Aij(nRowBlocks, nColBlocks)
29 {
30  row_offsets.MakeRef(const_cast< Array<int>& >(offsets));
31  col_offsets.MakeRef(const_cast< Array<int>& >(offsets));
32  Aij = (SparseMatrix *)NULL;
33 }
34 
35 BlockMatrix::BlockMatrix(const Array<int> & row_offsets_,
36  const Array<int> & col_offsets_):
37  AbstractSparseMatrix(row_offsets_.Last(), col_offsets_.Last()),
38  owns_blocks(false),
39  nRowBlocks(row_offsets_.Size()-1),
40  nColBlocks(col_offsets_.Size()-1),
41  Aij(nRowBlocks, nColBlocks)
42 {
43  row_offsets.MakeRef(const_cast< Array<int>& >(row_offsets_));
44  col_offsets.MakeRef(const_cast< Array<int>& >(col_offsets_));
45  Aij = (SparseMatrix *)NULL;
46 }
47 
48 
50 {
51  if (owns_blocks)
52  {
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 
61 void BlockMatrix::SetBlock(int i, int j, SparseMatrix * mat)
62 {
63 #ifdef MFEM_DEBUG
64  if (nRowBlocks <= i || nColBlocks <= j)
65  {
66  mfem_error("BlockMatrix::SetBlock #0");
67  }
68 
69  if (mat->Height() != row_offsets[i+1] - row_offsets[i])
70  {
71  mfem_error("BlockMatrix::SetBlock #1");
72  }
73 
74  if (mat->Width() != col_offsets[j+1] - col_offsets[j])
75  {
76  mfem_error("BlockMatrix::SetBlock #2");
77  }
78 #endif
79  Aij(i,j) = mat;
80 }
81 
83 {
84 #ifdef MFEM_DEBUG
85  if (nRowBlocks <= i || nColBlocks <= j)
86  {
87  mfem_error("BlockMatrix::Block #0");
88  }
89 
90  if (IsZeroBlock(i,j))
91  {
92  mfem_error("BlockMatrix::Block #1");
93  }
94 #endif
95  return *Aij(i,j);
96 }
97 
98 const SparseMatrix & BlockMatrix::GetBlock(int i, int j) const
99 {
100 #ifdef MFEM_DEBUG
101  if (nRowBlocks <= i || nColBlocks <= j)
102  {
103  mfem_error("BlockMatrix::Block const #0");
104  }
105 
106  if (IsZeroBlock(i,j))
107  {
108  mfem_error("BlockMatrix::Block const #1");
109  }
110 #endif
111 
112  return *Aij(i,j);
113 }
114 
115 
117 {
118  int nnz_elem = 0;
119  for (int jcol = 0; jcol != nColBlocks; ++jcol)
120  {
121  for (int irow = 0; irow != nRowBlocks; ++irow)
122  {
123  if (Aij(irow,jcol))
124  {
125  nnz_elem+= Aij(irow,jcol)->NumNonZeroElems();
126  }
127  }
128  }
129  return nnz_elem;
130 }
131 
132 
133 double& BlockMatrix::Elem (int i, int j)
134 {
135  int iloc, iblock;
136  int jloc, jblock;
137 
138  findGlobalRow(i, iblock, iloc);
139  findGlobalCol(j, jblock, jloc);
140 
141  if (IsZeroBlock(iblock, jblock))
142  {
143  mfem_error("BlockMatrix::Elem");
144  }
145 
146  return Aij(iblock, jblock)->Elem(iloc, jloc);
147 }
148 
149 const double& BlockMatrix::Elem (int i, int j) const
150 {
151  static const double zero = 0.0;
152  int iloc, iblock;
153  int jloc, jblock;
154 
155  findGlobalRow(i, iblock, iloc);
156  findGlobalCol(j, jblock, jloc);
157 
158  if (IsZeroBlock(iblock, jblock))
159  {
160  return zero;
161  }
162  return static_cast<const SparseMatrix *>(Aij(iblock, jblock))->Elem(iloc, jloc);
163 }
164 
165 int BlockMatrix::RowSize(const int i) const
166 {
167  int rowsize = 0;
168 
169  int iblock, iloc;
170  findGlobalRow(i, iblock, iloc);
171 
172  for (int jblock = 0; jblock < nColBlocks; ++jblock)
173  {
174  if (Aij(iblock,jblock) != NULL)
175  {
176  rowsize += Aij(iblock,jblock)->RowSize(iloc);
177  }
178  }
179 
180  return rowsize;
181 }
182 
183 int BlockMatrix::GetRow(const int row, Array<int> &cols, Vector &srow) const
184 {
185  int iblock, iloc, rowsize;
186  findGlobalRow(row, iblock, iloc);
187  rowsize = RowSize(row);
188  cols.SetSize(rowsize);
189  srow.SetSize(rowsize);
190 
191  Array<int> bcols;
192  Vector bsrow;
193 
194  int * it_cols = cols.GetData();
195  double *it_srow = srow.GetData();
196 
197  for (int jblock = 0; jblock < nColBlocks; ++jblock)
198  {
199  if (Aij(iblock,jblock) != NULL)
200  {
201  Aij(iblock,jblock)->GetRow(iloc, bcols, bsrow);
202  for (int i = 0; i < bcols.Size(); ++i)
203  {
204  *(it_cols++) = bcols[i] + col_offsets[jblock];
205  *(it_srow++) = bsrow(i);
206  }
207  }
208  }
209 
210  return 0;
211 }
212 
214 {
215  // Find the block to which the dof belongs and its local number
216  int idx, iiblock;
217  for (iiblock = 0; iiblock < nRowBlocks; ++iiblock)
218  {
219  idx = rc - row_offsets[iiblock];
220  if (idx < 0 ) { break; }
221  }
222  iiblock--;
223  idx = rc - row_offsets[iiblock];
224 
225  // Asserts
226  MFEM_ASSERT(nRowBlocks == nColBlocks,
227  "BlockMatrix::EliminateRowCol: nRowBlocks != nColBlocks");
228 
229  MFEM_ASSERT(row_offsets[iiblock] == col_offsets[iiblock],
230  "BlockMatrix::EliminateRowCol: row_offsets["
231  << iiblock << "] != col_offsets["<<iiblock<<"]");
232 
233  MFEM_ASSERT(Aij(iiblock, iiblock),
234  "BlockMatrix::EliminateRowCol: Null diagonal block");
235 
236  // Apply the constraint idx to the iiblock
237  for (int jjblock = 0; jjblock < nRowBlocks; ++jjblock)
238  {
239  if (iiblock == jjblock) { continue; }
240  if (Aij(iiblock,jjblock)) { Aij(iiblock,jjblock)->EliminateRow(idx); }
241  }
242  for (int jjblock = 0; jjblock < nRowBlocks; ++jjblock)
243  {
244  if (iiblock == jjblock) { continue; }
245  if (Aij(jjblock,iiblock)) { Aij(jjblock,iiblock)->EliminateCol(idx); }
246  }
247  Aij(iiblock, iiblock)->EliminateRowCol(idx,dpolicy);
248 }
249 
251  Vector & rhs)
252 {
253  if (nRowBlocks != nColBlocks)
254  {
255  mfem_error("BlockMatrix::EliminateRowCol: nRowBlocks != nColBlocks");
256  }
257 
258  for (int iiblock = 0; iiblock < nRowBlocks; ++iiblock)
259  {
260  if (row_offsets[iiblock] != col_offsets[iiblock])
261  {
262  mfem::out << "BlockMatrix::EliminateRowCol: row_offsets["
263  << iiblock << "] != col_offsets["<<iiblock<<"]\n";
264  mfem_error();
265  }
266  }
267 
268  // We also have to do the same for each Aij
269  Array<int> block_dofs;
270  Vector block_sol, block_rhs;
271 
272  for (int iiblock = 0; iiblock < nRowBlocks; ++iiblock)
273  {
274  int dsize = row_offsets[iiblock+1] - row_offsets[iiblock];
275  block_dofs.MakeRef(ess_bc_dofs.GetData()+row_offsets[iiblock], dsize);
276  block_sol.SetDataAndSize(sol.GetData()+row_offsets[iiblock], dsize);
277  block_rhs.SetDataAndSize(rhs.GetData()+row_offsets[iiblock], dsize);
278 
279  if (Aij(iiblock, iiblock))
280  {
281  for (int i = 0; i < block_dofs.Size(); ++i)
282  {
283  if (block_dofs[i])
284  {
285  Aij(iiblock, iiblock)->EliminateRowCol(i,block_sol(i), block_rhs);
286  }
287  }
288  }
289  else
290  {
291  for (int i = 0; i < block_dofs.Size(); ++i)
292  {
293  if (block_dofs[i])
294  {
295  mfem_error("BlockMatrix::EliminateRowCol: Null diagonal block \n");
296  }
297  }
298  }
299 
300  for (int jjblock = 0; jjblock < nRowBlocks; ++jjblock)
301  {
302  if (jjblock != iiblock && Aij(iiblock, jjblock))
303  {
304  for (int i = 0; i < block_dofs.Size(); ++i)
305  {
306  if (block_dofs[i])
307  {
308  Aij(iiblock, jjblock)->EliminateRow(i);
309  }
310  }
311  }
312  if (jjblock != iiblock && Aij(jjblock, iiblock))
313  {
314  block_rhs.SetDataAndSize(rhs.GetData()+row_offsets[jjblock],
315  row_offsets[jjblock+1] - row_offsets[jjblock]);
316  Aij(jjblock, iiblock)->EliminateCols(block_dofs, &block_sol, &block_rhs);
317  }
318  }
319  }
320 }
321 
323  DiagonalPolicy dpolicy)
324 {
325  MFEM_VERIFY(Ae,
326  "BlockMatrix::EliminateRowCols: Elimination matrix pointer is null");
327  MFEM_VERIFY(nRowBlocks == nColBlocks,
328  "BlockMatrix::EliminateRowCols supported only for"
329  "nRowBlocks = nColBlocks");
330 
331  std::vector<Array<int>> cols(nRowBlocks);
332  std::vector<Array<int>> rows(nRowBlocks);
333  SparseMatrix * tmp = nullptr;
334 
335  for (int k = 0; k < vdofs.Size(); k++)
336  {
337  int vdof = (vdofs[k]) >=0 ? vdofs[k] : -1 - vdofs[k];
338  // find block
339  int iblock, dof;
340  findGlobalCol(vdof,iblock,dof);
341  cols[iblock].Append(dof);
342  tmp = &GetBlock(iblock,iblock);
343  if (tmp)
344  {
345  tmp->EliminateRowCol(dof,Ae->GetBlock(iblock,iblock), dpolicy);
346  }
347  }
348 
349  // Eliminate col from off-diagonal blocks
350  for (int j = 0; j<nColBlocks; j++)
351  {
352  if (!cols[j].Size()) { continue; }
353  int blocksize = col_offsets[j+1] - col_offsets[j];
354  Array<int> colmarker(blocksize); colmarker = 0;
355  for (int i = 0; i < cols[j].Size(); i++) { colmarker[cols[j][i]] = 1; }
356  for (int i = 0; i<nRowBlocks; i++)
357  {
358  if (i == j) { continue; }
359  tmp = &GetBlock(i,j);
360  if (tmp) { tmp->EliminateCols(colmarker,Ae->GetBlock(i,j)); }
361  for (int k = 0; k < cols[j].Size(); k++)
362  {
363  tmp = &GetBlock(j,i);
364  if (tmp) { tmp->EliminateRow(cols[j][k]); }
365  }
366  }
367  }
368 }
369 
370 void BlockMatrix::EliminateZeroRows(const double threshold)
371 {
372  MFEM_VERIFY(nRowBlocks == nColBlocks, "not a square matrix");
373 
374  for (int iblock = 0; iblock < nRowBlocks; ++iblock)
375  {
376  if (Aij(iblock,iblock))
377  {
378  double norm;
379  for (int i = 0; i < Aij(iblock, iblock)->NumRows(); ++i)
380  {
381  norm = 0.;
382  for (int jblock = 0; jblock < nColBlocks; ++jblock)
383  if (Aij(iblock,jblock))
384  {
385  norm += Aij(iblock,jblock)->GetRowNorml1(i);
386  }
387 
388  if (norm <= threshold)
389  {
390  for (int jblock = 0; jblock < nColBlocks; ++jblock)
391  {
392  if (Aij(iblock,jblock))
393  {
394  Aij(iblock,jblock)->EliminateRow(
395  i, (iblock==jblock) ? DIAG_ONE : DIAG_ZERO);
396  }
397  }
398  }
399  }
400  }
401  else
402  {
403  double norm;
404  for (int i = 0; i < row_offsets[iblock+1] - row_offsets[iblock]; ++i)
405  {
406  norm = 0.;
407  for (int jblock = 0; jblock < nColBlocks; ++jblock)
408  {
409  if (Aij(iblock,jblock))
410  {
411  norm += Aij(iblock,jblock)->GetRowNorml1(i);
412  }
413  }
414 
415  MFEM_VERIFY(!(norm <= threshold), "diagonal block is NULL:"
416  " iblock = " << iblock << ", i = " << i << ", norm = "
417  << norm);
418  }
419  }
420  }
421 }
422 
423 void BlockMatrix::Finalize(int skip_zeros, bool fix_empty_rows)
424 {
425  for (int iblock = 0; iblock < nRowBlocks; ++iblock)
426  {
427  for (int jblock = 0; jblock < nColBlocks; ++jblock)
428  {
429  if (!Aij(iblock,jblock)) { continue; }
430  if (!Aij(iblock,jblock)->Finalized())
431  {
432  Aij(iblock,jblock)->Finalize(skip_zeros, fix_empty_rows);
433  }
434  }
435  }
436 }
437 
438 void BlockMatrix::Mult(const Vector & x, Vector & y) const
439 {
440  if (x.GetData() == y.GetData())
441  {
442  mfem_error("Error: x and y can't point to the same data \n");
443  }
444 
445  MFEM_ASSERT(width == x.Size(), "Input vector size (" << x.Size()
446  << ") must match matrix width (" << width << ")");
447  MFEM_ASSERT(height == y.Size(), "Output vector size (" << y.Size()
448  << ") must match matrix height (" << height << ")");
449  y = 0.;
450  AddMult(x, y, 1.0);
451 }
452 
453 void BlockMatrix::AddMult(const Vector & x, Vector & y, const double val) const
454 {
455  if (x.GetData() == y.GetData())
456  {
457  mfem_error("Error: x and y can't point to the same data \n");
458  }
459 
460  Vector xblockview, yblockview;
461 
462  for (int iblock = 0; iblock != nRowBlocks; ++iblock)
463  {
464  yblockview.SetDataAndSize(y.GetData() + row_offsets[iblock],
465  row_offsets[iblock+1] - row_offsets[iblock]);
466 
467  for (int jblock = 0; jblock != nColBlocks; ++jblock)
468  {
469  if (Aij(iblock, jblock) != NULL)
470  {
471  xblockview.SetDataAndSize(
472  x.GetData() + col_offsets[jblock],
473  col_offsets[jblock+1] - col_offsets[jblock]);
474 
475  Aij(iblock, jblock)->AddMult(xblockview, yblockview, val);
476  }
477  }
478  }
479 }
480 
481 void BlockMatrix::MultTranspose(const Vector & x, Vector & y) const
482 {
483  if (x.GetData() == y.GetData())
484  {
485  mfem_error("Error: x and y can't point to the same data \n");
486  }
487 
488  y = 0.;
489  AddMultTranspose(x, y, 1.0);
490 }
491 
493  const double val) const
494 {
495  if (x.GetData() == y.GetData())
496  {
497  mfem_error("Error: x and y can't point to the same data \n");
498  }
499 
500  Vector xblockview, yblockview;
501 
502  for (int iblock = 0; iblock != nColBlocks; ++iblock)
503  {
504  yblockview.SetDataAndSize(y.GetData() + col_offsets[iblock],
505  col_offsets[iblock+1] - col_offsets[iblock]);
506 
507  for (int jblock = 0; jblock != nRowBlocks; ++jblock)
508  {
509  if (Aij(jblock, iblock) != NULL)
510  {
511  xblockview.SetDataAndSize(
512  x.GetData() + row_offsets[jblock],
513  row_offsets[jblock+1] - row_offsets[jblock]);
514 
515  Aij(jblock, iblock)->AddMultTranspose(xblockview, yblockview, val);
516  }
517  }
518  }
519 }
520 
521 void BlockMatrix::PartMult(const Array<int> &rows, const Vector &x,
522  Vector &y) const
523 {
524  Array<int> cols;
525  Vector srow;
526  for (int i = 0; i<rows.Size(); i++)
527  {
528  int dof = (rows[i]>=0) ? rows[i] : -1-rows[i];
529  GetRow(dof,cols,srow);
530 
531  double s=0.0;
532  for (int k = 0; k <cols.Size(); k++)
533  {
534  s += srow[k] * x[cols[k]];
535  }
536  y[dof] = s;
537  }
538 }
539 void BlockMatrix::PartAddMult(const Array<int> &rows, const Vector &x,
540  Vector &y,
541  const double a) const
542 {
543  Array<int> cols;
544  Vector srow;
545  for (int i = 0; i<rows.Size(); i++)
546  {
547  int dof = (rows[i]>=0) ? rows[i] : -1-rows[i];
548  GetRow(dof,cols,srow);
549 
550  double s=0.0;
551  for (int k = 0; k <cols.Size(); k++)
552  {
553  s += srow[k] * x[cols[k]];
554  }
555  y[dof] += a * s;
556  }
557 }
558 
559 
561 {
562  int nnz = NumNonZeroElems();
563 
564  int * i_amono = Memory<int>(row_offsets[nRowBlocks]+2);
565  int * j_amono = Memory<int>(nnz);
566  double * data = Memory<double>(nnz);
567 
568  for (int i = 0; i < row_offsets[nRowBlocks]+2; i++)
569  {
570  i_amono[i] = 0;
571  }
572 
573  int * i_amono_construction = i_amono+1;
574 
575  for (int iblock = 0; iblock != nRowBlocks; ++iblock)
576  {
577  for (int irow(row_offsets[iblock]); irow < row_offsets[iblock+1]; ++irow)
578  {
579  int local_row = irow - row_offsets[iblock];
580  int ind = i_amono_construction[irow];
581  for (int jblock = 0; jblock < nColBlocks; ++jblock)
582  {
583  if (Aij(iblock,jblock) != NULL)
584  ind += Aij(iblock, jblock)->GetI()[local_row+1]
585  - Aij(iblock, jblock)->GetI()[local_row];
586  }
587  i_amono_construction[irow+1] = ind;
588  }
589  }
590 
591  // Fill in the jarray and copy the data
592  for (int iblock = 0; iblock != nRowBlocks; ++iblock)
593  {
594  for (int jblock = 0; jblock != nColBlocks; ++jblock)
595  {
596  if (Aij(iblock,jblock) != NULL)
597  {
598  int nrow = row_offsets[iblock+1]-row_offsets[iblock];
599  int * i_aij = Aij(iblock, jblock)->GetI();
600  int * j_aij = Aij(iblock, jblock)->GetJ();
601  double * data_aij = Aij(iblock, jblock)->GetData();
602  int *i_it = i_amono_construction+row_offsets[iblock];
603 
604  int loc_start_index = 0;
605  int loc_end_index = 0;
606  int glob_start_index = 0;
607 
608  int shift(col_offsets[jblock]);
609  for (int * i_it_aij(i_aij+1); i_it_aij != i_aij+nrow+1; ++i_it_aij)
610  {
611  glob_start_index = *i_it;
612 
613 #ifdef MFEM_DEBUG
614  if (glob_start_index > nnz)
615  {
616  mfem::out<<"glob_start_index = " << glob_start_index << "\n";
617  mfem::out<<"Block:" << iblock << " " << jblock << "\n";
618  mfem::out<<std::endl;
619  }
620 #endif
621  loc_end_index = *(i_it_aij);
622  for (int cnt = 0; cnt < loc_end_index-loc_start_index; cnt++)
623  {
624  data[glob_start_index+cnt] = data_aij[loc_start_index+cnt];
625  j_amono[glob_start_index+cnt] = j_aij[loc_start_index+cnt] + shift;
626  }
627 
628  *i_it += loc_end_index-loc_start_index;
629  ++i_it;
630  loc_start_index = loc_end_index;
631  }
632  }
633  }
634  }
635 
636  return new SparseMatrix(i_amono, j_amono, data, row_offsets[nRowBlocks],
637  col_offsets[nColBlocks]);
638 }
639 
640 void BlockMatrix::PrintMatlab(std::ostream & os) const
641 {
642 
643  Vector row_data;
644  Array<int> row_ind;
645  int nnz_elem = NumNonZeroElems();
646  os<<"% size " << row_offsets.Last() << " " << col_offsets.Last() << "\n";
647  os<<"% Non Zeros " << nnz_elem << "\n";
648  int i, j;
649  std::ios::fmtflags old_fmt = os.flags();
650  os.setf(std::ios::scientific);
651  std::streamsize old_prec = os.precision(14);
652  for (i = 0; i < row_offsets.Last(); i++)
653  {
654  GetRow(i, row_ind, row_data);
655  for (j = 0; j < row_ind.Size(); j++)
656  {
657  os << i+1 << " " << row_ind[j]+1 << " " << row_data[j] << std::endl;
658  }
659  }
660  // Write a zero entry at (m,n) to make sure MATLAB doesn't shrink the matrix
661  os << row_offsets.Last() << " " << col_offsets.Last () << " 0.0\n";
662 
663  os.precision(old_prec);
664  os.flags(old_fmt);
665 }
666 
668 {
669  BlockMatrix * At = new BlockMatrix(A.ColOffsets(), A.RowOffsets());
670  At->owns_blocks = 1;
671 
672  for (int irowAt = 0; irowAt < At->NumRowBlocks(); ++irowAt)
673  {
674  for (int jcolAt = 0; jcolAt < At->NumColBlocks(); ++jcolAt)
675  {
676  if (!A.IsZeroBlock(jcolAt, irowAt))
677  {
678  At->SetBlock(irowAt, jcolAt, Transpose(A.GetBlock(jcolAt, irowAt)));
679  }
680  }
681  }
682  return At;
683 }
684 
685 BlockMatrix * Mult(const BlockMatrix & A, const BlockMatrix & B)
686 {
687  BlockMatrix * C= new BlockMatrix(A.RowOffsets(), B.ColOffsets());
688  C->owns_blocks = 1;
689  Array<SparseMatrix *> CijPieces(A.NumColBlocks());
690 
691  for (int irowC = 0; irowC < A.NumRowBlocks(); ++irowC)
692  {
693  for (int jcolC = 0; jcolC < B.NumColBlocks(); ++jcolC)
694  {
695  CijPieces.SetSize(0, static_cast<SparseMatrix *>(NULL));
696  for (int k = 0; k < A.NumColBlocks(); ++k)
697  {
698  if (!A.IsZeroBlock(irowC, k) && !B.IsZeroBlock(k, jcolC))
699  {
700  CijPieces.Append(Mult(A.GetBlock(irowC, k), B.GetBlock(k, jcolC)));
701  }
702  }
703 
704  if (CijPieces.Size() > 1)
705  {
706  C->SetBlock(irowC, jcolC, Add(CijPieces));
707  for (SparseMatrix ** it = CijPieces.GetData();
708  it != CijPieces.GetData()+CijPieces.Size(); ++it)
709  {
710  delete *it;
711  }
712  }
713  else if (CijPieces.Size() == 1)
714  {
715  C->SetBlock(irowC, jcolC, CijPieces[0]);
716  }
717  }
718  }
719 
720  return C;
721 }
722 
723 }
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
Partial matrix vector multiplication of (*this) with x involving only the rows given by rows...
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:517
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
void PartAddMult(const Array< int > &rows, const Vector &x, Vector &y, const double a=1.0) const
Partial matrix vector multiplication of (*this) with x involving only the rows given by rows...
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
Abstract data type for sparse matrices.
Definition: matrix.hpp:73
T * GetData()
Returns the data.
Definition: array.hpp:115
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
virtual void PrintMatlab(std::ostream &os=mfem::out) const
Export the monolithic matrix to file.
BlockMatrix(const Array< int > &offsets)
Constructor for square block matrices.
Definition: blockmatrix.cpp:23
virtual int NumNonZeroElems() const
Returns the total number of non zeros in the matrix.
double * GetData()
Return the element data, i.e. the array A.
Definition: sparsemat.hpp:232
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2320
void SetBlock(int i, int j, SparseMatrix *mat)
Set A(i,j) = mat.
Definition: blockmatrix.cpp:61
int NumRowBlocks() const
Return the number of row blocks.
Definition: blockmatrix.hpp:46
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:91
Data type sparse matrix.
Definition: sparsemat.hpp:50
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:154
virtual void MultTranspose(const Vector &x, Vector &y) const
MatrixTranspose-Vector Multiplication y = A&#39;*x.
Array< int > & ColOffsets()
Return the columns offsets for block starts.
Definition: blockmatrix.hpp:59
int NumColBlocks() const
Return the number of column blocks.
Definition: blockmatrix.hpp:48
Array< int > & RowOffsets()
Return the row offsets for block starts.
Definition: blockmatrix.hpp:57
Set the diagonal value to one.
Definition: operator.hpp:50
virtual ~BlockMatrix()
Destructor.
Definition: blockmatrix.cpp:49
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:206
virtual int GetRow(const int row, Array< int > &cols, Vector &srow) const
Gets the columns indexes and values for row row.
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
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
void EliminateRowCol(int rc, const double sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate row rc and column rc and modify the rhs using sol.
Definition: sparsemat.cpp:1854
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
int RowSize(const int i) const
Return the number of non zeros in row i.
virtual void Mult(const Vector &x, Vector &y) const
Matrix-Vector Multiplication y = A*x.
void EliminateRowCol(int rc, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate the row and column rc from the matrix.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:154
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
void EliminateRow(int row, const double sol, Vector &rhs)
Eliminates a column from the transpose matrix.
Definition: sparsemat.cpp:1693
double a
Definition: lissajous.cpp:41
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition: operator.hpp:47
void EliminateCols(const Array< int > &cols, const Vector *x=NULL, Vector *b=NULL)
Eliminate all columns i for which cols[i] != 0.
Definition: sparsemat.cpp:1781
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:82
virtual void AddMult(const Vector &x, Vector &y, const double val=1.) const
Matrix-Vector Multiplication y = y + val*A*x.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
T & Last()
Return the last element in the array.
Definition: array.hpp:792
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:872
Vector data type.
Definition: vector.hpp:58
RefCoord s[3]
SparseMatrix * CreateMonolithic() const
Returns a monolithic CSR matrix that represents this operator.
virtual void AddMultTranspose(const Vector &x, Vector &y, const double val=1.) const
MatrixTranspose-Vector Multiplication y = y + val*A&#39;*x.
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 "almost" zero l1...
void EliminateRowCols(const Array< int > &vdofs, BlockMatrix *Ae, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate the rows and columns corresponding to the entries in vdofs + save the eliminated entries in...
Set the diagonal value to zero.
Definition: operator.hpp:49
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
int IsZeroBlock(int i, int j) const
Check if block (i,j) is a zero block.
Definition: blockmatrix.hpp:55