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