FEI Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
fei_Aztec_LSVector.cpp
Go to the documentation of this file.
1/*
2// @HEADER
3// ************************************************************************
4// FEI: Finite Element Interface to Linear Solvers
5// Copyright (2005) Sandia Corporation.
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
8// U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Alan Williams (william@sandia.gov)
38//
39// ************************************************************************
40// @HEADER
41*/
42
43
45#include <fei_iostream.hpp>
46
47#ifdef HAVE_FEI_AZTECOO
48
49#include <assert.h>
50#include <cmath> // needed for declaration of sqrt, abs
51#include <unistd.h>
52#include <stdio.h>
53
54#include <fei_mpi.h>
55
56#ifndef FEI_SER
57
58#define AZTEC_MPI AZTEC_MPI
59#define AZ_MPI AZ_MPI
60#ifndef MPI
61#define MPI MPI
62#endif
63
64#endif
65
66#include <az_blas_wrappers.h>
67#include <az_aztec.h>
68#include <fei_SharedPtr.hpp>
69#include <fei_Aztec_Map.hpp>
71
72namespace fei_trilinos {
73
74//=============================================================================
76 : amap_(map)
77{
78// Aztec_LSVector::Aztec_LSVector -- construct a zero filled distributed vector
79// object.
80
81 int tmp1 = data_org[AZ_N_internal];
82 int tmp2 = data_org[AZ_N_border];
83 length_ = tmp1 + tmp2 + data_org[AZ_N_external];
84
85 localCoeffs_ = new double[length_];
86
87 this->put(0.0);
88}
89
91Aztec_LSVector::~Aztec_LSVector(){
92 delete [] localCoeffs_;
93 localCoeffs_ = NULL;
94
95 length_ = 0;
96}
97
99Aztec_LSVector::Aztec_LSVector(const Aztec_LSVector& source)
100 : amap_(source.amap_)
101{
102
105 length_ = source.length_;
106 localCoeffs_ = new double[length_];
107
108// Since virtual dispatching will not occur in a constructor,
109// specify call explicitly for clarity.
110
111 Aztec_LSVector::operator=(source);
112}
113
114
116double Aztec_LSVector::dotProd (const Aztec_LSVector& y) const {
117
120 int N_update = amap_->localSize();
121
122 double *pv = (double*)localCoeffs_;
123 double *py = (double*)y.startPointer();
124
125 double dot = AZ_gdot(N_update, pv, py, amap_->getProcConfig());
126
127 return dot;
128}
129
130
132void Aztec_LSVector::put(double scalar) {
133
136 for (int i = 0; i < length_; i++)
137 localCoeffs_[i] = scalar;
138}
139
141void Aztec_LSVector::scale (double s) {
142
145 int N_update = amap_->localSize();
146 int one = 1;
147
148 DSCAL_F77(&N_update, &s, localCoeffs_, &one);
149
150 return;
151}
152
154void Aztec_LSVector::addVec (double s, const Aztec_LSVector& c) {
155
158 int N_update = amap_->localSize();
159 int one = 1;
160
161 double *pv = localCoeffs_;
162 double *pc = (double*)c.startPointer();
163
164 DAXPY_F77(&N_update,&s,pc,&one,pv,&one);
165
166 return;
167}
168
170double Aztec_LSVector::norm (void) const {
171
174 int N_update = amap_->localSize();
175
176 return(AZ_gvector_norm(N_update, 2,localCoeffs_, amap_->getProcConfig()));
177}
178
180double Aztec_LSVector::norm1 (void) const {
181
184 int N_update = amap_->localSize();
185
186 return(AZ_gvector_norm(N_update, 1,localCoeffs_, amap_->getProcConfig()));
187}
188
190double& Aztec_LSVector::operator [] (int index) {
191
192// Aztec_LSVector::operator [] --- return non-const reference
193
194 int offset = amap_->localOffset();
195
196 return(localCoeffs_[index - offset]);
197}
198
200const double& Aztec_LSVector::operator [] (int index) const {
201
202// Aztec_LSVector::operator [] --- return const reference
203
204 int offset = amap_->localOffset();
205
206 return(localCoeffs_[index - offset]);
207}
208
210Aztec_LSVector* Aztec_LSVector::newVector() const {
211
215 Aztec_LSVector* p = new Aztec_LSVector(*this);
216
217 return p;
218}
219
221Aztec_LSVector& Aztec_LSVector::operator= (const Aztec_LSVector& rhs) {
222
223// check for special case of a=a
224 if (this != &rhs) {
225 assign(rhs);
226 }
227
228 return(*this);
229}
230
232void Aztec_LSVector::assign(const Aztec_LSVector& rhs) {
233
234 if ((amap_->globalSize() != rhs.amap_->globalSize()) ||
235 (amap_->localSize() != rhs.amap_->localSize()) ) {
236 fei::console_out() << "Aztec_LSVector::assign: ERROR, incompatible maps."
237 << " Aborting assignment." << FEI_ENDL;
238 return;
239 }
240
241 int N_update = amap_->localSize();
242 double *pr = (double*)rhs.startPointer();
243
244 for(int i=0; i<N_update; ++i) {
245 localCoeffs_[i] = pr[i];
246 }
247
248 return;
249}
250
252bool Aztec_LSVector::readFromFile(const char *fileName) {
253//
254//For now, this function will just use a simple (not very scalable)
255//strategy of having all processors read their own piece of the vector
256//from the file.
257//
258//Important note: This function assumes that the equation numbers in
259//the file are 1-based, and converts them to 0-based.
260//
261 int globalSize = amap_->globalSize();
262
263 int i=0, nn, nnz;
264 double value;
265 bool binaryData;
266
267 char line[128];
268
269 if (fileName == NULL) {
270 fei::console_out() << "Aztec_LSVector::readFromFile: ERROR, NULL fileName." << FEI_ENDL;
271 return(false);
272 }
273
274 if (strstr(fileName, ".txt") != NULL) {
275 binaryData = false;
276 }
277 else {
278 fei::console_out() << "Aztec_LSVector::readFromFile: fileName doesn't contain "
279 << "'.txt', assuming binary data..." << FEI_ENDL;
280 binaryData = true;
281 }
282
283 FILE *file = fopen(fileName,"r");
284 if (!file){
285 fei::console_out() << "Aztec_LSVector: Error opening " << fileName << FEI_ENDL;
286 return false;
287 }
288
289 if (binaryData) {
290 if (fread((char *)&nn,sizeof(int),1,file) != 1)
291 throw "fei_Aztec_LSVector.cpp: I/O error.";
292 if (fread((char *)&nnz,sizeof(int),1,file) != 1)
293 throw "fei_Aztec_LSVector.cpp: I/O error.";
294 }
295 else {
296 do {
297 if (fgets(line,128,file) == NULL)
298 throw "fei_Aztec_LSVector.cpp: I/O error.";
299 } while(strchr(line,'%'));
300 sscanf(line,"%d",&nn);
301 }
302 if (nn != globalSize) {
303 fei::console_out() << "ERROR in Aztec_LSVector::readFromFile." << FEI_ENDL;
304 fei::console_out() << "Vector in file has wrong dimension." << FEI_ENDL;
305 fei::console_out() << "amap_->globalSize():" << globalSize << " file n:" << nn << FEI_ENDL;
306 return(false);
307 }
308
309 int start = amap_->localOffset();
310 int end = start + amap_->localSize() - 1;
311
312 while (!feof(file)) {
313 if(binaryData) {
314 if (fread((char *)&i,sizeof(int),1,file) != 1)
315 throw "fei_Aztec_LSVector.cpp: I/O error.";
316 if (fread((char *)&value,sizeof(double),1,file) != 1)
317 throw "fei_Aztec_LSVector.cpp: I/O error.";
318 }
319 else {
320 if (fgets(line,128,file) == NULL)
321 throw "fei_Aztec_LSVector.cpp: I/O error.";
322 sscanf(line,"%d %le",&i,&value);
323 }
324 if(feof(file))break;
325
326 if ((start <= i) && (i <= end)) {
327 localCoeffs_[i - start] = value;
328 }
329 } //end of 'while(!feof)
330
331 return true;
332}
333
335bool Aztec_LSVector::writeToFile(const char *fileName) const {
336 int i,p;
337 int N_update = amap_->localSize();
338 int start = amap_->localOffset();
339 int numProcs = amap_->getProcConfig()[AZ_N_procs];
340 int localRank = amap_->getProcConfig()[AZ_node];
341 int masterRank = 0;
342 MPI_Comm thisComm = amap_->getCommunicator();
343
344 for(p=0; p<numProcs; p++){
345
346 //A barrier inside the loop so each processor waits its turn.
347 MPI_Barrier(thisComm);
348
349 if (p == localRank){
350 FILE *file = NULL;
351
352 if (masterRank == localRank){
353 //This is the master processor, open a new file.
354 file = fopen(fileName,"w");
355
356 //Write the vector dimension n into the file.
357 fprintf(file,"%d\n",amap_->globalSize());
358 }
359 else {
360 //This is not the master proc, open file for appending
361 file = fopen(fileName,"a");
362 }
363
364 //Now loop over the local portion of the vector.
365 for(i=0; i<N_update; i++) {
366 fprintf(file,"%d %20.13e\n",start + i, localCoeffs_[i]);
367 }
368
369 fclose(file);
370 }
371 }
372
373 return(true);
374}
375
376}//namespace fei_trilinos
377
378#endif
379//HAVE_FEI_AZTECOO
void PREFIX DAXPY_F77(const int *n, const double *alpha, const double x[], const int *incx, double y[], const int *incy)
void PREFIX DSCAL_F77(const int *n, const double *alpha, double *x, const int *incx)
Aztec_LSVector(fei::SharedPtr< Aztec_Map > map, int *data_org)
#define FEI_ENDL
#define MPI_Barrier(a)
Definition fei_mpi.h:61
#define MPI_Comm
Definition fei_mpi.h:56
std::ostream & console_out()
int numProcs(MPI_Comm comm)