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