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