FEI Version of the Day
fei_LinProbMgr_EpetraBasic.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_SparseRowGraph.hpp>
00046 
00047 #ifdef HAVE_FEI_EPETRA
00048 
00049 #include <fei_LinProbMgr_EpetraBasic.hpp>
00050 
00051 LinProbMgr_EpetraBasic::LinProbMgr_EpetraBasic(MPI_Comm comm)
00052  : comm_(comm),
00053    ownedRows_(),
00054    epetra_comm_(),
00055    epetra_rowmap_(),
00056    fei_srgraph_(),
00057    crsgraph_(),
00058    A_(),
00059    numVectors_(1),
00060    x_(),
00061    b_()
00062 {
00063 #ifdef FEI_SER
00064   fei::SharedPtr<Epetra_Comm> ecomm(new Epetra_SerialComm);
00065 #else
00066   fei::SharedPtr<Epetra_Comm> ecomm(new Epetra_MpiComm(comm));
00067 #endif
00068 
00069   epetra_comm_ = ecomm;
00070 }
00071 
00072 LinProbMgr_EpetraBasic::~LinProbMgr_EpetraBasic()
00073 {
00074 }
00075 
00076 
00077 void LinProbMgr_EpetraBasic
00078 ::setRowDistribution(const std::vector<int>& ownedGlobalRows)
00079 {
00080   if (ownedGlobalRows == ownedRows_) {
00081     return;
00082   }
00083 
00084   if (!ownedRows_.empty()) {
00085     throw std::runtime_error("setRowDistribution called multiple times with different distributions. not allowed.");
00086   }
00087 
00088   int* rows = const_cast<int*>(&ownedGlobalRows[0]);
00089   epetra_rowmap_.reset(new Epetra_Map(-1, ownedGlobalRows.size(),
00090                                       rows, 0, //indexBase
00091                                       *epetra_comm_));
00092 
00093   x_.reset(new Epetra_MultiVector(*epetra_rowmap_, numVectors_));
00094 
00095   b_.reset(new Epetra_MultiVector(*epetra_rowmap_, numVectors_));
00096 }
00097 
00098 void LinProbMgr_EpetraBasic
00099 ::setMatrixGraph(fei::SharedPtr<fei::SparseRowGraph> matrixGraph)
00100 {
00101   if (fei_srgraph_.get() != NULL) {
00102     if (*fei_srgraph_ != *matrixGraph) {
00103       throw std::runtime_error("setMatrixGraph called multiple times with different graphs. not allowed.");
00104     }
00105     return;
00106   }
00107   else {
00108     fei_srgraph_ = matrixGraph;
00109   }
00110 
00111   if (epetra_rowmap_.get() == NULL) {
00112     setRowDistribution(matrixGraph->rowNumbers);
00113   }
00114 
00115   if ((int)fei_srgraph_->rowNumbers.size() != epetra_rowmap_->NumMyElements()) {
00116     throw std::runtime_error("setMatrixGraph: num-rows not consistent with value from setRowDistribution");
00117   }
00118 
00119   //We'll create and populate a Epetra_CrsGraph object.
00120 
00121   std::vector<int>& rowNumbers = fei_srgraph_->rowNumbers;
00122   std::vector<int>& rowOffsets = fei_srgraph_->rowOffsets;
00123 
00124   //First create an array of num-indices-per-row.
00125   std::vector<int> numIndicesPerRow; numIndicesPerRow.reserve(rowNumbers.size());
00126 
00127   int err;
00128   unsigned i;
00129   for(i=0; i<numIndicesPerRow.size(); ++i) {
00130     numIndicesPerRow.push_back(rowOffsets[i+1] - rowOffsets[i]);
00131   }
00132 
00133   bool static_profile = true;
00134 
00135   crsgraph_.reset(new Epetra_CrsGraph(Copy, *epetra_rowmap_,
00136                   &numIndicesPerRow[0], static_profile));
00137 
00138   //Now put in all the column-indices
00139 
00140   std::vector<int>& colIndices = fei_srgraph_->packedColumnIndices;
00141   for(i=0; i<rowNumbers.size(); ++i) {
00142     int offset = rowOffsets[i];
00143     err = crsgraph_->InsertGlobalIndices(rowNumbers[i], numIndicesPerRow[i],
00144                                         &colIndices[offset]);
00145     if (err != 0) {
00146       throw std::runtime_error("setMatrixGraph: err from Epetra_CrsGraph::InsertGlobalIndices.");
00147     }
00148   }
00149 
00150   err = crsgraph_->FillComplete();
00151   if (err != 0) {
00152     throw std::runtime_error("setMatrixGraph: err from Epetra_CrsGraph::FillComplete.");
00153   }
00154  
00155   //and finally, create a matrix.
00156   A_.reset(new Epetra_CrsMatrix(Copy, *crsgraph_));
00157 }
00158 
00159 void LinProbMgr_EpetraBasic::setMatrixValues(double scalar)
00160 {
00161   int err = A_->PutScalar(scalar);
00162   if (err != 0) {
00163     throw std::runtime_error("error in Epetra_CrsMatrix->PutScalar");
00164   }
00165 }
00166 
00167 void LinProbMgr_EpetraBasic::setVectorValues(double scalar,
00168                                              bool soln_vector)
00169 {
00170   int err = soln_vector ?
00171     x_->PutScalar(scalar) : b_->PutScalar(scalar);
00172   if (err != 0) {
00173     throw std::runtime_error("error in Epetra_MultiVector->PutScalar");
00174   }
00175 }
00176 
00177 int LinProbMgr_EpetraBasic::getLocalNumRows()
00178 {
00179   if (epetra_rowmap_.get() == NULL) return(-1);
00180 
00181   return(epetra_rowmap_->NumMyElements());
00182 }
00183 
00184 int LinProbMgr_EpetraBasic::getRowLength(int row)
00185 {
00186   if (A_.get() == NULL) return(-1);
00187 
00188   return( A_->NumGlobalEntries(row) );
00189 }
00190 
00191 int LinProbMgr_EpetraBasic::copyOutMatrixRow(int row, int len,
00192                                              double* coefs, int* indices)
00193 {
00194   int dummy;
00195   return( A_->ExtractGlobalRowCopy(row, len, dummy, coefs, indices) );
00196 }
00197 
00198 int LinProbMgr_EpetraBasic::insertMatrixValues(int numRows, const int* rows,
00199                                                int numCols, const int* cols,
00200                                                const double* const* values,
00201                                                bool sum_into)
00202 {
00203   int* nc_cols = const_cast<int*>(cols);
00204   double** nc_values = const_cast<double**>(values);
00205   int err = 0;
00206   if (sum_into) {
00207     for(int i=0; i<numRows; ++i) {
00208       err = A_->SumIntoGlobalValues(rows[i], numCols, nc_values[i], nc_cols);
00209       if (err < 0) {
00210         return(err);
00211       }
00212     }
00213   }
00214   else {
00215     for(int i=0; i<numRows; ++i) {
00216       err = A_->ReplaceGlobalValues(rows[i], numCols, nc_values[i], nc_cols);
00217       if (err < 0) {
00218         return(err);
00219       }
00220     }
00221   }
00222   return(err);
00223 }
00224 
00225 int LinProbMgr_EpetraBasic::insertVectorValues(int numValues,
00226                                                const int* globalIndices,
00227                                                const double* values,
00228                                                bool sum_into,
00229                                                bool soln_vector,
00230                                                int vectorIndex)
00231 {
00232   double* localvaluesptr = soln_vector ?
00233     x_->Pointers()[vectorIndex] : b_->Pointers()[vectorIndex];
00234 
00235   int min_my_gid = epetra_rowmap_->MinMyGID();
00236   int returnValue = 0;
00237 
00238   if (sum_into) {
00239     for(int i=0; i<numValues; ++i) {
00240       int offset = globalIndices[i] - min_my_gid;
00241       if (offset < 0) {
00242         returnValue = 1;
00243         continue;
00244       }
00245       localvaluesptr[offset] += values[i];
00246     }
00247   }
00248   else {
00249     for(int i=0; i<numValues; ++i) {
00250       int offset = globalIndices[i] - min_my_gid;
00251       if (offset < 0) {
00252         returnValue = 1;
00253         continue;
00254       }
00255       localvaluesptr[offset] = values[i];
00256     }
00257   }
00258 
00259   return(returnValue);
00260 }
00261 
00262 int LinProbMgr_EpetraBasic::copyOutVectorValues(int numValues,
00263                                                 const int* globalIndices,
00264                                                 double* values,
00265                                                 bool soln_vector,
00266                                                 int vectorIndex)
00267 {
00268   double* localvaluesptr = soln_vector ?
00269     x_->Pointers()[vectorIndex] : b_->Pointers()[vectorIndex];
00270 
00271   int min_my_gid = epetra_rowmap_->MinMyGID();
00272 
00273   for(int i=0; i<numValues; ++i) {
00274     int offset = globalIndices[i] - min_my_gid;
00275     values[i] = localvaluesptr[offset];
00276   }
00277   return(0);
00278 }
00279 
00280 double* LinProbMgr_EpetraBasic::getLocalVectorValuesPtr(bool soln_vector,
00281                                                         int vectorIndex)
00282 {
00283   double* localvaluesptr = soln_vector ?
00284     x_->Pointers()[vectorIndex] : b_->Pointers()[vectorIndex];
00285 
00286   return(localvaluesptr);
00287 }
00288 
00289 int LinProbMgr_EpetraBasic::globalAssemble()
00290 {
00291   if (!A_->Filled()) {
00292     //assumes the matrix is square!
00293     int err = A_->FillComplete();
00294     if (err != 0) {
00295       return(err);
00296     }
00297   }
00298 
00299   if (!A_->StorageOptimized()) {
00300     A_->OptimizeStorage();
00301   }
00302 
00303   return(0);
00304 }
00305 
00306 fei::SharedPtr<Epetra_CrsMatrix>
00307 LinProbMgr_EpetraBasic::get_A_matrix()
00308 {
00309   return( A_ );
00310 }
00311 
00312 fei::SharedPtr<Epetra_MultiVector>
00313 LinProbMgr_EpetraBasic::get_rhs_vector()
00314 {
00315   return( b_ );
00316 }
00317 
00318 fei::SharedPtr<Epetra_MultiVector>
00319 LinProbMgr_EpetraBasic::get_solution_vector()
00320 {
00321   return( x_ );
00322 }
00323 
00324 //HAVE_FEI_EPETRA
00325 #endif
00326 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends