FEI Version of the Day
fei_Matrix_Impl.hpp
00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 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 #ifndef _fei_Matrix_Impl_hpp_
00010 #define _fei_Matrix_Impl_hpp_
00011 
00012 #include <fei_fwd.hpp>
00013 #include <fei_mpi.h>
00014 #include <fei_defs.h>
00015 #include "fei_iostream.hpp"
00016 #include "fei_fstream.hpp"
00017 #include "fei_sstream.hpp"
00018 #include <fei_ArrayUtils.hpp>
00019 #include <fei_MatrixTraits.hpp>
00020 #include <fei_MatrixTraits_LinProbMgr.hpp>
00021 #include <fei_MatrixTraits_LinSysCore.hpp>
00022 #include <fei_MatrixTraits_FEData.hpp>
00023 #include <fei_MatrixTraits_FillableMat.hpp>
00024 #include <fei_FillableMat.hpp>
00025 
00026 #include <snl_fei_FEMatrixTraits.hpp>
00027 #include <snl_fei_FEMatrixTraits_FED.hpp>
00028 #include <snl_fei_BlockMatrixTraits.hpp>
00029 #include <fei_Matrix.hpp>
00030 #include <fei_MatrixGraph.hpp>
00031 #include <fei_Matrix_core.hpp>
00032 #include <snl_fei_Utils.hpp>
00033 
00034 #undef fei_file
00035 #define fei_file "fei_Matrix_Impl.hpp"
00036 #include <fei_ErrMacros.hpp>
00037 
00038 namespace fei {
00039 
00052   template<typename T>
00053   class Matrix_Impl : public fei::Matrix, public fei::Matrix_core {
00054   public:
00056     Matrix_Impl(fei::SharedPtr<T> matrix,
00057                fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00058                 int numLocalEqns,
00059                 bool zeroSharedRows=true);
00060 
00062     virtual ~Matrix_Impl();
00063 
00067     const char* typeName()
00068       {
00069         if (haveBlockMatrix()) {
00070           return(snl_fei::BlockMatrixTraits<T>::typeName());
00071         }
00072         else if (haveFEMatrix()) {
00073           return(snl_fei::FEMatrixTraits<T>::typeName());
00074         }
00075         else {
00076           return(fei::MatrixTraits<T>::typeName());
00077         }
00078       }
00079 
00082     int parameters(const fei::ParameterSet& paramset);
00083 
00087     fei::SharedPtr<T> getMatrix() { return( matrix_ ); }
00088     const fei::SharedPtr<T> getMatrix() const { return( matrix_ ); }
00089 
00090     fei::SharedPtr<fei::MatrixGraph> getMatrixGraph() const
00091       {return( Matrix_core::getMatrixGraph() ); }
00092 
00094     void setMatrixGraph(fei::SharedPtr<fei::MatrixGraph> matrixGraph);
00095 
00098     int getGlobalNumRows() const
00099       {
00100         if ((int)(globalOffsets().size()) < numProcs()+1) return(-1);
00101         return(globalOffsets()[numProcs()]);
00102       }
00103 
00106     int getLocalNumRows() const
00107       {
00108         return(lastLocalOffset() - firstLocalOffset() + 1);
00109       }
00110 
00112     int putScalar(double scalar);
00113 
00119     int getRowLength(int row, int& length) const;
00120 
00130     int copyOutRow(int row, int len, double* coefs, int* indices) const;
00131 
00143     int sumIn(int numRows, const int* rows,
00144               int numCols, const int* cols,
00145               const double* const* values,
00146               int format=0);
00147 
00159     int copyIn(int numRows, const int* rows,
00160                int numCols, const int* cols,
00161                const double* const* values,
00162                int format=0);
00163 
00179     int sumInFieldData(int fieldID,
00180                        int idType,
00181                        int rowID,
00182                        int colID,
00183                        const double* const* data,
00184                        int format=0);
00185 
00203     int sumInFieldData(int fieldID,
00204                        int idType,
00205                        int rowID,
00206                        int colID,
00207                        const double* data,
00208                        int format=0);
00209 
00219     int sumIn(int blockID, int connectivityID,
00220               const double* const* values,
00221               int format=0);
00222 
00227     int globalAssemble();
00228 
00231     int multiply(fei::Vector* x,
00232                  fei::Vector* y);
00233 
00234     void setCommSizes();
00235 
00241     int gatherFromOverlap(bool accumulate = true);
00242 
00244     int writeToFile(const char* filename,
00245                     bool matrixMarketFormat=true);
00246 
00249     int writeToStream(FEI_OSTREAM& ostrm,
00250                       bool matrixMarketFormat=true);
00251 
00252     bool usingBlockEntryStorage()
00253       {
00254         return(haveBlockMatrix());
00255       }
00256 
00258     int giveToUnderlyingMatrix(int numRows, const int* rows,
00259                                int numCols, const int* cols,
00260                                const double* const* values,
00261                                bool sumInto,
00262                                int format);
00263 
00265     int giveToUnderlyingBlockMatrix(int row,
00266                                     int rowDim,
00267                                     int numCols,
00268                                     const int* cols,
00269                                     const int* LDAs,
00270                                     const int* colDims,
00271                                     const double* const* values,
00272                                     bool sumInto);
00273 
00274     void markState();
00275 
00276     bool changedSinceMark();
00277 
00278     double* getBeginPointer()
00279       {
00280         return fei::MatrixTraits<T>::getBeginPointer(matrix_.get());
00281       }
00282 
00283     int getOffset(int row, int col)
00284       {
00285         fei::SharedPtr<fei::MatrixGraph> mgraph = getMatrixGraph();
00286         fei::SharedPtr<fei::VectorSpace> rowspace = mgraph->getRowSpace();
00287         fei::SharedPtr<fei::VectorSpace> colspace = mgraph->getColSpace();
00288         int row_index, col_index;
00289         int nodeType = 0;//fix this!!! hard-coded 0
00290         rowspace->getGlobalIndex(nodeType, row, row_index);
00291         colspace->getGlobalIndex(nodeType, col, col_index);
00292         
00293         return fei::MatrixTraits<T>::getOffset(matrix_.get(),row_index,col_index);
00294       }
00295 
00296   private:
00297     int giveToMatrix(int numRows, const int* rows,
00298                      int numCols, const int* cols,
00299                      const double* const* values,
00300                      bool sumInto,
00301                      int format);
00302  
00303     int giveToBlockMatrix(int numRows, const int* rows,
00304                           int numCols, const int* cols,
00305                           const double* const* values,
00306                           bool sumInto);
00307 
00308     fei::SharedPtr<T> matrix_;
00309     bool globalAssembleCalled_;
00310     bool changedSinceMark_;
00311     std::string dbgprefix_;
00312   };//class Matrix_Impl
00313 }//namespace fei
00314 
00315 //----------------------------------------------------------------------------
00316 template<typename T>
00317 inline int fei::Matrix_Impl<T>::globalAssemble()
00318 {
00319   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00320     FEI_OSTREAM& os = *output_stream_;
00321     os << dbgprefix_<<"globalAssemble"<<FEI_ENDL;
00322   }
00323 
00324   globalAssembleCalled_ = true;
00325 
00326   if (haveBlockMatrix()) {
00327     return( snl_fei::BlockMatrixTraits<T>::globalAssemble(matrix_.get()) );
00328   }
00329   return( fei::MatrixTraits<T>::globalAssemble(matrix_.get()) );
00330 }
00331 
00332 //----------------------------------------------------------------------------
00333 template<typename T>
00334 inline void fei::Matrix_Impl<T>::setMatrixGraph(fei::SharedPtr<fei::MatrixGraph> matrixGraph)
00335 {
00336   Matrix_core::setMatrixGraph(matrixGraph);
00337 }
00338 
00339 //----------------------------------------------------------------------------
00340 template<typename T>
00341 inline void fei::Matrix_Impl<T>::markState()
00342 {
00343 if(changedSinceMark_)
00344   changedSinceMark_ = false;
00345 }
00346 
00347 //----------------------------------------------------------------------------
00348 template<typename T>
00349 inline bool fei::Matrix_Impl<T>::changedSinceMark()
00350 {
00351   return(changedSinceMark_);
00352 }
00353 
00354 //----------------------------------------------------------------------------
00355 template<typename T>
00356 int fei::Matrix_Impl<T>::giveToUnderlyingMatrix(int numRows, const int* rows,
00357                                                int numCols, const int* cols,
00358                                                const double* const* values,
00359                                                bool sumInto,
00360                                                int format)
00361 {
00362   if (format != FEI_DENSE_ROW) {
00363     ERReturn(-1);
00364   }
00365 
00366   if (output_level_ > fei::BRIEF_LOGS && output_stream_ != 0) {
00367     FEI_OSTREAM& os = *output_stream_;
00368     for(int i=0; i<numRows; ++i) {
00369       os << dbgprefix_<<"giveToUnderlyingMatrix ";
00370       for(int j=0; j<numCols; ++j) {
00371         os << "("<<rows[i]<<","<<cols[j]<<","<<values[i][j]<<") ";
00372       }
00373       os << FEI_ENDL;
00374     }
00375   }
00376 
00377   int err = fei::MatrixTraits<T>::putValuesIn(matrix_.get(), numRows, rows,
00378                                         numCols, cols, values, sumInto);
00379   if (err != 0) {
00380     return(err);
00381   }
00382 
00383   changedSinceMark_ = true;
00384   return(0);
00385 }
00386 
00387 //----------------------------------------------------------------------------
00388 template<typename T>
00389 int fei::Matrix_Impl<T>::giveToUnderlyingBlockMatrix(int row,
00390                                                     int rowDim,
00391                                                     int numCols,
00392                                                     const int* cols,
00393                                                     const int* LDAs,
00394                                                     const int* colDims,
00395                                                     const double* const* values,
00396                                                     bool sumInto)
00397 {
00398   if (sumInto) {
00399     if ( snl_fei::BlockMatrixTraits<T>::sumIn(matrix_.get(),
00400                                          row, rowDim, numCols, cols,
00401                                          LDAs, colDims, values) != 0) {
00402       ERReturn(-1);
00403     }
00404   }
00405   else {
00406     if ( snl_fei::BlockMatrixTraits<T>::copyIn(matrix_.get(),
00407                                           row, rowDim, numCols, cols,
00408                                           LDAs, colDims, values) != 0) {
00409       ERReturn(-1);
00410     }
00411   }
00412 
00413   changedSinceMark_ = true;
00414 
00415   return(0);
00416 }
00417 
00418 #include <fei_macros.hpp>
00419 #include <fei_chk_mpi.hpp>
00420 
00421 #include <fei_TemplateUtils.hpp>
00422 
00423 #include <fei_Pattern.hpp>
00424 #include <fei_VectorSpace.hpp>
00425 
00426 #include <fei_ConnectivityBlock.hpp>
00427 
00428 #include <snl_fei_PointBlockMap.hpp>
00429 
00430 #include <fei_Record.hpp>
00431 
00432 #include <snl_fei_Utils.hpp>
00433 
00434 //----------------------------------------------------------------------------
00435 template<typename T>
00436 fei::Matrix_Impl<T>::Matrix_Impl(fei::SharedPtr<T> matrix,
00437                            fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00438                                 int numLocalEqns,
00439                                 bool zeroSharedRows)
00440   : Matrix_core(matrixGraph, numLocalEqns),
00441     matrix_(matrix),
00442     globalAssembleCalled_(false),
00443     changedSinceMark_(true),
00444     dbgprefix_("MatImpl: ")
00445 {
00446   if (strcmp(snl_fei::FEMatrixTraits<T>::typeName(), "unsupported")) {
00447     setFEMatrix(true);
00448   }
00449   else {
00450     setFEMatrix(false);
00451   }
00452 
00453   if (strcmp(snl_fei::BlockMatrixTraits<T>::typeName(), "unsupported")) {
00454     setBlockMatrix(true);
00455   }
00456   else {
00457     setBlockMatrix(false);
00458   }
00459 
00460   if (zeroSharedRows && matrixGraph->getGlobalNumSlaveConstraints() == 0) {
00461     std::vector<double> zeros;
00462     fei::SharedPtr<fei::SparseRowGraph> srg = matrixGraph->getRemotelyOwnedGraphRows();
00463     if (srg.get() != NULL) {
00464       for(size_t row=0; row<srg->rowNumbers.size(); ++row) {
00465         int rowLength = srg->rowOffsets[row+1] - srg->rowOffsets[row];
00466         if (rowLength == 0) continue;
00467         zeros.resize(rowLength, 0.0);
00468         const double* zerosPtr = &zeros[0];
00469         const int* cols = &srg->packedColumnIndices[srg->rowOffsets[row]];
00470         sumIn(1, &srg->rowNumbers[row], rowLength, cols, &zerosPtr);
00471       }
00472       setCommSizes();
00473       std::map<int,FillableMat*>& remote = getRemotelyOwnedMatrices();
00474       for(std::map<int,FillableMat*>::iterator iter=remote.begin(), end=remote.end(); iter!=end; ++iter)
00475       {
00476         iter->second->clear();
00477       }
00478     }
00479   }
00480 }
00481 
00482 //----------------------------------------------------------------------------
00483 template<typename T>
00484 fei::Matrix_Impl<T>::~Matrix_Impl()
00485 {
00486 }
00487 
00488 //----------------------------------------------------------------------------
00489 template<typename T>
00490 int fei::Matrix_Impl<T>::parameters(const fei::ParameterSet& paramset)
00491 {
00492   Matrix_core::parameters(paramset);
00493   return 0;
00494 }
00495 
00496 //----------------------------------------------------------------------------
00497 template<typename T>
00498 int fei::Matrix_Impl<T>::getRowLength(int row, int& length) const
00499 {
00500   if (haveBlockMatrix()) {
00501     return( snl_fei::BlockMatrixTraits<T>::getPointRowLength(matrix_.get(), row, length) );
00502   }
00503   else {
00504     int code = fei::MatrixTraits<T>::getRowLength(matrix_.get(), row, length);
00505     if (code < 0) {
00506       //matrix row probably not local. maybe it's shared.
00507       int proc = getOwnerProc(row);
00508       const FillableMat* remote_mat = getRemotelyOwnedMatrix(proc);
00509       if (remote_mat->hasRow(row)) {
00510         const CSVec* row_entries = remote_mat->getRow(row);
00511         length = row_entries->size();
00512       }
00513       else {
00514         length = 0;
00515       }
00516     }
00517   }
00518   return 0;
00519 }
00520 
00521 //----------------------------------------------------------------------------
00522 template<typename T>
00523 int fei::Matrix_Impl<T>::putScalar(double scalar)
00524 {
00525   if (output_level_ > fei::BRIEF_LOGS && output_stream_ != 0) {
00526     FEI_OSTREAM& os = *output_stream_;
00527     os << dbgprefix_<<"putScalar("<<scalar<<")"<<FEI_ENDL;
00528   }
00529 
00530   if (haveFEMatrix()) {
00531     if (scalar != 0.0) return(-1);
00532     CHK_ERR( snl_fei::FEMatrixTraits<T>::reset(matrix_.get()) );
00533   }
00534   else if (haveBlockMatrix()) {
00535     if (globalAssembleCalled_ == true) {
00536       CHK_ERR( snl_fei::BlockMatrixTraits<T>::putScalar(matrix_.get(), scalar) );
00537     }
00538   }
00539   else {
00540     CHK_ERR( fei::MatrixTraits<T>::setValues(matrix_.get(), scalar) );
00541   }
00542 
00543   putScalar_remotelyOwned(scalar);
00544 
00545   changedSinceMark_ = true;
00546 
00547   return(0);
00548 }
00549 
00550 //----------------------------------------------------------------------------
00551 template<typename T>
00552 int fei::Matrix_Impl<T>::copyOutRow(int row, int len,
00553                                    double* coefs, int* indices) const
00554 {
00555   if (len < 1) {
00556     return 0;
00557   }
00558 
00559   if (haveBlockMatrix()) {
00560     int dummy;
00561     return( snl_fei::BlockMatrixTraits<T>::copyOutPointRow(matrix_.get(), firstLocalOffset(),
00562                                                   row, len,
00563                                                   coefs, indices, dummy));
00564   }
00565   else {
00566     int code = fei::MatrixTraits<T>::copyOutRow(matrix_.get(), row, len,
00567                                         coefs, indices);
00568     if (code < 0) {
00569       int proc = getOwnerProc(row);
00570       const FillableMat* remote_mat = getRemotelyOwnedMatrix(proc);
00571       if (remote_mat->hasRow(row)) {
00572         const CSVec* row_entries = remote_mat->getRow(row);
00573         const std::vector<int>& row_indices = row_entries->indices();
00574         const std::vector<double>& row_coefs = row_entries->coefs();
00575         for(size_t i=0; i<row_indices.size(); ++i) {
00576           indices[i] = row_indices[i];
00577           coefs[i] = row_coefs[i];
00578         }
00579       }
00580     }
00581   }
00582   return 0;
00583 }
00584 
00585 //----------------------------------------------------------------------------
00586 template<typename T>
00587 int fei::Matrix_Impl<T>::sumIn(int numRows, const int* rows,
00588                               int numCols, const int* cols,
00589                               const double* const* values,
00590                               int format)
00591 {
00592   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00593     FEI_OSTREAM& os = *output_stream_;
00594     os << dbgprefix_<<"sumIn"<<FEI_ENDL;
00595     if (output_level_ >= fei::FULL_LOGS) {
00596       for(int i=0; i<numRows; ++i) {
00597         for(int j=0; j<numCols; ++j) {
00598           os << "("<<rows[i]<<","<<cols[j]<<","<<values[i][j]<<") ";
00599         }
00600         os << FEI_ENDL;
00601       }
00602     }
00603   }
00604 
00605   return( giveToMatrix( numRows, rows, numCols, cols, values, true, format) );
00606 }
00607 
00608 //----------------------------------------------------------------------------
00609 template<typename T>
00610 int fei::Matrix_Impl<T>::copyIn(int numRows, const int* rows,
00611                                int numCols, const int* cols,
00612                                const double* const* values,
00613                                int format)
00614 {
00615   if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
00616     FEI_OSTREAM& os = *output_stream_;
00617     os << "copyIn"<<FEI_ENDL;
00618     if (output_level_ >= fei::FULL_LOGS) {
00619       for(int i=0; i<numRows; ++i) {
00620         for(int j=0; j<numCols; ++j) {
00621           os << "("<<rows[i]<<","<<cols[j]<<","<<values[i][j]<<") ";
00622         }
00623         os << FEI_ENDL;
00624       }
00625     }
00626   }
00627   return( giveToMatrix( numRows, rows, numCols, cols, values, false, format) );
00628 }
00629 
00630 //----------------------------------------------------------------------------
00631 template<typename T>
00632 int fei::Matrix_Impl<T>::sumInFieldData(int fieldID,
00633                                        int idType,
00634                                        int rowID,
00635                                        int colID,
00636                                        const double* const* data,
00637                                        int format)
00638 {
00639   fei::SharedPtr<fei::VectorSpace> vspace = vecSpace();
00640 
00641   int fieldSize = vspace->getFieldSize(fieldID);
00642 
00643   if (fieldSize <= 0) ERReturn(-1);
00644 
00645   work_indices_.resize(fieldSize*2);
00646   int* indicesPtr = &work_indices_[0];
00647   int i;
00648 
00649   CHK_ERR( vspace->getGlobalIndices(1, &rowID, idType, fieldID, indicesPtr));
00650   for(i=1; i<fieldSize; ++i) {
00651     indicesPtr[i] = indicesPtr[i-1]+1;
00652   }
00653 
00654   CHK_ERR( vspace->getGlobalIndices(1, &colID, idType, fieldID,
00655                                         &(indicesPtr[fieldSize])) );
00656   for(i=fieldSize+1; i<2*fieldSize; ++i) {
00657     indicesPtr[i] = indicesPtr[i-1]+1;
00658   }
00659 
00660   CHK_ERR( sumIn(fieldSize, indicesPtr, fieldSize, &(indicesPtr[fieldSize]),
00661                  data, format) );
00662 
00663   return(0);
00664 }
00665 
00666 //----------------------------------------------------------------------------
00667 template<typename T>
00668 int fei::Matrix_Impl<T>::sumInFieldData(int fieldID,
00669                                        int idType,
00670                                        int rowID,
00671                                        int colID,
00672                                        const double* data,
00673                                        int format)
00674 {
00675   fei::SharedPtr<fei::VectorSpace> vspace = vecSpace();
00676 
00677   int fieldSize = vspace->getFieldSize(fieldID);
00678 
00679   if (fieldSize <= 0) ERReturn(-1);
00680 
00681   work_data2D_.resize(fieldSize);
00682 
00683   const double** data2DPtr = &work_data2D_[0];
00684   for(int i=0; i<fieldSize; ++i) {
00685     data2DPtr[i] = &(data[i*fieldSize]);
00686   }
00687 
00688   CHK_ERR( sumInFieldData(fieldID, idType, rowID, colID, data2DPtr, format) );
00689 
00690   return(0);
00691 }
00692 
00693 //----------------------------------------------------------------------------
00694 template<typename T>
00695 int fei::Matrix_Impl<T>::sumIn(int blockID, int connectivityID,
00696                               const double* const* values, int format)
00697 {
00698   if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
00699     FEI_OSTREAM& os = *output_stream_;
00700     os << dbgprefix_<<"sumIn blkID=" << blockID
00701        << ", connID=" << connectivityID << FEI_ENDL;
00702   }
00703 
00704   fei::SharedPtr<fei::MatrixGraph> mgraph = getMatrixGraph();
00705 
00706   if (mgraph.get() == NULL) ERReturn(-1);
00707 
00708   const fei::ConnectivityBlock* cblock = mgraph->getConnectivityBlock(blockID);
00709   if (cblock==NULL) {
00710     FEI_OSTRINGSTREAM osstr;
00711     osstr << "fei::Matrix_Impl::sumIn ERROR, unable to "
00712           << "look up connectivity-block with ID "<<blockID;
00713     throw std::runtime_error(osstr.str());
00714   }
00715 
00716   bool symmetric = cblock->isSymmetric();
00717   const fei::Pattern* pattern = cblock->getRowPattern();
00718   const fei::Pattern* colpattern = symmetric ? NULL : cblock->getColPattern();
00719   const int* rowConn = cblock->getRowConnectivity(connectivityID);
00720 
00721   fei::SharedPtr<fei::VectorSpace> rspace = mgraph->getRowSpace();
00722   rspace->getGlobalIndicesL(pattern, rowConn, work_indices2_);
00723 
00724   int numRowIndices = work_indices2_.size();
00725   int* rowIndices = &work_indices2_[0];
00726 
00727   if (haveFEMatrix() || haveBlockMatrix()) {
00728     FieldDofMap<int>& fdofmap = rspace->getFieldDofMap();
00729     const std::map<int,int>& connIDs = cblock->getConnectivityIDs();
00730     std::map<int,int>::const_iterator
00731       iter = connIDs.find(connectivityID);
00732     if (iter == connIDs.end()) ERReturn(-1);
00733     int connOffset = iter->second;
00734     int numIDs = pattern->getNumIDs();
00735     const int* indicesPerID = pattern->getNumIndicesPerID();
00736     const int* fieldsPerID = pattern->getNumFieldsPerID();
00737     const int* fieldIDs = pattern->getFieldIDs();
00738 
00739     int numDofs = 0;
00740     for(int ii=0; ii<numIDs; ++ii) {
00741       numDofs += indicesPerID[ii];
00742     }
00743 
00744     const int* numIndicesPerID = pattern->getNumIndicesPerID();
00745     work_indices_.resize(numIDs+numDofs);
00746     int i, *nodeNumbers = &work_indices_[0];
00747     int* dof_ids = nodeNumbers+numIDs;
00748 
00749     int nodeType = 0;
00750     snl_fei::RecordCollection* records = NULL;
00751     rspace->getRecordCollection(nodeType, records);
00752     int foffset = 0;
00753     int doffset = 0;
00754     for(i=0; i<numIDs; ++i) {
00755       fei::Record<int>* rec = records->getRecordWithLocalID(rowConn[i]);
00756       nodeNumbers[i] = rec->getNumber();
00757       for(int ii=0; ii<fieldsPerID[i]; ++ii) {
00758         int fieldSize = rspace->getFieldSize(fieldIDs[foffset]);
00759         int dof_id = fdofmap.get_dof_id(fieldIDs[foffset++], 0);
00760         for(int j=0; j<fieldSize; ++j) {
00761           dof_ids[doffset++] = dof_id + j;
00762         }
00763       }
00764     }
00765 
00766     if (haveFEMatrix()) {
00767       CHK_ERR( snl_fei::FEMatrixTraits<T>::sumInElemMatrix(matrix_.get(), blockID, connOffset,
00768                                                   numIDs, nodeNumbers,
00769                                                   numIndicesPerID, dof_ids, values) );
00770       changedSinceMark_ = true;
00771     }
00772 
00773     if (haveBlockMatrix()) {
00774       if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL &&
00775           format != FEI_BLOCK_DIAGONAL_ROW) {
00776         fei::console_out() << "fei::Matrix_Impl::sumIn ERROR, for block-matrix, format must"
00777                  << " be FEI_DENSE_ROW or FEI_DENSE_COL."<<FEI_ENDL;
00778         ERReturn(-1);
00779       }
00780 
00781       int numPtIndices = pattern->getNumIndices();
00782       int* ptIndices = &work_indices2_[0];
00783 
00784       int numPtColIndices = symmetric ? numPtIndices : colpattern->getNumIndices();
00785 
00786       int len = numPtIndices*numPtColIndices;
00787 
00788       if (format == FEI_BLOCK_DIAGONAL_ROW) {
00789         len = 0;
00790         for(i=0; i<numIDs; ++i) {
00791           len += numIndicesPerID[i]*numIndicesPerID[i];
00792         }
00793       }
00794 
00795       double* ccvalues = new double[len];
00796       //Ouch! Data copy! Remember to optimize this later...
00797       if (format == FEI_BLOCK_DIAGONAL_ROW) {
00798         snl_fei::copy2DBlockDiagToColumnContig(numIDs, numIndicesPerID,
00799                                                values, format, ccvalues);
00800       }
00801       else {
00802         snl_fei::copy2DToColumnContig(numPtIndices, numPtColIndices, values, format,
00803                                       ccvalues);
00804       }
00805 
00806       int ptRowOffset = 0;
00807       for(i=0; i<numIDs; ++i) {
00808 
00809         fei::Record<int>* record = records->getRecordWithLocalID(rowConn[i]);
00810         if (record->getOwnerProc() == localProc()) {
00811 
00812           int numColIDs = numIDs;
00813           int* colNodeNums = nodeNumbers;
00814           const int* colDims = numIndicesPerID;
00815           int LDA = numPtColIndices;
00816           if (format == FEI_BLOCK_DIAGONAL_ROW) {
00817             numColIDs = 1;
00818             colNodeNums = &(nodeNumbers[i]);
00819             colDims = &(numIndicesPerID[i]);
00820             LDA = numIndicesPerID[i];
00821           }
00822 
00823           CHK_ERR( snl_fei::BlockMatrixTraits<T>::sumIn(matrix_.get(),
00824                                                nodeNumbers[i],
00825                                                numIndicesPerID[i],
00826                                                numColIDs, colNodeNums,
00827                                                colDims, LDA,
00828                                                &(ccvalues[ptRowOffset])) );
00829           changedSinceMark_ = true;
00830         }
00831         else {
00832           if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
00833             FEI_OSTREAM& os = *output_stream_;
00834             for(int ii=0; ii<numIndicesPerID[i]; ++ii) {
00835               os << "#   remote pt-row " << ptIndices[ptRowOffset]+ii <<" ";
00836               for(int jj=0; jj<numPtIndices; ++jj) {
00837                 os << "("<<ptIndices[jj]<<","<<values[ptRowOffset+ii][jj]<<") ";
00838               }
00839               os << FEI_ENDL;
00840             }
00841           }
00842 
00843           for(int ii=0; ii<numIndicesPerID[i]; ++ii) {
00844             int p=eqnComm_->getOwnerProc(ptIndices[ptRowOffset]+ii);
00845             FillableMat* remote_mat = getRemotelyOwnedMatrix(p);
00846             remote_mat->sumInRow(ptIndices[ptRowOffset]+ii,
00847                                       ptIndices,
00848                                       values[ptRowOffset+ii],
00849                                       numPtIndices);
00850             changedSinceMark_ = true;
00851           }
00852         }
00853 
00854         ptRowOffset += numIndicesPerID[i];
00855       }
00856 
00857       delete [] ccvalues;
00858     }
00859 
00860     return(0);
00861   }
00862 
00863   int numColIndices = symmetric ? numRowIndices : colpattern->getNumIndices();
00864   int* colIndices = rowIndices;
00865   const int* colConn = NULL;
00866 
00867   if (!symmetric) {
00868     colConn = cblock->getColConnectivity(connectivityID);
00869     mgraph->getColSpace()->getGlobalIndicesL(colpattern,
00870                                             colConn, work_indices_);
00871     colIndices = &work_indices_[0];
00872   }
00873 
00874   if (symmetric) {
00875     if (format == FEI_DENSE_ROW || format == FEI_DENSE_COL) {
00876       CHK_ERR( sumIn(numRowIndices, rowIndices, numRowIndices, rowIndices,
00877                      values, format) );
00878     }
00879     else if (format == FEI_BLOCK_DIAGONAL_ROW) {
00880       int totalnumfields = pattern->getTotalNumFields();
00881       const int* fieldIDs = pattern->getFieldIDs();
00882       fei::SharedPtr<fei::VectorSpace> rowspace = mgraph->getRowSpace();
00883       fei::VectorSpace* rowspaceptr = rowspace.get();
00884       int ioffset = 0;
00885       for(int i=0; i<totalnumfields; ++i) {
00886         int fieldsize = rowspaceptr->getFieldSize(fieldIDs[i]);
00887         if (ioffset+fieldsize > numRowIndices) {
00888           FEI_OSTRINGSTREAM osstr;
00889           osstr<<"snl_fei::sumIn, format=FEI_BLOCK_DIAGONAL_ROW, block-sizes"
00890                << " not consistent with total num-indices."<<FEI_ENDL;
00891           throw std::runtime_error(osstr.str());
00892         }
00893 
00894         CHK_ERR( sumIn(fieldsize, &(rowIndices[ioffset]),
00895                        fieldsize, &(rowIndices[ioffset]),
00896                        &(values[ioffset]), FEI_DENSE_ROW) );
00897         ioffset += fieldsize;
00898         changedSinceMark_ = true;
00899       }
00900     }
00901     else {
00902       FEI_OSTRINGSTREAM osstr;
00903       osstr << "fei::Matrix_Impl::sumIn, format="<<format<<" not supported."
00904             << FEI_ENDL;
00905       throw std::runtime_error(osstr.str());
00906     }
00907   }
00908   else {
00909     if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL) {
00910       FEI_OSTRINGSTREAM osstr;
00911       osstr << "fei::Matrix_Impl::sumIn, format="<<format<<" not valid with"
00912             << " un-symmetric matrix contributions."<<FEI_ENDL;
00913       throw std::runtime_error(osstr.str());
00914     }
00915 
00916     CHK_ERR( sumIn(numRowIndices, rowIndices, numColIndices, colIndices,
00917                    values, format) );
00918     changedSinceMark_ = true;
00919   }
00920 
00921   return(0);
00922 }
00923 
00924 //----------------------------------------------------------------------------
00925 template<typename T>
00926 int fei::Matrix_Impl<T>::multiply(fei::Vector* x,
00927                                  fei::Vector* y)
00928 {
00929   return( fei::MatrixTraits<T>::matvec(matrix_.get(), x, y) );
00930 }
00931 
00932 //----------------------------------------------------------------------------
00933 template<typename T>
00934 void fei::Matrix_Impl<T>::setCommSizes()
00935 {
00936   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00937     (*output_stream_) << dbgprefix_<<"setCommSizes"<<FEI_ENDL;
00938   }
00939 
00940   Matrix_core::setCommSizes();
00941 }
00942 
00943 //----------------------------------------------------------------------------
00944 template<typename T>
00945 int fei::Matrix_Impl<T>::gatherFromOverlap(bool accumulate)
00946 {
00947   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00948     (*output_stream_) << dbgprefix_<<"gatherFromOverlap"<<FEI_ENDL;
00949   }
00950 
00951   CHK_ERR( Matrix_core::gatherFromOverlap(accumulate) );
00952 
00953   return(0);
00954 }
00955 
00956 //----------------------------------------------------------------------------
00957 template<typename T>
00958 int fei::Matrix_Impl<T>::giveToMatrix(int numRows, const int* rows,
00959                                      int numCols, const int* cols,
00960                                      const double* const* values,
00961                                      bool sumInto,
00962                                      int format)
00963 {
00964   if (numRows == 0 || numCols == 0) {
00965     return(0);
00966   }
00967 
00968   if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL) {
00969     return(-1);
00970   }
00971 
00972   const double** myvalues = const_cast<const double**>(values);
00973   if (format != FEI_DENSE_ROW) {
00974     copyTransposeToWorkArrays(numRows, numCols, values,
00975                               work_data1D_, work_data2D_);
00976     myvalues = &work_data2D_[0];
00977   }
00978 
00979   if ((int)work_ints_.size() < numRows) {
00980     work_ints_.resize(numRows);
00981   }
00982 
00983   if (haveBlockMatrix()) {
00984     return( giveToBlockMatrix(numRows, rows, numCols, cols,
00985                               myvalues, sumInto) );
00986   }
00987 
00988   int i; 
00989   int numRemote = 0;
00990   int* workIntPtr = &work_ints_[0];
00991   for(i=0; i<numRows; ++i) {
00992     int row = rows[i];
00993     if (row < firstLocalOffset() || row > lastLocalOffset()) {
00994       ++numRemote;
00995       workIntPtr[i] = 1;
00996     }
00997     else {
00998       workIntPtr[i] = 0;
00999     }
01000   }
01001 
01002   if (numRemote < 1) {
01003     int err = giveToUnderlyingMatrix(numRows, rows, numCols, cols, myvalues,
01004                                     sumInto, FEI_DENSE_ROW);
01005     changedSinceMark_ = true;
01006     if (err != 0) {
01007       FEI_OSTRINGSTREAM osstr;
01008       osstr << "fei::Matrix_Impl::giveToMatrix ERROR: err="<<err
01009         << " returned from giveToUnderlyingMatrix.";
01010       throw std::runtime_error(osstr.str());
01011     }
01012     return(0);
01013   }
01014 
01015   for(i=0; i<numRows; ++i) {
01016     int row = rows[i];
01017     const double*const rowvalues = myvalues[i];
01018 
01019     if (workIntPtr[i] > 0) {
01020       int proc = eqnComm_->getOwnerProc(row);
01021       FillableMat* remote_mat = getRemotelyOwnedMatrix(proc);
01022 
01023       if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
01024         FEI_OSTREAM& os = *output_stream_;
01025         os << dbgprefix_<<" remote_mat["<<proc<<"]: ";
01026         for(int jj=0; jj<numCols; ++jj) {
01027           os << "("<<row<<","<<cols[jj]<<","<<rowvalues[jj]<<") ";
01028         }
01029         os << FEI_ENDL;
01030       }
01031 
01032       if (sumInto) {
01033         remote_mat->sumInRow(row, cols, rowvalues, numCols);
01034       }
01035       else {
01036         remote_mat->putRow(row, cols, rowvalues, numCols);
01037       }
01038     }
01039     else {
01040       CHK_ERR( giveToUnderlyingMatrix(1, &row, numCols, cols, &rowvalues,
01041                                       sumInto, 0) );
01042     }
01043   }
01044 
01045   changedSinceMark_ = true;
01046 
01047   return(0);
01048 }
01049 
01050 //----------------------------------------------------------------------------
01051 template<typename T>
01052 int fei::Matrix_Impl<T>::giveToBlockMatrix(int numRows, const int* rows,
01053                                           int numCols, const int* cols,
01054                                           const double* const* values,
01055                                           bool sumInto)
01056 {
01057   if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
01058     FEI_OSTREAM& os = *output_stream_;
01059     os << "#  giveToBlockMatrix numRows: " << numRows
01060        << ", numCols: " << numCols << FEI_ENDL;
01061   }
01062 
01063   if (numRows == 0 || numCols == 0) {
01064     return(0);
01065   }
01066 
01067   snl_fei::PointBlockMap* pointBlockMap = vecSpace()->getPointBlockMap();
01068 
01069   if (sumInto && numProcs() == 1) {
01070     std::vector<int> temp((numRows+numCols)*2);
01071     int* tempPtr = &temp[0];
01072     int* blkRows = tempPtr;
01073     int* blkRowOffsets = tempPtr+numRows;
01074     int* blkCols = blkRowOffsets+numRows;
01075     int* blkColOffsets = blkCols+numCols;
01076 
01077     CHK_ERR( convertPtToBlk(numRows, rows, numCols, cols,
01078                             blkRows, blkRowOffsets,
01079                             blkCols, blkColOffsets) );
01080 
01081     std::vector<int> blockRows, blockRowSizes;
01082     std::vector<int> blockCols, blockColSizes;
01083     for(int i=0; i<numRows; ++i) fei::sortedListInsert(blkRows[i], blockRows);
01084     for(int i=0; i<numCols; ++i) fei::sortedListInsert(blkCols[i], blockCols);
01085 
01086     int rowSizeTotal = 0, colSizeTotal = 0;
01087 
01088     for(size_t i=0; i<blockRows.size(); ++i) {
01089       int size = pointBlockMap->getBlkEqnSize(blockRows[i]);
01090       blockRowSizes.push_back(size);
01091       rowSizeTotal += size;
01092     }
01093     for(size_t i=0; i<blockCols.size(); ++i) {
01094       int size = pointBlockMap->getBlkEqnSize(blockCols[i]);
01095       blockColSizes.push_back(size);
01096       colSizeTotal += size;
01097     }
01098     std::vector<double> coefs_1d(rowSizeTotal*colSizeTotal, 0.0);
01099     double* coefs1dPtr = &coefs_1d[0];
01100     std::vector<double*> coefs_2d(blockRows.size()*blockCols.size());
01101     double** coefs2dPtr = &coefs_2d[0];
01102 
01103     int blkCounter = 0;
01104     int offset = 0;
01105     for(size_t i=0; i<blockRows.size(); ++i) {
01106       for(size_t j=0; j<blockCols.size(); ++j) {
01107         coefs2dPtr[blkCounter++] = &(coefs1dPtr[offset]);
01108         offset += blockRowSizes[i]*blockColSizes[j];
01109       }
01110     }
01111 
01112     for(int i=0; i<numRows; ++i) {
01113       int rowind = fei::binarySearch(blkRows[i], blockRows);
01114       int rowsize = blockRowSizes[rowind];
01115 
01116       for(int j=0; j<numCols; ++j) {
01117         int colind = fei::binarySearch(blkCols[j], blockCols);
01118         int pos = blkColOffsets[j]*rowsize + blkRowOffsets[i];
01119 
01120         coefs2dPtr[rowind*blockCols.size()+colind][pos] += values[i][j];
01121       }
01122     }
01123 
01124     for(size_t i=0; i<blockRows.size(); ++i) {
01125       CHK_ERR( giveToUnderlyingBlockMatrix(blockRows[i],
01126                                            blockRowSizes[i],
01127                                            blockCols.size(),
01128                                            &blockCols[0],
01129                                            &blockColSizes[0],
01130                                            &blockColSizes[0],
01131                                            &coefs2dPtr[i*blockCols.size()],
01132                                            true) );
01133     }
01134     changedSinceMark_ = true;
01135 
01136     return(0);
01137   }
01138 
01139   int maxBlkEqnSize = pointBlockMap->getMaxBlkEqnSize();
01140   int coefBlkLen = maxBlkEqnSize*maxBlkEqnSize*2;
01141 
01142   for(int i=0; i<numRows; ++i) {
01143     int row = rows[i];
01144 
01145     if (row < firstLocalOffset() || row > lastLocalOffset()) {
01146       int proc = eqnComm_->getOwnerProc(row);
01147       FillableMat* remote_mat = getRemotelyOwnedMatrix(proc);
01148       if (sumInto) {
01149         remote_mat->sumInRow(row, cols, values[i], numCols);
01150       }
01151       else {
01152         remote_mat->putRow(row, cols, values[i], numCols);
01153       }
01154       changedSinceMark_ = true;
01155       continue;
01156     }
01157 
01158     int blockRow = pointBlockMap->eqnToBlkEqn(row);
01159     int blockRowSize = pointBlockMap->getBlkEqnSize(blockRow);
01160     int blockRowOffset = pointBlockMap->getBlkEqnOffset(blockRow, row);
01161 
01162     int blockRowLength = 0;
01163     CHK_ERR( snl_fei::BlockMatrixTraits<T>::getRowLength(matrix_.get(), blockRow,
01164                                                 blockRowLength) );
01165 
01166     std::vector<int> blkCols(blockRowLength);
01167     int* blkCols_ptr = &blkCols[0];
01168     std::vector<int> blkColDims(blockRowLength);
01169     int* blkColDims_ptr = &blkColDims[0];
01170     std::vector<double> coefs_1D(blockRowLength*coefBlkLen);
01171     double* coefs_1D_ptr = &coefs_1D[0];
01172     int coefsLen = coefs_1D.size();
01173     std::vector<double*> coefs_2D(blockRowLength);
01174     double** coefs_2D_ptr = &coefs_2D[0];
01175 
01176     std::vector<int> LDAs(blockRowLength);
01177     std::fill(LDAs.begin(), LDAs.end(), blockRowSize);
01178     std::fill(coefs_1D.begin(), coefs_1D.end(), 0.0);
01179 
01180     int checkRowLen = 0;
01181     CHK_ERR( snl_fei::BlockMatrixTraits<T>::copyOutRow(matrix_.get(),
01182                                               blockRow, blockRowLength,
01183                                               blockRowSize,
01184                                               blkCols_ptr,
01185                                               blkColDims_ptr,
01186                                               coefs_1D_ptr,
01187                                               coefsLen,
01188                                               checkRowLen) );
01189     int coefs_1D_offset = 0;
01190     for(int j=0; j<checkRowLen; ++j) {
01191       coefs_2D_ptr[j] = &(coefs_1D_ptr[coefs_1D_offset]);
01192       coefs_1D_offset += blockRowSize*blkColDims_ptr[j];
01193     }
01194 
01195     for(int j=0; j<numCols; ++j) {
01196       int blockCol, blkOffset;
01197       CHK_ERR( pointBlockMap->getPtEqnInfo(cols[j], blockCol, blkOffset) );
01198 
01199       for(int jj=0; jj<blockRowLength; ++jj) {
01200 
01201         if (blockCol == blkCols_ptr[jj]) {
01202           if (sumInto) {
01203             coefs_2D_ptr[jj][blkOffset*blockRowSize+blockRowOffset] += values[i][j];
01204           }
01205           else {
01206             coefs_2D_ptr[jj][blkOffset*blockRowSize+blockRowOffset] = values[i][j];
01207           }
01208 
01209           break;
01210         }
01211       }
01212     }
01213 
01214     //Now put the block-row back into the matrix
01215     CHK_ERR( giveToUnderlyingBlockMatrix(blockRow, blockRowSize,
01216                                          blockRowLength, blkCols_ptr,
01217                                          &LDAs[0],
01218                                          blkColDims_ptr,
01219                                          coefs_2D_ptr,
01220                                          false) );
01221 
01222     changedSinceMark_ = true;
01223   }
01224 
01225   return(0);
01226 }
01227 
01228 //----------------------------------------------------------------------------
01229 template<typename T>
01230 int fei::Matrix_Impl<T>::writeToFile(const char* filename,
01231                                     bool matrixMarketFormat)
01232 {
01233   fei::SharedPtr<fei::MatrixGraph> mgraph = getMatrixGraph();
01234   if (mgraph.get() == NULL) {
01235     return(-1);
01236   }
01237 
01238   if (haveFEMatrix()) {
01239     return(0);//temporary no-op...
01240   }
01241 
01242   if (haveBlockMatrix()) {
01243     CHK_ERR( snl_fei::BlockMatrixTraits<T>::globalAssemble(matrix_.get()) );
01244   }
01245 
01246   static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
01247 
01248   fei::SharedPtr<fei::VectorSpace> vspace = mgraph->getRowSpace();
01249 
01250   int globalNNZ = 0;
01251   int globalNumRows = vspace->getGlobalNumIndices();
01252   int localNumRows = vspace->getNumIndices_Owned();
01253 
01254   fei::SharedPtr<fei::VectorSpace> cspace = mgraph->getColSpace();
01255   int globalNumCols = globalNumRows;
01256   if (cspace.get() != NULL) {
01257     globalNumCols = cspace->getGlobalNumIndices();
01258   }
01259 
01260   std::vector<int> indices_owned;
01261   int localNNZ = 0;
01262   CHK_ERR( vspace->getIndices_Owned(indices_owned) );
01263   int* rowsPtr = &indices_owned[0];
01264   for(int i=0; i<localNumRows; ++i) {
01265     int len;
01266     CHK_ERR( getRowLength(rowsPtr[i], len) );
01267     localNNZ += len;
01268   }
01269 
01270   CHK_MPI( GlobalSum(getCommunicator(), localNNZ, globalNNZ) );
01271 
01272   for(int p=0; p<numProcs(); ++p) {
01273     fei::Barrier(getCommunicator());
01274     if (p != localProc()) continue;
01275 
01276     FEI_OFSTREAM* outFile = NULL;
01277     if (p==0) {
01278       outFile = new FEI_OFSTREAM(filename, IOS_OUT);
01279       FEI_OFSTREAM& ofs = *outFile;
01280       if (matrixMarketFormat) {
01281         ofs << mmbanner << FEI_ENDL;
01282         ofs <<globalNumRows << " " <<globalNumCols << " " <<globalNNZ <<FEI_ENDL;
01283       }
01284       else {
01285         ofs <<globalNumRows << " " <<globalNumCols <<FEI_ENDL;
01286       }
01287     }
01288     else outFile = new FEI_OFSTREAM(filename, IOS_APP);
01289 
01290     outFile->setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
01291     outFile->precision(13);
01292     FEI_OFSTREAM& ofs = *outFile;
01293 
01294     int rowLength;
01295 
01296     for(int i=firstLocalOffset(); i<=lastLocalOffset(); ++i) {
01297       CHK_ERR( getRowLength(i, rowLength) );
01298 
01299       work_indices_.resize(rowLength);
01300       work_data1D_.resize(rowLength);
01301 
01302       int* indPtr = &work_indices_[0];
01303       double* coefPtr = &work_data1D_[0];
01304 
01305       CHK_ERR( copyOutRow(i, rowLength, coefPtr, indPtr) );
01306 
01307       fei::insertion_sort_with_companions<double>(rowLength,
01308                                             indPtr, coefPtr);
01309 
01310       for(int j=0; j<rowLength; ++j) {
01311         if (matrixMarketFormat) {
01312           ofs << i+1 <<" "<<indPtr[j]+1<<" "<<coefPtr[j]<<FEI_ENDL;
01313         }
01314         else {
01315           ofs << i <<" "<<indPtr[j]<<" "<<coefPtr[j]<<FEI_ENDL;
01316         }
01317       }
01318     }
01319 
01320     delete outFile;
01321   }
01322 
01323   return(0);
01324 }
01325 
01326 //----------------------------------------------------------------------------
01327 template<typename T>
01328 int fei::Matrix_Impl<T>::writeToStream(FEI_OSTREAM& ostrm,
01329                                       bool matrixMarketFormat)
01330 {
01331   fei::SharedPtr<fei::MatrixGraph> mgraph = getMatrixGraph();
01332   if (mgraph.get() == NULL) {
01333     return(-1);
01334   }
01335 
01336   if (haveFEMatrix()) {
01337     return(0);//temporary no-op...
01338   }
01339 
01340   if (haveBlockMatrix()) {
01341     CHK_ERR( snl_fei::BlockMatrixTraits<T>::globalAssemble(matrix_.get()) );
01342   }
01343 
01344   static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
01345 
01346   fei::SharedPtr<fei::VectorSpace> vspace = mgraph->getRowSpace();
01347 
01348   int globalNNZ = 0;
01349   int localNumRows = vspace->getNumIndices_Owned();
01350   std::vector<int> indices_owned;
01351   int localNNZ = 0;
01352   CHK_ERR( vspace->getIndices_Owned(indices_owned));
01353   int* rowsPtr = &indices_owned[0];
01354   for(int i=0; i<localNumRows; ++i) {
01355     int len;
01356     CHK_ERR( getRowLength(rowsPtr[i], len) );
01357     localNNZ += len;
01358   }
01359 
01360   CHK_MPI( fei::GlobalSum(getCommunicator(), localNNZ, globalNNZ) );
01361 
01362   IOS_FMTFLAGS oldf = ostrm.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
01363 
01364   for(int p=0; p<numProcs(); ++p) {
01365     fei::Barrier(getCommunicator());
01366     if (p != localProc()) continue;
01367 
01368     if (p==0) {
01369       int globalSize = globalOffsets()[numProcs()]-1;
01370       if (matrixMarketFormat) {
01371         ostrm << mmbanner << FEI_ENDL;
01372         ostrm << globalSize << " " << globalSize << " " << globalNNZ << FEI_ENDL;
01373       }
01374       else {
01375         ostrm << globalSize << " " << globalSize << FEI_ENDL;
01376       }
01377     }
01378 
01379     int rowLength;
01380 
01381     for(int i=firstLocalOffset(); i<=lastLocalOffset(); ++i) {
01382       CHK_ERR( getRowLength(i, rowLength) );
01383 
01384       work_indices_.resize(rowLength);
01385       work_data1D_.resize(rowLength);
01386 
01387       int* indPtr = &work_indices_[0];
01388       double* coefPtr = &work_data1D_[0];
01389 
01390       CHK_ERR( copyOutRow(i, rowLength, coefPtr, indPtr) );
01391 
01392       for(int j=0; j<rowLength; ++j) {
01393         ostrm << i << " " << indPtr[j] << " " << coefPtr[j] << FEI_ENDL;
01394       }
01395     }
01396   }
01397 
01398   ostrm.setf(oldf, IOS_FLOATFIELD);
01399 
01400   return(0);
01401 }
01402 
01403 #endif // _fei_Matrix_Impl_hpp_
01404 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends