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