MOAB: Mesh Oriented datABase  (version 5.2.1)
HypreParMatrix.cpp
Go to the documentation of this file.
00001 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
00002 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
00003 // reserved. See file COPYRIGHT for details.
00004 //
00005 // This file is part of the MFEM library. For more information and source code
00006 // availability see http://mfem.org.
00007 //
00008 // MFEM is free software; you can redistribute it and/or modify it under the
00009 // terms of the GNU Lesser General Public License (as published by the Free
00010 // Software Foundation) version 2.1 dated February 1999.
00011 
00012 #include "HypreParVector.hpp"
00013 #include "HypreParMatrix.hpp"
00014 
00015 #ifdef MOAB_HAVE_MPI
00016 
00017 #include "hypre_parcsr.hpp"
00018 
00019 #include <fstream>
00020 #include <iostream>
00021 #include <cmath>
00022 #include <cstdlib>
00023 #include <cassert>
00024 
00025 using namespace std;
00026 
00027 #define MOAB_DEBUG
00028 
00029 namespace moab
00030 {
00031 #define moab_hypre_assert( a, b ) \
00032     {                             \
00033     }
00034 #define moab_hypre_assert_t( a, b )                         \
00035     {                                                       \
00036         if( !a )                                            \
00037         {                                                   \
00038             std::cout << "HYPRE Error: " << b << std::endl; \
00039             exit( -1 );                                     \
00040         }                                                   \
00041     }
00042 
00043 template < typename TargetT, typename SourceT >
00044 static TargetT* DuplicateAs( const SourceT* array, int size, bool cplusplus = true )
00045 {
00046     TargetT* target_array = cplusplus ? new TargetT[size]
00047                                       /*     */
00048                                       : hypre_TAlloc( TargetT, size );
00049 
00050     for( int i = 0; i < size; i++ )
00051     {
00052         target_array[i] = array[i];
00053     }
00054 
00055     return target_array;
00056 }
00057 
00058 void HypreParMatrix::Init()
00059 {
00060     A           = NULL;
00061     A_parcsr    = NULL;
00062     ParCSROwner = 1;
00063 }
00064 
00065 HypreParMatrix::HypreParMatrix( moab::ParallelComm* p_comm ) : pcomm( p_comm )
00066 {
00067     Init();
00068     height = width = 0;
00069 }
00070 
00071 // Square block-diagonal constructor
00072 HypreParMatrix::HypreParMatrix( moab::ParallelComm* p_comm, HYPRE_Int glob_size, HYPRE_Int* row_starts,
00073                                 HYPRE_Int nnz_pr )
00074     : pcomm( p_comm )
00075 {
00076     Init();
00077     resize( glob_size, row_starts, nnz_pr );
00078 }
00079 
00080 void HypreParMatrix::resize( HYPRE_Int glob_size, HYPRE_Int* row_starts, HYPRE_Int nnz_pr_diag,
00081                              HYPRE_Int nnz_pr_offdiag )
00082 {
00083     /* Create the matrix.
00084         Note that this is a square matrix, so we indicate the row partition
00085         size twice (since number of rows = number of cols) */
00086     HYPRE_IJMatrixCreate( pcomm->comm(), row_starts[0], row_starts[1], row_starts[0], row_starts[1], &A );
00087     /* Choose a parallel csr format storage (see the User's Manual) */
00088     HYPRE_IJMatrixSetObjectType( A, HYPRE_PARCSR );
00089 
00090     int locsize = row_starts[1] - row_starts[0] + 1;
00091     int num_procs, rank;
00092     MPI_Comm_size( pcomm->comm(), &num_procs );
00093     MPI_Comm_rank( pcomm->comm(), &rank );
00094 
00095     if( nnz_pr_diag != 0 || nnz_pr_offdiag != 0 )
00096     {
00097         if( num_procs > 1 )
00098         {
00099             std::vector< HYPRE_Int > m_nnz_pr_diag( locsize, nnz_pr_diag );
00100             std::vector< HYPRE_Int > m_onz_pr_diag( locsize, nnz_pr_offdiag );
00101             HYPRE_IJMatrixSetDiagOffdSizes( A, &m_nnz_pr_diag[0], &m_onz_pr_diag[0] );
00102         }
00103         else
00104         {
00105             std::vector< HYPRE_Int > m_nnz_pr_diag( locsize, nnz_pr_diag );
00106             HYPRE_IJMatrixSetRowSizes( A, &m_nnz_pr_diag[0] );
00107         }
00108     }
00109 
00110     /* Initialize before setting coefficients */
00111     HYPRE_IJMatrixInitialize( A );
00112     /* Get the parcsr matrix object to use */
00113     HYPRE_IJMatrixGetObject( A, (void**)&A_parcsr );
00114 
00115     /* define an interpreter for the ParCSR interface */
00116     HYPRE_MatvecFunctions matvec_fn;
00117     interpreter = hypre_CTAlloc( mv_InterfaceInterpreter, 1 );
00118     HYPRE_ParCSRSetupInterpreter( interpreter );
00119     HYPRE_ParCSRSetupMatvec( &matvec_fn );
00120     gnrows = glob_size;
00121     gncols = glob_size;
00122     height = GetNumRows();
00123     width  = GetNumCols();
00124 }
00125 
00126 void HypreParMatrix::resize( HYPRE_Int global_num_rows, HYPRE_Int global_num_cols, HYPRE_Int* row_starts,
00127                              HYPRE_Int* col_starts, HYPRE_Int* nnz_pr_diag, HYPRE_Int* onz_pr_diag,
00128                              HYPRE_Int nnz_pr_offdiag )
00129 {
00130     /* Create the matrix.
00131         Note that this is a square matrix, so we indicate the row partition
00132         size twice (since number of rows = number of cols) */
00133     HYPRE_IJMatrixCreate( pcomm->comm(), row_starts[0], row_starts[1], col_starts[0], col_starts[1], &A );
00134     /* Choose a parallel csr format storage (see the User's Manual) */
00135     HYPRE_IJMatrixSetObjectType( A, HYPRE_PARCSR );
00136 
00137     int num_procs;
00138     MPI_Comm_size( pcomm->comm(), &num_procs );
00139 
00140     /* Set the matrix pre-allocation data */
00141     if( num_procs > 1 )
00142     {
00143         if( nnz_pr_diag == NULL && onz_pr_diag == NULL )
00144         {
00145             std::cout << "Parameter nnz_pr_diag and onz_pr_diag cannot be NULL.\n";
00146             moab_hypre_assert_t( nnz_pr_diag, true );
00147         }
00148 
00149         HYPRE_IJMatrixSetDiagOffdSizes( A, nnz_pr_diag, onz_pr_diag );
00150     }
00151     else
00152     {
00153         HYPRE_IJMatrixSetRowSizes( A, nnz_pr_diag );
00154     }
00155 
00156     if( nnz_pr_offdiag ) { HYPRE_IJMatrixSetMaxOffProcElmts( A, nnz_pr_offdiag ); }
00157 
00158     /* Initialize before setting coefficients */
00159     HYPRE_IJMatrixInitialize( A );
00160     /* Get the parcsr matrix object to use */
00161     HYPRE_IJMatrixGetObject( A, (void**)&A_parcsr );
00162 
00163     /* define an interpreter for the ParCSR interface */
00164     HYPRE_MatvecFunctions matvec_fn;
00165     interpreter = hypre_CTAlloc( mv_InterfaceInterpreter, 1 );
00166     HYPRE_ParCSRSetupInterpreter( interpreter );
00167     HYPRE_ParCSRSetupMatvec( &matvec_fn );
00168     gnrows = global_num_rows;
00169     gncols = global_num_cols;
00170     height = GetNumRows();
00171     width  = GetNumCols();
00172 }
00173 
00174 // Rectangular block-diagonal constructor
00175 HypreParMatrix::HypreParMatrix( moab::ParallelComm* p_comm, HYPRE_Int global_num_rows, HYPRE_Int global_num_cols,
00176                                 HYPRE_Int* row_starts, HYPRE_Int* col_starts, HYPRE_Int nnz_pr_diag,
00177                                 HYPRE_Int onz_pr_diag, HYPRE_Int nnz_pr_offdiag )
00178     : pcomm( p_comm )
00179 {
00180     Init();
00181     std::vector< HYPRE_Int > m_nnz_pr_diag( row_starts[1] - row_starts[0], nnz_pr_diag );
00182     std::vector< HYPRE_Int > m_onz_pr_diag( row_starts[1] - row_starts[0], onz_pr_diag );
00183     resize( global_num_rows, global_num_cols, row_starts, col_starts, &m_nnz_pr_diag[0], &m_onz_pr_diag[0],
00184             nnz_pr_offdiag );
00185 }
00186 
00187 void HypreParMatrix::MakeRef( const HypreParMatrix& master )
00188 {
00189     Destroy();
00190     Init();
00191     A           = master.A;
00192     A_parcsr    = master.A_parcsr;
00193     ParCSROwner = 0;
00194     height      = master.GetNumRows();
00195     width       = master.GetNumCols();
00196 }
00197 
00198 HYPRE_IJMatrix HypreParMatrix::StealData()
00199 {
00200     // Only safe when (diagOwner == -1 && offdOwner == -1 && colMapOwner == -1)
00201     // Otherwise, there may be memory leaks or hypre may destroy arrays allocated
00202     // with operator new.
00203     moab_hypre_assert( diagOwner == -1 && offdOwner == -1 && colMapOwner == -1, "" );
00204     moab_hypre_assert( ParCSROwner, "" );
00205     HYPRE_IJMatrix R = A;
00206     A                = NULL;
00207     ParCSROwner      = 0;
00208     Destroy();
00209     Init();
00210     return R;
00211 }
00212 
00213 void HypreParMatrix::StealData( HypreParMatrix& X )
00214 {
00215     // Only safe when (diagOwner == -1 && offdOwner == -1 && colMapOwner == -1)
00216     // Otherwise, there may be memory leaks or hypre may destroy arrays allocated
00217     // with operator new.
00218     moab_hypre_assert( diagOwner == -1 && offdOwner == -1 && colMapOwner == -1, "" );
00219     moab_hypre_assert( ParCSROwner, "" );
00220     moab_hypre_assert( X.diagOwner == -1 && X.offdOwner == -1 && X.colMapOwner == -1, "" );
00221     moab_hypre_assert( X.ParCSROwner, "" );
00222     A             = X.A;
00223     A_parcsr      = X.A_parcsr;
00224     ParCSROwner   = 1;
00225     X.A           = NULL;
00226     X.A_parcsr    = NULL;
00227     X.ParCSROwner = 0;
00228     X.Destroy();
00229     X.Init();
00230 }
00231 
00232 void HypreParMatrix::GetDiag( Eigen::VectorXd& diag ) const
00233 {
00234     int size = this->height;
00235     diag.resize( size );
00236 
00237     for( int j = 0; j < size; j++ )
00238     {
00239         diag( j ) = A_parcsr->diag->data[A_parcsr->diag->i[j]];
00240         moab_hypre_assert( A_parcsr->diag->j[A_parcsr->diag->i[j]] == j,
00241                            "the first entry in each row must be the diagonal one" );
00242     }
00243 }
00244 
00245 static void MakeWrapper( const hypre_CSRMatrix* mat, Eigen::SparseMatrix< double, Eigen::RowMajor >& wrapper )
00246 {
00247     HYPRE_Int nr = hypre_CSRMatrixNumRows( mat );
00248     HYPRE_Int nc = hypre_CSRMatrixNumCols( mat );
00249     Eigen::SparseMatrix< double, Eigen::RowMajor > tmp( nr, nc );
00250     typedef Eigen::Triplet< double > T;
00251     HYPRE_Int nnz = hypre_CSRMatrixNumNonzeros( mat );
00252     int *rindx, *cindx;
00253     double* dindx;
00254 #ifndef HYPRE_BIGINT
00255     rindx = hypre_CSRMatrixI( mat );
00256     cindx = hypre_CSRMatrixJ( mat );
00257     dindx = hypre_CSRMatrixData( mat );
00258 #else
00259     rindx = DuplicateAs< int >( hypre_CSRMatrixI( mat ), nr + 1 );
00260     cindx = DuplicateAs< int >( hypre_CSRMatrixJ( mat ), nnz );
00261     dindx = hypre_CSRMatrixData( mat );
00262 #endif
00263     std::vector< T > tripletList( nnz );
00264 
00265     for( int i = 0; i < nnz; ++i )
00266     {
00267         tripletList.push_back( T( rindx[i], cindx[i], dindx[i] ) );
00268     }
00269 
00270     tmp.setFromTriplets( tripletList.begin(), tripletList.end() );
00271     wrapper.swap( tmp );
00272 }
00273 
00274 void HypreParMatrix::GetDiag( Eigen::SparseMatrix< double, Eigen::RowMajor >& diag ) const
00275 {
00276     MakeWrapper( A_parcsr->diag, diag );
00277 }
00278 
00279 void HypreParMatrix::GetOffd( Eigen::SparseMatrix< double, Eigen::RowMajor >& offd, HYPRE_Int*& cmap ) const
00280 {
00281     MakeWrapper( A_parcsr->offd, offd );
00282     cmap = A_parcsr->col_map_offd;
00283 }
00284 
00285 // void HypreParMatrix::GetBlocks(Eigen::Array<HypreParMatrix*,Eigen::Dynamic,Eigen::Dynamic>
00286 // &blocks,
00287 //                                bool interleaved_rows,
00288 //                                bool interleaved_cols) const
00289 // {
00290 //    int nr = blocks.rows();
00291 //    int nc = blocks.cols();
00292 
00293 //    hypre_ParCSRMatrix **hypre_blocks = new hypre_ParCSRMatrix*[nr * nc];
00294 //    internal::hypre_ParCSRMatrixSplit(A, nr, nc, hypre_blocks,
00295 //                                      interleaved_rows, interleaved_cols);
00296 
00297 //    for (int i = 0; i < nr; i++)
00298 //    {
00299 //       for (int j = 0; j < nc; j++)
00300 //       {
00301 //          blocks[i][j] = new HypreParMatrix(hypre_blocks[i*nc + j]);
00302 //       }
00303 //    }
00304 
00305 //    delete [] hypre_blocks;
00306 // }
00307 
00308 // HypreParMatrix * HypreParMatrix::Transpose()
00309 // {
00310 //    hypre_ParCSRMatrix * At;
00311 //    hypre_ParCSRMatrixTranspose(A_parcsr, &At, 1);
00312 //    hypre_ParCSRMatrixSetNumNonzeros(At);
00313 
00314 //    hypre_MatvecCommPkgCreate(At);
00315 
00316 //    return new HypreParMatrix(At);
00317 // }
00318 
00319 HYPRE_Int HypreParMatrix::Mult( HypreParVector& x, HypreParVector& y, double a, double b )
00320 {
00321     return hypre_ParCSRMatrixMatvec( a, A_parcsr, x.x_par, b, y.x_par );
00322 }
00323 
00324 void HypreParMatrix::Mult( double a, const HypreParVector& x, double b, HypreParVector& y ) const
00325 {
00326     hypre_ParCSRMatrixMatvec( a, A_parcsr, x.x_par, b, y.x_par );
00327 }
00328 
00329 void HypreParMatrix::MultTranspose( double a, const HypreParVector& x, double b, HypreParVector& y ) const
00330 {
00331     hypre_ParCSRMatrixMatvecT( a, A_parcsr, y.x_par, b, x.x_par );
00332 }
00333 
00334 HYPRE_Int HypreParMatrix::Mult( HYPRE_ParVector x, HYPRE_ParVector y, double a, double b )
00335 {
00336     return hypre_ParCSRMatrixMatvec( a, A_parcsr, (hypre_ParVector*)x, b, (hypre_ParVector*)y );
00337 }
00338 
00339 HYPRE_Int HypreParMatrix::MultTranspose( HypreParVector& x, HypreParVector& y, double a, double b )
00340 {
00341     return hypre_ParCSRMatrixMatvecT( a, A_parcsr, x.x_par, b, y.x_par );
00342 }
00343 
00344 // HypreParMatrix* HypreParMatrix::LeftDiagMult(const Eigen::SparseMatrix<double, Eigen::RowMajor>
00345 // &D,
00346 //                                              HYPRE_Int* row_starts) const
00347 // {
00348 //    const bool assumed_partition = HYPRE_AssumedPartitionCheck();
00349 //    const bool same_rows = (D.rows() == hypre_CSRMatrixNumRows(A_parcsr->diag));
00350 
00351 //    int part_size;
00352 //    if (assumed_partition)
00353 //    {
00354 //       part_size = 2;
00355 //    }
00356 //    else
00357 //    {
00358 //       MPI_Comm_size(GetComm(), &part_size);
00359 //       part_size++;
00360 //    }
00361 
00362 //    HYPRE_Int global_num_rows;
00363 //    if (same_rows)
00364 //    {
00365 //       row_starts = hypre_ParCSRMatrixRowStarts(A_parcsr);
00366 //       global_num_rows = hypre_ParCSRMatrixGlobalNumRows(A_parcsr);
00367 //    }
00368 //    else
00369 //    {
00370 //       // MFEM_VERIFY(row_starts != NULL, "the number of rows in D and A is not "
00371 //       //             "the same; row_starts must be given (not NULL)");
00372 //       // Here, when assumed_partition is true we use row_starts[2], so
00373 //       // row_starts must come from the GetDofOffsets/GetTrueDofOffsets methods
00374 //       // of ParFiniteElementSpace (HYPRE's partitions have only 2 entries).
00375 //       global_num_rows =
00376 //          assumed_partition ? row_starts[2] : row_starts[part_size-1];
00377 //    }
00378 
00379 //    HYPRE_Int *col_starts = hypre_ParCSRMatrixColStarts(A_parcsr);
00380 //    HYPRE_Int *col_map_offd;
00381 
00382 //    // get the diag and offd blocks as SparseMatrix wrappers
00383 //    Eigen::SparseMatrix<double, Eigen::RowMajor> A_diag, A_offd;
00384 //    GetDiag(A_diag);
00385 //    GetOffd(A_offd, col_map_offd);
00386 
00387 //    // multiply the blocks with D and create a new HypreParMatrix
00388 //    Eigen::SparseMatrix<double, Eigen::RowMajor> DA_diag = (D * A_diag);
00389 //    Eigen::SparseMatrix<double, Eigen::RowMajor> DA_offd = (D * A_offd);
00390 
00391 //    HypreParMatrix* DA =
00392 //       new HypreParMatrix(GetComm(),
00393 //                          global_num_rows, hypre_ParCSRMatrixGlobalNumCols(A),
00394 //                          DuplicateAs<HYPRE_Int>(row_starts, part_size, false),
00395 //                          DuplicateAs<HYPRE_Int>(col_starts, part_size, false),
00396 //                          &DA_diag, &DA_offd,
00397 //                          DuplicateAs<HYPRE_Int>(col_map_offd, A_offd.cols()));
00398 
00399 //    // When HYPRE_BIGINT is defined, we want DA_{diag,offd} to delete their I and
00400 //    // J arrays but not their data arrays; when HYPRE_BIGINT is not defined, we
00401 //    // don't want DA_{diag,offd} to delete anything.
00402 // // #ifndef HYPRE_BIGINT
00403 // //    DA_diag.LoseData();
00404 // //    DA_offd.LoseData();
00405 // // #else
00406 // //    DA_diag.SetDataOwner(false);
00407 // //    DA_offd.SetDataOwner(false);
00408 // // #endif
00409 
00410 //    // delete DA_diag;
00411 //    // delete DA_offd;
00412 
00413 //    hypre_ParCSRMatrixSetRowStartsOwner(DA->A, 1);
00414 //    hypre_ParCSRMatrixSetColStartsOwner(DA->A, 1);
00415 
00416 //    DA->diagOwner = DA->offdOwner = 3;
00417 //    DA->colMapOwner = 1;
00418 
00419 //    return DA;
00420 // }
00421 
00422 void HypreParMatrix::ScaleRows( const Eigen::VectorXd& diag )
00423 {
00424     if( hypre_CSRMatrixNumRows( A_parcsr->diag ) != hypre_CSRMatrixNumRows( A_parcsr->offd ) )
00425     { MB_SET_ERR_RET( "Row does not match" ); }
00426 
00427     if( hypre_CSRMatrixNumRows( A_parcsr->diag ) != diag.size() )
00428     { MB_SET_ERR_RET( "Note the Eigen::VectorXd diag is not of compatible dimensions with A\n" ); }
00429 
00430     int size           = this->height;
00431     double* Adiag_data = hypre_CSRMatrixData( A_parcsr->diag );
00432     HYPRE_Int* Adiag_i = hypre_CSRMatrixI( A_parcsr->diag );
00433     double* Aoffd_data = hypre_CSRMatrixData( A_parcsr->offd );
00434     HYPRE_Int* Aoffd_i = hypre_CSRMatrixI( A_parcsr->offd );
00435     double val;
00436     HYPRE_Int jj;
00437 
00438     for( int i( 0 ); i < size; ++i )
00439     {
00440         val = diag[i];
00441 
00442         for( jj = Adiag_i[i]; jj < Adiag_i[i + 1]; ++jj )
00443         {
00444             Adiag_data[jj] *= val;
00445         }
00446 
00447         for( jj = Aoffd_i[i]; jj < Aoffd_i[i + 1]; ++jj )
00448         {
00449             Aoffd_data[jj] *= val;
00450         }
00451     }
00452 }
00453 
00454 void HypreParMatrix::InvScaleRows( const Eigen::VectorXd& diag )
00455 {
00456     if( hypre_CSRMatrixNumRows( A_parcsr->diag ) != hypre_CSRMatrixNumRows( A_parcsr->offd ) )
00457     { MB_SET_ERR_RET( "Row does not match" ); }
00458 
00459     if( hypre_CSRMatrixNumRows( A_parcsr->diag ) != diag.size() )
00460     { MB_SET_ERR_RET( "Note the Eigen::VectorXd diag is not of compatible dimensions with A_parcsr\n" ); }
00461 
00462     int size           = this->height;
00463     double* Adiag_data = hypre_CSRMatrixData( A_parcsr->diag );
00464     HYPRE_Int* Adiag_i = hypre_CSRMatrixI( A_parcsr->diag );
00465     double* Aoffd_data = hypre_CSRMatrixData( A_parcsr->offd );
00466     HYPRE_Int* Aoffd_i = hypre_CSRMatrixI( A_parcsr->offd );
00467     double val;
00468     HYPRE_Int jj;
00469 
00470     for( int i( 0 ); i < size; ++i )
00471     {
00472 #ifdef MOAB_DEBUG
00473 
00474         if( 0.0 == diag( i ) ) { MB_SET_ERR_RET( "HypreParMatrix::InvDiagScale : Division by 0" ); }
00475 
00476 #endif
00477         val = 1. / diag( i );
00478 
00479         for( jj = Adiag_i[i]; jj < Adiag_i[i + 1]; ++jj )
00480         {
00481             Adiag_data[jj] *= val;
00482         }
00483 
00484         for( jj = Aoffd_i[i]; jj < Aoffd_i[i + 1]; ++jj )
00485         {
00486             Aoffd_data[jj] *= val;
00487         }
00488     }
00489 }
00490 
00491 void HypreParMatrix::operator*=( double s )
00492 {
00493     if( hypre_CSRMatrixNumRows( A_parcsr->diag ) != hypre_CSRMatrixNumRows( A_parcsr->offd ) )
00494     { MB_SET_ERR_RET( "Row does not match" ); }
00495 
00496     HYPRE_Int size = hypre_CSRMatrixNumRows( A_parcsr->diag );
00497     HYPRE_Int jj;
00498     double* Adiag_data = hypre_CSRMatrixData( A_parcsr->diag );
00499     HYPRE_Int* Adiag_i = hypre_CSRMatrixI( A_parcsr->diag );
00500 
00501     for( jj = 0; jj < Adiag_i[size]; ++jj )
00502     {
00503         Adiag_data[jj] *= s;
00504     }
00505 
00506     double* Aoffd_data = hypre_CSRMatrixData( A_parcsr->offd );
00507     HYPRE_Int* Aoffd_i = hypre_CSRMatrixI( A_parcsr->offd );
00508 
00509     for( jj = 0; jj < Aoffd_i[size]; ++jj )
00510     {
00511         Aoffd_data[jj] *= s;
00512     }
00513 }
00514 
00515 HYPRE_Int HypreParMatrix::GetValues( const HYPRE_Int nrows, HYPRE_Int* ncols, HYPRE_Int* rows, HYPRE_Int* cols,
00516                                      HYPRE_Complex* values )
00517 {
00518     return HYPRE_IJMatrixGetValues( A, nrows, ncols, rows, cols, values );
00519 }
00520 
00521 HYPRE_Int HypreParMatrix::SetValues( const HYPRE_Int nrows, HYPRE_Int* ncols, const HYPRE_Int* rows,
00522                                      const HYPRE_Int* cols, const HYPRE_Complex* values )
00523 {
00524     return HYPRE_IJMatrixSetValues( A, nrows, ncols, rows, cols, values );
00525 }
00526 
00527 HYPRE_Int HypreParMatrix::AddToValues( const HYPRE_Int nrows, HYPRE_Int* ncols, const HYPRE_Int* rows,
00528                                        const HYPRE_Int* cols, const HYPRE_Complex* values )
00529 {
00530     return HYPRE_IJMatrixAddToValues( A, nrows, ncols, rows, cols, values );
00531 }
00532 
00533 HYPRE_Int HypreParMatrix::verbosity( const HYPRE_Int level )
00534 {
00535     return HYPRE_IJMatrixSetPrintLevel( A, level );
00536 }
00537 
00538 HYPRE_Int HypreParMatrix::FinalizeAssembly()
00539 {
00540     return HYPRE_IJMatrixAssemble( A );
00541 }
00542 
00543 static void get_sorted_rows_cols( const std::vector< HYPRE_Int >& rows_cols, std::vector< HYPRE_Int >& hypre_sorted )
00544 {
00545     hypre_sorted.resize( rows_cols.size() );
00546     std::copy( rows_cols.begin(), rows_cols.end(), hypre_sorted.begin() );
00547     std::sort( hypre_sorted.begin(), hypre_sorted.end() );
00548 }
00549 
00550 void HypreParMatrix::Threshold( double threshold )
00551 {
00552     int ierr = 0;
00553     int num_procs;
00554     hypre_CSRMatrix* csr_A;
00555     hypre_CSRMatrix* csr_A_wo_z;
00556     hypre_ParCSRMatrix* parcsr_A_ptr;
00557     int* row_starts = NULL;
00558     int* col_starts = NULL;
00559     int row_start   = -1;
00560     int row_end     = -1;
00561     int col_start   = -1;
00562     int col_end     = -1;
00563     MPI_Comm_size( pcomm->comm(), &num_procs );
00564     ierr += hypre_ParCSRMatrixGetLocalRange( A_parcsr, &row_start, &row_end, &col_start, &col_end );
00565     row_starts   = hypre_ParCSRMatrixRowStarts( A_parcsr );
00566     col_starts   = hypre_ParCSRMatrixColStarts( A_parcsr );
00567     parcsr_A_ptr = hypre_ParCSRMatrixCreate( pcomm->comm(), row_starts[num_procs], col_starts[num_procs], row_starts,
00568                                              col_starts, 0, 0, 0 );
00569     csr_A        = hypre_MergeDiagAndOffd( A_parcsr );
00570     csr_A_wo_z   = hypre_CSRMatrixDeleteZeros( csr_A, threshold );
00571 
00572     /* hypre_CSRMatrixDeleteZeros will return a NULL pointer rather than a usable
00573        CSR matrix if it finds no non-zeros */
00574     if( csr_A_wo_z == NULL ) { csr_A_wo_z = csr_A; }
00575     else
00576     {
00577         ierr += hypre_CSRMatrixDestroy( csr_A );
00578     }
00579 
00580     ierr += GenerateDiagAndOffd( csr_A_wo_z, parcsr_A_ptr, col_start, col_end );
00581     ierr += hypre_CSRMatrixDestroy( csr_A_wo_z );
00582     ierr += hypre_ParCSRMatrixDestroy( A_parcsr );
00583     hypre_IJMatrixObject( A ) = parcsr_A_ptr;
00584     /* Get the parcsr matrix object to use */
00585     HYPRE_IJMatrixGetObject( A, (void**)A_parcsr );
00586 }
00587 
00588 void HypreParMatrix::EliminateRowsCols( const std::vector< HYPRE_Int >& rows_cols, const HypreParVector& X,
00589                                         HypreParVector& B )
00590 {
00591     std::vector< HYPRE_Int > rc_sorted;
00592     get_sorted_rows_cols( rows_cols, rc_sorted );
00593     internal::hypre_ParCSRMatrixEliminateAXB( A_parcsr, rc_sorted.size(), rc_sorted.data(), X.x_par, B.x_par );
00594 }
00595 
00596 HypreParMatrix* HypreParMatrix::EliminateRowsCols( const std::vector< HYPRE_Int >& rows_cols )
00597 {
00598     std::vector< HYPRE_Int > rc_sorted;
00599     get_sorted_rows_cols( rows_cols, rc_sorted );
00600     hypre_ParCSRMatrix* Ae;
00601     internal::hypre_ParCSRMatrixEliminateAAe( A_parcsr, &Ae, rc_sorted.size(), rc_sorted.data() );
00602     HypreParMatrix* tmpMat            = new HypreParMatrix( pcomm, M(), N(), RowPart(), ColPart(), 0, 0 );
00603     hypre_IJMatrixObject( tmpMat->A ) = Ae;
00604     /* Get the parcsr matrix object to use */
00605     HYPRE_IJMatrixGetObject( tmpMat->A, (void**)tmpMat->A_parcsr );
00606     return tmpMat;
00607 }
00608 
00609 void HypreParMatrix::Print( const char* fname, HYPRE_Int offi, HYPRE_Int offj )
00610 {
00611     hypre_ParCSRMatrixPrintIJ( A_parcsr, offi, offj, fname );
00612 }
00613 
00614 void HypreParMatrix::Read( const char* fname )
00615 {
00616     Destroy();
00617     Init();
00618     HYPRE_Int base_i, base_j;
00619     hypre_ParCSRMatrixReadIJ( pcomm->comm(), fname, &base_i, &base_j, &A_parcsr );
00620     hypre_ParCSRMatrixSetNumNonzeros( A_parcsr );
00621     hypre_MatvecCommPkgCreate( A_parcsr );
00622     height = GetNumRows();
00623     width  = GetNumCols();
00624 }
00625 
00626 void HypreParMatrix::Destroy()
00627 {
00628     if( interpreter ) { hypre_TFree( interpreter ); }
00629 
00630     if( A == NULL ) { return; }
00631 
00632     else
00633     {
00634         HYPRE_IJMatrixDestroy( A );
00635     }
00636 }
00637 
00638 HypreParMatrix* ParMult( HypreParMatrix* A, HypreParMatrix* B )
00639 {
00640     hypre_ParCSRMatrix* ab;
00641     ab = hypre_ParMatmul( A->A_parcsr, B->A_parcsr );
00642     hypre_MatvecCommPkgCreate( ab );
00643     HypreParMatrix* tmpMat = new HypreParMatrix( A->pcomm, A->M(), A->N(), A->RowPart(), A->ColPart(), 0, 0 );
00644     hypre_IJMatrixObject( tmpMat->A ) = ab;
00645     /* Get the parcsr matrix object to use */
00646     HYPRE_IJMatrixGetObject( tmpMat->A, (void**)tmpMat->A_parcsr );
00647     return tmpMat;
00648 }
00649 
00650 HypreParMatrix* RAP( HypreParMatrix* A, HypreParMatrix* P )
00651 {
00652     HYPRE_Int P_owns_its_col_starts = hypre_ParCSRMatrixOwnsColStarts( (hypre_ParCSRMatrix*)( P->A_parcsr ) );
00653     hypre_ParCSRMatrix* rap;
00654     hypre_BoomerAMGBuildCoarseOperator( P->A_parcsr, A->A_parcsr, P->A_parcsr, &rap );
00655     hypre_ParCSRMatrixSetNumNonzeros( rap );
00656     // hypre_MatvecCommPkgCreate(rap);
00657     /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
00658        from P (even if it does not own them)! */
00659     hypre_ParCSRMatrixSetRowStartsOwner( rap, 0 );
00660     hypre_ParCSRMatrixSetColStartsOwner( rap, 0 );
00661 
00662     if( P_owns_its_col_starts ) { hypre_ParCSRMatrixSetColStartsOwner( P->A_parcsr, 1 ); }
00663 
00664     HypreParMatrix* tmpMat = new HypreParMatrix( A->pcomm, A->M(), A->N(), A->RowPart(), A->ColPart(), 0, 0 );
00665     hypre_IJMatrixObject( tmpMat->A ) = rap;
00666     /* Get the parcsr matrix object to use */
00667     HYPRE_IJMatrixGetObject( tmpMat->A, (void**)tmpMat->A_parcsr );
00668     return tmpMat;
00669 }
00670 
00671 HypreParMatrix* RAP( HypreParMatrix* Rt, HypreParMatrix* A, HypreParMatrix* P )
00672 {
00673     HYPRE_Int P_owns_its_col_starts  = hypre_ParCSRMatrixOwnsColStarts( (hypre_ParCSRMatrix*)( P->A_parcsr ) );
00674     HYPRE_Int Rt_owns_its_col_starts = hypre_ParCSRMatrixOwnsColStarts( (hypre_ParCSRMatrix*)( Rt->A_parcsr ) );
00675     hypre_ParCSRMatrix* rap;
00676     hypre_BoomerAMGBuildCoarseOperator( Rt->A_parcsr, A->A_parcsr, P->A_parcsr, &rap );
00677     hypre_ParCSRMatrixSetNumNonzeros( rap );
00678 
00679     // hypre_MatvecCommPkgCreate(rap);
00680     if( !P_owns_its_col_starts )
00681     {
00682         /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
00683            from P (even if it does not own them)! */
00684         hypre_ParCSRMatrixSetColStartsOwner( rap, 0 );
00685     }
00686 
00687     if( !Rt_owns_its_col_starts )
00688     {
00689         /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
00690            from P (even if it does not own them)! */
00691         hypre_ParCSRMatrixSetRowStartsOwner( rap, 0 );
00692     }
00693 
00694     HypreParMatrix* tmpMat = new HypreParMatrix( A->pcomm, A->M(), A->N(), A->RowPart(), A->ColPart(), 0, 0 );
00695     hypre_IJMatrixObject( tmpMat->A ) = rap;
00696     /* Get the parcsr matrix object to use */
00697     HYPRE_IJMatrixGetObject( tmpMat->A, (void**)tmpMat->A_parcsr );
00698     return tmpMat;
00699 }
00700 
00701 void EliminateBC( HypreParMatrix& A, HypreParMatrix& Ae, const std::vector< int >& ess_dof_list,
00702                   const HypreParVector& X, HypreParVector& B )
00703 {
00704     // B -= Ae*X
00705     Ae.Mult( -1.0, X, 1.0, B );
00706     hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( (hypre_ParCSRMatrix*)A.A_parcsr );
00707     double* data            = hypre_CSRMatrixData( A_diag );
00708     HYPRE_Int* I            = hypre_CSRMatrixI( A_diag );
00709 #ifdef MOAB_DEBUG
00710     HYPRE_Int* J            = hypre_CSRMatrixJ( A_diag );
00711     hypre_CSRMatrix* A_offd = hypre_ParCSRMatrixOffd( (hypre_ParCSRMatrix*)A.A_parcsr );
00712     HYPRE_Int* I_offd       = hypre_CSRMatrixI( A_offd );
00713     double* data_offd       = hypre_CSRMatrixData( A_offd );
00714 #endif
00715     std::vector< HYPRE_Complex > bdata( ess_dof_list.size() ), xdata( ess_dof_list.size() );
00716     B.GetValues( ess_dof_list.size(), ess_dof_list.data(), bdata.data() );
00717     X.GetValues( ess_dof_list.size(), ess_dof_list.data(), xdata.data() );
00718 
00719     for( size_t i = 0; i < ess_dof_list.size(); i++ )
00720     {
00721         int r    = ess_dof_list[i];
00722         bdata[r] = data[I[r]] * xdata[r];
00723 #ifdef MOAB_DEBUG
00724 
00725         // Check that in the rows specified by the ess_dof_list, the matrix A has
00726         // only one entry -- the diagonal.
00727         // if (I[r+1] != I[r]+1 || J[I[r]] != r || I_offd[r] != I_offd[r+1])
00728         if( J[I[r]] != r ) { MB_SET_ERR_RET( "the diagonal entry must be the first entry in the row!" ); }
00729 
00730         for( int j = I[r] + 1; j < I[r + 1]; j++ )
00731         {
00732             if( data[j] != 0.0 ) { MB_SET_ERR_RET( "all off-diagonal entries must be zero!" ); }
00733         }
00734 
00735         for( int j = I_offd[r]; j < I_offd[r + 1]; j++ )
00736         {
00737             if( data_offd[j] != 0.0 ) { MB_SET_ERR_RET( "all off-diagonal entries must be zero!" ); }
00738         }
00739 
00740 #endif
00741     }
00742 }
00743 
00744 // Taubin or "lambda-mu" scheme, which alternates between positive and
00745 // negative step sizes to approximate low-pass filter effect.
00746 
00747 int ParCSRRelax_Taubin( hypre_ParCSRMatrix* A,  // matrix to relax with
00748                         hypre_ParVector* f,     // right-hand side
00749                         double lambda, double mu, int N, double max_eig,
00750                         hypre_ParVector* u,  // initial/updated approximation
00751                         hypre_ParVector* r   // another temp vector
00752 )
00753 {
00754     hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( A );
00755     HYPRE_Int num_rows      = hypre_CSRMatrixNumRows( A_diag );
00756     double* u_data          = hypre_VectorData( hypre_ParVectorLocalVector( u ) );
00757     double* r_data          = hypre_VectorData( hypre_ParVectorLocalVector( r ) );
00758 
00759     for( int i = 0; i < N; i++ )
00760     {
00761         // get residual: r = f - A*u
00762         hypre_ParVectorCopy( f, r );
00763         hypre_ParCSRMatrixMatvec( -1.0, A, u, 1.0, r );
00764         double coef;
00765         ( 0 == ( i % 2 ) ) ? coef = lambda : coef = mu;
00766 
00767         for( HYPRE_Int j = 0; j < num_rows; j++ )
00768         {
00769             u_data[j] += coef * r_data[j] / max_eig;
00770         }
00771     }
00772 
00773     return 0;
00774 }
00775 
00776 // FIR scheme, which uses Chebyshev polynomials and a window function
00777 // to approximate a low-pass step filter.
00778 
00779 int ParCSRRelax_FIR( hypre_ParCSRMatrix* A,  // matrix to relax with
00780                      hypre_ParVector* f,     // right-hand side
00781                      double max_eig, int poly_order, double* fir_coeffs,
00782                      hypre_ParVector* u,   // initial/updated approximation
00783                      hypre_ParVector* x0,  // temporaries
00784                      hypre_ParVector* x1, hypre_ParVector* x2, hypre_ParVector* x3 )
00785 
00786 {
00787     hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( A );
00788     HYPRE_Int num_rows      = hypre_CSRMatrixNumRows( A_diag );
00789     double* u_data          = hypre_VectorData( hypre_ParVectorLocalVector( u ) );
00790     double* x0_data         = hypre_VectorData( hypre_ParVectorLocalVector( x0 ) );
00791     double* x1_data         = hypre_VectorData( hypre_ParVectorLocalVector( x1 ) );
00792     double* x2_data         = hypre_VectorData( hypre_ParVectorLocalVector( x2 ) );
00793     double* x3_data         = hypre_VectorData( hypre_ParVectorLocalVector( x3 ) );
00794     hypre_ParVectorCopy( u, x0 );
00795     // x1 = f -A*x0/max_eig
00796     hypre_ParVectorCopy( f, x1 );
00797     hypre_ParCSRMatrixMatvec( -1.0, A, x0, 1.0, x1 );
00798 
00799     for( HYPRE_Int i = 0; i < num_rows; i++ )
00800     {
00801         x1_data[i] /= -max_eig;
00802     }
00803 
00804     // x1 = x0 -x1
00805     for( HYPRE_Int i = 0; i < num_rows; i++ )
00806     {
00807         x1_data[i] = x0_data[i] - x1_data[i];
00808     }
00809 
00810     // x3 = f0*x0 +f1*x1
00811     for( HYPRE_Int i = 0; i < num_rows; i++ )
00812     {
00813         x3_data[i] = fir_coeffs[0] * x0_data[i] + fir_coeffs[1] * x1_data[i];
00814     }
00815 
00816     for( int n = 2; n <= poly_order; n++ )
00817     {
00818         // x2 = f - A*x1/max_eig
00819         hypre_ParVectorCopy( f, x2 );
00820         hypre_ParCSRMatrixMatvec( -1.0, A, x1, 1.0, x2 );
00821 
00822         for( HYPRE_Int i = 0; i < num_rows; i++ )
00823         {
00824             x2_data[i] /= -max_eig;
00825         }
00826 
00827         // x2 = (x1-x0) +(x1-2*x2)
00828         // x3 = x3 +f[n]*x2
00829         // x0 = x1
00830         // x1 = x2
00831 
00832         for( HYPRE_Int i = 0; i < num_rows; i++ )
00833         {
00834             x2_data[i] = ( x1_data[i] - x0_data[i] ) + ( x1_data[i] - 2 * x2_data[i] );
00835             x3_data[i] += fir_coeffs[n] * x2_data[i];
00836             x0_data[i] = x1_data[i];
00837             x1_data[i] = x2_data[i];
00838         }
00839     }
00840 
00841     for( HYPRE_Int i = 0; i < num_rows; i++ )
00842     {
00843         u_data[i] = x3_data[i];
00844     }
00845 
00846     return 0;
00847 }
00848 
00849 }  // namespace moab
00850 
00851 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines