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