12 #include "../general/array.hpp" 13 #include "../general/globals.hpp" 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)
39 nRowBlocks(row_offsets_.Size()-1),
40 nColBlocks(col_offsets_.Size()-1),
41 row_offsets(const_cast<
Array<int>& >(row_offsets_).GetData(),
43 col_offsets(const_cast<
Array<int>& >(col_offsets_).GetData(),
45 Aij(nRowBlocks, nColBlocks)
56 it != Aij.GetRow(0)+(Aij.NumRows()*Aij.NumCols()); ++it)
66 if (nRowBlocks <= i || nColBlocks <= j)
71 if (mat->
Height() != row_offsets[i+1] - row_offsets[i])
76 if (mat->
Width() != col_offsets[j+1] - col_offsets[j])
87 if (nRowBlocks <= i || nColBlocks <= j)
103 if (nRowBlocks <= i || nColBlocks <= j)
121 for (
int jcol = 0; jcol != nColBlocks; ++jcol)
123 for (
int irow = 0; irow != nRowBlocks; ++irow)
127 nnz_elem+= Aij(irow,jcol)->NumNonZeroElems();
140 findGlobalRow(i, iblock, iloc);
141 findGlobalCol(j, jblock, jloc);
148 return Aij(iblock, jblock)->Elem(iloc, jloc);
153 static const double zero = 0.0;
157 findGlobalRow(i, iblock, iloc);
158 findGlobalCol(j, jblock, jloc);
164 return static_cast<const SparseMatrix *
>(Aij(iblock, jblock))->
Elem(iloc, jloc);
172 findGlobalRow(i, iblock, iloc);
174 for (
int jblock = 0; jblock < nColBlocks; ++jblock)
176 if (Aij(iblock,jblock) != NULL)
178 rowsize += Aij(iblock,jblock)->RowSize(iloc);
187 int iblock, iloc, rowsize;
188 findGlobalRow(row, iblock, iloc);
196 int * it_cols = cols.
GetData();
197 double *it_srow = srow.
GetData();
199 for (
int jblock = 0; jblock < nColBlocks; ++jblock)
201 if (Aij(iblock,jblock) != NULL)
203 Aij(iblock,jblock)->GetRow(iloc, bcols, bsrow);
204 for (
int i = 0; i < bcols.
Size(); ++i)
206 *(it_cols++) = bcols[i] + col_offsets[jblock];
207 *(it_srow++) = bsrow(i);
219 for (iiblock = 0; iiblock < nRowBlocks; ++iiblock)
221 idx = rc - row_offsets[iiblock];
222 if (idx < 0 ) {
break; }
225 idx = rc - row_offsets[iiblock];
228 MFEM_ASSERT(nRowBlocks == nColBlocks,
229 "BlockMatrix::EliminateRowCol: nRowBlocks != nColBlocks");
231 MFEM_ASSERT(row_offsets[iiblock] == col_offsets[iiblock],
232 "BlockMatrix::EliminateRowCol: row_offsets[" 233 << iiblock <<
"] != col_offsets["<<iiblock<<
"]");
235 MFEM_ASSERT(Aij(iiblock, iiblock),
236 "BlockMatrix::EliminateRowCol: Null diagonal block");
239 for (
int jjblock = 0; jjblock < nRowBlocks; ++jjblock)
241 if (iiblock == jjblock) {
continue; }
242 if (Aij(iiblock,jjblock)) { Aij(iiblock,jjblock)->EliminateRow(idx); }
244 for (
int jjblock = 0; jjblock < nRowBlocks; ++jjblock)
246 if (iiblock == jjblock) {
continue; }
247 if (Aij(jjblock,iiblock)) { Aij(jjblock,iiblock)->EliminateCol(idx); }
249 Aij(iiblock, iiblock)->EliminateRowCol(idx,dpolicy);
255 if (nRowBlocks != nColBlocks)
257 mfem_error(
"BlockMatrix::EliminateRowCol: nRowBlocks != nColBlocks");
260 for (
int iiblock = 0; iiblock < nRowBlocks; ++iiblock)
262 if (row_offsets[iiblock] != col_offsets[iiblock])
264 mfem::out <<
"BlockMatrix::EliminateRowCol: row_offsets[" 265 << iiblock <<
"] != col_offsets["<<iiblock<<
"]\n";
272 Vector block_sol, block_rhs;
274 for (
int iiblock = 0; iiblock < nRowBlocks; ++iiblock)
276 int dsize = row_offsets[iiblock+1] - row_offsets[iiblock];
277 block_dofs.
MakeRef(ess_bc_dofs.
GetData()+row_offsets[iiblock], dsize);
281 if (Aij(iiblock, iiblock))
283 for (
int i = 0; i < block_dofs.
Size(); ++i)
287 Aij(iiblock, iiblock)->EliminateRowCol(i,block_sol(i), block_rhs);
293 for (
int i = 0; i < block_dofs.
Size(); ++i)
297 mfem_error(
"BlockMatrix::EliminateRowCol: Null diagonal block \n");
302 for (
int jjblock = 0; jjblock < nRowBlocks; ++jjblock)
304 if (jjblock != iiblock && Aij(iiblock, jjblock))
306 for (
int i = 0; i < block_dofs.
Size(); ++i)
310 Aij(iiblock, jjblock)->EliminateRow(i);
314 if (jjblock != iiblock && Aij(jjblock, iiblock))
317 row_offsets[jjblock+1] - row_offsets[jjblock]);
318 Aij(jjblock, iiblock)->EliminateCols(block_dofs, &block_sol, &block_rhs);
328 "BlockMatrix::EliminateRowCols: Elimination matrix pointer is null");
329 MFEM_VERIFY(nRowBlocks == nColBlocks,
330 "BlockMatrix::EliminateRowCols supported only for" 331 "nRowBlocks = nColBlocks");
333 std::vector<Array<int>> cols(nRowBlocks);
334 std::vector<Array<int>> rows(nRowBlocks);
337 for (
int k = 0; k < vdofs.
Size(); k++)
339 int vdof = (vdofs[k]) >=0 ? vdofs[k] : -1 - vdofs[k];
342 findGlobalCol(vdof,iblock,dof);
343 cols[iblock].Append(dof);
352 for (
int j = 0; j<nColBlocks; j++)
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++)
360 if (i == j) {
continue; }
363 for (
int k = 0; k < cols[j].Size(); k++)
374 MFEM_VERIFY(nRowBlocks == nColBlocks,
"not a square matrix");
376 for (
int iblock = 0; iblock < nRowBlocks; ++iblock)
378 if (Aij(iblock,iblock))
381 for (
int i = 0; i < Aij(iblock, iblock)->NumRows(); ++i)
384 for (
int jblock = 0; jblock < nColBlocks; ++jblock)
385 if (Aij(iblock,jblock))
387 norm += Aij(iblock,jblock)->GetRowNorml1(i);
390 if (norm <= threshold)
392 for (
int jblock = 0; jblock < nColBlocks; ++jblock)
394 if (Aij(iblock,jblock))
396 Aij(iblock,jblock)->EliminateRow(
406 for (
int i = 0; i < row_offsets[iblock+1] - row_offsets[iblock]; ++i)
409 for (
int jblock = 0; jblock < nColBlocks; ++jblock)
411 if (Aij(iblock,jblock))
413 norm += Aij(iblock,jblock)->GetRowNorml1(i);
417 MFEM_VERIFY(!(norm <= threshold),
"diagonal block is NULL:" 418 " iblock = " << iblock <<
", i = " << i <<
", norm = " 427 for (
int iblock = 0; iblock < nRowBlocks; ++iblock)
429 for (
int jblock = 0; jblock < nColBlocks; ++jblock)
431 if (!Aij(iblock,jblock)) {
continue; }
432 if (!Aij(iblock,jblock)->Finalized())
434 Aij(iblock,jblock)->Finalize(skip_zeros, fix_empty_rows);
444 mfem_error(
"Error: x and y can't point to the same data \n");
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 <<
")");
459 mfem_error(
"Error: x and y can't point to the same data \n");
462 Vector xblockview, yblockview;
464 for (
int iblock = 0; iblock != nRowBlocks; ++iblock)
467 row_offsets[iblock+1] - row_offsets[iblock]);
469 for (
int jblock = 0; jblock != nColBlocks; ++jblock)
471 if (Aij(iblock, jblock) != NULL)
474 x.
GetData() + col_offsets[jblock],
475 col_offsets[jblock+1] - col_offsets[jblock]);
477 Aij(iblock, jblock)->AddMult(xblockview, yblockview, val);
487 mfem_error(
"Error: x and y can't point to the same data \n");
495 const double val)
const 499 mfem_error(
"Error: x and y can't point to the same data \n");
502 Vector xblockview, yblockview;
504 for (
int iblock = 0; iblock != nColBlocks; ++iblock)
507 col_offsets[iblock+1] - col_offsets[iblock]);
509 for (
int jblock = 0; jblock != nRowBlocks; ++jblock)
511 if (Aij(jblock, iblock) != NULL)
514 x.
GetData() + row_offsets[jblock],
515 row_offsets[jblock+1] - row_offsets[jblock]);
517 Aij(jblock, iblock)->AddMultTranspose(xblockview, yblockview, val);
528 for (
int i = 0; i<rows.
Size(); i++)
530 int dof = (rows[i]>=0) ? rows[i] : -1-rows[i];
534 for (
int k = 0; k <cols.
Size(); k++)
536 s += srow[k] * x[cols[k]];
543 const double a)
const 547 for (
int i = 0; i<rows.
Size(); i++)
549 int dof = (rows[i]>=0) ? rows[i] : -1-rows[i];
553 for (
int k = 0; k <cols.
Size(); k++)
555 s += srow[k] * x[cols[k]];
566 int * i_amono =
Memory<int>(row_offsets[nRowBlocks]+2);
570 for (
int i = 0; i < row_offsets[nRowBlocks]+2; i++)
575 int * i_amono_construction = i_amono+1;
577 for (
int iblock = 0; iblock != nRowBlocks; ++iblock)
579 for (
int irow(row_offsets[iblock]); irow < row_offsets[iblock+1]; ++irow)
581 int local_row = irow - row_offsets[iblock];
582 int ind = i_amono_construction[irow];
583 for (
int jblock = 0; jblock < nColBlocks; ++jblock)
585 if (Aij(iblock,jblock) != NULL)
586 ind += Aij(iblock, jblock)->GetI()[local_row+1]
587 - Aij(iblock, jblock)->GetI()[local_row];
589 i_amono_construction[irow+1] = ind;
594 for (
int iblock = 0; iblock != nRowBlocks; ++iblock)
596 for (
int jblock = 0; jblock != nColBlocks; ++jblock)
598 if (Aij(iblock,jblock) != NULL)
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];
606 int loc_start_index = 0;
607 int loc_end_index = 0;
608 int glob_start_index = 0;
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)
613 glob_start_index = *i_it;
616 if (glob_start_index > nnz)
618 mfem::out<<
"glob_start_index = " << glob_start_index <<
"\n";
619 mfem::out<<
"Block:" << iblock <<
" " << jblock <<
"\n";
623 loc_end_index = *(i_it_aij);
624 for (
int cnt = 0; cnt < loc_end_index-loc_start_index; cnt++)
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;
630 *i_it += loc_end_index-loc_start_index;
632 loc_start_index = loc_end_index;
638 return new SparseMatrix(i_amono, j_amono, data, row_offsets[nRowBlocks],
639 col_offsets[nColBlocks]);
648 os<<
"% size " << row_offsets.
Last() <<
" " << col_offsets.
Last() <<
"\n";
649 os<<
"% Non Zeros " << nnz_elem <<
"\n";
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++)
656 GetRow(i, row_ind, row_data);
657 for (j = 0; j < row_ind.
Size(); j++)
659 os << i+1 <<
" " << row_ind[j]+1 <<
" " << row_data[j] << std::endl;
663 os << row_offsets.
Last() <<
" " << col_offsets.
Last () <<
" 0.0\n";
665 os.precision(old_prec);
674 for (
int irowAt = 0; irowAt < At->
NumRowBlocks(); ++irowAt)
676 for (
int jcolAt = 0; jcolAt < At->
NumColBlocks(); ++jcolAt)
697 CijPieces.
SetSize(0, static_cast<SparseMatrix *>(NULL));
706 if (CijPieces.Size() > 1)
710 it != CijPieces.GetData()+CijPieces.Size(); ++it)
715 else if (CijPieces.Size() == 1)
717 C->
SetBlock(irowC, jcolC, CijPieces[0]);
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.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
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().
Abstract data type for sparse matrices.
T * GetData()
Returns the data.
int Size() const
Returns the size of the vector.
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.
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.
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
void SetBlock(int i, int j, SparseMatrix *mat)
Set A(i,j) = mat.
int NumRowBlocks() const
Return the number of row blocks.
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
virtual void Finalize(int skip_zeros=1)
Finalize all the submatrices.
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
virtual void MultTranspose(const Vector &x, Vector &y) const
MatrixTranspose-Vector Multiplication y = A'*x.
Array< int > & ColOffsets()
Return the columns offsets for block starts.
int NumColBlocks() const
Return the number of column blocks.
Array< int > & RowOffsets()
Return the row offsets for block starts.
Set the diagonal value to one.
virtual ~BlockMatrix()
Destructor.
double * GetData() const
Return a pointer to the beginning of the Vector data.
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...
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
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.
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
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().
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
int height
Dimension of the output / number of rows in the matrix.
void EliminateRow(int row, const double sol, Vector &rhs)
Eliminates a column from the transpose matrix.
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
void EliminateCols(const Array< int > &cols, const Vector *x=NULL, Vector *b=NULL)
Eliminate all columns i for which cols[i] != 0.
SparseMatrix & GetBlock(int i, int j)
Return a reference to block (i,j). Reference may be invalid if Aij(i,j) == NULL.
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.
T & Last()
Return the last element in the array.
void MakeRef(T *, int)
Make this Array a reference to a pointer.
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'*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.
int width
Dimension of the input / number of columns in the matrix.
int IsZeroBlock(int i, int j) const
Check if block (i,j) is a zero block.