 |
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 <fstream>
18 #include <iostream>
19 #include <cmath>
20 #include <cstdlib>
21 #include <cassert>
22
23 using namespace std;
24
25 namespace moab
26 {
27 #define moab_hypre_assert( a, b ) \
28 { \
29 }
30 #define moab_hypre_assert_t( a, b ) \
31 { \
32 if( !a ) \
33 { \
34 std::cout << "HYPRE Error: " << b << std::endl; \
35 exit( -1 ); \
36 }
37
38 HypreParVector::HypreParVector( moab::ParallelComm* p_comm ) : pcomm( p_comm )
39 {
40 x = NULL;
41 initialized = size = gsize = rstart = rend = 0;
42 }
43
44 HypreParVector::HypreParVector( moab::ParallelComm* p_comm,
45 HYPRE_Int glob_size,
46 HYPRE_Int p_irstart,
47 HYPRE_Int p_irend )
48 : rstart( p_irstart ), rend( p_irend ), pcomm( p_comm )
49 {
50 HYPRE_IJVectorCreate( pcomm->comm(), rstart, rend, &x );
51 HYPRE_IJVectorSetObjectType( x, HYPRE_PARCSR );
52 HYPRE_IJVectorInitialize( x );
53 HYPRE_IJVectorAssemble( x );
54 HYPRE_IJVectorGetObject( x, (void**)&x_par );
55 size = rstart - rend;
56 gsize = glob_size;
57 own_ParVector = 1;
58 initialized = 1;
59 }
60
61 HypreParVector::HypreParVector( const HypreParVector& y )
62 {
63 pcomm = y.pcomm;
64 rstart = y.rstart;
65 rend = y.rend;
66 size = y.size;
67 gsize = y.gsize;
68 HYPRE_IJVectorCreate( pcomm->comm(), y.rstart, y.rend, &x );
69 HYPRE_IJVectorSetObjectType( x, HYPRE_PARCSR );
70 HYPRE_IJVectorInitialize( x );
71 HYPRE_IJVectorAssemble( x );
72 HYPRE_IJVectorGetObject( x, (void**)&x_par );
73 HYPRE_Complex* array = NULL;
74 HYPRE_IJVectorGetValues( y.x, y.size, NULL, array );
75 HYPRE_IJVectorSetValues( x, size, NULL, array );
76 array = NULL;
77 own_ParVector = 1;
78 initialized = 1;
79 }
80
81 HypreParVector::HypreParVector( HypreParMatrix& A, int tr )
82 {
83 pcomm = A.GetParallelCommunicator();
84 int* part;
85
86 if( tr )
87 {
88 part = A.ColPart();
89 }
90 else
91 part = A.RowPart();
92
93 HYPRE_IJVectorCreate( pcomm->comm(), part[0], part[1], &x );
94 HYPRE_IJVectorSetObjectType( x, HYPRE_PARCSR );
95 HYPRE_IJVectorInitialize( x );
96 HYPRE_IJVectorAssemble( x );
97 HYPRE_IJVectorGetObject( x, (void**)&x_par );
98 // if (!tr)
99 // {
100 // x_par = hypre_ParVectorInDomainOf(static_cast<hypre_ParCSRMatrix*>(A.A_parcsr));
101 // }
102 // else
103 // {
104 // x_par = hypre_ParVectorInRangeOf(static_cast<hypre_ParCSRMatrix*>(A.A_parcsr));
105 // }
106 // comm = hypre_ParVectorComm(x);
107 rstart = part[0];
108 rend = part[1];
109 size = rstart - rend;
110 own_ParVector = 1;
111 initialized = 1;
112 }
113
114 HypreParVector::operator HYPRE_IJVector() const
115 {
116 return x;
117 }
118
119 // #ifndef HYPRE_PAR_VECTOR_STRUCT
120 HypreParVector::operator HYPRE_ParVector() const
121 {
122 return x_par;
123 }
124 // #endif
125
126 HypreParVector& HypreParVector::operator=( double d )
127 {
128 hypre_ParVectorSetConstantValues( x_par, d );
129 return *this;
130 }
131
132 HypreParVector& HypreParVector::operator=( const HypreParVector& y )
133 {
134 #ifndef NDEBUG
135
136 if( this->GlobalSize() != y.GlobalSize() )
137 { // || local_size != y.local_size) {
138 MB_SET_ERR_RET_VAL( "HypreParVector::operator failed. Incompatible vector sizes", *this );
139 }
140
141 #endif
142 pcomm = y.pcomm;
143 rstart = y.rstart;
144 rend = y.rend;
145 size = y.size;
146
147 if( x )
148 {
149 HYPRE_IJVectorDestroy( x );
150 }
151
152 x = y.x;
153 x_par = y.x_par;
154 own_ParVector = 0;
155 initialized = y.initialized;
156 return *this;
157 }
158
159 HYPRE_Int HypreParVector::resize( HYPRE_Int /*glob_size*/, HYPRE_Int p_irstart, HYPRE_Int p_irend )
160 {
161 if( initialized || x != NULL ) MB_SET_ERR_RET_VAL( "Vector is already initialized and partitioned", -1 );
162
163 HYPRE_IJVectorCreate( this->pcomm->comm(), p_irstart, p_irend, &x );
164 HYPRE_IJVectorSetObjectType( x, HYPRE_PARCSR );
165 HYPRE_IJVectorInitialize( x );
166 HYPRE_IJVectorAssemble( x );
167 HYPRE_IJVectorGetObject( x, (void**)&x_par );
168 rstart = p_irstart;
169 rend = p_irend;
170 size = rstart - rend;
171 own_ParVector = 1;
172 initialized = 1;
173 return 0;
174 }
175
176 HYPRE_Int HypreParVector::SetData( HYPRE_Complex* p_data, HYPRE_Int* p_col )
177 {
178 return HYPRE_IJVectorSetValues( x, size, p_col, p_data );
179 }
180
181 HYPRE_Int HypreParVector::AddData( HYPRE_Complex* p_data, HYPRE_Int* p_col )
182 {
183 return HYPRE_IJVectorAddToValues( x, size, p_col, p_data );
184 }
185
186 HYPRE_Int HypreParVector::GetValues( const int ndata, const HYPRE_Int* indices, HYPRE_Complex* const _data ) const
187 {
188 return HYPRE_IJVectorGetValues( x, ndata, indices, _data );
189 }
190
191 HYPRE_Int HypreParVector::SetValues( const int ndata, const HYPRE_Int* indices, const HYPRE_Complex* const _data )
192 {
193 return HYPRE_IJVectorSetValues( x, ndata, indices, _data );
194 }
195
196 HYPRE_Int HypreParVector::AddValues( const int ndata, const HYPRE_Int* indices, const HYPRE_Complex* const _data )
197 {
198 return HYPRE_IJVectorAddToValues( x, ndata, indices, _data );
199 }
200
201 HYPRE_Int HypreParVector::GetValue( HYPRE_Int index, HYPRE_Complex* const _data ) const
202 {
203 return HYPRE_IJVectorGetValues( x, 1, &index, _data );
204 }
205
206 HYPRE_Int HypreParVector::SetValue( const HYPRE_Int index, const HYPRE_Complex _data )
207 {
208 return HYPRE_IJVectorSetValues( x, 1, &index, &_data );
209 }
210
211 HYPRE_Int HypreParVector::AddValue( const HYPRE_Int index, const HYPRE_Complex _data )
212 {
213 return HYPRE_IJVectorAddToValues( x, 1, &index, &_data );
214 }
215
216 HYPRE_Int HypreParVector::FinalizeAssembly()
217 {
218 return HYPRE_IJVectorAssemble( x );
219 }
220
221 HYPRE_Int HypreParVector::verbosity( const HYPRE_Int level )
222 {
223 return HYPRE_IJVectorSetPrintLevel( x, level );
224 }
225
226 void HypreParVector::Print( const char* fname ) const
227 {
228 HYPRE_IJVectorPrint( x, fname );
229 }
230
231 HypreParVector::~HypreParVector()
232 {
233 if( own_ParVector )
234 {
235 HYPRE_IJVectorDestroy( x );
236 }
237 }
238
239 double HypreParVector::InnerProduct( HypreParVector& x, HypreParVector& y )
240 {
241 return hypre_ParVectorInnerProd( x.x_par, y.x_par );
242 }
243
244 } // namespace moab
245
246 #endif
1.9.1.