 |
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 #include "HypreParVector.hpp"
13 #include "HypreParMatrix.hpp"
14
15 #ifdef MOAB_HAVE_MPI
16
17 #include "hypre_parcsr.hpp"
18
19 #include <fstream>
20 #include <iostream>
21 #include <cmath>
22 #include <cstdlib>
23 #include <cassert>
24
25 using namespace std;
26
27 #define MOAB_DEBUG
28
29 namespace moab
30 {
31 #define moab_hypre_assert( a, b ) \
32 { \
33 }
34 #define moab_hypre_assert_t( a, b ) \
35 { \
36 if( !a ) \
37 { \
38 std::cout << "HYPRE Error: " << b << std::endl; \
39 exit( -1 ); \
40 } \
41 }
42
43 template < typename TargetT, typename SourceT >
44 static TargetT* DuplicateAs( const SourceT* array, int size, bool cplusplus = true )
45 {
46 TargetT* target_array = cplusplus ? new TargetT[size]
47 /* */
48 : hypre_TAlloc( TargetT, size );
49
50 for( int i = 0; i < size; i++ )
51 {
52 target_array[i] = array[i];
53 }
54
55 return target_array;
56 }
57
58 void HypreParMatrix::Init()
59 {
60 A = NULL;
61 A_parcsr = NULL;
62 ParCSROwner = 1;
63 }
64
65 HypreParMatrix::HypreParMatrix( moab::ParallelComm* p_comm ) : pcomm( p_comm )
66 {
67 Init();
68 height = width = 0;
69 }
70
71 // Square block-diagonal constructor
72 HypreParMatrix::HypreParMatrix( moab::ParallelComm* p_comm,
73 HYPRE_Int glob_size,
74 HYPRE_Int* row_starts,
75 HYPRE_Int nnz_pr )
76 : pcomm( p_comm )
77 {
78 Init();
79 resize( glob_size, row_starts, nnz_pr );
80 }
81
82 void HypreParMatrix::resize( HYPRE_Int glob_size,
83 HYPRE_Int* row_starts,
84 HYPRE_Int nnz_pr_diag,
85 HYPRE_Int nnz_pr_offdiag )
86 {
87 /* Create the matrix.
88 Note that this is a square matrix, so we indicate the row partition
89 size twice (since number of rows = number of cols) */
90 HYPRE_IJMatrixCreate( pcomm->comm(), row_starts[0], row_starts[1], row_starts[0], row_starts[1], &A );
91 /* Choose a parallel csr format storage (see the User's Manual) */
92 HYPRE_IJMatrixSetObjectType( A, HYPRE_PARCSR );
93
94 int locsize = row_starts[1] - row_starts[0] + 1;
95 int num_procs, rank;
96 MPI_Comm_size( pcomm->comm(), &num_procs );
97 MPI_Comm_rank( pcomm->comm(), &rank );
98
99 if( nnz_pr_diag != 0 || nnz_pr_offdiag != 0 )
100 {
101 if( num_procs > 1 )
102 {
103 std::vector< HYPRE_Int > m_nnz_pr_diag( locsize, nnz_pr_diag );
104 std::vector< HYPRE_Int > m_onz_pr_diag( locsize, nnz_pr_offdiag );
105 HYPRE_IJMatrixSetDiagOffdSizes( A, &m_nnz_pr_diag[0], &m_onz_pr_diag[0] );
106 }
107 else
108 {
109 std::vector< HYPRE_Int > m_nnz_pr_diag( locsize, nnz_pr_diag );
110 HYPRE_IJMatrixSetRowSizes( A, &m_nnz_pr_diag[0] );
111 }
112 }
113
114 /* Initialize before setting coefficients */
115 HYPRE_IJMatrixInitialize( A );
116 /* Get the parcsr matrix object to use */
117 HYPRE_IJMatrixGetObject( A, (void**)&A_parcsr );
118
119 /* define an interpreter for the ParCSR interface */
120 HYPRE_MatvecFunctions matvec_fn;
121 interpreter = hypre_CTAlloc( mv_InterfaceInterpreter, 1 );
122 HYPRE_ParCSRSetupInterpreter( interpreter );
123 HYPRE_ParCSRSetupMatvec( &matvec_fn );
124 gnrows = glob_size;
125 gncols = glob_size;
126 height = GetNumRows();
127 width = GetNumCols();
128 }
129
130 void HypreParMatrix::resize( HYPRE_Int global_num_rows,
131 HYPRE_Int global_num_cols,
132 HYPRE_Int* row_starts,
133 HYPRE_Int* col_starts,
134 HYPRE_Int* nnz_pr_diag,
135 HYPRE_Int* onz_pr_diag,
136 HYPRE_Int nnz_pr_offdiag )
137 {
138 /* Create the matrix.
139 Note that this is a square matrix, so we indicate the row partition
140 size twice (since number of rows = number of cols) */
141 HYPRE_IJMatrixCreate( pcomm->comm(), row_starts[0], row_starts[1], col_starts[0], col_starts[1], &A );
142 /* Choose a parallel csr format storage (see the User's Manual) */
143 HYPRE_IJMatrixSetObjectType( A, HYPRE_PARCSR );
144
145 int num_procs;
146 MPI_Comm_size( pcomm->comm(), &num_procs );
147
148 /* Set the matrix pre-allocation data */
149 if( num_procs > 1 )
150 {
151 if( nnz_pr_diag == NULL && onz_pr_diag == NULL )
152 {
153 std::cout << "Parameter nnz_pr_diag and onz_pr_diag cannot be NULL.\n";
154 moab_hypre_assert_t( nnz_pr_diag, true );
155 }
156
157 HYPRE_IJMatrixSetDiagOffdSizes( A, nnz_pr_diag, onz_pr_diag );
158 }
159 else
160 {
161 HYPRE_IJMatrixSetRowSizes( A, nnz_pr_diag );
162 }
163
164 if( nnz_pr_offdiag )
165 {
166 HYPRE_IJMatrixSetMaxOffProcElmts( A, nnz_pr_offdiag );
167 }
168
169 /* Initialize before setting coefficients */
170 HYPRE_IJMatrixInitialize( A );
171 /* Get the parcsr matrix object to use */
172 HYPRE_IJMatrixGetObject( A, (void**)&A_parcsr );
173
174 /* define an interpreter for the ParCSR interface */
175 HYPRE_MatvecFunctions matvec_fn;
176 interpreter = hypre_CTAlloc( mv_InterfaceInterpreter, 1 );
177 HYPRE_ParCSRSetupInterpreter( interpreter );
178 HYPRE_ParCSRSetupMatvec( &matvec_fn );
179 gnrows = global_num_rows;
180 gncols = global_num_cols;
181 height = GetNumRows();
182 width = GetNumCols();
183 }
184
185 // Rectangular block-diagonal constructor
186 HypreParMatrix::HypreParMatrix( moab::ParallelComm* p_comm,
187 HYPRE_Int global_num_rows,
188 HYPRE_Int global_num_cols,
189 HYPRE_Int* row_starts,
190 HYPRE_Int* col_starts,
191 HYPRE_Int nnz_pr_diag,
192 HYPRE_Int onz_pr_diag,
193 HYPRE_Int nnz_pr_offdiag )
194 : pcomm( p_comm )
195 {
196 Init();
197 std::vector< HYPRE_Int > m_nnz_pr_diag( row_starts[1] - row_starts[0], nnz_pr_diag );
198 std::vector< HYPRE_Int > m_onz_pr_diag( row_starts[1] - row_starts[0], onz_pr_diag );
199 resize( global_num_rows, global_num_cols, row_starts, col_starts, &m_nnz_pr_diag[0], &m_onz_pr_diag[0],
200 nnz_pr_offdiag );
201 }
202
203 void HypreParMatrix::MakeRef( const HypreParMatrix& master )
204 {
205 Destroy();
206 Init();
207 A = master.A;
208 A_parcsr = master.A_parcsr;
209 ParCSROwner = 0;
210 height = master.GetNumRows();
211 width = master.GetNumCols();
212 }
213
214 HYPRE_IJMatrix HypreParMatrix::StealData()
215 {
216 // Only safe when (diagOwner == -1 && offdOwner == -1 && colMapOwner == -1)
217 // Otherwise, there may be memory leaks or hypre may destroy arrays allocated
218 // with operator new.
219 moab_hypre_assert( diagOwner == -1 && offdOwner == -1 && colMapOwner == -1, "" );
220 moab_hypre_assert( ParCSROwner, "" );
221 HYPRE_IJMatrix R = A;
222 A = NULL;
223 ParCSROwner = 0;
224 Destroy();
225 Init();
226 return R;
227 }
228
229 void HypreParMatrix::StealData( HypreParMatrix& X )
230 {
231 // Only safe when (diagOwner == -1 && offdOwner == -1 && colMapOwner == -1)
232 // Otherwise, there may be memory leaks or hypre may destroy arrays allocated
233 // with operator new.
234 moab_hypre_assert( diagOwner == -1 && offdOwner == -1 && colMapOwner == -1, "" );
235 moab_hypre_assert( ParCSROwner, "" );
236 moab_hypre_assert( X.diagOwner == -1 && X.offdOwner == -1 && X.colMapOwner == -1, "" );
237 moab_hypre_assert( X.ParCSROwner, "" );
238 A = X.A;
239 A_parcsr = X.A_parcsr;
240 ParCSROwner = 1;
241 X.A = NULL;
242 X.A_parcsr = NULL;
243 X.ParCSROwner = 0;
244 X.Destroy();
245 X.Init();
246 }
247
248 void HypreParMatrix::GetDiag( Eigen::VectorXd& diag ) const
249 {
250 int size = this->height;
251 diag.resize( size );
252
253 for( int j = 0; j < size; j++ )
254 {
255 diag( j ) = A_parcsr->diag->data[A_parcsr->diag->i[j]];
256 moab_hypre_assert( A_parcsr->diag->j[A_parcsr->diag->i[j]] == j,
257 "the first entry in each row must be the diagonal one" );
258 }
259 }
260
261 static void MakeWrapper( const hypre_CSRMatrix* mat, Eigen::SparseMatrix< double, Eigen::RowMajor >& wrapper )
262 {
263 HYPRE_Int nr = hypre_CSRMatrixNumRows( mat );
264 HYPRE_Int nc = hypre_CSRMatrixNumCols( mat );
265 Eigen::SparseMatrix< double, Eigen::RowMajor > tmp( nr, nc );
266 typedef Eigen::Triplet< double > T;
267 HYPRE_Int nnz = hypre_CSRMatrixNumNonzeros( mat );
268 int *rindx, *cindx;
269 double* dindx;
270 #ifndef HYPRE_BIGINT
271 rindx = hypre_CSRMatrixI( mat );
272 cindx = hypre_CSRMatrixJ( mat );
273 dindx = hypre_CSRMatrixData( mat );
274 #else
275 rindx = DuplicateAs< int >( hypre_CSRMatrixI( mat ), nr + 1 );
276 cindx = DuplicateAs< int >( hypre_CSRMatrixJ( mat ), nnz );
277 dindx = hypre_CSRMatrixData( mat );
278 #endif
279 std::vector< T > tripletList( nnz );
280
281 for( int i = 0; i < nnz; ++i )
282 {
283 tripletList.push_back( T( rindx[i], cindx[i], dindx[i] ) );
284 }
285
286 tmp.setFromTriplets( tripletList.begin(), tripletList.end() );
287 wrapper.swap( tmp );
288 }
289
290 void HypreParMatrix::GetDiag( Eigen::SparseMatrix< double, Eigen::RowMajor >& diag ) const
291 {
292 MakeWrapper( A_parcsr->diag, diag );
293 }
294
295 void HypreParMatrix::GetOffd( Eigen::SparseMatrix< double, Eigen::RowMajor >& offd, HYPRE_Int*& cmap ) const
296 {
297 MakeWrapper( A_parcsr->offd, offd );
298 cmap = A_parcsr->col_map_offd;
299 }
300
301 // void HypreParMatrix::GetBlocks(Eigen::Array<HypreParMatrix*,Eigen::Dynamic,Eigen::Dynamic>
302 // &blocks,
303 // bool interleaved_rows,
304 // bool interleaved_cols) const
305 // {
306 // int nr = blocks.rows();
307 // int nc = blocks.cols();
308
309 // hypre_ParCSRMatrix **hypre_blocks = new hypre_ParCSRMatrix*[nr * nc];
310 // internal::hypre_ParCSRMatrixSplit(A, nr, nc, hypre_blocks,
311 // interleaved_rows, interleaved_cols);
312
313 // for (int i = 0; i < nr; i++)
314 // {
315 // for (int j = 0; j < nc; j++)
316 // {
317 // blocks[i][j] = new HypreParMatrix(hypre_blocks[i*nc + j]);
318 // }
319 // }
320
321 // delete [] hypre_blocks;
322 // }
323
324 // HypreParMatrix * HypreParMatrix::Transpose()
325 // {
326 // hypre_ParCSRMatrix * At;
327 // hypre_ParCSRMatrixTranspose(A_parcsr, &At, 1);
328 // hypre_ParCSRMatrixSetNumNonzeros(At);
329
330 // hypre_MatvecCommPkgCreate(At);
331
332 // return new HypreParMatrix(At);
333 // }
334
335 HYPRE_Int HypreParMatrix::Mult( HypreParVector& x, HypreParVector& y, double a, double b )
336 {
337 return hypre_ParCSRMatrixMatvec( a, A_parcsr, x.x_par, b, y.x_par );
338 }
339
340 void HypreParMatrix::Mult( double a, const HypreParVector& x, double b, HypreParVector& y ) const
341 {
342 hypre_ParCSRMatrixMatvec( a, A_parcsr, x.x_par, b, y.x_par );
343 }
344
345 void HypreParMatrix::MultTranspose( double a, const HypreParVector& x, double b, HypreParVector& y ) const
346 {
347 hypre_ParCSRMatrixMatvecT( a, A_parcsr, y.x_par, b, x.x_par );
348 }
349
350 HYPRE_Int HypreParMatrix::Mult( HYPRE_ParVector x, HYPRE_ParVector y, double a, double b )
351 {
352 return hypre_ParCSRMatrixMatvec( a, A_parcsr, (hypre_ParVector*)x, b, (hypre_ParVector*)y );
353 }
354
355 HYPRE_Int HypreParMatrix::MultTranspose( HypreParVector& x, HypreParVector& y, double a, double b )
356 {
357 return hypre_ParCSRMatrixMatvecT( a, A_parcsr, x.x_par, b, y.x_par );
358 }
359
360 // HypreParMatrix* HypreParMatrix::LeftDiagMult(const Eigen::SparseMatrix<double, Eigen::RowMajor>
361 // &D,
362 // HYPRE_Int* row_starts) const
363 // {
364 // const bool assumed_partition = HYPRE_AssumedPartitionCheck();
365 // const bool same_rows = (D.rows() == hypre_CSRMatrixNumRows(A_parcsr->diag));
366
367 // int part_size;
368 // if (assumed_partition)
369 // {
370 // part_size = 2;
371 // }
372 // else
373 // {
374 // MPI_Comm_size(GetComm(), &part_size);
375 // part_size++;
376 // }
377
378 // HYPRE_Int global_num_rows;
379 // if (same_rows)
380 // {
381 // row_starts = hypre_ParCSRMatrixRowStarts(A_parcsr);
382 // global_num_rows = hypre_ParCSRMatrixGlobalNumRows(A_parcsr);
383 // }
384 // else
385 // {
386 // // MFEM_VERIFY(row_starts != NULL, "the number of rows in D and A is not "
387 // // "the same; row_starts must be given (not NULL)");
388 // // Here, when assumed_partition is true we use row_starts[2], so
389 // // row_starts must come from the GetDofOffsets/GetTrueDofOffsets methods
390 // // of ParFiniteElementSpace (HYPRE's partitions have only 2 entries).
391 // global_num_rows =
392 // assumed_partition ? row_starts[2] : row_starts[part_size-1];
393 // }
394
395 // HYPRE_Int *col_starts = hypre_ParCSRMatrixColStarts(A_parcsr);
396 // HYPRE_Int *col_map_offd;
397
398 // // get the diag and offd blocks as SparseMatrix wrappers
399 // Eigen::SparseMatrix<double, Eigen::RowMajor> A_diag, A_offd;
400 // GetDiag(A_diag);
401 // GetOffd(A_offd, col_map_offd);
402
403 // // multiply the blocks with D and create a new HypreParMatrix
404 // Eigen::SparseMatrix<double, Eigen::RowMajor> DA_diag = (D * A_diag);
405 // Eigen::SparseMatrix<double, Eigen::RowMajor> DA_offd = (D * A_offd);
406
407 // HypreParMatrix* DA =
408 // new HypreParMatrix(GetComm(),
409 // global_num_rows, hypre_ParCSRMatrixGlobalNumCols(A),
410 // DuplicateAs<HYPRE_Int>(row_starts, part_size, false),
411 // DuplicateAs<HYPRE_Int>(col_starts, part_size, false),
412 // &DA_diag, &DA_offd,
413 // DuplicateAs<HYPRE_Int>(col_map_offd, A_offd.cols()));
414
415 // // When HYPRE_BIGINT is defined, we want DA_{diag,offd} to delete their I and
416 // // J arrays but not their data arrays; when HYPRE_BIGINT is not defined, we
417 // // don't want DA_{diag,offd} to delete anything.
418 // // #ifndef HYPRE_BIGINT
419 // // DA_diag.LoseData();
420 // // DA_offd.LoseData();
421 // // #else
422 // // DA_diag.SetDataOwner(false);
423 // // DA_offd.SetDataOwner(false);
424 // // #endif
425
426 // // delete DA_diag;
427 // // delete DA_offd;
428
429 // hypre_ParCSRMatrixSetRowStartsOwner(DA->A, 1);
430 // hypre_ParCSRMatrixSetColStartsOwner(DA->A, 1);
431
432 // DA->diagOwner = DA->offdOwner = 3;
433 // DA->colMapOwner = 1;
434
435 // return DA;
436 // }
437
438 void HypreParMatrix::ScaleRows( const Eigen::VectorXd& diag )
439 {
440 if( hypre_CSRMatrixNumRows( A_parcsr->diag ) != hypre_CSRMatrixNumRows( A_parcsr->offd ) )
441 {
442 MB_SET_ERR_RET( "Row does not match" );
443 }
444
445 if( hypre_CSRMatrixNumRows( A_parcsr->diag ) != diag.size() )
446 {
447 MB_SET_ERR_RET( "Note the Eigen::VectorXd diag is not of compatible dimensions with A\n" );
448 }
449
450 int size = this->height;
451 double* Adiag_data = hypre_CSRMatrixData( A_parcsr->diag );
452 HYPRE_Int* Adiag_i = hypre_CSRMatrixI( A_parcsr->diag );
453 double* Aoffd_data = hypre_CSRMatrixData( A_parcsr->offd );
454 HYPRE_Int* Aoffd_i = hypre_CSRMatrixI( A_parcsr->offd );
455 double val;
456 HYPRE_Int jj;
457
458 for( int i( 0 ); i < size; ++i )
459 {
460 val = diag[i];
461
462 for( jj = Adiag_i[i]; jj < Adiag_i[i + 1]; ++jj )
463 {
464 Adiag_data[jj] *= val;
465 }
466
467 for( jj = Aoffd_i[i]; jj < Aoffd_i[i + 1]; ++jj )
468 {
469 Aoffd_data[jj] *= val;
470 }
471 }
472 }
473
474 void HypreParMatrix::InvScaleRows( const Eigen::VectorXd& diag )
475 {
476 if( hypre_CSRMatrixNumRows( A_parcsr->diag ) != hypre_CSRMatrixNumRows( A_parcsr->offd ) )
477 {
478 MB_SET_ERR_RET( "Row does not match" );
479 }
480
481 if( hypre_CSRMatrixNumRows( A_parcsr->diag ) != diag.size() )
482 {
483 MB_SET_ERR_RET( "Note the Eigen::VectorXd diag is not of compatible dimensions with A_parcsr\n" );
484 }
485
486 int size = this->height;
487 double* Adiag_data = hypre_CSRMatrixData( A_parcsr->diag );
488 HYPRE_Int* Adiag_i = hypre_CSRMatrixI( A_parcsr->diag );
489 double* Aoffd_data = hypre_CSRMatrixData( A_parcsr->offd );
490 HYPRE_Int* Aoffd_i = hypre_CSRMatrixI( A_parcsr->offd );
491 double val;
492 HYPRE_Int jj;
493
494 for( int i( 0 ); i < size; ++i )
495 {
496 #ifdef MOAB_DEBUG
497
498 if( 0.0 == diag( i ) )
499 {
500 MB_SET_ERR_RET( "HypreParMatrix::InvDiagScale : Division by 0" );
501 }
502
503 #endif
504 val = 1. / diag( i );
505
506 for( jj = Adiag_i[i]; jj < Adiag_i[i + 1]; ++jj )
507 {
508 Adiag_data[jj] *= val;
509 }
510
511 for( jj = Aoffd_i[i]; jj < Aoffd_i[i + 1]; ++jj )
512 {
513 Aoffd_data[jj] *= val;
514 }
515 }
516 }
517
518 void HypreParMatrix::operator*=( double s )
519 {
520 if( hypre_CSRMatrixNumRows( A_parcsr->diag ) != hypre_CSRMatrixNumRows( A_parcsr->offd ) )
521 {
522 MB_SET_ERR_RET( "Row does not match" );
523 }
524
525 HYPRE_Int size = hypre_CSRMatrixNumRows( A_parcsr->diag );
526 HYPRE_Int jj;
527 double* Adiag_data = hypre_CSRMatrixData( A_parcsr->diag );
528 HYPRE_Int* Adiag_i = hypre_CSRMatrixI( A_parcsr->diag );
529
530 for( jj = 0; jj < Adiag_i[size]; ++jj )
531 {
532 Adiag_data[jj] *= s;
533 }
534
535 double* Aoffd_data = hypre_CSRMatrixData( A_parcsr->offd );
536 HYPRE_Int* Aoffd_i = hypre_CSRMatrixI( A_parcsr->offd );
537
538 for( jj = 0; jj < Aoffd_i[size]; ++jj )
539 {
540 Aoffd_data[jj] *= s;
541 }
542 }
543
544 HYPRE_Int HypreParMatrix::GetValues( const HYPRE_Int nrows,
545 HYPRE_Int* ncols,
546 HYPRE_Int* rows,
547 HYPRE_Int* cols,
548 HYPRE_Complex* values )
549 {
550 return HYPRE_IJMatrixGetValues( A, nrows, ncols, rows, cols, values );
551 }
552
553 HYPRE_Int HypreParMatrix::SetValues( const HYPRE_Int nrows,
554 HYPRE_Int* ncols,
555 const HYPRE_Int* rows,
556 const HYPRE_Int* cols,
557 const HYPRE_Complex* values )
558 {
559 return HYPRE_IJMatrixSetValues( A, nrows, ncols, rows, cols, values );
560 }
561
562 HYPRE_Int HypreParMatrix::AddToValues( const HYPRE_Int nrows,
563 HYPRE_Int* ncols,
564 const HYPRE_Int* rows,
565 const HYPRE_Int* cols,
566 const HYPRE_Complex* values )
567 {
568 return HYPRE_IJMatrixAddToValues( A, nrows, ncols, rows, cols, values );
569 }
570
571 HYPRE_Int HypreParMatrix::verbosity( const HYPRE_Int level )
572 {
573 return HYPRE_IJMatrixSetPrintLevel( A, level );
574 }
575
576 HYPRE_Int HypreParMatrix::FinalizeAssembly()
577 {
578 return HYPRE_IJMatrixAssemble( A );
579 }
580
581 static void get_sorted_rows_cols( const std::vector< HYPRE_Int >& rows_cols, std::vector< HYPRE_Int >& hypre_sorted )
582 {
583 hypre_sorted.resize( rows_cols.size() );
584 std::copy( rows_cols.begin(), rows_cols.end(), hypre_sorted.begin() );
585 std::sort( hypre_sorted.begin(), hypre_sorted.end() );
586 }
587
588 void HypreParMatrix::Threshold( double threshold )
589 {
590 int ierr = 0;
591 int num_procs;
592 hypre_CSRMatrix* csr_A;
593 hypre_CSRMatrix* csr_A_wo_z;
594 hypre_ParCSRMatrix* parcsr_A_ptr;
595 int* row_starts = NULL;
596 int* col_starts = NULL;
597 int row_start = -1;
598 int row_end = -1;
599 int col_start = -1;
600 int col_end = -1;
601 MPI_Comm_size( pcomm->comm(), &num_procs );
602 ierr += hypre_ParCSRMatrixGetLocalRange( A_parcsr, &row_start, &row_end, &col_start, &col_end );
603 row_starts = hypre_ParCSRMatrixRowStarts( A_parcsr );
604 col_starts = hypre_ParCSRMatrixColStarts( A_parcsr );
605 parcsr_A_ptr = hypre_ParCSRMatrixCreate( pcomm->comm(), row_starts[num_procs], col_starts[num_procs], row_starts,
606 col_starts, 0, 0, 0 );
607 csr_A = hypre_MergeDiagAndOffd( A_parcsr );
608 csr_A_wo_z = hypre_CSRMatrixDeleteZeros( csr_A, threshold );
609
610 /* hypre_CSRMatrixDeleteZeros will return a NULL pointer rather than a usable
611 CSR matrix if it finds no non-zeros */
612 if( csr_A_wo_z == NULL )
613 {
614 csr_A_wo_z = csr_A;
615 }
616 else
617 {
618 ierr += hypre_CSRMatrixDestroy( csr_A );
619 }
620
621 ierr += GenerateDiagAndOffd( csr_A_wo_z, parcsr_A_ptr, col_start, col_end );
622 ierr += hypre_CSRMatrixDestroy( csr_A_wo_z );
623 ierr += hypre_ParCSRMatrixDestroy( A_parcsr );
624 hypre_IJMatrixObject( A ) = parcsr_A_ptr;
625 /* Get the parcsr matrix object to use */
626 HYPRE_IJMatrixGetObject( A, (void**)A_parcsr );
627 }
628
629 void HypreParMatrix::EliminateRowsCols( const std::vector< HYPRE_Int >& rows_cols,
630 const HypreParVector& X,
631 HypreParVector& B )
632 {
633 std::vector< HYPRE_Int > rc_sorted;
634 get_sorted_rows_cols( rows_cols, rc_sorted );
635 internal::hypre_ParCSRMatrixEliminateAXB( A_parcsr, rc_sorted.size(), rc_sorted.data(), X.x_par, B.x_par );
636 }
637
638 HypreParMatrix* HypreParMatrix::EliminateRowsCols( const std::vector< HYPRE_Int >& rows_cols )
639 {
640 std::vector< HYPRE_Int > rc_sorted;
641 get_sorted_rows_cols( rows_cols, rc_sorted );
642 hypre_ParCSRMatrix* Ae;
643 internal::hypre_ParCSRMatrixEliminateAAe( A_parcsr, &Ae, rc_sorted.size(), rc_sorted.data() );
644 HypreParMatrix* tmpMat = new HypreParMatrix( pcomm, M(), N(), RowPart(), ColPart(), 0, 0 );
645 hypre_IJMatrixObject( tmpMat->A ) = Ae;
646 /* Get the parcsr matrix object to use */
647 HYPRE_IJMatrixGetObject( tmpMat->A, (void**)tmpMat->A_parcsr );
648 return tmpMat;
649 }
650
651 void HypreParMatrix::Print( const char* fname, HYPRE_Int offi, HYPRE_Int offj )
652 {
653 hypre_ParCSRMatrixPrintIJ( A_parcsr, offi, offj, fname );
654 }
655
656 void HypreParMatrix::Read( const char* fname )
657 {
658 Destroy();
659 Init();
660 HYPRE_Int base_i, base_j;
661 hypre_ParCSRMatrixReadIJ( pcomm->comm(), fname, &base_i, &base_j, &A_parcsr );
662 hypre_ParCSRMatrixSetNumNonzeros( A_parcsr );
663 hypre_MatvecCommPkgCreate( A_parcsr );
664 height = GetNumRows();
665 width = GetNumCols();
666 }
667
668 void HypreParMatrix::Destroy()
669 {
670 if( interpreter )
671 {
672 hypre_TFree( interpreter );
673 }
674
675 if( A == NULL )
676 {
677 return;
678 }
679
680 else
681 {
682 HYPRE_IJMatrixDestroy( A );
683 }
684 }
685
686 HypreParMatrix* ParMult( HypreParMatrix* A, HypreParMatrix* B )
687 {
688 hypre_ParCSRMatrix* ab;
689 ab = hypre_ParMatmul( A->A_parcsr, B->A_parcsr );
690 hypre_MatvecCommPkgCreate( ab );
691 HypreParMatrix* tmpMat = new HypreParMatrix( A->pcomm, A->M(), A->N(), A->RowPart(), A->ColPart(), 0, 0 );
692 hypre_IJMatrixObject( tmpMat->A ) = ab;
693 /* Get the parcsr matrix object to use */
694 HYPRE_IJMatrixGetObject( tmpMat->A, (void**)tmpMat->A_parcsr );
695 return tmpMat;
696 }
697
698 HypreParMatrix* RAP( HypreParMatrix* A, HypreParMatrix* P )
699 {
700 HYPRE_Int P_owns_its_col_starts = hypre_ParCSRMatrixOwnsColStarts( (hypre_ParCSRMatrix*)( P->A_parcsr ) );
701 hypre_ParCSRMatrix* rap;
702 hypre_BoomerAMGBuildCoarseOperator( P->A_parcsr, A->A_parcsr, P->A_parcsr, &rap );
703 hypre_ParCSRMatrixSetNumNonzeros( rap );
704 // hypre_MatvecCommPkgCreate(rap);
705 /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
706 from P (even if it does not own them)! */
707 hypre_ParCSRMatrixSetRowStartsOwner( rap, 0 );
708 hypre_ParCSRMatrixSetColStartsOwner( rap, 0 );
709
710 if( P_owns_its_col_starts )
711 {
712 hypre_ParCSRMatrixSetColStartsOwner( P->A_parcsr, 1 );
713 }
714
715 HypreParMatrix* tmpMat = new HypreParMatrix( A->pcomm, A->M(), A->N(), A->RowPart(), A->ColPart(), 0, 0 );
716 hypre_IJMatrixObject( tmpMat->A ) = rap;
717 /* Get the parcsr matrix object to use */
718 HYPRE_IJMatrixGetObject( tmpMat->A, (void**)tmpMat->A_parcsr );
719 return tmpMat;
720 }
721
722 HypreParMatrix* RAP( HypreParMatrix* Rt, HypreParMatrix* A, HypreParMatrix* P )
723 {
724 HYPRE_Int P_owns_its_col_starts = hypre_ParCSRMatrixOwnsColStarts( (hypre_ParCSRMatrix*)( P->A_parcsr ) );
725 HYPRE_Int Rt_owns_its_col_starts = hypre_ParCSRMatrixOwnsColStarts( (hypre_ParCSRMatrix*)( Rt->A_parcsr ) );
726 hypre_ParCSRMatrix* rap;
727 hypre_BoomerAMGBuildCoarseOperator( Rt->A_parcsr, A->A_parcsr, P->A_parcsr, &rap );
728 hypre_ParCSRMatrixSetNumNonzeros( rap );
729
730 // hypre_MatvecCommPkgCreate(rap);
731 if( !P_owns_its_col_starts )
732 {
733 /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
734 from P (even if it does not own them)! */
735 hypre_ParCSRMatrixSetColStartsOwner( rap, 0 );
736 }
737
738 if( !Rt_owns_its_col_starts )
739 {
740 /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
741 from P (even if it does not own them)! */
742 hypre_ParCSRMatrixSetRowStartsOwner( rap, 0 );
743 }
744
745 HypreParMatrix* tmpMat = new HypreParMatrix( A->pcomm, A->M(), A->N(), A->RowPart(), A->ColPart(), 0, 0 );
746 hypre_IJMatrixObject( tmpMat->A ) = rap;
747 /* Get the parcsr matrix object to use */
748 HYPRE_IJMatrixGetObject( tmpMat->A, (void**)tmpMat->A_parcsr );
749 return tmpMat;
750 }
751
752 void EliminateBC( HypreParMatrix& A,
753 HypreParMatrix& Ae,
754 const std::vector< int >& ess_dof_list,
755 const HypreParVector& X,
756 HypreParVector& B )
757 {
758 // B -= Ae*X
759 Ae.Mult( -1.0, X, 1.0, B );
760 hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( (hypre_ParCSRMatrix*)A.A_parcsr );
761 double* data = hypre_CSRMatrixData( A_diag );
762 HYPRE_Int* I = hypre_CSRMatrixI( A_diag );
763 #ifdef MOAB_DEBUG
764 HYPRE_Int* J = hypre_CSRMatrixJ( A_diag );
765 hypre_CSRMatrix* A_offd = hypre_ParCSRMatrixOffd( (hypre_ParCSRMatrix*)A.A_parcsr );
766 HYPRE_Int* I_offd = hypre_CSRMatrixI( A_offd );
767 double* data_offd = hypre_CSRMatrixData( A_offd );
768 #endif
769 std::vector< HYPRE_Complex > bdata( ess_dof_list.size() ), xdata( ess_dof_list.size() );
770 B.GetValues( ess_dof_list.size(), ess_dof_list.data(), bdata.data() );
771 X.GetValues( ess_dof_list.size(), ess_dof_list.data(), xdata.data() );
772
773 for( size_t i = 0; i < ess_dof_list.size(); i++ )
774 {
775 int r = ess_dof_list[i];
776 bdata[r] = data[I[r]] * xdata[r];
777 #ifdef MOAB_DEBUG
778
779 // Check that in the rows specified by the ess_dof_list, the matrix A has
780 // only one entry -- the diagonal.
781 // if (I[r+1] != I[r]+1 || J[I[r]] != r || I_offd[r] != I_offd[r+1])
782 if( J[I[r]] != r )
783 {
784 MB_SET_ERR_RET( "the diagonal entry must be the first entry in the row!" );
785 }
786
787 for( int j = I[r] + 1; j < I[r + 1]; j++ )
788 {
789 if( data[j] != 0.0 )
790 {
791 MB_SET_ERR_RET( "all off-diagonal entries must be zero!" );
792 }
793 }
794
795 for( int j = I_offd[r]; j < I_offd[r + 1]; j++ )
796 {
797 if( data_offd[j] != 0.0 )
798 {
799 MB_SET_ERR_RET( "all off-diagonal entries must be zero!" );
800 }
801 }
802
803 #endif
804 }
805 }
806
807 // Taubin or "lambda-mu" scheme, which alternates between positive and
808 // negative step sizes to approximate low-pass filter effect.
809
810 int ParCSRRelax_Taubin( hypre_ParCSRMatrix* A, // matrix to relax with
811 hypre_ParVector* f, // right-hand side
812 double lambda,
813 double mu,
814 int N,
815 double max_eig,
816 hypre_ParVector* u, // initial/updated approximation
817 hypre_ParVector* r // another temp vector
818 )
819 {
820 hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( A );
821 HYPRE_Int num_rows = hypre_CSRMatrixNumRows( A_diag );
822 double* u_data = hypre_VectorData( hypre_ParVectorLocalVector( u ) );
823 double* r_data = hypre_VectorData( hypre_ParVectorLocalVector( r ) );
824
825 for( int i = 0; i < N; i++ )
826 {
827 // get residual: r = f - A*u
828 hypre_ParVectorCopy( f, r );
829 hypre_ParCSRMatrixMatvec( -1.0, A, u, 1.0, r );
830 double coef;
831 ( 0 == ( i % 2 ) ) ? coef = lambda : coef = mu;
832
833 for( HYPRE_Int j = 0; j < num_rows; j++ )
834 {
835 u_data[j] += coef * r_data[j] / max_eig;
836 }
837 }
838
839 return 0;
840 }
841
842 // FIR scheme, which uses Chebyshev polynomials and a window function
843 // to approximate a low-pass step filter.
844
845 int ParCSRRelax_FIR( hypre_ParCSRMatrix* A, // matrix to relax with
846 hypre_ParVector* f, // right-hand side
847 double max_eig,
848 int poly_order,
849 double* fir_coeffs,
850 hypre_ParVector* u, // initial/updated approximation
851 hypre_ParVector* x0, // temporaries
852 hypre_ParVector* x1,
853 hypre_ParVector* x2,
854 hypre_ParVector* x3 )
855
856 {
857 hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( A );
858 HYPRE_Int num_rows = hypre_CSRMatrixNumRows( A_diag );
859 double* u_data = hypre_VectorData( hypre_ParVectorLocalVector( u ) );
860 double* x0_data = hypre_VectorData( hypre_ParVectorLocalVector( x0 ) );
861 double* x1_data = hypre_VectorData( hypre_ParVectorLocalVector( x1 ) );
862 double* x2_data = hypre_VectorData( hypre_ParVectorLocalVector( x2 ) );
863 double* x3_data = hypre_VectorData( hypre_ParVectorLocalVector( x3 ) );
864 hypre_ParVectorCopy( u, x0 );
865 // x1 = f -A*x0/max_eig
866 hypre_ParVectorCopy( f, x1 );
867 hypre_ParCSRMatrixMatvec( -1.0, A, x0, 1.0, x1 );
868
869 for( HYPRE_Int i = 0; i < num_rows; i++ )
870 {
871 x1_data[i] /= -max_eig;
872 }
873
874 // x1 = x0 -x1
875 for( HYPRE_Int i = 0; i < num_rows; i++ )
876 {
877 x1_data[i] = x0_data[i] - x1_data[i];
878 }
879
880 // x3 = f0*x0 +f1*x1
881 for( HYPRE_Int i = 0; i < num_rows; i++ )
882 {
883 x3_data[i] = fir_coeffs[0] * x0_data[i] + fir_coeffs[1] * x1_data[i];
884 }
885
886 for( int n = 2; n <= poly_order; n++ )
887 {
888 // x2 = f - A*x1/max_eig
889 hypre_ParVectorCopy( f, x2 );
890 hypre_ParCSRMatrixMatvec( -1.0, A, x1, 1.0, x2 );
891
892 for( HYPRE_Int i = 0; i < num_rows; i++ )
893 {
894 x2_data[i] /= -max_eig;
895 }
896
897 // x2 = (x1-x0) +(x1-2*x2)
898 // x3 = x3 +f[n]*x2
899 // x0 = x1
900 // x1 = x2
901
902 for( HYPRE_Int i = 0; i < num_rows; i++ )
903 {
904 x2_data[i] = ( x1_data[i] - x0_data[i] ) + ( x1_data[i] - 2 * x2_data[i] );
905 x3_data[i] += fir_coeffs[n] * x2_data[i];
906 x0_data[i] = x1_data[i];
907 x1_data[i] = x2_data[i];
908 }
909 }
910
911 for( HYPRE_Int i = 0; i < num_rows; i++ )
912 {
913 u_data[i] = x3_data[i];
914 }
915
916 return 0;
917 }
918
919 } // namespace moab
920
921 #endif
1.9.1.