MFEM  v3.2
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 "matrix.hpp"
14 #include "sparsemat.hpp"
15 #include "blockvector.hpp"
16 #include "blockmatrix.hpp"
17 
18 namespace mfem
19 {
20 
22  AbstractSparseMatrix(offsets.Last()),
23  owns_blocks(false),
24  nRowBlocks(offsets.Size()-1),
25  nColBlocks(offsets.Size()-1),
26  row_offsets(const_cast< Array<int>& >(offsets).GetData(), offsets.Size()),
27  col_offsets(const_cast< Array<int>& >(offsets).GetData(), offsets.Size()),
28  Aij(nRowBlocks, nColBlocks)
29 {
30  Aij = (SparseMatrix *)NULL;
31 }
32 
33 BlockMatrix::BlockMatrix(const Array<int> & row_offsets_,
34  const Array<int> & col_offsets_):
35  AbstractSparseMatrix(row_offsets_.Last(), col_offsets_.Last()),
36  owns_blocks(false),
37  nRowBlocks(row_offsets_.Size()-1),
38  nColBlocks(col_offsets_.Size()-1),
39  row_offsets(const_cast< Array<int>& >(row_offsets_).GetData(),
40  row_offsets_.Size()),
41  col_offsets(const_cast< Array<int>& >(col_offsets_).GetData(),
42  col_offsets_.Size()),
43  Aij(nRowBlocks, nColBlocks)
44 {
45  Aij = (SparseMatrix *)NULL;
46 }
47 
48 
50 {
51  if (owns_blocks)
52  for (SparseMatrix ** it = Aij.GetRow(0);
53  it != Aij.GetRow(0)+(Aij.NumRows()*Aij.NumCols()); ++it)
54  {
55  delete *it;
56  }
57 }
58 
59 void BlockMatrix::SetBlock(int i, int j, SparseMatrix * mat)
60 {
61 #ifdef MFEM_DEBUG
62  if (nRowBlocks <= i || nColBlocks <= j)
63  {
64  mfem_error("BlockMatrix::SetBlock #0");
65  }
66 
67  if (mat->Height() != row_offsets[i+1] - row_offsets[i])
68  {
69  mfem_error("BlockMatrix::SetBlock #1");
70  }
71 
72  if (mat->Width() != col_offsets[j+1] - col_offsets[j])
73  {
74  mfem_error("BlockMatrix::SetBlock #2");
75  }
76 #endif
77  Aij(i,j) = mat;
78 }
79 
81 {
82 #ifdef MFEM_DEBUG
83  if (nRowBlocks <= i || nColBlocks <= j)
84  {
85  mfem_error("BlockMatrix::Block #0");
86  }
87 
88  if (IsZeroBlock(i,j))
89  {
90  mfem_error("BlockMatrix::Block #1");
91  }
92 #endif
93  return *Aij(i,j);
94 }
95 
96 const SparseMatrix & BlockMatrix::GetBlock(int i, int j) const
97 {
98 #ifdef MFEM_DEBUG
99  if (nRowBlocks <= i || nColBlocks <= j)
100  {
101  mfem_error("BlockMatrix::Block const #0");
102  }
103 
104  if (IsZeroBlock(i,j))
105  {
106  mfem_error("BlockMatrix::Block const #1");
107  }
108 #endif
109 
110  return *Aij(i,j);
111 }
112 
113 
115 {
116  int nnz_elem = 0;
117  for (int jcol = 0; jcol != nColBlocks; ++jcol)
118  for (int irow = 0; irow != nRowBlocks; ++irow)
119  {
120  if (Aij(irow,jcol))
121  {
122  nnz_elem+= Aij(irow,jcol)->NumNonZeroElems();
123  }
124  }
125  return nnz_elem;
126 }
127 
128 
129 double& BlockMatrix::Elem (int i, int j)
130 {
131  int iloc, iblock;
132  int jloc, jblock;
133 
134  findGlobalRow(i, iblock, iloc);
135  findGlobalCol(j, jblock, jloc);
136 
137  if (IsZeroBlock(i, j))
138  {
139  mfem_error("BlockMatrix::Elem");
140  }
141 
142  return Aij(iblock, jblock)->Elem(iloc, jloc);
143 }
144 
145 const double& BlockMatrix::Elem (int i, int j) const
146 {
147  int iloc, iblock;
148  int jloc, jblock;
149 
150  findGlobalRow(i, iblock, iloc);
151  findGlobalCol(j, jblock, jloc);
152 
153  if (IsZeroBlock(i, j))
154  {
155  mfem_error("BlockMatrix::Elem");
156  }
157 
158  return Aij(iblock, jblock)->Elem(iloc, jloc);
159 }
160 
161 int BlockMatrix::RowSize(const int i) const
162 {
163  int rowsize = 0;
164 
165  int iblock, iloc;
166  findGlobalRow(i, iblock, iloc);
167 
168  for (int jblock = 0; jblock < nColBlocks; ++jblock)
169  if (Aij(iblock,jblock) != NULL)
170  {
171  rowsize += Aij(iblock,jblock)->RowSize(iloc);
172  }
173 
174  return rowsize;
175 }
176 
177 int BlockMatrix::GetRow(const int row, Array<int> &cols, Vector &srow) const
178 {
179  int iblock, iloc, rowsize;
180  findGlobalRow(row, iblock, iloc);
181  rowsize = RowSize(row);
182  cols.SetSize(rowsize);
183  srow.SetSize(rowsize);
184 
185  Array<int> bcols;
186  Vector bsrow;
187 
188  int * it_cols = cols.GetData();
189  double *it_srow = srow.GetData();
190 
191  for (int jblock = 0; jblock < nColBlocks; ++jblock)
192  if (Aij(iblock,jblock) != NULL)
193  {
194  Aij(iblock,jblock)->GetRow(iloc, bcols, bsrow);
195  for (int i = 0; i < bcols.Size(); ++i)
196  {
197  *(it_cols++) = bcols[i] + col_offsets[jblock];
198  *(it_srow++) = bsrow(i);
199  }
200  }
201 
202  return 0;
203 }
204 
206  Vector & rhs)
207 {
208  if (nRowBlocks != nColBlocks)
209  {
210  mfem_error("BlockMatrix::EliminateRowCol: nRowBlocks != nColBlocks");
211  }
212 
213  for (int iiblock = 0; iiblock < nRowBlocks; ++iiblock)
214  if (row_offsets[iiblock] != col_offsets[iiblock])
215  {
216  std::cout << "BlockMatrix::EliminateRowCol: row_offests["
217  << iiblock << "] != col_offsets["<<iiblock<<"]\n";
218  mfem_error();
219  }
220 
221  // We also have to do the same for each Aij
222  Array<int> block_dofs;
223  Vector block_sol, block_rhs;
224 
225  for (int iiblock = 0; iiblock < nRowBlocks; ++iiblock)
226  {
227  int dsize = row_offsets[iiblock+1] - row_offsets[iiblock];
228  block_dofs.MakeRef(ess_bc_dofs.GetData()+row_offsets[iiblock], dsize);
229  block_sol.SetDataAndSize(sol.GetData()+row_offsets[iiblock], dsize);
230  block_rhs.SetDataAndSize(rhs.GetData()+row_offsets[iiblock], dsize);
231 
232  if (Aij(iiblock, iiblock))
233  {
234  for (int i = 0; i < block_dofs.Size(); ++i)
235  if (block_dofs[i])
236  {
237  Aij(iiblock, iiblock)->EliminateRowCol(i,block_sol(i), block_rhs);
238  }
239  }
240  else
241  {
242  for (int i = 0; i < block_dofs.Size(); ++i)
243  if (block_dofs[i])
244  {
245  mfem_error("BlockMatrix::EliminateRowCol: Null diagonal block \n");
246  }
247  }
248 
249  for (int jjblock = 0; jjblock < nRowBlocks; ++jjblock)
250  {
251  if (jjblock != iiblock && Aij(iiblock, jjblock))
252  {
253  for (int i = 0; i < block_dofs.Size(); ++i)
254  if (block_dofs[i])
255  {
256  Aij(iiblock, jjblock)->EliminateRow(i);
257  }
258  }
259  if (jjblock != iiblock && Aij(jjblock, iiblock))
260  {
261  block_rhs.SetDataAndSize(rhs.GetData()+row_offsets[jjblock],
262  row_offsets[jjblock+1] - row_offsets[jjblock]);
263  Aij(jjblock, iiblock)->EliminateCols(block_dofs, &block_sol, &block_rhs);
264  }
265  }
266  }
267 }
268 
270 {
271  if (nRowBlocks != nColBlocks)
272  {
273  mfem_error("BlockMatrix::EliminateZeroRows() #1");
274  }
275 
276  for (int iblock = 0; iblock < nRowBlocks; ++iblock)
277  {
278  if (Aij(iblock,iblock))
279  {
280  double norm;
281  for (int i = 0; i < Aij(iblock, iblock)->NumRows(); ++i)
282  {
283  norm = 0.;
284  for (int jblock = 0; jblock < nColBlocks; ++jblock)
285  if (Aij(iblock,jblock))
286  {
287  norm += Aij(iblock,jblock)->GetRowNorml1(i);
288  }
289 
290  if (norm < 1e-12)
291  {
292  for (int jblock = 0; jblock < nColBlocks; ++jblock)
293  if (Aij(iblock,jblock))
294  {
295  Aij(iblock,jblock)->EliminateRow(i, iblock==jblock);
296  }
297  }
298  }
299  }
300  else
301  {
302  double norm;
303  for (int i = 0; i < row_offsets[iblock+1] - row_offsets[iblock]; ++i)
304  {
305  norm = 0.;
306  for (int jblock = 0; jblock < nColBlocks; ++jblock)
307  if (Aij(iblock,jblock))
308  {
309  norm += Aij(iblock,jblock)->GetRowNorml1(i);
310  }
311 
312  if (norm < 1e-12)
313  {
314  std::cout<<"i = " << i << "\n";
315  std::cout<<"norm = " << norm << "\n";
316  mfem_error("BlockMatrix::EliminateZeroRows() #2");
317  }
318  }
319  }
320  }
321 }
322 
323 void BlockMatrix::Mult(const Vector & x, Vector & y) const
324 {
325  if (x.GetData() == y.GetData())
326  {
327  mfem_error("Error: x and y can't point to the same datas \n");
328  }
329 
330  y = 0.;
331  AddMult(x, y, 1.0);
332 }
333 
334 void BlockMatrix::AddMult(const Vector & x, Vector & y, const double val) const
335 {
336  if (x.GetData() == y.GetData())
337  {
338  mfem_error("Error: x and y can't point to the same datas \n");
339  }
340 
341  Vector xblockview, yblockview;
342 
343  for (int iblock = 0; iblock != nRowBlocks; ++iblock)
344  {
345  yblockview.SetDataAndSize(y.GetData() + row_offsets[iblock],
346  row_offsets[iblock+1] - row_offsets[iblock]);
347 
348  for (int jblock = 0; jblock != nColBlocks; ++jblock)
349  {
350  if (Aij(iblock, jblock) != NULL)
351  {
352  xblockview.SetDataAndSize(
353  x.GetData() + col_offsets[jblock],
354  col_offsets[jblock+1] - col_offsets[jblock]);
355 
356  Aij(iblock, jblock)->AddMult(xblockview, yblockview, val);
357  }
358  }
359  }
360 }
361 
362 void BlockMatrix::MultTranspose(const Vector & x, Vector & y) const
363 {
364  if (x.GetData() == y.GetData())
365  {
366  mfem_error("Error: x and y can't point to the same datas \n");
367  }
368 
369  y = 0.;
370  AddMultTranspose(x, y, 1.0);
371 }
372 
374  const double val) const
375 {
376  if (x.GetData() == y.GetData())
377  {
378  mfem_error("Error: x and y can't point to the same datas \n");
379  }
380 
381  Vector xblockview, yblockview;
382 
383  for (int iblock = 0; iblock != nColBlocks; ++iblock)
384  {
385  yblockview.SetDataAndSize(y.GetData() + col_offsets[iblock],
386  col_offsets[iblock+1] - col_offsets[iblock]);
387 
388  for (int jblock = 0; jblock != nRowBlocks; ++jblock)
389  {
390  if (Aij(jblock, iblock) != NULL)
391  {
392  xblockview.SetDataAndSize(
393  x.GetData() + row_offsets[jblock],
394  row_offsets[jblock+1] - row_offsets[jblock]);
395 
396  Aij(jblock, iblock)->AddMultTranspose(xblockview, yblockview, val);
397  }
398  }
399  }
400 }
401 
403 {
404  int nnz = NumNonZeroElems();
405 
406  int * i_amono = new int[ row_offsets[nRowBlocks]+2 ];
407  int * j_amono = new int[ nnz ];
408  double * data = new double[ nnz ];
409 
410  for (int i = 0; i < row_offsets[nRowBlocks]+2; i++)
411  {
412  i_amono[i] = 0;
413  }
414 
415  int * i_amono_construction = i_amono+1;
416 
417  int * i_it(i_amono_construction);
418 
419  for (int iblock = 0; iblock != nRowBlocks; ++iblock)
420  {
421  for (int irow(row_offsets[iblock]); irow < row_offsets[iblock+1]; ++irow)
422  {
423  int local_row = irow - row_offsets[iblock];
424  int ind = i_amono_construction[irow];
425  for (int jblock = 0; jblock < nColBlocks; ++jblock)
426  {
427  if (Aij(iblock,jblock) != NULL)
428  ind += Aij(iblock, jblock)->GetI()[local_row+1]
429  - Aij(iblock, jblock)->GetI()[local_row];
430  }
431  i_amono_construction[irow+1] = ind;
432  }
433  }
434 
435  // Fill in the jarray and copy the data
436  for (int iblock = 0; iblock != nRowBlocks; ++iblock)
437  {
438  for (int jblock = 0; jblock != nColBlocks; ++jblock)
439  {
440  if (Aij(iblock,jblock) != NULL)
441  {
442  int nrow = row_offsets[iblock+1]-row_offsets[iblock];
443  int * i_aij = Aij(iblock, jblock)->GetI();
444  int * j_aij = Aij(iblock, jblock)->GetJ();
445  double * data_aij = Aij(iblock, jblock)->GetData();
446  i_it = i_amono_construction+row_offsets[iblock];
447 
448  int loc_start_index = 0;
449  int loc_end_index = 0;
450  int glob_start_index = 0;
451 
452  int shift(col_offsets[jblock]);
453  for (int * i_it_aij(i_aij+1); i_it_aij != i_aij+nrow+1; ++i_it_aij)
454  {
455  glob_start_index = *i_it;
456 
457 #ifdef MFEM_DEBUG
458  if (glob_start_index > nnz)
459  {
460  std::cout<<"glob_start_index = " << glob_start_index << "\n";
461  std::cout<<"Block:" << iblock << " " << jblock << "\n";
462  std::cout<<std::endl;
463  }
464 #endif
465  loc_end_index = *(i_it_aij);
466  for (int cnt = 0; cnt < loc_end_index-loc_start_index; cnt++)
467  {
468  data[glob_start_index+cnt] = data_aij[loc_start_index+cnt];
469  j_amono[glob_start_index+cnt] = j_aij[loc_start_index+cnt] + shift;
470  }
471 
472  *i_it += loc_end_index-loc_start_index;
473  ++i_it;
474  loc_start_index = loc_end_index;
475  }
476  }
477  }
478  }
479 
480  return new SparseMatrix(i_amono, j_amono, data, row_offsets[nRowBlocks],
481  col_offsets[nColBlocks]);
482 }
483 
484 void BlockMatrix::PrintMatlab(std::ostream & os) const
485 {
486 
487  Vector row_data;
488  Array<int> row_ind;
489  int nnz_elem = NumNonZeroElems();
490  os<<"% size " << row_offsets.Last() << " " << col_offsets.Last() << "\n";
491  os<<"% Non Zeros " << nnz_elem << "\n";
492  int i, j;
493  std::ios::fmtflags old_fmt = os.flags();
494  os.setf(std::ios::scientific);
495  std::streamsize old_prec = os.precision(14);
496  for (i = 0; i < row_offsets.Last(); i++)
497  {
498  GetRow(i, row_ind, row_data);
499  for (j = 0; j < row_ind.Size(); j++)
500  {
501  os << i+1 << " " << row_ind[j]+1 << " " << row_data[j] << std::endl;
502  }
503  }
504 
505  os.precision(old_prec);
506  os.flags(old_fmt);
507 }
508 
510 {
511  BlockMatrix * At = new BlockMatrix(A.ColOffsets(), A.RowOffsets());
512  At->owns_blocks = 1;
513 
514  for (int irowAt = 0; irowAt < At->NumRowBlocks(); ++irowAt)
515  for (int jcolAt = 0; jcolAt < At->NumColBlocks(); ++jcolAt)
516  if (!A.IsZeroBlock(jcolAt, irowAt))
517  {
518  At->SetBlock(irowAt, jcolAt, Transpose(A.GetBlock(jcolAt, irowAt)));
519  }
520  return At;
521 }
522 
523 BlockMatrix * Mult(const BlockMatrix & A, const BlockMatrix & B)
524 {
525  BlockMatrix * C= new BlockMatrix(A.RowOffsets(), B.ColOffsets());
526  C->owns_blocks = 1;
527  Array<SparseMatrix *> CijPieces(A.NumColBlocks());
528 
529  for (int irowC = 0; irowC < A.NumRowBlocks(); ++irowC)
530  for (int jcolC = 0; jcolC < B.NumColBlocks(); ++jcolC)
531  {
532  CijPieces.SetSize(0, static_cast<SparseMatrix *>(NULL));
533  for (int k = 0; k < A.NumColBlocks(); ++k)
534  if (!A.IsZeroBlock(irowC, k) && !B.IsZeroBlock(k, jcolC))
535  {
536  CijPieces.Append(Mult(A.GetBlock(irowC, k), B.GetBlock(k, jcolC)));
537  }
538 
539  if (CijPieces.Size() > 1)
540  {
541  C->SetBlock(irowC, jcolC, Add(CijPieces));
542  for (SparseMatrix ** it = CijPieces.GetData();
543  it != CijPieces.GetData()+CijPieces.Size(); ++it)
544  {
545  delete *it;
546  }
547  }
548  else if (CijPieces.Size() == 1)
549  {
550  C->SetBlock(irowC, jcolC, CijPieces[0]);
551  }
552  }
553 
554  return C;
555 }
556 
557 }
int Size() const
Logical size of the array.
Definition: array.hpp:109
void PrintMatlab(std::ostream &os=std::cout) const
Export the monolithic matrix to file.
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 if the new size is different.
Definition: vector.hpp:263
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:451
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols.
Definition: operator.hpp:41
Abstract data type for sparse matrices.
Definition: matrix.hpp:69
T * GetData()
Returns the data.
Definition: array.hpp:91
BlockMatrix(const Array< int > &offsets)
Constructor for square block matrices.
Definition: blockmatrix.cpp:21
double * GetData() const
Definition: vector.hpp:90
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
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2809
void SetBlock(int i, int j, SparseMatrix *mat)
Set A(i,j) = mat.
Definition: blockmatrix.cpp:59
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
virtual void EliminateZeroRows()
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:35
int IsZeroBlock(int i, int j) const
Check if block (i,j) is a zero block.
Definition: blockmatrix.hpp:51
int NumRowBlocks() const
Return the number of row blocks.
Definition: blockmatrix.hpp:42
Array< int > & ColOffsets()
Return the columns offsets for block starts.
Definition: blockmatrix.hpp:55
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:391
Array< int > & RowOffsets()
Return the row offsets for block starts.
Definition: blockmatrix.hpp:53
void EliminateRowCol(Array< int > &ess_bc_dofs, Vector &sol, Vector &rhs)
Symmetric elimination of the marked degree of freedom.
virtual int NumNonZeroElems() const
Returns the total number of non zeros in the matrix.
virtual ~BlockMatrix()
Destructor.
Definition: blockmatrix.cpp:49
double * GetData() const
Return element data.
Definition: sparsemat.hpp:121
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 mfem_error(const char *msg)
Definition: error.cpp:23
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:331
int NumColBlocks() const
Return the number of column blocks.
Definition: blockmatrix.hpp:44
void SetDataAndSize(double *d, int s)
Definition: vector.hpp:71
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:80
T & Last()
Return the last element in the array.
Definition: array.hpp:411
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:474
Vector data type.
Definition: vector.hpp:33
int RowSize(const int i) const
Return the number of non zeros in row i.