fei_Aztec_Vector.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_Aztec_Map.hpp>
00034 #include <fei_Aztec_Vector.hpp>
00035 
00036 namespace fei_trilinos {
00037 
00038 //=============================================================================
00039 Aztec_Vector::Aztec_Vector(const Aztec_Map& map, int* data_org)
00040  : amap_(map)
00041 {
00042 //  Aztec_Vector::Aztec_Vector -- construct a zero filled distributed vector 
00043 //  object.
00044 
00045     int tmp1 = data_org[AZ_N_internal];
00046     int tmp2 = data_org[AZ_N_border];
00047     length_ = tmp1 + tmp2 + data_org[AZ_N_external];
00048 
00049     localCoeffs_ = new double[length_];
00050 
00051     this->put(0.0);
00052 }
00053 
00055 Aztec_Vector::~Aztec_Vector(){
00056     delete [] localCoeffs_;
00057     localCoeffs_ = NULL;
00058 
00059     length_ = 0;
00060 }
00061 
00063 Aztec_Vector::Aztec_Vector(const Aztec_Vector& source)
00064  : amap_(source.amap_)
00065 {
00066 
00069    length_ = source.length_;
00070    localCoeffs_ = new double[length_];
00071 
00072 //  Since virtual dispatching will not occur in a constructor, 
00073 //  specify call explicitly for clarity.
00074 
00075    Aztec_Vector::operator=(source);
00076 }
00077 
00078 
00080 double Aztec_Vector::dotProd (const Aztec_Vector& y) const {
00081 
00084    int N_update = amap_.localSize();
00085 
00086    double *pv = (double*)localCoeffs_;
00087    double *py = (double*)y.startPointer();
00088 
00089    double dot = AZ_gdot(N_update, pv, py, amap_.getProcConfig());
00090 
00091    return dot;
00092 }
00093 
00094 
00096 void Aztec_Vector::put(double scalar) {
00097 
00100    for (int i = 0; i < length_; i++)
00101       localCoeffs_[i] = scalar;
00102 }
00103 
00105 void Aztec_Vector::scale (double s) {
00106 
00109    int N_update = amap_.localSize();
00110    int one = 1;
00111 
00112    DSCAL_F77(&N_update, &s, localCoeffs_, &one);
00113 
00114    return;
00115 }
00116 
00118 void Aztec_Vector::linComb(const Aztec_Vector& b, double s,
00119                            const Aztec_Vector& c) {
00120 
00123    int N_update = amap_.localSize();
00124 
00125    double *pv = localCoeffs_;
00126    const double *pb = (double*)b.startPointer();
00127    const double *pc = (double*)c.startPointer();
00128    for (int i = 0; i < N_update; i++)
00129       pv[i] = pb[i] + s * pc[i];
00130 
00131    return;
00132 }
00133 
00135 void Aztec_Vector::addVec (double s, const Aztec_Vector& c) {
00136 
00139    int N_update = amap_.localSize();
00140    int one = 1;
00141 
00142    double *pv = localCoeffs_;
00143    double *pc = (double*)c.startPointer();
00144 
00145    DAXPY_F77(&N_update,&s,pc,&one,pv,&one);
00146 
00147    return;
00148 }
00149 
00151 double Aztec_Vector::norm (void) const  {
00152 
00155    int N_update = amap_.localSize();
00156 
00157    return(AZ_gvector_norm(N_update, 2,localCoeffs_, amap_.getProcConfig()));
00158 }
00159 
00161 double Aztec_Vector::norm1 (void) const  {
00162 
00165    int N_update = amap_.localSize();
00166 
00167    return(AZ_gvector_norm(N_update, 1,localCoeffs_, amap_.getProcConfig()));
00168 }
00169  
00171 double& Aztec_Vector::operator [] (int index)  {
00172 
00173 //  Aztec_Vector::operator [] --- return non-const reference 
00174 
00175    int offset = amap_.localOffset();
00176 
00177    return(localCoeffs_[index - offset]);
00178 }
00179 
00181 const double& Aztec_Vector::operator [] (int index) const  {
00182 
00183 // Aztec_Vector::operator [] --- return const reference 
00184 
00185    int offset = amap_.localOffset();
00186 
00187    return(localCoeffs_[index - offset]);
00188 }
00189 
00191 Aztec_Vector* Aztec_Vector::newVector() const  {
00192 
00196     Aztec_Vector* p = new Aztec_Vector(*this);
00197 
00198     return p;
00199 }
00200 
00202 Aztec_Vector& Aztec_Vector::operator= (const Aztec_Vector& rhs) {
00203 
00204 //    check for special case of a=a
00205    if (this != &rhs) {
00206       assign(rhs);
00207    }
00208         
00209    return(*this);
00210 }
00211 
00213 void Aztec_Vector::assign(const Aztec_Vector& rhs) {
00214 
00215    if ((amap_.globalSize() != rhs.amap_.globalSize()) ||
00216       (amap_.localSize() != rhs.amap_.localSize()) ) {
00217       FEI_CERR << "Aztec_Vector::assign: ERROR, incompatible maps."
00218            << " Aborting assignment." << FEI_ENDL;
00219       return;
00220    }
00221 
00222    int N_update = amap_.localSize();
00223    double *pr = (double*)rhs.startPointer();
00224 
00225    for(int i=0; i<N_update; ++i) {
00226      localCoeffs_[i] = pr[i];
00227    }
00228 
00229    return;    
00230 }
00231 
00233 bool Aztec_Vector::readFromFile(const char *fileName) {
00234 //
00235 //For now, this function will just use a simple (not very scalable)
00236 //strategy of having all processors read their own piece of the vector
00237 //from the file.
00238 //
00239 //Important note: This function assumes that the equation numbers in
00240 //the file are 1-based, and converts them to 0-based.
00241 //
00242    int globalSize = amap_.globalSize();
00243 
00244    int i=0, nn, nnz;
00245    double value;
00246    bool binaryData;
00247 
00248    char line[128];
00249 
00250    if (fileName == NULL) {
00251       FEI_CERR << "Aztec_Vector::readFromFile: ERROR, NULL fileName." << FEI_ENDL;
00252       return(false);
00253    }
00254 
00255    if (strstr(fileName, ".txt") != NULL) {
00256       binaryData = false;
00257    }
00258    else {
00259       FEI_CERR << "Aztec_Vector::readFromFile: fileName doesn't contain "
00260            << "'.txt', assuming binary data..." << FEI_ENDL;
00261       binaryData = true;
00262    }
00263 
00264    FILE *file = fopen(fileName,"r");
00265    if (!file){
00266       FEI_CERR << "Aztec_Vector: Error opening " << fileName << FEI_ENDL;
00267       return false;
00268    }
00269 
00270    if (binaryData) {
00271       fread((char *)&nn,sizeof(int),1,file);
00272       fread((char *)&nnz,sizeof(int),1,file);
00273    }
00274    else {
00275       do {
00276          fgets(line,128,file);
00277       } while(strchr(line,'%'));
00278       sscanf(line,"%d",&nn);
00279    }
00280    if (nn != globalSize) {
00281       FEI_CERR << "ERROR in Aztec_Vector::readFromFile." << FEI_ENDL;
00282       FEI_CERR << "Vector in file has wrong dimension." << FEI_ENDL;
00283       FEI_CERR << "amap_.globalSize():" << globalSize << " file n:" << nn << FEI_ENDL;
00284       return(false);
00285    }
00286 
00287    int start = amap_.localOffset();
00288    int end = start + amap_.localSize() - 1;
00289 
00290    while (!feof(file)) {
00291       if(binaryData) {
00292          fread((char *)&i,sizeof(int),1,file);
00293          fread((char *)&value,sizeof(double),1,file);
00294       }
00295       else {
00296          fgets(line,128,file);
00297          sscanf(line,"%d %le",&i,&value);
00298       }
00299       if(feof(file))break;
00300 
00301       if ((start <= i) && (i <= end)) {
00302          localCoeffs_[i - start] = value;
00303       }
00304    } //end of 'while(!feof)
00305 
00306    return true;
00307 }   
00308 
00310 bool Aztec_Vector::writeToFile(const char *fileName) const {
00311    int i,p;
00312    int N_update = amap_.localSize();
00313    int start = amap_.localOffset();
00314    int numProcs = amap_.getProcConfig()[AZ_N_procs];
00315    int localRank = amap_.getProcConfig()[AZ_node];
00316    int masterRank = 0;
00317    MPI_Comm thisComm = amap_.getCommunicator();
00318 
00319    for(p=0; p<numProcs; p++){
00320 
00321       //A barrier inside the loop so each processor waits its turn.
00322       MPI_Barrier(thisComm);
00323 
00324       if (p == localRank){
00325          FILE *file = NULL;
00326 
00327          if (masterRank == localRank){
00328             //This is the master processor, open a new file.
00329             file = fopen(fileName,"w");
00330 
00331             //Write the vector dimension n into the file.
00332             fprintf(file,"%d\n",amap_.globalSize());
00333          }
00334          else {
00335             //This is not the master proc, open file for appending
00336             file = fopen(fileName,"a");
00337          }
00338 
00339          //Now loop over the local portion of the vector.
00340          for(i=0; i<N_update; i++) {
00341             fprintf(file,"%d %20.13e\n",start + i, localCoeffs_[i]);
00342          }
00343 
00344          fclose(file);
00345       }
00346    }
00347 
00348    return(true);
00349 }
00350 
00351 }//namespace fei_trilinos
00352 
00353 #endif
00354 //HAVE_FEI_AZTECOO

Generated on Tue Jul 13 09:27:44 2010 for FEI by  doxygen 1.4.7