fei_Matrix_Local.cpp

00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2007 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_ParameterSet.hpp>
00010 #include "fei_Matrix_Local.hpp"
00011 #include "fei_Matrix_core.hpp"
00012 #include "fei_sstream.hpp"
00013 #include "fei_fstream.hpp"
00014 
00015 namespace fei {
00016 
00017 Matrix_Local::Matrix_Local(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00018                fei::SharedPtr<fei::SparseRowGraph> sparseRowGraph)
00019  : matrixGraph_(matrixGraph),
00020    sparseRowGraph_(sparseRowGraph),
00021    coefs_(sparseRowGraph->packedColumnIndices.size(), 0.0),
00022    stateChanged_(false),
00023    work_data1D_(),
00024    work_data2D_()
00025 {
00026 }
00027 
00028 Matrix_Local::~Matrix_Local()
00029 {
00030 }
00031 
00032 //----------------------------------------------------------------------------
00033 fei::SharedPtr<fei::Matrix>
00034 Matrix_Local::create_Matrix_Local(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00035                                 bool blockEntry)
00036 {
00037   fei::SharedPtr<fei::SparseRowGraph> srg =
00038     matrixGraph->createGraph(blockEntry, true);
00039   fei::SharedPtr<fei::Matrix> mat(new fei::Matrix_Local(matrixGraph, srg));
00040   return(mat);
00041 }
00042 
00043 //----------------------------------------------------------------------------
00044 const char*
00045 Matrix_Local::typeName()
00046 { return( "fei::Matrix_Local" ); }
00047 
00048 //----------------------------------------------------------------------------
00049 int
00050 Matrix_Local::parameters(const fei::ParameterSet& /*paramset*/)
00051 {
00052   return(0);
00053 }
00054 
00055 //----------------------------------------------------------------------------
00056 int
00057 Matrix_Local::parameters(int /*numParams*/, const char* const* /*paramStrings*/)
00058 {
00059   return(0);
00060 }
00061 
00062 fei::SharedPtr<fei::MatrixGraph>
00063 Matrix_Local::getMatrixGraph()
00064 { return( matrixGraph_ ); }
00065 
00066 void
00067 Matrix_Local::setMatrixGraph(fei::SharedPtr<fei::MatrixGraph> matrixGraph)
00068 { matrixGraph_ = matrixGraph; }
00069 
00070 int
00071 Matrix_Local::getGlobalNumRows() const
00072 { return( sparseRowGraph_->rowNumbers.size() ); }
00073 
00074 int
00075 Matrix_Local::getLocalNumRows() const
00076 { return( getGlobalNumRows() ); }
00077 
00078 int
00079 Matrix_Local::getRowIndex(int rowNumber) const
00080 {
00081   int* rows = &(sparseRowGraph_->rowNumbers[0]);
00082   int numRows = getLocalNumRows();
00083   return( fei::binarySearch(rowNumber, rows, numRows) );
00084 }
00085 
00086 int
00087 Matrix_Local::getRowLength(int row, int& length) const
00088 {
00089   int idx = getRowIndex(row);
00090   if (idx < 0) return(idx);
00091 
00092   length = sparseRowGraph_->rowOffsets[idx+1] -
00093            sparseRowGraph_->rowOffsets[idx];
00094   return(0);
00095 }
00096 
00097 int
00098 Matrix_Local::putScalar(double scalar)
00099 {
00100   for(unsigned i=0; i<coefs_.size(); ++i) coefs_[i] = scalar;
00101   stateChanged_ = true;
00102   return(0);
00103 }
00104 
00105 int
00106 Matrix_Local::copyOutRow(int row, int len, double* coefs, int* indices) const
00107 {
00108   int idx = getRowIndex(row);
00109   if (idx < 0) return(idx);
00110 
00111   int offset = sparseRowGraph_->rowOffsets[idx];
00112   int length = sparseRowGraph_->rowOffsets[idx+1]-offset;
00113   if (length > len) length = len;
00114 
00115   for(int i=0; i<length; ++i) {
00116     indices[i] = sparseRowGraph_->packedColumnIndices[offset+i];
00117     coefs[i] = coefs_[offset+i];
00118   }
00119 
00120   return(0);
00121 }
00122 
00123 int
00124 Matrix_Local::giveToMatrix(int numRows, const int* rows,
00125                       int numCols, const int* cols,
00126                       const double* const* values,
00127                       bool sumInto,
00128                       int format)
00129 {
00130   if (numRows == 0 || numCols == 0) {
00131     return(0);
00132   }
00133 
00134   if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL) {
00135     return(-1);
00136   }
00137 
00138   const double** myvalues = const_cast<const double**>(values);
00139   if (format != FEI_DENSE_ROW) {
00140     fei::Matrix_core::copyTransposeToWorkArrays(numRows, numCols, values,
00141                               work_data1D_, work_data2D_);
00142     myvalues = &work_data2D_[0];
00143   }
00144 
00145   for(int i=0; i<numRows; ++i) {
00146     int idx = getRowIndex(rows[i]);
00147     if (idx < 0) {
00148       throw std::runtime_error("fei::Matrix_Local::sumIn ERROR, row not found.");
00149     }
00150 
00151     int offset = sparseRowGraph_->rowOffsets[idx];
00152     int len = sparseRowGraph_->rowOffsets[idx+1] - offset;
00153 
00154     int* colInds = &(sparseRowGraph_->packedColumnIndices[offset]);
00155     double* coefs   = &(coefs_[offset]);
00156 
00157     for(int j=0; j<numCols; ++j) {
00158       int idx2 = fei::binarySearch(cols[j], colInds, len);
00159       if (idx2 < 0) {
00160         throw std::runtime_error("fei::Matrix_Local::sumIn ERROR, col not found.");
00161       }
00162 
00163       if (sumInto) {
00164         coefs[idx2] += myvalues[i][j];
00165       }
00166       else {
00167         coefs[idx2] = myvalues[i][j];
00168       }
00169     }
00170   }
00171 
00172   stateChanged_ = true;
00173   return(0);
00174 }
00175 
00176 int
00177 Matrix_Local::sumIn(int numRows, const int* rows,
00178                     int numCols, const int* cols,
00179                     const double* const* values,
00180                     int format)
00181 {
00182   return( giveToMatrix(numRows, rows, numCols, cols, values,
00183                        true, format) );
00184 }
00185 
00186 int
00187 Matrix_Local::copyIn(int numRows, const int* rows,
00188                        int numCols, const int* cols,
00189                        const double* const* values,
00190                       int format)
00191 {
00192   return( giveToMatrix(numRows, rows, numCols, cols, values,
00193                        false, format) );
00194 }
00195 
00196 int
00197 Matrix_Local::sumInFieldData(int fieldID,
00198                                int idType,
00199                                int rowID,
00200                                int colID,
00201                                const double* const* data,
00202                                int format)
00203 {
00204   fei::SharedPtr<fei::VectorSpace> rspace = matrixGraph_->getRowSpace();
00205   fei::SharedPtr<fei::VectorSpace> cspace = matrixGraph_->getColSpace();
00206 
00207   int fieldSize = (int)rspace->getFieldSize(fieldID);
00208   std::vector<int> indices(2*fieldSize);
00209 
00210   rspace->getGlobalIndex(idType, rowID, fieldID, indices[0]);
00211   for(int i=1; i<fieldSize; ++i) {
00212     indices[i] = indices[0]+i;
00213   }
00214 
00215   cspace->getGlobalIndex(idType, colID, fieldID, indices[fieldSize]);
00216   for(int i=1; i<fieldSize; ++i) {
00217     indices[fieldSize+i] = indices[fieldSize]+i;
00218   }
00219 
00220   return( giveToMatrix(fieldSize, &indices[0], fieldSize, &indices[fieldSize],
00221                        data, true, format) );
00222 }
00223 
00224 int
00225 Matrix_Local::sumInFieldData(int fieldID,
00226                                int idType,
00227                                int rowID,
00228                                int colID,
00229                                const double* data,
00230                                int format)
00231 {
00232   fei::SharedPtr<fei::VectorSpace> rspace = matrixGraph_->getRowSpace();
00233 
00234   int fieldSize = (int)rspace->getFieldSize(fieldID);
00235   std::vector<const double*> data2D(fieldSize);
00236 
00237   int offset = 0;
00238   for(int i=0; i<fieldSize; ++i) {
00239     data2D[i] = &data[offset];
00240     offset += fieldSize;
00241   }
00242 
00243   return( sumInFieldData(fieldID, idType, rowID, colID,
00244                          &data2D[0], format) );
00245 }
00246 
00247 int
00248 Matrix_Local::sumIn(int blockID, int connectivityID,
00249                     const double* const* values,
00250                     int format)
00251 {
00252   int numIndices = matrixGraph_->getConnectivityNumIndices(blockID);
00253   std::vector<int> indices(numIndices);
00254 
00255   matrixGraph_->getConnectivityIndices(blockID, connectivityID,
00256                                        numIndices, &indices[0], numIndices);
00257 
00258   return( giveToMatrix(numIndices, &indices[0], numIndices, &indices[0],
00259                        values, true, format) );
00260 }
00261 
00262 int
00263 Matrix_Local::globalAssemble()
00264 { return(0); }
00265 
00266 int
00267 Matrix_Local::multiply(fei::Vector* x,
00268                        fei::Vector* y)
00269 {
00270   FEI_COUT << "fei::Matrix_Local::multiply NOT IMPLEMENTED."<<FEI_ENDL;
00271   return(-1);
00272 }
00273 
00274 int
00275 Matrix_Local::gatherFromOverlap(bool accumulate)
00276 {
00277   (void)accumulate;
00278   return(0);
00279 }
00280 
00281 int
00282 Matrix_Local::writeToFile(const char* filename,
00283                           bool matrixMarketFormat)
00284 {
00285   fei::SharedPtr<fei::VectorSpace> vspace = matrixGraph_->getRowSpace();
00286 
00287   MPI_Comm comm = vspace->getCommunicator();
00288 
00289   int localProc = fei::localProc(comm);
00290 
00291   FEI_OSTRINGSTREAM osstr;
00292   osstr << filename << "." << localProc << ".mtx";
00293   std::string fullname = osstr.str();
00294 
00295   FEI_OFSTREAM ofstr(fullname.c_str(), IOS_OUT);
00296 
00297   return( writeToStream(ofstr, matrixMarketFormat) );
00298 }
00299 
00300 int
00301 Matrix_Local::writeToStream(FEI_OSTREAM& ostrm,
00302                             bool matrixMarketFormat)
00303 {
00304   static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
00305 
00306   fei::SharedPtr<fei::VectorSpace> rspace = matrixGraph_->getRowSpace();
00307   fei::SharedPtr<fei::VectorSpace> cspace = matrixGraph_->getColSpace();
00308 
00309   int numRows = getLocalNumRows();
00310   int numCols = cspace->getEqnNumbers().size();
00311   int nnz = coefs_.size();
00312 
00313   if (matrixMarketFormat) {
00314     ostrm << mmbanner << FEI_ENDL;
00315     ostrm << numRows << " " << numCols << " " << nnz << FEI_ENDL;
00316   }
00317   else {
00318     ostrm << numRows << " " << numCols << " "<< FEI_ENDL;
00319   }
00320 
00321   std::vector<int>& rowNumbers = sparseRowGraph_->rowNumbers;
00322   std::vector<int>& rowOffsets = sparseRowGraph_->rowOffsets;
00323   std::vector<int>& colIndices = sparseRowGraph_->packedColumnIndices;
00324 
00325   ostrm.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
00326   ostrm.precision(13);
00327 
00328   int offset = 0;
00329   for(unsigned i=0; i<rowNumbers.size(); ++i) {
00330     int rowlen = rowOffsets[i+1]-rowOffsets[i];
00331 
00332     for(int j=0; j<rowlen; ++j) {
00333       if (matrixMarketFormat) {
00334         ostrm << rowNumbers[i]+1 << " " << colIndices[offset]+1
00335            << " " << coefs_[offset] << FEI_ENDL;
00336       }
00337       else {
00338         ostrm << rowNumbers[i] << " " << colIndices[offset]
00339            << " " << coefs_[offset] << FEI_ENDL;
00340       }
00341       ++offset;
00342     }
00343   }
00344 
00345   return(0);
00346 }
00347 
00348 bool
00349 Matrix_Local::usingBlockEntryStorage()
00350 { return( false ); }
00351 
00352 void
00353 Matrix_Local::markState()
00354 {
00355   stateChanged_ = false;
00356 }
00357 
00358 bool
00359 Matrix_Local::changedSinceMark()
00360 { return(stateChanged_); }
00361 
00362 const std::vector<int>&
00363 Matrix_Local::getRowNumbers() const
00364 { return( sparseRowGraph_->rowNumbers ); }
00365 
00366 const std::vector<int>&
00367 Matrix_Local::getRowOffsets() const
00368 { return( sparseRowGraph_->rowOffsets ); }
00369 
00370 const std::vector<int>&
00371 Matrix_Local::getColumnIndices() const
00372 { return( sparseRowGraph_->packedColumnIndices ); }
00373 
00374 const std::vector<double>&
00375 Matrix_Local::getCoefs() const
00376 { return( coefs_ ); }
00377 
00378 }//namespace fei
00379 

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