FEI Version of the Day
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() const
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 void
00275 Matrix_Local::setCommSizes()
00276 {
00277 }
00278 
00279 int
00280 Matrix_Local::gatherFromOverlap(bool accumulate)
00281 {
00282   (void)accumulate;
00283   return(0);
00284 }
00285 
00286 int
00287 Matrix_Local::writeToFile(const char* filename,
00288                           bool matrixMarketFormat)
00289 {
00290   fei::SharedPtr<fei::VectorSpace> vspace = matrixGraph_->getRowSpace();
00291 
00292   MPI_Comm comm = vspace->getCommunicator();
00293 
00294   int localProc = fei::localProc(comm);
00295 
00296   FEI_OSTRINGSTREAM osstr;
00297   osstr << filename << "." << localProc << ".mtx";
00298   std::string fullname = osstr.str();
00299 
00300   FEI_OFSTREAM ofstr(fullname.c_str(), IOS_OUT);
00301 
00302   return( writeToStream(ofstr, matrixMarketFormat) );
00303 }
00304 
00305 int
00306 Matrix_Local::writeToStream(FEI_OSTREAM& ostrm,
00307                             bool matrixMarketFormat)
00308 {
00309   static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
00310 
00311   fei::SharedPtr<fei::VectorSpace> rspace = matrixGraph_->getRowSpace();
00312   fei::SharedPtr<fei::VectorSpace> cspace = matrixGraph_->getColSpace();
00313 
00314   int numRows = getLocalNumRows();
00315   int numCols = cspace->getEqnNumbers().size();
00316   int nnz = coefs_.size();
00317 
00318   if (matrixMarketFormat) {
00319     ostrm << mmbanner << FEI_ENDL;
00320     ostrm << numRows << " " << numCols << " " << nnz << FEI_ENDL;
00321   }
00322   else {
00323     ostrm << numRows << " " << numCols << " "<< FEI_ENDL;
00324   }
00325 
00326   std::vector<int>& rowNumbers = sparseRowGraph_->rowNumbers;
00327   std::vector<int>& rowOffsets = sparseRowGraph_->rowOffsets;
00328   std::vector<int>& colIndices = sparseRowGraph_->packedColumnIndices;
00329 
00330   ostrm.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
00331   ostrm.precision(13);
00332 
00333   int offset = 0;
00334   for(unsigned i=0; i<rowNumbers.size(); ++i) {
00335     int rowlen = rowOffsets[i+1]-rowOffsets[i];
00336 
00337     for(int j=0; j<rowlen; ++j) {
00338       if (matrixMarketFormat) {
00339         ostrm << rowNumbers[i]+1 << " " << colIndices[offset]+1
00340            << " " << coefs_[offset] << FEI_ENDL;
00341       }
00342       else {
00343         ostrm << rowNumbers[i] << " " << colIndices[offset]
00344            << " " << coefs_[offset] << FEI_ENDL;
00345       }
00346       ++offset;
00347     }
00348   }
00349 
00350   return(0);
00351 }
00352 
00353 bool
00354 Matrix_Local::usingBlockEntryStorage()
00355 { return( false ); }
00356 
00357 void
00358 Matrix_Local::markState()
00359 {
00360   stateChanged_ = false;
00361 }
00362 
00363 bool
00364 Matrix_Local::changedSinceMark()
00365 { return(stateChanged_); }
00366 
00367 const std::vector<int>&
00368 Matrix_Local::getRowNumbers() const
00369 { return( sparseRowGraph_->rowNumbers ); }
00370 
00371 const std::vector<int>&
00372 Matrix_Local::getRowOffsets() const
00373 { return( sparseRowGraph_->rowOffsets ); }
00374 
00375 const std::vector<int>&
00376 Matrix_Local::getColumnIndices() const
00377 { return( sparseRowGraph_->packedColumnIndices ); }
00378 
00379 const std::vector<double>&
00380 Matrix_Local::getCoefs() const
00381 { return( coefs_ ); }
00382 
00383 }//namespace fei
00384 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends