fei_LinProbMgr_EpetraBasic.cpp

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

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