 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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 /// This HYPRE library interface has been taken originally from MFEM and modified
13 /// to suit the needs for the MOAB library.
14 /// Modified by: Vijay Mahadevan
15
16 #ifndef MOAB_HYPREPARMATRIX
17 #define MOAB_HYPREPARMATRIX
18
19 #include "moab/MOABConfig.h"
20 #include "moab/Core.hpp"
21
22 #ifdef MOAB_HAVE_EIGEN3
23 #include <Eigen/Core>
24 #include <Eigen/Sparse>
25 #else
26 #error Configure with Eigen3 enabled
27 #endif
28
29 #ifdef MOAB_HAVE_MPI
30 #include "moab/ParallelComm.hpp"
31
32 // hypre header files
33 // #include "HYPRE.h"
34 // #include "HYPRE_parcsr_ls.h"
35 // #include "HYPRE_MatvecFunctions.h"
36
37 #include "seq_mv.h"
38 #include "HypreParVector.hpp"
39 #include "_hypre_parcsr_mv.h"
40 #include "_hypre_parcsr_ls.h"
41 #include "interpreter.h"
42 #include "HYPRE_MatvecFunctions.h"
43
44 #include "temp_multivector.h"
45
46 #ifdef HYPRE_COMPLEX
47 #error "MOAB does not work with HYPRE's complex numbers support"
48 #endif
49
50 #include "hypre_parcsr.hpp"
51
52 namespace moab
53 {
54
55 namespace internal
56 {
57
58 // Convert a HYPRE_Int to int
59 inline int to_int( HYPRE_Int i )
60 {
61 #ifdef HYPRE_BIGINT
62 MFEM_ASSERT( HYPRE_Int( int( i ) ) == i, "overflow converting HYPRE_Int to int" );
63 #endif
64 return int( i );
65 }
66
67 } // namespace internal
68
69 /// Wrapper for hypre's ParCSR matrix class
70 class HypreParMatrix
71 {
72 private:
73 /// The actual object
74 HYPRE_IJMatrix A;
75 hypre_ParCSRMatrix* A_parcsr;
76
77 mv_InterfaceInterpreter* interpreter;
78
79 // Does the object own the pointer A?
80 char ParCSROwner;
81
82 // Store the dimensions of the matrix
83 int height, gnrows;
84 int width, gncols;
85
86 moab::ParallelComm* pcomm;
87 char initialized;
88
89 // Initialize with defaults. Does not initialize inherited members.
90 void Init();
91
92 // Delete all owned data. Does not perform re-initialization with defaults.
93 void Destroy();
94
95 friend class HypreSolver;
96
97 public:
98 typedef Eigen::VectorXd Vector;
99 typedef Eigen::SparseMatrix< double, Eigen::RowMajor > MOABSparseMatrix;
100 // typedef template Eigen::ArrayXX<T> Array2D<T>;
101
102 public:
103 /// An empty matrix to be used as a reference to an existing matrix
104 HypreParMatrix( moab::ParallelComm* p_comm );
105
106 /// Converts hypre's format to HypreParMatrix
107 HypreParMatrix( HYPRE_IJMatrix a )
108 {
109 Init();
110 A = a;
111 height = GetNumRows();
112 width = GetNumCols();
113 }
114
115 /** Creates block-diagonal square parallel matrix. Diagonal is given by diag
116 which must be in CSR format (finalized). The new HypreParMatrix does not
117 take ownership of any of the input arrays. */
118 HypreParMatrix( moab::ParallelComm* p_comm, HYPRE_Int glob_size, HYPRE_Int* row_starts, HYPRE_Int nnz_pr_diag = 0 );
119
120 /** Creates block-diagonal rectangular parallel matrix. Diagonal is given by
121 diag which must be in CSR format (finalized). The new HypreParMatrix does
122 not take ownership of any of the input arrays. */
123 HypreParMatrix( moab::ParallelComm* p_comm,
124 HYPRE_Int global_num_rows,
125 HYPRE_Int global_num_cols,
126 HYPRE_Int* row_starts,
127 HYPRE_Int* col_starts,
128 HYPRE_Int nnz_pr_diag = 0,
129 HYPRE_Int onz_pr_diag = 0,
130 HYPRE_Int nnz_pr_offdiag = 0 );
131
132 /** Creates a general parallel matrix from a local CSR matrix on each
133 processor described by the I, J and data arrays. The local matrix should
134 be of size (local) nrows by (global) glob_ncols. The new parallel matrix
135 contains copies of all input arrays (so they can be deleted). */
136 // HypreParMatrix(MPI_Comm comm, int nrows, HYPRE_Int glob_nrows,
137 // HYPRE_Int glob_ncols, int *I, HYPRE_Int *J,
138 // double *data, HYPRE_Int *rows, HYPRE_Int *cols);
139
140 /// Make this HypreParMatrix a reference to 'master'
141 void MakeRef( const HypreParMatrix& master );
142
143 /// MPI communicator
144 moab::ParallelComm* GetParallelCommunicator() const
145 {
146 return pcomm;
147 }
148
149 /// Typecasting to hypre's HYPRE_IJMatrix*
150 operator HYPRE_IJMatrix()
151 {
152 return A;
153 }
154 /// Typecasting to hypre's HYPRE_ParCSRMatrix, a.k.a. void *
155 operator HYPRE_ParCSRMatrix()
156 {
157 return (HYPRE_ParCSRMatrix)A_parcsr;
158 }
159
160 /// Changes the ownership of the matrix
161 HYPRE_IJMatrix StealData();
162 /// Changes the ownership of the matrix to A
163 void StealData( HypreParMatrix& A );
164
165 /** If the HypreParMatrix does not own the row-starts array, make a copy of
166 it that the HypreParMatrix will own. If the col-starts array is the same
167 as the row-starts array, col-starts is also replaced. */
168 // void CopyRowStarts();
169 /** If the HypreParMatrix does not own the col-starts array, make a copy of
170 it that the HypreParMatrix will own. If the row-starts array is the same
171 as the col-starts array, row-starts is also replaced. */
172 // void CopyColStarts();
173
174 /// Returns the global number of nonzeros
175 inline HYPRE_Int NNZ()
176 {
177 return A_parcsr->num_nonzeros;
178 }
179 /// Returns the row partitioning
180 inline HYPRE_Int* RowPart()
181 {
182 return A->row_partitioning;
183 }
184 /// Returns the column partitioning
185 inline HYPRE_Int* ColPart()
186 {
187 return A->col_partitioning;
188 }
189 /// Returns the global number of rows
190 inline HYPRE_Int M()
191 {
192 return A->global_num_rows;
193 }
194 /// Returns the global number of columns
195 inline HYPRE_Int N()
196 {
197 return A->global_num_cols;
198 }
199 /// Returns the global number of rows
200 inline HYPRE_Int FirstRow()
201 {
202 return A->global_first_row;
203 }
204 /// Returns the global number of columns
205 inline HYPRE_Int FirstCol()
206 {
207 return A->global_first_col;
208 }
209
210 /// Get the local diagonal of the matrix.
211 void GetDiag( Vector& diag ) const;
212 /// Get the local diagonal block. NOTE: 'diag' will not own any data.
213 void GetDiag( MOABSparseMatrix& diag ) const;
214 /// Get the local off-diagonal block. NOTE: 'offd' will not own any data.
215 void GetOffd( MOABSparseMatrix& offd, HYPRE_Int*& cmap ) const;
216
217 /** Split the matrix into M x N equally sized blocks of parallel matrices.
218 The size of 'blocks' must already be set to M x N. */
219 // void GetBlocks(Eigen::Array<HypreParMatrix*,Eigen::Dynamic,Eigen::Dynamic> &blocks,
220 // bool interleaved_rows = false,
221 // bool interleaved_cols = false) const;
222
223 /// Returns the transpose of *this
224 // HypreParMatrix * Transpose();
225
226 /** Destroy and resize block-diagonal square parallel matrix. Diagonal is given by diag
227 which must be in CSR format (finalized). The new HypreParMatrix does not
228 take ownership of any of the input arrays. */
229 void resize( HYPRE_Int glob_size, HYPRE_Int* row_starts, HYPRE_Int nnz_pr_diag = 0, HYPRE_Int nnz_pr_offdiag = 0 );
230
231 /** Destroy and resize block-diagonal rectangular parallel matrix. Diagonal is given by
232 diag which must be in CSR format (finalized). The new HypreParMatrix does
233 not take ownership of any of the input arrays. */
234 void resize( HYPRE_Int global_num_rows,
235 HYPRE_Int global_num_cols,
236 HYPRE_Int* row_starts,
237 HYPRE_Int* col_starts,
238 HYPRE_Int* nnz_pr_diag = NULL,
239 HYPRE_Int* onz_pr_diag = NULL,
240 HYPRE_Int nnz_pr_offdiag = 0 );
241
242 /// Returns the number of rows in the diagonal block of the ParCSRMatrix
243 int GetNumRows() const
244 {
245 return internal::to_int( hypre_CSRMatrixNumRows( hypre_ParCSRMatrixDiag( A_parcsr ) ) );
246 }
247
248 /// Returns the number of columns in the diagonal block of the ParCSRMatrix
249 int GetNumCols() const
250 {
251 return internal::to_int( hypre_CSRMatrixNumCols( hypre_ParCSRMatrixDiag( A_parcsr ) ) );
252 }
253
254 /// Computes y = alpha * A * x + beta * y
255 HYPRE_Int Mult( HypreParVector& x, HypreParVector& y, double alpha = 1.0, double beta = 0.0 );
256 /// Computes y = alpha * A * x + beta * y
257 HYPRE_Int Mult( HYPRE_ParVector x, HYPRE_ParVector y, double alpha = 1.0, double beta = 0.0 );
258 /// Computes y = alpha * A^t * x + beta * y
259 HYPRE_Int MultTranspose( HypreParVector& x, HypreParVector& y, double alpha = 1.0, double beta = 0.0 );
260
261 void Mult( double a, const HypreParVector& x, double b, HypreParVector& y ) const;
262 void MultTranspose( double a, const HypreParVector& x, double b, HypreParVector& y ) const;
263
264 // virtual void Mult(const Vector &x, Vector &y) const
265 // { Mult(1.0, x, 0.0, y); }
266 // virtual void MultTranspose(const Vector &x, Vector &y) const
267 // { MultTranspose(1.0, x, 0.0, y); }
268
269 /** The "Boolean" analog of y = alpha * A * x + beta * y, where elements in
270 the sparsity pattern of the matrix are treated as "true". */
271 // void BooleanMult(int alpha, int *x, int beta, int *y)
272 // {
273 // internal::hypre_ParCSRMatrixBooleanMatvec(A_parcsr, alpha, x, beta, y);
274 // }
275
276 /** Multiply A on the left by a block-diagonal parallel matrix D. Return
277 a new parallel matrix, D*A. If D has a different number of rows than A,
278 D's row starts array needs to be given (as returned by the methods
279 GetDofOffsets/GetTrueDofOffsets of ParFiniteElementSpace). The new matrix
280 D*A uses copies of the row-, column-starts arrays, so "this" matrix and
281 "row_starts" can be deleted.
282 NOTE: this operation is local and does not require communication. */
283 // HypreParMatrix* LeftDiagMult(const MOABSparseMatrix &D,
284 // HYPRE_Int* row_starts = NULL) const;
285
286 /// Scale the local row i by s(i).
287 void ScaleRows( const Vector& s );
288 /// Scale the local row i by 1./s(i)
289 void InvScaleRows( const Vector& s );
290 /// Scale all entries by s: A_scaled = s*A.
291 void operator*=( double s );
292
293 /// Remove values smaller in absolute value than some threshold
294 void Threshold( double threshold = 0.0 );
295
296 /// If a row contains only zeros, set its diagonal to 1.
297 void EliminateZeroRows()
298 {
299 hypre_ParCSRMatrixFixZeroRows( A_parcsr );
300 }
301
302 /** Eliminate rows and columns from the matrix, and rows from the vector B.
303 Modify B with the BC values in X. */
304 void EliminateRowsCols( const std::vector< HYPRE_Int >& rows_cols, const HypreParVector& X, HypreParVector& B );
305
306 /** Eliminate rows and columns from the matrix and store the eliminated
307 elements in a new matrix Ae (returned), so that the modified matrix and
308 Ae sum to the original matrix. */
309 HypreParMatrix* EliminateRowsCols( const std::vector< HYPRE_Int >& rows_cols );
310
311 HYPRE_Int GetValues( const HYPRE_Int nrows,
312 HYPRE_Int* ncols,
313 HYPRE_Int* rows,
314 HYPRE_Int* cols,
315 HYPRE_Complex* values );
316 HYPRE_Int SetValues( const HYPRE_Int nrows,
317 HYPRE_Int* ncols,
318 const HYPRE_Int* rows,
319 const HYPRE_Int* cols,
320 const HYPRE_Complex* values );
321 HYPRE_Int AddToValues( const HYPRE_Int nrows,
322 HYPRE_Int* ncols,
323 const HYPRE_Int* rows,
324 const HYPRE_Int* cols,
325 const HYPRE_Complex* values );
326
327 HYPRE_Int FinalizeAssembly();
328
329 HYPRE_Int verbosity( const HYPRE_Int level );
330
331 /// Prints the locally owned rows in parallel
332 void Print( const char* fname, HYPRE_Int offi = 0, HYPRE_Int offj = 0 );
333 /// Reads the matrix from a file
334 void Read( const char* fname );
335
336 /// Calls hypre's destroy function
337 virtual ~HypreParMatrix()
338 {
339 Destroy();
340 }
341
342 /// Returns the matrix A * B
343 friend HypreParMatrix* ParMult( HypreParMatrix* A, HypreParMatrix* B );
344
345 /// Returns the matrix P^t * A * P
346 friend HypreParMatrix* RAP( HypreParMatrix* A, HypreParMatrix* P );
347 /// Returns the matrix Rt^t * A * P
348 friend HypreParMatrix* RAP( HypreParMatrix* Rt, HypreParMatrix* A, HypreParMatrix* P );
349
350 /** Eliminate essential BC specified by 'ess_dof_list' from the solution X to
351 the r.h.s. B. Here A is a matrix with eliminated BC, while Ae is such that
352 (A+Ae) is the original (Neumann) matrix before elimination. */
353 friend void EliminateBC( HypreParMatrix& A,
354 HypreParMatrix& Ae,
355 const std::vector< int >& ess_dof_list,
356 const HypreParVector& X,
357 HypreParVector& B );
358
359 friend class HypreSolver;
360 };
361
362 } // namespace moab
363
364 #endif // MOAB_HAVE_MPI
365
366 #endif
1.9.1.