FEI Version of the Day
fei_Aztec_LSVector.cpp
00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 #include <fei_trilinos_macros.hpp>
00010 #include <fei_iostream.hpp>
00011 
00012 #ifdef HAVE_FEI_AZTECOO
00013 
00014 #include <assert.h>
00015 #include <cmath>         // needed for declaration of sqrt, abs
00016 #include <unistd.h>
00017 #include <stdio.h>
00018 
00019 #include <fei_mpi.h>
00020 
00021 #ifndef FEI_SER
00022 
00023 #define AZTEC_MPI AZTEC_MPI
00024 #define AZ_MPI AZ_MPI
00025 #ifndef MPI
00026 #define MPI MPI
00027 #endif
00028 
00029 #endif
00030 
00031 #include <az_blas_wrappers.h>
00032 #include <az_aztec.h>
00033 #include <fei_SharedPtr.hpp>
00034 #include <fei_Aztec_Map.hpp>
00035 #include <fei_Aztec_LSVector.hpp>
00036 
00037 namespace fei_trilinos {
00038 
00039 //=============================================================================
00040 Aztec_LSVector::Aztec_LSVector(fei::SharedPtr<Aztec_Map> map, int* data_org)
00041  : amap_(map)
00042 {
00043 //  Aztec_LSVector::Aztec_LSVector -- construct a zero filled distributed vector 
00044 //  object.
00045 
00046     int tmp1 = data_org[AZ_N_internal];
00047     int tmp2 = data_org[AZ_N_border];
00048     length_ = tmp1 + tmp2 + data_org[AZ_N_external];
00049 
00050     localCoeffs_ = new double[length_];
00051 
00052     this->put(0.0);
00053 }
00054 
00056 Aztec_LSVector::~Aztec_LSVector(){
00057     delete [] localCoeffs_;
00058     localCoeffs_ = NULL;
00059 
00060     length_ = 0;
00061 }
00062 
00064 Aztec_LSVector::Aztec_LSVector(const Aztec_LSVector& source)
00065  : amap_(source.amap_)
00066 {
00067 
00070    length_ = source.length_;
00071    localCoeffs_ = new double[length_];
00072 
00073 //  Since virtual dispatching will not occur in a constructor, 
00074 //  specify call explicitly for clarity.
00075 
00076    Aztec_LSVector::operator=(source);
00077 }
00078 
00079 
00081 double Aztec_LSVector::dotProd (const Aztec_LSVector& y) const {
00082 
00085    int N_update = amap_->localSize();
00086 
00087    double *pv = (double*)localCoeffs_;
00088    double *py = (double*)y.startPointer();
00089 
00090    double dot = AZ_gdot(N_update, pv, py, amap_->getProcConfig());
00091 
00092    return dot;
00093 }
00094 
00095 
00097 void Aztec_LSVector::put(double scalar) {
00098 
00101    for (int i = 0; i < length_; i++)
00102       localCoeffs_[i] = scalar;
00103 }
00104 
00106 void Aztec_LSVector::scale (double s) {
00107 
00110    int N_update = amap_->localSize();
00111    int one = 1;
00112 
00113    DSCAL_F77(&N_update, &s, localCoeffs_, &one);
00114 
00115    return;
00116 }
00117 
00119 void Aztec_LSVector::addVec (double s, const Aztec_LSVector& c) {
00120 
00123    int N_update = amap_->localSize();
00124    int one = 1;
00125 
00126    double *pv = localCoeffs_;
00127    double *pc = (double*)c.startPointer();
00128 
00129    DAXPY_F77(&N_update,&s,pc,&one,pv,&one);
00130 
00131    return;
00132 }
00133 
00135 double Aztec_LSVector::norm (void) const  {
00136 
00139    int N_update = amap_->localSize();
00140 
00141    return(AZ_gvector_norm(N_update, 2,localCoeffs_, amap_->getProcConfig()));
00142 }
00143 
00145 double Aztec_LSVector::norm1 (void) const  {
00146 
00149    int N_update = amap_->localSize();
00150 
00151    return(AZ_gvector_norm(N_update, 1,localCoeffs_, amap_->getProcConfig()));
00152 }
00153  
00155 double& Aztec_LSVector::operator [] (int index)  {
00156 
00157 //  Aztec_LSVector::operator [] --- return non-const reference 
00158 
00159    int offset = amap_->localOffset();
00160 
00161    return(localCoeffs_[index - offset]);
00162 }
00163 
00165 const double& Aztec_LSVector::operator [] (int index) const  {
00166 
00167 // Aztec_LSVector::operator [] --- return const reference 
00168 
00169    int offset = amap_->localOffset();
00170 
00171    return(localCoeffs_[index - offset]);
00172 }
00173 
00175 Aztec_LSVector* Aztec_LSVector::newVector() const  {
00176 
00180     Aztec_LSVector* p = new Aztec_LSVector(*this);
00181 
00182     return p;
00183 }
00184 
00186 Aztec_LSVector& Aztec_LSVector::operator= (const Aztec_LSVector& rhs) {
00187 
00188 //    check for special case of a=a
00189    if (this != &rhs) {
00190       assign(rhs);
00191    }
00192         
00193    return(*this);
00194 }
00195 
00197 void Aztec_LSVector::assign(const Aztec_LSVector& rhs) {
00198 
00199    if ((amap_->globalSize() != rhs.amap_->globalSize()) ||
00200       (amap_->localSize() != rhs.amap_->localSize()) ) {
00201       FEI_CERR << "Aztec_LSVector::assign: ERROR, incompatible maps."
00202            << " Aborting assignment." << FEI_ENDL;
00203       return;
00204    }
00205 
00206    int N_update = amap_->localSize();
00207    double *pr = (double*)rhs.startPointer();
00208 
00209    for(int i=0; i<N_update; ++i) {
00210      localCoeffs_[i] = pr[i];
00211    }
00212 
00213    return;    
00214 }
00215 
00217 bool Aztec_LSVector::readFromFile(const char *fileName) {
00218 //
00219 //For now, this function will just use a simple (not very scalable)
00220 //strategy of having all processors read their own piece of the vector
00221 //from the file.
00222 //
00223 //Important note: This function assumes that the equation numbers in
00224 //the file are 1-based, and converts them to 0-based.
00225 //
00226    int globalSize = amap_->globalSize();
00227 
00228    int i=0, nn, nnz;
00229    double value;
00230    bool binaryData;
00231 
00232    char line[128];
00233 
00234    if (fileName == NULL) {
00235       FEI_CERR << "Aztec_LSVector::readFromFile: ERROR, NULL fileName." << FEI_ENDL;
00236       return(false);
00237    }
00238 
00239    if (strstr(fileName, ".txt") != NULL) {
00240       binaryData = false;
00241    }
00242    else {
00243       FEI_CERR << "Aztec_LSVector::readFromFile: fileName doesn't contain "
00244            << "'.txt', assuming binary data..." << FEI_ENDL;
00245       binaryData = true;
00246    }
00247 
00248    FILE *file = fopen(fileName,"r");
00249    if (!file){
00250       FEI_CERR << "Aztec_LSVector: Error opening " << fileName << FEI_ENDL;
00251       return false;
00252    }
00253 
00254    if (binaryData) {
00255       fread((char *)&nn,sizeof(int),1,file);
00256       fread((char *)&nnz,sizeof(int),1,file);
00257    }
00258    else {
00259       do {
00260          fgets(line,128,file);
00261       } while(strchr(line,'%'));
00262       sscanf(line,"%d",&nn);
00263    }
00264    if (nn != globalSize) {
00265       FEI_CERR << "ERROR in Aztec_LSVector::readFromFile." << FEI_ENDL;
00266       FEI_CERR << "Vector in file has wrong dimension." << FEI_ENDL;
00267       FEI_CERR << "amap_->globalSize():" << globalSize << " file n:" << nn << FEI_ENDL;
00268       return(false);
00269    }
00270 
00271    int start = amap_->localOffset();
00272    int end = start + amap_->localSize() - 1;
00273 
00274    while (!feof(file)) {
00275       if(binaryData) {
00276          fread((char *)&i,sizeof(int),1,file);
00277          fread((char *)&value,sizeof(double),1,file);
00278       }
00279       else {
00280          fgets(line,128,file);
00281          sscanf(line,"%d %le",&i,&value);
00282       }
00283       if(feof(file))break;
00284 
00285       if ((start <= i) && (i <= end)) {
00286          localCoeffs_[i - start] = value;
00287       }
00288    } //end of 'while(!feof)
00289 
00290    return true;
00291 }   
00292 
00294 bool Aztec_LSVector::writeToFile(const char *fileName) const {
00295    int i,p;
00296    int N_update = amap_->localSize();
00297    int start = amap_->localOffset();
00298    int numProcs = amap_->getProcConfig()[AZ_N_procs];
00299    int localRank = amap_->getProcConfig()[AZ_node];
00300    int masterRank = 0;
00301    MPI_Comm thisComm = amap_->getCommunicator();
00302 
00303    for(p=0; p<numProcs; p++){
00304 
00305       //A barrier inside the loop so each processor waits its turn.
00306       MPI_Barrier(thisComm);
00307 
00308       if (p == localRank){
00309          FILE *file = NULL;
00310 
00311          if (masterRank == localRank){
00312             //This is the master processor, open a new file.
00313             file = fopen(fileName,"w");
00314 
00315             //Write the vector dimension n into the file.
00316             fprintf(file,"%d\n",amap_->globalSize());
00317          }
00318          else {
00319             //This is not the master proc, open file for appending
00320             file = fopen(fileName,"a");
00321          }
00322 
00323          //Now loop over the local portion of the vector.
00324          for(i=0; i<N_update; i++) {
00325             fprintf(file,"%d %20.13e\n",start + i, localCoeffs_[i]);
00326          }
00327 
00328          fclose(file);
00329       }
00330    }
00331 
00332    return(true);
00333 }
00334 
00335 }//namespace fei_trilinos
00336 
00337 #endif
00338 //HAVE_FEI_AZTECOO
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends