FEI Version of the Day
fei_Aztec_LSVector.cpp
00001 /*
00002 // @HEADER
00003 // ************************************************************************
00004 //             FEI: Finite Element Interface to Linear Solvers
00005 //                  Copyright (2005) Sandia Corporation.
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
00008 // U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Alan Williams (william@sandia.gov) 
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 */
00042 
00043 
00044 #include <fei_trilinos_macros.hpp>
00045 #include <fei_iostream.hpp>
00046 
00047 #ifdef HAVE_FEI_AZTECOO
00048 
00049 #include <assert.h>
00050 #include <cmath>         // needed for declaration of sqrt, abs
00051 #include <unistd.h>
00052 #include <stdio.h>
00053 
00054 #include <fei_mpi.h>
00055 
00056 #ifndef FEI_SER
00057 
00058 #define AZTEC_MPI AZTEC_MPI
00059 #define AZ_MPI AZ_MPI
00060 #ifndef MPI
00061 #define MPI MPI
00062 #endif
00063 
00064 #endif
00065 
00066 #include <az_blas_wrappers.h>
00067 #include <az_aztec.h>
00068 #include <fei_SharedPtr.hpp>
00069 #include <fei_Aztec_Map.hpp>
00070 #include <fei_Aztec_LSVector.hpp>
00071 
00072 namespace fei_trilinos {
00073 
00074 //=============================================================================
00075 Aztec_LSVector::Aztec_LSVector(fei::SharedPtr<Aztec_Map> map, int* data_org)
00076  : amap_(map)
00077 {
00078 //  Aztec_LSVector::Aztec_LSVector -- construct a zero filled distributed vector 
00079 //  object.
00080 
00081     int tmp1 = data_org[AZ_N_internal];
00082     int tmp2 = data_org[AZ_N_border];
00083     length_ = tmp1 + tmp2 + data_org[AZ_N_external];
00084 
00085     localCoeffs_ = new double[length_];
00086 
00087     this->put(0.0);
00088 }
00089 
00091 Aztec_LSVector::~Aztec_LSVector(){
00092     delete [] localCoeffs_;
00093     localCoeffs_ = NULL;
00094 
00095     length_ = 0;
00096 }
00097 
00099 Aztec_LSVector::Aztec_LSVector(const Aztec_LSVector& source)
00100  : amap_(source.amap_)
00101 {
00102 
00105    length_ = source.length_;
00106    localCoeffs_ = new double[length_];
00107 
00108 //  Since virtual dispatching will not occur in a constructor, 
00109 //  specify call explicitly for clarity.
00110 
00111    Aztec_LSVector::operator=(source);
00112 }
00113 
00114 
00116 double Aztec_LSVector::dotProd (const Aztec_LSVector& y) const {
00117 
00120    int N_update = amap_->localSize();
00121 
00122    double *pv = (double*)localCoeffs_;
00123    double *py = (double*)y.startPointer();
00124 
00125    double dot = AZ_gdot(N_update, pv, py, amap_->getProcConfig());
00126 
00127    return dot;
00128 }
00129 
00130 
00132 void Aztec_LSVector::put(double scalar) {
00133 
00136    for (int i = 0; i < length_; i++)
00137       localCoeffs_[i] = scalar;
00138 }
00139 
00141 void Aztec_LSVector::scale (double s) {
00142 
00145    int N_update = amap_->localSize();
00146    int one = 1;
00147 
00148    DSCAL_F77(&N_update, &s, localCoeffs_, &one);
00149 
00150    return;
00151 }
00152 
00154 void Aztec_LSVector::addVec (double s, const Aztec_LSVector& c) {
00155 
00158    int N_update = amap_->localSize();
00159    int one = 1;
00160 
00161    double *pv = localCoeffs_;
00162    double *pc = (double*)c.startPointer();
00163 
00164    DAXPY_F77(&N_update,&s,pc,&one,pv,&one);
00165 
00166    return;
00167 }
00168 
00170 double Aztec_LSVector::norm (void) const  {
00171 
00174    int N_update = amap_->localSize();
00175 
00176    return(AZ_gvector_norm(N_update, 2,localCoeffs_, amap_->getProcConfig()));
00177 }
00178 
00180 double Aztec_LSVector::norm1 (void) const  {
00181 
00184    int N_update = amap_->localSize();
00185 
00186    return(AZ_gvector_norm(N_update, 1,localCoeffs_, amap_->getProcConfig()));
00187 }
00188  
00190 double& Aztec_LSVector::operator [] (int index)  {
00191 
00192 //  Aztec_LSVector::operator [] --- return non-const reference 
00193 
00194    int offset = amap_->localOffset();
00195 
00196    return(localCoeffs_[index - offset]);
00197 }
00198 
00200 const double& Aztec_LSVector::operator [] (int index) const  {
00201 
00202 // Aztec_LSVector::operator [] --- return const reference 
00203 
00204    int offset = amap_->localOffset();
00205 
00206    return(localCoeffs_[index - offset]);
00207 }
00208 
00210 Aztec_LSVector* Aztec_LSVector::newVector() const  {
00211 
00215     Aztec_LSVector* p = new Aztec_LSVector(*this);
00216 
00217     return p;
00218 }
00219 
00221 Aztec_LSVector& Aztec_LSVector::operator= (const Aztec_LSVector& rhs) {
00222 
00223 //    check for special case of a=a
00224    if (this != &rhs) {
00225       assign(rhs);
00226    }
00227         
00228    return(*this);
00229 }
00230 
00232 void Aztec_LSVector::assign(const Aztec_LSVector& rhs) {
00233 
00234    if ((amap_->globalSize() != rhs.amap_->globalSize()) ||
00235       (amap_->localSize() != rhs.amap_->localSize()) ) {
00236       fei::console_out() << "Aztec_LSVector::assign: ERROR, incompatible maps."
00237            << " Aborting assignment." << FEI_ENDL;
00238       return;
00239    }
00240 
00241    int N_update = amap_->localSize();
00242    double *pr = (double*)rhs.startPointer();
00243 
00244    for(int i=0; i<N_update; ++i) {
00245      localCoeffs_[i] = pr[i];
00246    }
00247 
00248    return;    
00249 }
00250 
00252 bool Aztec_LSVector::readFromFile(const char *fileName) {
00253 //
00254 //For now, this function will just use a simple (not very scalable)
00255 //strategy of having all processors read their own piece of the vector
00256 //from the file.
00257 //
00258 //Important note: This function assumes that the equation numbers in
00259 //the file are 1-based, and converts them to 0-based.
00260 //
00261    int globalSize = amap_->globalSize();
00262 
00263    int i=0, nn, nnz;
00264    double value;
00265    bool binaryData;
00266 
00267    char line[128];
00268 
00269    if (fileName == NULL) {
00270       fei::console_out() << "Aztec_LSVector::readFromFile: ERROR, NULL fileName." << FEI_ENDL;
00271       return(false);
00272    }
00273 
00274    if (strstr(fileName, ".txt") != NULL) {
00275       binaryData = false;
00276    }
00277    else {
00278       fei::console_out() << "Aztec_LSVector::readFromFile: fileName doesn't contain "
00279            << "'.txt', assuming binary data..." << FEI_ENDL;
00280       binaryData = true;
00281    }
00282 
00283    FILE *file = fopen(fileName,"r");
00284    if (!file){
00285       fei::console_out() << "Aztec_LSVector: Error opening " << fileName << FEI_ENDL;
00286       return false;
00287    }
00288 
00289    if (binaryData) {
00290       fread((char *)&nn,sizeof(int),1,file);
00291       fread((char *)&nnz,sizeof(int),1,file);
00292    }
00293    else {
00294       do {
00295          fgets(line,128,file);
00296       } while(strchr(line,'%'));
00297       sscanf(line,"%d",&nn);
00298    }
00299    if (nn != globalSize) {
00300       fei::console_out() << "ERROR in Aztec_LSVector::readFromFile." << FEI_ENDL;
00301       fei::console_out() << "Vector in file has wrong dimension." << FEI_ENDL;
00302       fei::console_out() << "amap_->globalSize():" << globalSize << " file n:" << nn << FEI_ENDL;
00303       return(false);
00304    }
00305 
00306    int start = amap_->localOffset();
00307    int end = start + amap_->localSize() - 1;
00308 
00309    while (!feof(file)) {
00310       if(binaryData) {
00311          fread((char *)&i,sizeof(int),1,file);
00312          fread((char *)&value,sizeof(double),1,file);
00313       }
00314       else {
00315          fgets(line,128,file);
00316          sscanf(line,"%d %le",&i,&value);
00317       }
00318       if(feof(file))break;
00319 
00320       if ((start <= i) && (i <= end)) {
00321          localCoeffs_[i - start] = value;
00322       }
00323    } //end of 'while(!feof)
00324 
00325    return true;
00326 }   
00327 
00329 bool Aztec_LSVector::writeToFile(const char *fileName) const {
00330    int i,p;
00331    int N_update = amap_->localSize();
00332    int start = amap_->localOffset();
00333    int numProcs = amap_->getProcConfig()[AZ_N_procs];
00334    int localRank = amap_->getProcConfig()[AZ_node];
00335    int masterRank = 0;
00336    MPI_Comm thisComm = amap_->getCommunicator();
00337 
00338    for(p=0; p<numProcs; p++){
00339 
00340       //A barrier inside the loop so each processor waits its turn.
00341       MPI_Barrier(thisComm);
00342 
00343       if (p == localRank){
00344          FILE *file = NULL;
00345 
00346          if (masterRank == localRank){
00347             //This is the master processor, open a new file.
00348             file = fopen(fileName,"w");
00349 
00350             //Write the vector dimension n into the file.
00351             fprintf(file,"%d\n",amap_->globalSize());
00352          }
00353          else {
00354             //This is not the master proc, open file for appending
00355             file = fopen(fileName,"a");
00356          }
00357 
00358          //Now loop over the local portion of the vector.
00359          for(i=0; i<N_update; i++) {
00360             fprintf(file,"%d %20.13e\n",start + i, localCoeffs_[i]);
00361          }
00362 
00363          fclose(file);
00364       }
00365    }
00366 
00367    return(true);
00368 }
00369 
00370 }//namespace fei_trilinos
00371 
00372 #endif
00373 //HAVE_FEI_AZTECOO
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends