00001
00002
00003
00004
00005 #include "PlayaEpetraMatrixFactory.hpp"
00006 #include "PlayaEpetraMatrix.hpp"
00007 #include "PlayaEpetraVector.hpp"
00008 #include "PlayaVectorSpaceDecl.hpp"
00009 #include "PlayaVectorDecl.hpp"
00010 #include "Teuchos_Array.hpp"
00011 #include "PlayaMPIComm.hpp"
00012 #include "PlayaLinearOperatorDecl.hpp"
00013
00014 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00015 #include "PlayaLinearOperatorImpl.hpp"
00016 #include "PlayaVectorImpl.hpp"
00017 #endif
00018
00019
00020 namespace Playa
00021 {
00022
00023 using namespace Teuchos;
00024
00025
00026 EpetraMatrixFactory::EpetraMatrixFactory(
00027 const RCP<const EpetraVectorSpace>& domain,
00028 const RCP<const EpetraVectorSpace>& range)
00029 : graph_(rcp(new Epetra_CrsGraph(Copy, *(range->epetraMap()), 0))),
00030 range_(range),
00031 domain_(domain)
00032 {}
00033
00034
00035 void EpetraMatrixFactory::finalize()
00036 {
00037 int ierr = graph_->FillComplete(*(domain_->epetraMap()), *(range_->epetraMap()));
00038
00039 TEUCHOS_TEST_FOR_EXCEPTION(ierr < 0, std::runtime_error,
00040 "EpetraMatrixFactory::finalize() failed during call "
00041 "to FillComplete(). Error code was " << ierr);
00042
00043 if (!graph_->StorageOptimized())
00044 {
00045 ierr = graph_->OptimizeStorage();
00046
00047 TEUCHOS_TEST_FOR_EXCEPTION(ierr < 0, std::runtime_error,
00048 "EpetraMatrixFactory::freezeValues() failed during call "
00049 "to OptimizeStorage(). Error code was " << ierr);
00050 }
00051 }
00052
00053 void EpetraMatrixFactory::initializeNonzerosInRow(int globalRowIndex,
00054 int nElemsToInsert,
00055 const int* globalColumnIndices)
00056 {
00057 int ierr = graph_->InsertGlobalIndices(globalRowIndex,
00058 nElemsToInsert,
00059 (int*) globalColumnIndices);
00060
00061 TEUCHOS_TEST_FOR_EXCEPTION(ierr < 0, std::runtime_error,
00062 "failed to add to row " << globalRowIndex
00063 << " in EpetraMatrixFactory::setRowValues() with nnz="
00064 << nElemsToInsert
00065 << ". Error code was " << ierr);
00066 }
00067
00068
00069 void EpetraMatrixFactory::initializeNonzeroBatch(int numRows,
00070 int rowBlockSize,
00071 const int* globalRowIndices,
00072 int numColumnsPerRow,
00073 const int* globalColumnIndices,
00074 const int* skipRow)
00075 {
00076 int numRowBlocks = numRows/rowBlockSize;
00077 int row = 0;
00078
00079 for (int rb=0; rb<numRowBlocks; rb++)
00080 {
00081 const int* cols = globalColumnIndices + rb*numColumnsPerRow;
00082 for (int r=0; r<rowBlockSize; r++, row++)
00083 {
00084 if (skipRow[row]) continue;
00085 graph_->InsertGlobalIndices(globalRowIndices[row],
00086 numColumnsPerRow, (int*) cols);
00087 }
00088 }
00089 }
00090
00091
00092 void EpetraMatrixFactory::configure(int lowestRow,
00093 const std::vector<int>& rowPtrs,
00094 const std::vector<int>& nnzPerRow,
00095 const std::vector<int>& data)
00096 {
00097
00098 graph_ = rcp(new Epetra_CrsGraph(Copy, *(range_->epetraMap()),
00099 (const int*) &(nnzPerRow[0]),
00100 true));
00101
00102 for (unsigned int i=0; i<rowPtrs.size(); i++)
00103 {
00104 graph_->InsertGlobalIndices(lowestRow + i, nnzPerRow[i],
00105 (int*) &(data[rowPtrs[i]]));
00106 }
00107
00108 finalize();
00109 }
00110
00111 const Epetra_CrsGraph& EpetraMatrixFactory::graph() const
00112 {
00113 return *(graph_.get());
00114 }
00115
00116
00117 LinearOperator<double> EpetraMatrixFactory::createMatrix() const
00118 {
00119 RCP<const VectorSpaceBase<double> > dp = epDomain();
00120 RCP<const VectorSpaceBase<double> > rp = epRange();
00121 VectorSpace<double> d = dp;
00122 VectorSpace<double> r = rp;
00123
00124 RCP<LinearOperatorBase<double> > A
00125 = rcp(new EpetraMatrix(graph(), d, r));
00126 return A;
00127 }
00128
00129 }