MOAB: Mesh Oriented datABase  (version 5.2.1)
HypreParVector.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 <fstream>
00018 #include <iostream>
00019 #include <cmath>
00020 #include <cstdlib>
00021 #include <cassert>
00022 
00023 using namespace std;
00024 
00025 namespace moab
00026 {
00027 #define moab_hypre_assert( a, b ) \
00028     {                             \
00029     }
00030 #define moab_hypre_assert_t( a, b )                         \
00031     {                                                       \
00032         if( !a )                                            \
00033         {                                                   \
00034             std::cout << "HYPRE Error: " << b << std::endl; \
00035             exit( -1 );                                     \
00036         }
00037 
00038 HypreParVector::HypreParVector( moab::ParallelComm* p_comm ) : pcomm( p_comm )
00039 {
00040     x           = NULL;
00041     initialized = size = gsize = rstart = rend = 0;
00042 }
00043 
00044 HypreParVector::HypreParVector( moab::ParallelComm* p_comm, HYPRE_Int glob_size, HYPRE_Int p_irstart,
00045                                 HYPRE_Int p_irend )
00046     : rstart( p_irstart ), rend( p_irend ), pcomm( p_comm )
00047 {
00048     HYPRE_IJVectorCreate( pcomm->comm(), rstart, rend, &x );
00049     HYPRE_IJVectorSetObjectType( x, HYPRE_PARCSR );
00050     HYPRE_IJVectorInitialize( x );
00051     HYPRE_IJVectorAssemble( x );
00052     HYPRE_IJVectorGetObject( x, (void**)&x_par );
00053     size          = rstart - rend;
00054     gsize         = glob_size;
00055     own_ParVector = 1;
00056     initialized   = 1;
00057 }
00058 
00059 HypreParVector::HypreParVector( const HypreParVector& y )
00060 {
00061     pcomm  = y.pcomm;
00062     rstart = y.rstart;
00063     rend   = y.rend;
00064     size   = y.size;
00065     gsize  = y.gsize;
00066     HYPRE_IJVectorCreate( pcomm->comm(), y.rstart, y.rend, &x );
00067     HYPRE_IJVectorSetObjectType( x, HYPRE_PARCSR );
00068     HYPRE_IJVectorInitialize( x );
00069     HYPRE_IJVectorAssemble( x );
00070     HYPRE_IJVectorGetObject( x, (void**)&x_par );
00071     HYPRE_Complex* array = NULL;
00072     HYPRE_IJVectorGetValues( y.x, y.size, NULL, array );
00073     HYPRE_IJVectorSetValues( x, size, NULL, array );
00074     array         = NULL;
00075     own_ParVector = 1;
00076     initialized   = 1;
00077 }
00078 
00079 HypreParVector::HypreParVector( HypreParMatrix& A, int tr )
00080 {
00081     pcomm = A.GetParallelCommunicator();
00082     int* part;
00083 
00084     if( tr ) { part = A.ColPart(); }
00085     else
00086         part = A.RowPart();
00087 
00088     HYPRE_IJVectorCreate( pcomm->comm(), part[0], part[1], &x );
00089     HYPRE_IJVectorSetObjectType( x, HYPRE_PARCSR );
00090     HYPRE_IJVectorInitialize( x );
00091     HYPRE_IJVectorAssemble( x );
00092     HYPRE_IJVectorGetObject( x, (void**)&x_par );
00093     // if (!tr)
00094     // {
00095     //    x_par = hypre_ParVectorInDomainOf(static_cast<hypre_ParCSRMatrix*>(A.A_parcsr));
00096     // }
00097     // else
00098     // {
00099     //    x_par = hypre_ParVectorInRangeOf(static_cast<hypre_ParCSRMatrix*>(A.A_parcsr));
00100     // }
00101     // comm = hypre_ParVectorComm(x);
00102     rstart        = part[0];
00103     rend          = part[1];
00104     size          = rstart - rend;
00105     own_ParVector = 1;
00106     initialized   = 1;
00107 }
00108 
00109 HypreParVector::operator HYPRE_IJVector() const
00110 {
00111     return x;
00112 }
00113 
00114 // #ifndef HYPRE_PAR_VECTOR_STRUCT
00115 HypreParVector::operator HYPRE_ParVector() const
00116 {
00117     return x_par;
00118 }
00119 // #endif
00120 
00121 HypreParVector& HypreParVector::operator=( double d )
00122 {
00123     hypre_ParVectorSetConstantValues( x_par, d );
00124     return *this;
00125 }
00126 
00127 HypreParVector& HypreParVector::operator=( const HypreParVector& y )
00128 {
00129 #ifndef NDEBUG
00130 
00131     if( this->GlobalSize() != y.GlobalSize() )
00132     {  // || local_size != y.local_size) {
00133         MB_SET_ERR_RET_VAL( "HypreParVector::operator failed. Incompatible vector sizes", *this );
00134     }
00135 
00136 #endif
00137     pcomm  = y.pcomm;
00138     rstart = y.rstart;
00139     rend   = y.rend;
00140     size   = y.size;
00141 
00142     if( x ) { HYPRE_IJVectorDestroy( x ); }
00143 
00144     x             = y.x;
00145     x_par         = y.x_par;
00146     own_ParVector = 0;
00147     initialized   = y.initialized;
00148     return *this;
00149 }
00150 
00151 HYPRE_Int HypreParVector::resize( HYPRE_Int /*glob_size*/, HYPRE_Int p_irstart, HYPRE_Int p_irend )
00152 {
00153     if( initialized || x != NULL ) MB_SET_ERR_RET_VAL( "Vector is already initialized and partitioned", -1 );
00154 
00155     HYPRE_IJVectorCreate( this->pcomm->comm(), p_irstart, p_irend, &x );
00156     HYPRE_IJVectorSetObjectType( x, HYPRE_PARCSR );
00157     HYPRE_IJVectorInitialize( x );
00158     HYPRE_IJVectorAssemble( x );
00159     HYPRE_IJVectorGetObject( x, (void**)&x_par );
00160     rstart        = p_irstart;
00161     rend          = p_irend;
00162     size          = rstart - rend;
00163     own_ParVector = 1;
00164     initialized   = 1;
00165     return 0;
00166 }
00167 
00168 HYPRE_Int HypreParVector::SetData( HYPRE_Complex* p_data, HYPRE_Int* p_col )
00169 {
00170     return HYPRE_IJVectorSetValues( x, size, p_col, p_data );
00171 }
00172 
00173 HYPRE_Int HypreParVector::AddData( HYPRE_Complex* p_data, HYPRE_Int* p_col )
00174 {
00175     return HYPRE_IJVectorAddToValues( x, size, p_col, p_data );
00176 }
00177 
00178 HYPRE_Int HypreParVector::GetValues( const int ndata, const HYPRE_Int* indices, HYPRE_Complex* const _data ) const
00179 {
00180     return HYPRE_IJVectorGetValues( x, ndata, indices, _data );
00181 }
00182 
00183 HYPRE_Int HypreParVector::SetValues( const int ndata, const HYPRE_Int* indices, const HYPRE_Complex* const _data )
00184 {
00185     return HYPRE_IJVectorSetValues( x, ndata, indices, _data );
00186 }
00187 
00188 HYPRE_Int HypreParVector::AddValues( const int ndata, const HYPRE_Int* indices, const HYPRE_Complex* const _data )
00189 {
00190     return HYPRE_IJVectorAddToValues( x, ndata, indices, _data );
00191 }
00192 
00193 HYPRE_Int HypreParVector::GetValue( HYPRE_Int index, HYPRE_Complex* const _data ) const
00194 {
00195     return HYPRE_IJVectorGetValues( x, 1, &index, _data );
00196 }
00197 
00198 HYPRE_Int HypreParVector::SetValue( const HYPRE_Int index, const HYPRE_Complex _data )
00199 {
00200     return HYPRE_IJVectorSetValues( x, 1, &index, &_data );
00201 }
00202 
00203 HYPRE_Int HypreParVector::AddValue( const HYPRE_Int index, const HYPRE_Complex _data )
00204 {
00205     return HYPRE_IJVectorAddToValues( x, 1, &index, &_data );
00206 }
00207 
00208 HYPRE_Int HypreParVector::FinalizeAssembly()
00209 {
00210     return HYPRE_IJVectorAssemble( x );
00211 }
00212 
00213 HYPRE_Int HypreParVector::verbosity( const HYPRE_Int level )
00214 {
00215     return HYPRE_IJVectorSetPrintLevel( x, level );
00216 }
00217 
00218 void HypreParVector::Print( const char* fname ) const
00219 {
00220     HYPRE_IJVectorPrint( x, fname );
00221 }
00222 
00223 HypreParVector::~HypreParVector()
00224 {
00225     if( own_ParVector ) { HYPRE_IJVectorDestroy( x ); }
00226 }
00227 
00228 double HypreParVector::InnerProduct( HypreParVector& x, HypreParVector& y )
00229 {
00230     return hypre_ParVectorInnerProd( x.x_par, y.x_par );
00231 }
00232 
00233 }  // namespace moab
00234 
00235 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines