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

Generated on Wed May 12 21:30:41 2010 for FEI by  doxygen 1.4.7