MFEM  v4.6.0
Finite element discretization library
blockoperator.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 
13 #include "../general/array.hpp"
14 #include "operator.hpp"
15 #include "blockvector.hpp"
16 #include "blockoperator.hpp"
17 
18 namespace mfem
19 {
20 
22  : Operator(offsets.Last()),
23  owns_blocks(0),
24  nRowBlocks(offsets.Size() - 1),
25  nColBlocks(offsets.Size() - 1),
26  row_offsets(0),
27  col_offsets(0),
28  op(nRowBlocks, nRowBlocks),
29  coef(nRowBlocks, nColBlocks)
30 {
31  op = static_cast<Operator *>(NULL);
32  row_offsets.MakeRef(offsets);
33  col_offsets.MakeRef(offsets);
34 }
35 
37  const Array<int> & col_offsets_)
38  : Operator(row_offsets_.Last(), col_offsets_.Last()),
39  owns_blocks(0),
40  nRowBlocks(row_offsets_.Size()-1),
41  nColBlocks(col_offsets_.Size()-1),
42  row_offsets(0),
43  col_offsets(0),
44  op(nRowBlocks, nColBlocks),
45  coef(nRowBlocks, nColBlocks)
46 {
47  op = static_cast<Operator *>(NULL);
48  row_offsets.MakeRef(row_offsets_);
49  col_offsets.MakeRef(col_offsets_);
50 }
51 
52 void BlockOperator::SetDiagonalBlock(int iblock, Operator *opt, double c)
53 {
54  SetBlock(iblock, iblock, opt, c);
55 }
56 
57 void BlockOperator::SetBlock(int iRow, int iCol, Operator *opt, double c)
58 {
59  if (owns_blocks && op(iRow, iCol))
60  {
61  delete op(iRow, iCol);
62  }
63  op(iRow, iCol) = opt;
64  coef(iRow, iCol) = c;
65 
66  MFEM_VERIFY(row_offsets[iRow+1] - row_offsets[iRow] == opt->NumRows() &&
67  col_offsets[iCol+1] - col_offsets[iCol] == opt->NumCols(),
68  "incompatible Operator dimensions");
69 }
70 
71 // Operator application
72 void BlockOperator::Mult (const Vector & x, Vector & y) const
73 {
74  MFEM_ASSERT(x.Size() == width, "incorrect input Vector size");
75  MFEM_ASSERT(y.Size() == height, "incorrect output Vector size");
76 
77  x.Read();
78  y.Write(); y = 0.0;
79 
80  xblock.Update(const_cast<Vector&>(x),col_offsets);
81  yblock.Update(y,row_offsets);
82 
83  for (int iRow=0; iRow < nRowBlocks; ++iRow)
84  {
85  tmp.SetSize(row_offsets[iRow+1] - row_offsets[iRow]);
86  for (int jCol=0; jCol < nColBlocks; ++jCol)
87  {
88  if (op(iRow,jCol))
89  {
90  op(iRow,jCol)->Mult(xblock.GetBlock(jCol), tmp);
91  yblock.GetBlock(iRow).Add(coef(iRow,jCol), tmp);
92  }
93  }
94  }
95 
96  for (int iRow=0; iRow < nRowBlocks; ++iRow)
97  {
98  yblock.GetBlock(iRow).SyncAliasMemory(y);
99  }
100 }
101 
102 // Action of the transpose operator
103 void BlockOperator::MultTranspose (const Vector & x, Vector & y) const
104 {
105  MFEM_ASSERT(x.Size() == height, "incorrect input Vector size");
106  MFEM_ASSERT(y.Size() == width, "incorrect output Vector size");
107 
108  x.Read();
109  y.Write(); y = 0.0;
110 
111  xblock.Update(const_cast<Vector&>(x),row_offsets);
112  yblock.Update(y,col_offsets);
113 
114  for (int iRow=0; iRow < nColBlocks; ++iRow)
115  {
116  tmp.SetSize(col_offsets[iRow+1] - col_offsets[iRow]);
117  for (int jCol=0; jCol < nRowBlocks; ++jCol)
118  {
119  if (op(jCol,iRow))
120  {
121  op(jCol,iRow)->MultTranspose(xblock.GetBlock(jCol), tmp);
122  yblock.GetBlock(iRow).Add(coef(jCol,iRow), tmp);
123  }
124  }
125  }
126 
127  for (int iRow=0; iRow < nColBlocks; ++iRow)
128  {
129  yblock.GetBlock(iRow).SyncAliasMemory(y);
130  }
131 }
132 
134 {
135  if (owns_blocks)
136  {
137  for (int iRow=0; iRow < nRowBlocks; ++iRow)
138  {
139  for (int jCol=0; jCol < nColBlocks; ++jCol)
140  {
141  delete op(iRow,jCol);
142  }
143  }
144  }
145 }
146 
147 //-----------------------------------------------------------------------
149  const Array<int> & offsets_):
150  Solver(offsets_.Last()),
151  owns_blocks(0),
152  nBlocks(offsets_.Size() - 1),
153  offsets(0),
154  ops(nBlocks)
155 {
156  ops = static_cast<Operator *>(NULL);
157  offsets.MakeRef(offsets_);
158 }
159 
161 {
162  MFEM_VERIFY(offsets[iblock+1] - offsets[iblock] == op->Height() &&
163  offsets[iblock+1] - offsets[iblock] == op->Width(),
164  "incompatible Operator dimensions");
165 
166  if (owns_blocks && ops[iblock])
167  {
168  delete ops[iblock];
169  }
170  ops[iblock] = op;
171 }
172 
173 // Operator application
175 {
176  MFEM_ASSERT(x.Size() == width, "incorrect input Vector size");
177  MFEM_ASSERT(y.Size() == height, "incorrect output Vector size");
178 
179  x.Read();
180  y.Write(); y = 0.0;
181 
182  xblock.Update(const_cast<Vector&>(x),offsets);
183  yblock.Update(y,offsets);
184 
185  for (int i=0; i<nBlocks; ++i)
186  {
187  if (ops[i])
188  {
189  ops[i]->Mult(xblock.GetBlock(i), yblock.GetBlock(i));
190  }
191  else
192  {
193  yblock.GetBlock(i) = xblock.GetBlock(i);
194  }
195  }
196 
197  for (int i=0; i<nBlocks; ++i)
198  {
199  yblock.GetBlock(i).SyncAliasMemory(y);
200  }
201 }
202 
203 // Action of the transpose operator
205  Vector & y) const
206 {
207  MFEM_ASSERT(x.Size() == height, "incorrect input Vector size");
208  MFEM_ASSERT(y.Size() == width, "incorrect output Vector size");
209 
210  x.Read();
211  y.Write(); y = 0.0;
212 
213  xblock.Update(const_cast<Vector&>(x),offsets);
214  yblock.Update(y,offsets);
215 
216  for (int i=0; i<nBlocks; ++i)
217  {
218  if (ops[i])
219  {
220  (ops[i])->MultTranspose(xblock.GetBlock(i), yblock.GetBlock(i));
221  }
222  else
223  {
224  yblock.GetBlock(i) = xblock.GetBlock(i);
225  }
226  }
227 
228  for (int i=0; i<nBlocks; ++i)
229  {
230  yblock.GetBlock(i).SyncAliasMemory(y);
231  }
232 }
233 
235 {
236  if (owns_blocks)
237  {
238  for (int i=0; i<nBlocks; ++i)
239  {
240  delete ops[i];
241  }
242  }
243 }
244 
246  const Array<int> & offsets_)
247  : Solver(offsets_.Last()),
248  owns_blocks(0),
249  nBlocks(offsets_.Size() - 1),
250  offsets(0),
251  ops(nBlocks, nBlocks)
252 {
253  ops = static_cast<Operator *>(NULL);
254  offsets.MakeRef(offsets_);
255 }
256 
258  Operator *op)
259 {
260  MFEM_VERIFY(offsets[iblock+1] - offsets[iblock] == op->Height() &&
261  offsets[iblock+1] - offsets[iblock] == op->Width(),
262  "incompatible Operator dimensions");
263 
264  SetBlock(iblock, iblock, op);
265 }
266 
268  Operator *op)
269 {
270  MFEM_VERIFY(iRow >= iCol,"cannot set block in upper triangle");
271  MFEM_VERIFY(offsets[iRow+1] - offsets[iRow] == op->NumRows() &&
272  offsets[iCol+1] - offsets[iCol] == op->NumCols(),
273  "incompatible Operator dimensions");
274 
275  ops(iRow, iCol) = op;
276 }
277 
278 // Operator application
280  Vector & y) const
281 {
282  MFEM_ASSERT(x.Size() == width, "incorrect input Vector size");
283  MFEM_ASSERT(y.Size() == height, "incorrect output Vector size");
284 
285  yblock.Update(y.GetData(),offsets);
286  xblock.Update(x.GetData(),offsets);
287 
288  y = 0.0;
289  for (int iRow=0; iRow < nBlocks; ++iRow)
290  {
291  tmp.SetSize(offsets[iRow+1] - offsets[iRow]);
292  tmp2.SetSize(offsets[iRow+1] - offsets[iRow]);
293  tmp2 = 0.0;
294  tmp2 += xblock.GetBlock(iRow);
295  for (int jCol=0; jCol < iRow; ++jCol)
296  {
297  if (ops(iRow,jCol))
298  {
299  ops(iRow,jCol)->Mult(yblock.GetBlock(jCol), tmp);
300  tmp2 -= tmp;
301  }
302  }
303  if (ops(iRow,iRow))
304  {
305  ops(iRow,iRow)->Mult(tmp2, yblock.GetBlock(iRow));
306  }
307  else
308  {
309  yblock.GetBlock(iRow) = tmp2;
310  }
311  }
312 }
313 
314 // Action of the transpose operator
316  Vector & y) const
317 {
318  MFEM_ASSERT(x.Size() == height, "incorrect input Vector size");
319  MFEM_ASSERT(y.Size() == width, "incorrect output Vector size");
320 
321  yblock.Update(y.GetData(),offsets);
322  xblock.Update(x.GetData(),offsets);
323 
324  y = 0.0;
325  for (int iRow=nBlocks-1; iRow >=0; --iRow)
326  {
327  tmp.SetSize(offsets[iRow+1] - offsets[iRow]);
328  tmp2.SetSize(offsets[iRow+1] - offsets[iRow]);
329  tmp2 = 0.0;
330  tmp2 += xblock.GetBlock(iRow);
331  for (int jCol=iRow+1; jCol < nBlocks; ++jCol)
332  {
333  if (ops(jCol,iRow))
334  {
335  ops(jCol,iRow)->MultTranspose(yblock.GetBlock(jCol), tmp);
336  tmp2 -= tmp;
337  }
338  }
339  if (ops(iRow,iRow))
340  {
341  ops(iRow,iRow)->MultTranspose(tmp2, yblock.GetBlock(iRow));
342  }
343  else
344  {
345  yblock.GetBlock(iRow) = tmp2;
346  }
347  }
348 }
349 
351 {
352  if (owns_blocks)
353  {
354  for (int iRow=0; iRow < nBlocks; ++iRow)
355  {
356  for (int jCol=0; jCol < nBlocks; ++jCol)
357  {
358  delete ops(jCol,iRow);
359  }
360  }
361  }
362 }
363 
364 }
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
void Update(double *data, const Array< int > &bOffsets)
Update method.
Definition: blockvector.cpp:95
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:453
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
int NumRows() const
Get the number of rows (size of output) of the Operator. Synonym with Height().
Definition: operator.hpp:69
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:461
void SetBlock(int iRow, int iCol, Operator *op)
Add a block opt in the block-entry (iblock, jblock).
BlockDiagonalPreconditioner(const Array< int > &offsets)
Constructor that specifies the block structure.
BlockOperator(const Array< int > &offsets)
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:206
void SetDiagonalBlock(int iblock, Operator *op, double c=1.0)
Add block op in the block-entry (iblock, iblock).
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
Definition: operator.hpp:75
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
void SetDiagonalBlock(int iblock, Operator *op)
Add block op in the block-entry (iblock, iblock).
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:248
void SyncAliasMemory(const Vector &v) const
Update the alias memory location of the vector to match v.
Definition: vector.hpp:238
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:872
Vector data type.
Definition: vector.hpp:58
BlockLowerTriangularPreconditioner(const Array< int > &offsets)
Base class for solvers.
Definition: operator.hpp:682
Abstract operator.
Definition: operator.hpp:24
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
void SetBlock(int iRow, int iCol, Operator *op, double c=1.0)
Add a block op in the block-entry (iblock, jblock).
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
Vector & GetBlock(int i)
Get the i-th vector in the block.
Definition: blockvector.hpp:93
virtual void Mult(const Vector &x, Vector &y) const
Operator application.