MFEM  v3.4
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
blockoperator.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 
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 *op, double c)
53 {
54  SetBlock(iblock, iblock, op, c);
55 }
56 
57 void BlockOperator::SetBlock(int iRow, int iCol, Operator *opt, double c)
58 {
59  op(iRow, iCol) = opt;
60  coef(iRow, iCol) = c;
61 
62  MFEM_VERIFY(row_offsets[iRow+1] - row_offsets[iRow] == opt->NumRows() &&
63  col_offsets[iCol+1] - col_offsets[iCol] == opt->NumCols(),
64  "incompatible Operator dimensions");
65 }
66 
67 // Operator application
68 void BlockOperator::Mult (const Vector & x, Vector & y) const
69 {
70  MFEM_ASSERT(x.Size() == width, "incorrect input Vector size");
71  MFEM_ASSERT(y.Size() == height, "incorrect output Vector size");
72 
73  yblock.Update(y.GetData(),row_offsets);
74  xblock.Update(x.GetData(),col_offsets);
75 
76  y = 0.0;
77  for (int iRow=0; iRow < nRowBlocks; ++iRow)
78  {
79  tmp.SetSize(row_offsets[iRow+1] - row_offsets[iRow]);
80  for (int jCol=0; jCol < nColBlocks; ++jCol)
81  {
82  if (op(iRow,jCol))
83  {
84  op(iRow,jCol)->Mult(xblock.GetBlock(jCol), tmp);
85  yblock.GetBlock(iRow).Add(coef(iRow,jCol), tmp);
86  }
87  }
88  }
89 }
90 
91 // Action of the transpose operator
92 void BlockOperator::MultTranspose (const Vector & x, Vector & y) const
93 {
94  MFEM_ASSERT(x.Size() == height, "incorrect input Vector size");
95  MFEM_ASSERT(y.Size() == width, "incorrect output Vector size");
96 
97  y = 0.0;
98 
99  xblock.Update(x.GetData(),row_offsets);
100  yblock.Update(y.GetData(),col_offsets);
101 
102  for (int iRow=0; iRow < nColBlocks; ++iRow)
103  {
104  tmp.SetSize(col_offsets[iRow+1] - col_offsets[iRow]);
105  for (int jCol=0; jCol < nRowBlocks; ++jCol)
106  {
107  if (op(jCol,iRow))
108  {
109  op(jCol,iRow)->MultTranspose(xblock.GetBlock(jCol), tmp);
110  yblock.GetBlock(iRow).Add(coef(jCol,iRow), tmp);
111  }
112  }
113  }
114 
115 }
116 
118 {
119  if (owns_blocks)
120  for (int iRow=0; iRow < nRowBlocks; ++iRow)
121  for (int jCol=0; jCol < nColBlocks; ++jCol)
122  {
123  delete op(jCol,iRow);
124  }
125 }
126 
127 //-----------------------------------------------------------------------
129  const Array<int> & offsets_):
130  Solver(offsets_.Last()),
131  owns_blocks(0),
132  nBlocks(offsets_.Size() - 1),
133  offsets(0),
134  op(nBlocks)
135 
136 {
137  op = static_cast<Operator *>(NULL);
138  offsets.MakeRef(offsets_);
139 }
140 
142 {
143  MFEM_VERIFY(offsets[iblock+1] - offsets[iblock] == opt->Height() &&
144  offsets[iblock+1] - offsets[iblock] == opt->Width(),
145  "incompatible Operator dimensions");
146 
147  op[iblock] = opt;
148 }
149 
150 // Operator application
152 {
153  MFEM_ASSERT(x.Size() == width, "incorrect input Vector size");
154  MFEM_ASSERT(y.Size() == height, "incorrect output Vector size");
155 
156  yblock.Update(y.GetData(), offsets);
157  xblock.Update(x.GetData(), offsets);
158 
159  for (int i=0; i<nBlocks; ++i)
160  if (op[i])
161  {
162  op[i]->Mult(xblock.GetBlock(i), yblock.GetBlock(i));
163  }
164  else
165  {
166  yblock.GetBlock(i) = xblock.GetBlock(i);
167  }
168 }
169 
170 // Action of the transpose operator
172  Vector & y) const
173 {
174  MFEM_ASSERT(x.Size() == height, "incorrect input Vector size");
175  MFEM_ASSERT(y.Size() == width, "incorrect output Vector size");
176 
177  yblock.Update(y.GetData(), offsets);
178  xblock.Update(x.GetData(), offsets);
179 
180  for (int i=0; i<nBlocks; ++i)
181  if (op[i])
182  {
183  (op[i])->MultTranspose(xblock.GetBlock(i), yblock.GetBlock(i));
184  }
185  else
186  {
187  yblock.GetBlock(i) = xblock.GetBlock(i);
188  }
189 }
190 
192 {
193  if (owns_blocks)
194  for (int i=0; i<nBlocks; ++i)
195  {
196  delete op[i];
197  }
198 }
199 
201  const Array<int> & offsets_)
202  : Solver(offsets_.Last()),
203  owns_blocks(0),
204  nBlocks(offsets_.Size() - 1),
205  offsets(0),
206  op(nBlocks, nBlocks)
207 {
208  op = static_cast<Operator *>(NULL);
209  offsets.MakeRef(offsets_);
210 }
211 
213  Operator *op)
214 {
215  MFEM_VERIFY(offsets[iblock+1] - offsets[iblock] == op->Height() &&
216  offsets[iblock+1] - offsets[iblock] == op->Width(),
217  "incompatible Operator dimensions");
218 
219  SetBlock(iblock, iblock, op);
220 }
221 
223  Operator *opt)
224 {
225  MFEM_VERIFY(iRow >= iCol,"cannot set block in upper triangle");
226  MFEM_VERIFY(offsets[iRow+1] - offsets[iRow] == opt->NumRows() &&
227  offsets[iCol+1] - offsets[iCol] == opt->NumCols(),
228  "incompatible Operator dimensions");
229 
230  op(iRow, iCol) = opt;
231 }
232 
233 // Operator application
235  Vector & y) const
236 {
237  MFEM_ASSERT(x.Size() == width, "incorrect input Vector size");
238  MFEM_ASSERT(y.Size() == height, "incorrect output Vector size");
239 
240  yblock.Update(y.GetData(),offsets);
241  xblock.Update(x.GetData(),offsets);
242 
243  y = 0.0;
244  for (int iRow=0; iRow < nBlocks; ++iRow)
245  {
246  tmp.SetSize(offsets[iRow+1] - offsets[iRow]);
247  tmp2.SetSize(offsets[iRow+1] - offsets[iRow]);
248  tmp2 = 0.0;
249  tmp2 += xblock.GetBlock(iRow);
250  for (int jCol=0; jCol < iRow; ++jCol)
251  {
252  if (op(iRow,jCol))
253  {
254  op(iRow,jCol)->Mult(yblock.GetBlock(jCol), tmp);
255  tmp2 -= tmp;
256  }
257  }
258  if (op(iRow,iRow))
259  {
260  op(iRow,iRow)->Mult(tmp2, yblock.GetBlock(iRow));
261  }
262  else
263  {
264  yblock.GetBlock(iRow) = tmp2;
265  }
266  }
267 }
268 
269 // Action of the transpose operator
271  Vector & y) const
272 {
273  MFEM_ASSERT(x.Size() == height, "incorrect input Vector size");
274  MFEM_ASSERT(y.Size() == width, "incorrect output Vector size");
275 
276  yblock.Update(y.GetData(),offsets);
277  xblock.Update(x.GetData(),offsets);
278 
279  y = 0.0;
280  for (int iRow=nBlocks-1; iRow >=0; --iRow)
281  {
282  tmp.SetSize(offsets[iRow+1] - offsets[iRow]);
283  tmp2.SetSize(offsets[iRow+1] - offsets[iRow]);
284  tmp2 = 0.0;
285  tmp2 += xblock.GetBlock(iRow);
286  for (int jCol=iRow+1; jCol < nBlocks; ++jCol)
287  {
288  if (op(jCol,iRow))
289  {
290  op(jCol,iRow)->MultTranspose(yblock.GetBlock(jCol), tmp);
291  tmp2 -= tmp;
292  }
293  }
294  if (op(iRow,iRow))
295  {
296  op(iRow,iRow)->MultTranspose(tmp2, yblock.GetBlock(iRow));
297  }
298  else
299  {
300  yblock.GetBlock(iRow) = tmp2;
301  }
302  }
303 }
304 
306 {
307  if (owns_blocks)
308  {
309  for (int iRow=0; iRow < nBlocks; ++iRow)
310  {
311  for (int jCol=0; jCol < nBlocks; ++jCol)
312  {
313  delete op(jCol,iRow);
314  }
315  }
316  }
317 }
318 
319 }
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:328
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
int Size() const
Returns the size of the vector.
Definition: vector.hpp:120
void Update(double *data, const Array< int > &bOffsets)
Update method.
Definition: blockvector.cpp:67
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:129
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
void SetBlock(int iRow, int iCol, Operator *op)
Add a block op in the block-entry (iblock, jblock).
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
BlockDiagonalPreconditioner(const Array< int > &offsets)
Constructor that specifies the block structure.
BlockOperator(const Array< int > &offsets)
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
void SetDiagonalBlock(int iblock, Operator *op, double c=1.0)
Add block op in the block-entry (iblock, iblock).
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
int NumRows() const
Get the number of rows (size of output) of the Operator. Synonym with Height().
Definition: operator.hpp:39
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
Definition: operator.hpp:45
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:24
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:201
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:720
Vector data type.
Definition: vector.hpp:48
BlockLowerTriangularPreconditioner(const Array< int > &offsets)
Base class for solvers.
Definition: operator.hpp:268
Abstract operator.
Definition: operator.hpp:21
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:25
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:77