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