FEI Version of the Day
fei_MatrixGraph_Impl2.cpp
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 #include <fei_sstream.hpp>
00045 
00046 #include <limits>
00047 #include <cmath>
00048 
00049 #include <fei_MatrixGraph_Impl2.hpp>
00050 
00051 #include <fei_utils.hpp>
00052 #include "fei_TemplateUtils.hpp"
00053 
00054 #include <fei_Pattern.hpp>
00055 #include <fei_LogManager.hpp>
00056 #include <fei_TemplateUtils.hpp>
00057 #include <fei_impl_utils.hpp>
00058 #include <snl_fei_Utils.hpp>
00059 #include <fei_FieldMask.hpp>
00060 #include <fei_Record.hpp>
00061 #include <snl_fei_RecordCollection.hpp>
00062 #include <fei_VectorSpace.hpp>
00063 #include <fei_ParameterSet.hpp>
00064 #include <fei_ostream_ops.hpp>
00065 #include <fei_Reducer.hpp>
00066 #include <fei_GraphReducer.hpp>
00067 #include <fei_ConnectivityBlock.hpp>
00068 #include <snl_fei_BlkSizeMsgHandler.hpp>
00069 #include <fei_Graph_Impl.hpp>
00070 #include <snl_fei_Constraint.hpp>
00071 
00072 #include <fei_EqnBuffer.hpp>
00073 #include <fei_EqnCommMgr.hpp>
00074 #include <SNL_FEI_Structure.hpp>
00075 
00076 #undef fei_file
00077 #define fei_file "fei_MatrixGraph.cpp"
00078 
00079 #include <fei_ErrMacros.hpp>
00080 
00081 namespace snl_fei {
00082 static unsigned getFieldSize(int fieldID,
00083                       fei::VectorSpace* space1,
00084                       fei::VectorSpace* space2)
00085 {
00086   unsigned fieldsize = 0;
00087   bool foundfield = false;
00088   if (space1 != NULL) {
00089     try {
00090       fieldsize = space1->getFieldSize(fieldID);
00091       foundfield = true;
00092     }
00093     catch (...) {
00094       foundfield = false;
00095     }
00096   }
00097 
00098   if (!foundfield && space2 != NULL) {
00099     try {
00100       fieldsize = space2->getFieldSize(fieldID);
00101     }
00102     catch (std::runtime_error& exc) {
00103       std::string msg("snl_fei::getFieldSize: ");
00104       msg += exc.what();
00105       throw std::runtime_error(msg);
00106     }
00107   }
00108 
00109   return(fieldsize);
00110 }
00111 
00112 }//namespace snl_fei
00113 
00114 //------------------------------------------------------------------------------
00115 fei::SharedPtr<fei::MatrixGraph>
00116 fei::MatrixGraph_Impl2::Factory::createMatrixGraph(fei::SharedPtr<fei::VectorSpace> rowSpace,
00117                                              fei::SharedPtr<fei::VectorSpace> columnSpace,
00118                                              const char* name)
00119 {
00120   fei::SharedPtr<fei::MatrixGraph> mgptr(new fei::MatrixGraph_Impl2(rowSpace,
00121                                                               columnSpace,
00122                                                               name));
00123 
00124   return(mgptr);
00125 }
00126 
00127 //=====Constructor==============================================================
00128 fei::MatrixGraph_Impl2::MatrixGraph_Impl2(fei::SharedPtr<fei::VectorSpace> rowSpace,
00129                               fei::SharedPtr<fei::VectorSpace> colSpace,
00130                               const char* name)
00131  : comm_(rowSpace->getCommunicator()),
00132    rowSpace_(rowSpace),
00133    colSpace_(colSpace),
00134    haveRowSpace_(false),
00135    haveColSpace_(false),
00136    symmetric_(false),
00137    remotelyOwnedGraphRows_(NULL),
00138    simpleProblem_(false),
00139    blockEntryGraph_(false),
00140    patterns_(),
00141    connectivityBlocks_(),
00142    arbitraryBlockCounter_(1),
00143    sparseBlocks_(),
00144    lagrangeConstraints_(),
00145    penaltyConstraints_(),
00146    slaveConstraints_(),
00147    ptEqualBlk_(false),
00148    newSlaveData_(false),
00149    localNumSlaves_(0),
00150    globalNumSlaves_(0),
00151    D_(NULL),
00152    g_(),
00153    g_nonzero_(false),
00154    reducer_(),
00155    name_(),
00156    dbgprefix_("MatGrph: "),
00157    tmpIntArray1_(),
00158    tmpIntArray2_()
00159 {
00160   localProc_ = fei::localProc(comm_);
00161   numProcs_ = fei::numProcs(comm_);
00162 
00163   if (rowSpace_.get() == NULL) {
00164     voidERReturn;
00165   }
00166   else haveRowSpace_ = true;
00167 
00168   if (colSpace_.get() != NULL) haveColSpace_ = true;
00169   else colSpace_ = rowSpace_;
00170 
00171   setName(name);
00172 }
00173 
00174 //------------------------------------------------------------------------------
00175 fei::MatrixGraph_Impl2::~MatrixGraph_Impl2()
00176 {
00177   fei::destroyValues(patterns_);
00178   patterns_.clear();
00179 
00180   fei::destroyValues(connectivityBlocks_);
00181   connectivityBlocks_.clear();
00182 
00183   int i, len = sparseBlocks_.size();
00184   for(i=0; i<len; ++i) {
00185     delete sparseBlocks_[i];
00186   }
00187 
00188   fei::destroyValues(lagrangeConstraints_);
00189   lagrangeConstraints_.clear();
00190 
00191   fei::destroyValues(penaltyConstraints_);
00192   penaltyConstraints_.clear();
00193 
00194   fei::destroyValues(slaveConstraints_);
00195   slaveConstraints_.clear();
00196 }
00197 
00198 //------------------------------------------------------------------------------
00199 void fei::MatrixGraph_Impl2::setParameters(const fei::ParameterSet& params)
00200 {
00201   const fei::Param* param = 0;
00202   fei::Param::ParamType ptype = fei::Param::BAD_TYPE;
00203 
00204   param = params.get("FEI_OUTPUT_LEVEL");
00205   ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
00206   if (ptype == fei::Param::STRING) {
00207     fei::LogManager& log_manager = fei::LogManager::getLogManager();
00208     log_manager.setOutputLevel(param->getStringValue().c_str());
00209     setOutputLevel(fei::utils::string_to_output_level(param->getStringValue()));
00210   }
00211 
00212   param = params.get("FEI_LOG_EQN");
00213   ptype =  param != NULL ? param->getType() : fei::Param::BAD_TYPE;  
00214   if (ptype == fei::Param::INT) {
00215     addLogEqn(param->getIntValue());
00216   }
00217 
00218   param = params.get("FEI_LOG_ID");
00219   ptype =  param != NULL ? param->getType() : fei::Param::BAD_TYPE;  
00220   if (ptype == fei::Param::INT) {
00221     addLogID(param->getIntValue());
00222   }
00223 
00224   param = params.get("name");
00225   ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
00226   if (ptype == fei::Param::STRING) {
00227     setName(param->getStringValue().c_str());
00228   }
00229 
00230   param = params.get("BLOCK_GRAPH");
00231   ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
00232   if (ptype != fei::Param::BAD_TYPE) {
00233     blockEntryGraph_ = true;
00234   }
00235 
00236   param = params.get("BLOCK_MATRIX");
00237   ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
00238   if (ptype != fei::Param::BAD_TYPE) {
00239     if (ptype == fei::Param::BOOL) {
00240       blockEntryGraph_ = param->getBoolValue();
00241     }
00242     else {
00243       blockEntryGraph_ = true;
00244     }
00245   }
00246 }
00247 
00248 //----------------------------------------------------------------------------
00249 void fei::MatrixGraph_Impl2::setRowSpace(fei::SharedPtr<fei::VectorSpace> rowSpace)
00250 {
00251   rowSpace_ = rowSpace;
00252   haveRowSpace_ = true;
00253   if (!haveColSpace_) symmetric_ = true;
00254 }
00255 
00256 //----------------------------------------------------------------------------
00257 void fei::MatrixGraph_Impl2::setColumnSpace(fei::SharedPtr<fei::VectorSpace> colSpace)
00258 {
00259   colSpace_ = colSpace;
00260   haveColSpace_ = true;
00261   if (colSpace_ == rowSpace_) symmetric_ = true;
00262   else symmetric_ = false;
00263 }
00264 
00265 //----------------------------------------------------------------------------
00266 int fei::MatrixGraph_Impl2::addPattern(fei::Pattern* pattern)
00267 {
00268   std::map<int,fei::Pattern*>::iterator
00269     p_iter = patterns_.begin(), p_iter_end = patterns_.end();
00270 
00271   bool matches_existing_pattern = false;
00272   int pattern_id = -1;
00273   while(p_iter != p_iter_end && !matches_existing_pattern) {
00274     if (*pattern == *(p_iter->second)) {
00275       matches_existing_pattern = true;
00276       pattern_id = p_iter->first;
00277     }
00278     else ++p_iter;
00279   }
00280 
00281   if (!matches_existing_pattern) {
00282     pattern_id = patterns_.size();
00283     patterns_.insert(std::make_pair(pattern_id, pattern));
00284   }
00285   else {
00286     delete pattern;
00287   }
00288 
00289   return pattern_id;
00290 }
00291 
00292 //----------------------------------------------------------------------------
00293 int fei::MatrixGraph_Impl2::definePattern(int numIDs,
00294                                     int idType)
00295 {
00296   snl_fei::RecordCollection* rec_coll = NULL;
00297   rowSpace_->getRecordCollection(idType, rec_coll);
00298 
00299   fei::Pattern* pattern = new fei::Pattern(numIDs, idType, rec_coll);
00300   return addPattern(pattern);
00301 }
00302 
00303 //----------------------------------------------------------------------------
00304 int fei::MatrixGraph_Impl2::definePattern(int numIDs,
00305                                       int idType,
00306                                       int fieldID)
00307 {
00308   unsigned fieldsize = 0;
00309   try {
00310     fieldsize = snl_fei::getFieldSize(fieldID, rowSpace_.get(), colSpace_.get());
00311   }
00312   catch (std::runtime_error& exc) {
00313     FEI_OSTRINGSTREAM osstr;
00314     osstr << "fei::MatrixGraph_Impl2::definePattern caught error: "<<exc.what();
00315     throw std::runtime_error(osstr.str());
00316   }
00317 
00318   snl_fei::RecordCollection* rec_coll = NULL;
00319   rowSpace_->getRecordCollection(idType, rec_coll);
00320 
00321   fei::Pattern* pattern =
00322     new fei::Pattern(numIDs, idType, rec_coll, fieldID, fieldsize);
00323   return addPattern(pattern);
00324 }
00325 
00326 //----------------------------------------------------------------------------
00327 int fei::MatrixGraph_Impl2::definePattern(int numIDs,
00328                                      int idType,
00329                                      const int* numFieldsPerID,
00330                                      const int* fieldIDs)
00331 { 
00332   std::vector<int> fieldSizes;
00333   try {
00334     int offset = 0;
00335     for(int i=0; i<numIDs; ++i) {
00336       for(int j=0; j<numFieldsPerID[i]; ++j) {
00337         fieldSizes.push_back(snl_fei::getFieldSize(fieldIDs[offset++],
00338                                               rowSpace_.get(), colSpace_.get()));
00339       }
00340     }
00341   }
00342   catch (std::runtime_error& exc) {
00343     FEI_OSTRINGSTREAM osstr;
00344     osstr << "fei::MatrixGraph_Impl2::definePattern caught error: "<<exc.what();
00345     throw std::runtime_error(osstr.str());
00346   }
00347 
00348   snl_fei::RecordCollection* rec_coll = NULL;
00349   rowSpace_->getRecordCollection(idType, rec_coll);
00350 
00351   fei::Pattern* pattern = new fei::Pattern(numIDs, idType, rec_coll,
00352                         numFieldsPerID, fieldIDs, &(fieldSizes[0]));
00353   return addPattern(pattern);
00354 }
00355 
00356 //----------------------------------------------------------------------------
00357 int fei::MatrixGraph_Impl2::definePattern(int numIDs,
00358                                       const int* idTypes,
00359                                       const int* numFieldsPerID,
00360                                       const int* fieldIDs)
00361 {
00362   std::vector<int> fieldSizes;
00363   try {
00364     int offset = 0;
00365     for(int i=0; i<numIDs; ++i) {
00366       for(int j=0; j<numFieldsPerID[i]; ++j) {
00367         fieldSizes.push_back(snl_fei::getFieldSize(fieldIDs[offset++],
00368                                               rowSpace_.get(), colSpace_.get()));
00369       }
00370     }
00371   }
00372   catch (std::runtime_error& exc) {
00373     FEI_OSTRINGSTREAM osstr;
00374     osstr << "fei::MatrixGraph_Impl2::definePattern caught error: "<<exc.what();
00375     throw std::runtime_error(osstr.str());
00376   }
00377 
00378   std::vector<snl_fei::RecordCollection*> recordCollections(numIDs);
00379   for(int i=0; i<numIDs; ++i) {
00380     rowSpace_->getRecordCollection(idTypes[i], recordCollections[i]);
00381   }
00382 
00383   fei::Pattern* pattern = new fei::Pattern(numIDs, idTypes, &recordCollections[0],
00384                         numFieldsPerID, fieldIDs, &(fieldSizes[0]));
00385   return addPattern(pattern);
00386 }
00387 
00388 //------------------------------------------------------------------------------
00389 fei::Pattern* fei::MatrixGraph_Impl2::getPattern(int patternID)
00390 {
00391   std::map<int,fei::Pattern*>::iterator
00392     p_iter = patterns_.find(patternID);
00393 
00394   if (p_iter == patterns_.end()) {
00395     return NULL;
00396   }
00397 
00398   fei::Pattern* pattern = (*p_iter).second;
00399   return pattern;
00400 }
00401 
00402 //------------------------------------------------------------------------------
00403 int fei::MatrixGraph_Impl2::initConnectivityBlock(int blockID,
00404                                             int numConnectivityLists,
00405                                             int patternID,
00406                                             bool diagonal)
00407 {
00408   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00409     (*output_stream_) <<dbgprefix_<< "initConnectivityBlock symm, blkID="
00410      << blockID<<", num="<<numConnectivityLists<<", patternID="
00411      <<patternID<<FEI_ENDL;
00412   }
00413 
00414   if (blockID < 0) {
00415     fei::console_out() << "fei::MatrixGraph_Impl2::initConnectivityBlock: blockID ("
00416          << blockID << ") must be non-negative." << FEI_ENDL;
00417     ERReturn(-1);
00418   }
00419 
00420   std::map<int,fei::Pattern*>::iterator
00421     p_iter = patterns_.find(patternID);
00422 
00423   if (p_iter == patterns_.end()) {
00424     ERReturn(-1);
00425   }
00426 
00427   fei::Pattern* pattern = (*p_iter).second;
00428 
00429   if (getConnectivityBlock(blockID) != NULL) ERReturn(-1);
00430 
00431   fei::ConnectivityBlock* cblock =
00432     new fei::ConnectivityBlock(blockID, pattern, numConnectivityLists);
00433 
00434   connectivityBlocks_.insert(std::pair<int,fei::ConnectivityBlock*>(blockID, cblock));
00435   if (diagonal) cblock->setIsDiagonal(diagonal);
00436 
00437   return(0);
00438 }
00439 
00440 //------------------------------------------------------------------------------
00441 int fei::MatrixGraph_Impl2::initConnectivityBlock(int numConnectivityLists,
00442                                                   int patternID,
00443                                                   bool diagonal)
00444 {
00445   int blockID = connectivityBlocks_.size();
00446   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00447     (*output_stream_) <<dbgprefix_<< "initConnectivityBlock symm, blkID="
00448      << blockID<<", num="<<numConnectivityLists<<", patternID="
00449      <<patternID<<FEI_ENDL;
00450   }
00451 
00452   std::map<int,fei::Pattern*>::iterator
00453     p_iter = patterns_.find(patternID);
00454 
00455   if (p_iter == patterns_.end()) {
00456     FEI_OSTRINGSTREAM osstr;
00457     osstr << "fei::MatrixGraph_Impl2::initConnectivityBlock: ERROR, patternID "
00458       << patternID << " not found.";
00459     throw std::runtime_error(osstr.str());
00460   }
00461 
00462   fei::Pattern* pattern = (*p_iter).second;
00463 
00464   fei::ConnectivityBlock* cblock =
00465     new fei::ConnectivityBlock(blockID, pattern, numConnectivityLists);
00466 
00467   connectivityBlocks_.insert(std::pair<int,fei::ConnectivityBlock*>(blockID, cblock));
00468   if (diagonal) cblock->setIsDiagonal(diagonal);
00469 
00470   return(blockID);
00471 }
00472 
00473 //------------------------------------------------------------------------------
00474 int fei::MatrixGraph_Impl2::initConnectivityBlock(int blockID,
00475                                             int numConnectivityLists,
00476                                             int rowPatternID,
00477                                             int colPatternID)
00478 {
00479   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00480     (*output_stream_) <<dbgprefix_<< "initConnectivityBlock nonsymm, blkID="
00481      << blockID<<", num="<<numConnectivityLists<<", row-patternID="
00482      <<rowPatternID<<", col-patternID="<<colPatternID<<FEI_ENDL;
00483   }
00484 
00485   if (blockID < 0) {
00486     FEI_OSTRINGSTREAM osstr;
00487     osstr << "fei::MatrixGraph_Impl2::initConnectivityBlock: blockID ("
00488           << blockID << ") must be non-negative." << FEI_ENDL;
00489     throw std::runtime_error(osstr.str());
00490   }
00491 
00492   std::map<int,fei::Pattern*>::iterator
00493     rp_iter = patterns_.find(rowPatternID);
00494   if (rp_iter == patterns_.end()) ERReturn(-1);
00495 
00496   fei::Pattern* rowpattern = (*rp_iter).second;
00497 
00498   std::map<int,fei::Pattern*>::iterator
00499     cp_iter = patterns_.find(colPatternID);
00500   if (cp_iter == patterns_.end()) ERReturn(-1);
00501 
00502   fei::Pattern* colpattern = (*cp_iter).second;
00503 
00504   if (getConnectivityBlock(blockID) != NULL) ERReturn(-1);
00505 
00506   fei::ConnectivityBlock* cblock =
00507     new fei::ConnectivityBlock(blockID, rowpattern,
00508                                    colpattern, numConnectivityLists);
00509 
00510   connectivityBlocks_.insert(std::pair<int,fei::ConnectivityBlock*>(blockID, cblock));
00511 
00512   return(0);
00513 }
00514 
00515 //------------------------------------------------------------------------------
00516 int fei::MatrixGraph_Impl2::initConnectivity(int blockID,
00517                                        int connectivityID,
00518                                        const int* connectedIdentifiers)
00519 {
00520   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00521     (*output_stream_) <<dbgprefix_<< "initConnectivity blkID="
00522      <<blockID<<", connID="<<connectivityID<<FEI_ENDL;
00523   }
00524 
00525   //for each connected identifier, get a pointer to its record from the
00526   //solution-space.
00527   //store those pointers in the appropriate connectivity-table, mapped to
00528   //the user-provided connectivityID.
00529 
00530   fei::ConnectivityBlock* connblk = getConnectivityBlock(blockID);
00531   if (connblk == NULL) {
00532     FEI_OSTRINGSTREAM osstr;
00533     osstr << "fei::MatrixGraph_Impl2::initConnectivity ERROR, blockID " << blockID
00534           << " doesn't correspond to an existing ConnectivityBlock.";
00535     throw std::runtime_error(osstr.str());
00536   }
00537 
00538   fei::Pattern* pattern = connblk->getRowPattern();
00539   int numIDs = pattern->getNumIDs();
00540 
00541   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00542     for(int ii=0; ii<numIDs; ++ii) {
00543       if (isLogID(connectedIdentifiers[ii])) {
00544         FEI_OSTREAM& os = *output_stream_;
00545         os << dbgprefix_<<"initConnectivity blockID " << blockID << ", connID "
00546            << connectivityID << ": ";
00547         for(int jj=0; jj<numIDs; ++jj) {
00548           os << connectedIdentifiers[jj] << " ";
00549         }
00550         os << FEI_ENDL;
00551       }
00552     }
00553   }
00554 
00555   if (rowSpace_.get() == NULL) ERReturn(-1);
00556 
00557   IndexType<int,int>& connectivityIDs = connblk->getNativeConnectivityIDs();
00558 
00559   int idOffset = -1;
00560   IndexType<int,int>::iterator
00561    iter = connectivityIDs.find(connectivityID);
00562   if (iter == connectivityIDs.end()) {
00563     idOffset = connectivityIDs.size();
00564     connectivityIDs.insert(iter, std::make_pair(connectivityID,idOffset));
00565   }
00566   else {
00567     idOffset = iter->second;
00568   }
00569 
00570   idOffset *= pattern->getNumIDs();
00571 
00572   int* rlist = &(connblk->getRowConnectivities()[idOffset]);
00573 
00574   CHK_ERR( getConnectivityRecords(pattern, rowSpace_.get(),
00575                                   connectedIdentifiers, rlist) );
00576 
00577   for(int i=0; i<numIDs; ++i) {
00578     if (pattern->getNumFieldsPerID()[i] > 0) {
00579       pattern->getRecordCollections()[i]->getRecordWithLocalID(rlist[i])->isInLocalSubdomain_ = true;
00580     }
00581   }
00582 
00583   return(0);
00584 }
00585 
00586 //------------------------------------------------------------------------------
00587 int fei::MatrixGraph_Impl2::initConnectivity(int idType,
00588                                            int numRows,
00589                                            const int* rowIDs,
00590                                            const int* rowOffsets,
00591                                            const int* packedColumnIDs)
00592 {
00593   fei::ConnectivityBlock* block = new fei::ConnectivityBlock(numRows,
00594                                                    rowIDs, rowOffsets);
00595 
00596   int max_row_len = 0;
00597   for(int i=0; i<numRows; ++i) {
00598     int row_len = rowOffsets[i+1]-rowOffsets[i];
00599     if (row_len > max_row_len) max_row_len = row_len;
00600   }
00601 
00602   int patternID = definePattern(max_row_len, idType);
00603   fei::Pattern* pattern = getPattern(patternID);
00604   block->setRowPattern(pattern);
00605 
00606   sparseBlocks_.push_back(block);
00607 
00608   int* row_records = &(block->getRowConnectivities()[0]);
00609 
00610   CHK_ERR( getConnectivityRecords(rowSpace_.get(),
00611                                   idType, numRows, rowIDs, row_records) );
00612 
00613   fei::VectorSpace* colSpace = rowSpace_.get();
00614   if (colSpace_.get() != NULL) {
00615     colSpace = colSpace_.get();
00616   }
00617 
00618   int* col_records = &(block->getColConnectivities()[0]);
00619 
00620   CHK_ERR( getConnectivityRecords(colSpace, idType, rowOffsets[numRows],
00621                                   packedColumnIDs, col_records));
00622 
00623   return(0);
00624 }
00625 
00626 //------------------------------------------------------------------------------
00627 int fei::MatrixGraph_Impl2::initConnectivity(int idType,
00628                                        int numRows,
00629                                        const int* rowIDs,
00630                                        const int* rowLengths,
00631                                        const int*const* columnIDs)
00632 {
00633   fei::ConnectivityBlock* block = new fei::ConnectivityBlock(numRows,
00634                                                    rowIDs, rowLengths, true);
00635 
00636   int max_row_len = 0;
00637   for(int i=0; i<numRows; ++i) {
00638     int row_len = rowLengths[i];
00639     if (row_len > max_row_len) max_row_len = row_len;
00640   }
00641 
00642   int patternID = definePattern(max_row_len, idType);
00643   fei::Pattern* pattern = getPattern(patternID);
00644   block->setRowPattern(pattern);
00645 
00646   sparseBlocks_.push_back(block);
00647 
00648   int* row_records = &(block->getRowConnectivities()[0]);
00649 
00650   CHK_ERR( getConnectivityRecords(rowSpace_.get(),
00651                                   idType, numRows, rowIDs, row_records) );
00652 
00653   fei::VectorSpace* colSpace = rowSpace_.get();
00654   if (colSpace_.get() != NULL) {
00655     colSpace = colSpace_.get();
00656   }
00657 
00658   int* col_records = &(block->getColConnectivities()[0]);
00659 
00660   int offset = 0;
00661   for(int i=0; i<numRows; ++i) {
00662     CHK_ERR( getConnectivityRecords(colSpace, idType, rowLengths[i],
00663                                     columnIDs[i], &(col_records[offset])));
00664     offset += rowLengths[i];
00665   }
00666 
00667   return(0);
00668 }
00669 
00670 //------------------------------------------------------------------------------
00671 int fei::MatrixGraph_Impl2::getConnectivityRecords(fei::VectorSpace* vecSpace,
00672                                                  int idType,
00673                                                  int numIDs,
00674                                                  const int* IDs,
00675                                                  int* records)
00676 {
00677   snl_fei::RecordCollection* collection = NULL;
00678   CHK_ERR( vecSpace->getRecordCollection(idType, collection) );
00679 
00680   for(int i=0; i<numIDs; ++i) {
00681     records[i] = collection->getLocalID(IDs[i]);
00682     if (records[i] == -1) {
00683       CHK_ERR( vecSpace->addDOFs(idType, 1, &IDs[i], &records[i]) );
00684     }
00685   }
00686 
00687   return(0);
00688 }
00689 
00690 //------------------------------------------------------------------------------
00691 int fei::MatrixGraph_Impl2::getConnectivityRecords(fei::VectorSpace* vecSpace,
00692                                                  int idType,
00693                                                  int fieldID,
00694                                                  int numIDs,
00695                                                  const int* IDs,
00696                                                  int* records)
00697 {
00698   snl_fei::RecordCollection* collection = NULL;
00699   CHK_ERR( vecSpace->getRecordCollection(idType, collection) );
00700 
00701   for(int i=0; i<numIDs; ++i) {
00702     records[i] = collection->getLocalID(IDs[i]);
00703     if (records[i] == -1) {
00704       CHK_ERR( vecSpace->addDOFs(fieldID, idType, 1, &IDs[i], &records[i]));
00705     }
00706   }
00707 
00708   return(0);
00709 }
00710 
00711 //------------------------------------------------------------------------------
00712 int fei::MatrixGraph_Impl2::getConnectivityRecords(fei::Pattern* pattern,
00713                                              fei::VectorSpace* vecSpace,
00714                                              const int* connectedIdentifiers,
00715                                              int* recordList)
00716 {
00717   fei::Pattern::PatternType pType = pattern->getPatternType();
00718   int i, numIDs = pattern->getNumIDs();
00719   const int* numFieldsPerID = pattern->getNumFieldsPerID();
00720   const int* fieldIDs = pattern->getFieldIDs();
00721 
00722   if (pType == fei::Pattern::GENERAL) {
00723     const int* idTypes = pattern->getIDTypes();
00724 
00725     int fieldOffset = 0;
00726     for(i=0; i<numIDs; ++i) {
00727       int id = connectedIdentifiers[i];
00728 
00729       for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
00730         CHK_ERR( vecSpace->addDOFs(fieldIDs[fieldOffset++],
00731                                                idTypes[i], 1, &id,
00732                                                &recordList[i]));
00733       }
00734     }
00735   }
00736   else if (pType == fei::Pattern::SINGLE_IDTYPE ||
00737            pType == fei::Pattern::SIMPLE ||
00738            pType == fei::Pattern::NO_FIELD) {
00739 
00740     int idType = pattern->getIDTypes()[0];
00741 
00742     switch(pType) {
00743     case fei::Pattern::SINGLE_IDTYPE:
00744       {
00745         int fieldOffset = 0;
00746         for(i=0; i<numIDs; ++i) {
00747           int id = connectedIdentifiers[i];
00748           for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
00749             CHK_ERR( vecSpace->addDOFs(fieldIDs[fieldOffset++],
00750                                                    idType, 1, &id,
00751                                                    &(recordList[i])));
00752           }
00753         }
00754       }
00755       break;
00756     case fei::Pattern::SIMPLE:
00757       {
00758         CHK_ERR( vecSpace->addDOFs(fieldIDs[0], idType,
00759                                    numIDs, connectedIdentifiers,
00760                                    recordList) );
00761       }
00762       break;
00763     case fei::Pattern::NO_FIELD:
00764       CHK_ERR( vecSpace->addDOFs(idType, numIDs,
00765                                  connectedIdentifiers, recordList));
00766       break;
00767     case fei::Pattern::GENERAL:
00768       //Include a stub for this case to avoid compiler warning...
00769       std::abort(); break;
00770     }
00771   }
00772   else {
00773     ERReturn(-1);
00774   }
00775 
00776   return(0);
00777 }
00778 
00779 //------------------------------------------------------------------------------
00780 int fei::MatrixGraph_Impl2::initConnectivity(int blockID,
00781                                        int connectivityID,
00782                                        const int* rowConnectedIdentifiers,
00783                                        const int* colConnectedIdentifiers)
00784 {
00785   //for each connected identifier, get a pointer to its record from the
00786   //solution-space.
00787   //store those pointers in the appropriate connectivity-table, mapped to
00788   //the user-provided connectivityID.
00789 
00790   fei::ConnectivityBlock* connblk = getConnectivityBlock(blockID);
00791   if (connblk == NULL) {
00792     FEI_OSTRINGSTREAM osstr;
00793     osstr << "fei::MatrixGraph_Impl2::initConnectivity ERROR, blockID " << blockID
00794           << " doesn't correspond to an existing ConnectivityBlock.";
00795     throw std::runtime_error(osstr.str());
00796   }
00797 
00798   fei::Pattern* pattern = connblk->getRowPattern();
00799   int numIDs = pattern->getNumIDs();
00800   fei::Pattern* colPattern = connblk->getColPattern();
00801   int numColIDs = colPattern->getNumIDs();
00802   if (rowSpace_.get() == NULL) {
00803     ERReturn(-1);
00804   }
00805   if (colSpace_.get() == NULL) {
00806     ERReturn(-1);
00807   }
00808 
00809   IndexType<int,int>& connectivityIDs = connblk->getNativeConnectivityIDs();
00810 
00811   int i, idOffset = -1;
00812   IndexType<int,int>::iterator
00813     iter = connectivityIDs.find(connectivityID);
00814   if (iter == connectivityIDs.end()) {
00815     idOffset = connectivityIDs.size();
00816     connectivityIDs.insert(iter, std::make_pair(connectivityID,idOffset)); 
00817   }
00818   else {
00819     idOffset = iter->second;
00820   }
00821   int* row_rlist =
00822     &(connblk->getRowConnectivities()[idOffset*numIDs]);
00823   int* col_rlist =
00824     &(connblk->getColConnectivities()[idOffset*numColIDs]);
00825 
00826   CHK_ERR( getConnectivityRecords(pattern, rowSpace_.get(),
00827                                   rowConnectedIdentifiers, row_rlist) );
00828 
00829   for(i=0; i<numIDs; ++i)
00830     pattern->getRecordCollections()[i]->getRecordWithLocalID(row_rlist[i])->isInLocalSubdomain_ = true;
00831   CHK_ERR( getConnectivityRecords(colPattern, colSpace_.get(),
00832                                   colConnectedIdentifiers, col_rlist) );
00833 
00834   for(i=0; i<numColIDs; ++i)
00835     colPattern->getRecordCollections()[i]->getRecordWithLocalID(col_rlist[i])->isInLocalSubdomain_ = true;
00836 
00837   return(0);
00838 }
00839 
00840 //------------------------------------------------------------------------------
00841 int fei::MatrixGraph_Impl2::initConnectivity(int patternID,
00842                                        const int* connectedIdentifiers)
00843 {
00844   std::map<int,fei::Pattern*>::iterator
00845     p_iter = patterns_.find(patternID);
00846   if (p_iter == patterns_.end()) ERReturn(-1);
00847 
00848   fei::Pattern* pattern = p_iter->second;
00849 
00850   int blockID = -arbitraryBlockCounter_++;
00851   fei::ConnectivityBlock* cblock = new fei::ConnectivityBlock(blockID, pattern, 1);
00852 
00853   connectivityBlocks_.insert(std::pair<int,fei::ConnectivityBlock*>(blockID, cblock));
00854 
00855   CHK_ERR( initConnectivity(blockID, 0, connectedIdentifiers) );
00856   return(0);
00857 }
00858 
00859 //------------------------------------------------------------------------------
00860 int fei::MatrixGraph_Impl2::initConnectivity(int rowPatternID,
00861                                        const int* rowConnectedIdentifiers,
00862                                        int colPatternID,
00863                                        const int* colConnectedIdentifiers)
00864 {
00865   std::map<int,fei::Pattern*>::iterator
00866     rp_iter = patterns_.find(rowPatternID);
00867   if (rp_iter == patterns_.end()) ERReturn(-1);
00868 
00869   fei::Pattern* rowpattern = (*rp_iter).second;
00870 
00871   std::map<int,fei::Pattern*>::iterator
00872     cp_iter = patterns_.find(colPatternID);
00873   if (cp_iter == patterns_.end()) ERReturn(-1);
00874 
00875   fei::Pattern* colpattern = (*cp_iter).second;
00876 
00877   int blockID = -arbitraryBlockCounter_++;
00878   fei::ConnectivityBlock* cblock = new fei::ConnectivityBlock(blockID, rowpattern,
00879                                                     colpattern, 1);
00880 
00881   connectivityBlocks_.insert(std::pair<int,fei::ConnectivityBlock*>(blockID, cblock));
00882 
00883   CHK_ERR( initConnectivity(blockID, 0,
00884                             rowConnectedIdentifiers,
00885                             colConnectedIdentifiers) );
00886   return(0);
00887 }
00888 
00889 //------------------------------------------------------------------------------
00890 int fei::MatrixGraph_Impl2::getPatternNumIndices(int patternID,
00891                                            int& numIndices)
00892 {
00893   std::map<int,fei::Pattern*>::iterator
00894     p_iter = patterns_.find(patternID);
00895   if (p_iter == patterns_.end()) ERReturn(-1);
00896   fei::Pattern* pattern = (*p_iter).second;
00897 
00898   numIndices = pattern->getNumIndices();
00899 
00900   return(0);
00901 }
00902 
00903 //------------------------------------------------------------------------------
00904 int fei::MatrixGraph_Impl2::getPatternIndices(int patternID,
00905                                         const int* IDs,
00906                                         std::vector<int>& indices)
00907 {
00908   std::map<int,fei::Pattern*>::iterator
00909     p_iter = patterns_.find(patternID);
00910   if (p_iter == patterns_.end()) ERReturn(-1);
00911 
00912   fei::Pattern* pattern = (*p_iter).second;
00913 
00914   indices.resize(pattern->getNumIndices());
00915 
00916   int numIDs                = pattern->getNumIDs();
00917   const int* idTypes        = pattern->getIDTypes();
00918   const int* numFieldsPerID = pattern->getNumFieldsPerID();
00919   const int* fieldIDs       = pattern->getFieldIDs();
00920 
00921   int offset = 0, fieldOffset = 0;
00922   for(int i=0; i<numIDs; ++i) {
00923     for(int j=0; j<numFieldsPerID[i]; ++j) {
00924       CHK_ERR( rowSpace_->getGlobalIndices(1, &(IDs[i]), idTypes[i],
00925                                                   fieldIDs[fieldOffset],
00926                                                   &(indices[offset])) );
00927       offset += snl_fei::getFieldSize(fieldIDs[fieldOffset++],
00928                                         rowSpace_.get(), colSpace_.get());
00929     }
00930   }
00931 
00932   return(0);
00933 }
00934 
00935 //------------------------------------------------------------------------------
00936 int fei::MatrixGraph_Impl2::initConnectivity(int idType,
00937                                        int fieldID,
00938                                        int numRows,
00939                                        const int* rowIDs,
00940                                        const int* rowOffsets,
00941                                        const int* packedColumnIDs)
00942 {
00943   fei::ConnectivityBlock* block = new fei::ConnectivityBlock(fieldID, numRows,
00944                                                    rowIDs, rowOffsets);
00945 
00946   sparseBlocks_.push_back(block);
00947 
00948   int* row_records = &(block->getRowConnectivities()[0]);
00949 
00950   CHK_ERR( getConnectivityRecords(rowSpace_.get(), idType, fieldID,
00951                                   numRows, rowIDs, row_records) );
00952 
00953   fei::VectorSpace* colSpace = rowSpace_.get();
00954   if (colSpace_.get() != NULL) {
00955     colSpace = colSpace_.get();
00956   }
00957 
00958   int* col_records = &(block->getColConnectivities()[0]);
00959 
00960   CHK_ERR( getConnectivityRecords(colSpace, idType, fieldID, rowOffsets[numRows],
00961                                   packedColumnIDs, col_records));
00962 
00963   return(0);
00964 }
00965 
00966 //------------------------------------------------------------------------------
00967 int fei::MatrixGraph_Impl2::initLagrangeConstraint(int constraintID,
00968                                                int constraintIDType,
00969                                                int numIDs,
00970                                                const int* idTypes,
00971                                                const int* IDs,
00972                                                const int* fieldIDs)
00973 {
00974   if (rowSpace_.get() == NULL) ERReturn(-1);
00975 
00976   ConstraintType* constraint = getLagrangeConstraint(constraintID);
00977 
00978   CHK_ERR( rowSpace_->addDOFs(constraintIDType, 1, &constraintID) );
00979 
00980   if (haveColSpace_) {
00981     if (colSpace_.get() == NULL) {
00982       ERReturn(-1);
00983     }
00984     CHK_ERR( colSpace_->addDOFs(constraintIDType,
00985                                             1, &constraintID) );
00986   }
00987 
00988   ConstraintType* newconstraint =
00989     new ConstraintType(constraintID, constraintIDType,
00990                       false, //isSlave
00991                       false,  //isPenalty
00992                       numIDs, idTypes, IDs, fieldIDs,
00993                       0, //slaveOffset
00994                       0, //offsetIntoSlaveField
00995                       NULL, //weights
00996                       0.0, //rhsValue
00997                       rowSpace_.get());
00998  
00999   if (constraint != NULL) {
01000     if (*constraint != *newconstraint) {
01001       delete constraint;
01002     }
01003     else {
01004       delete newconstraint;
01005       return(0);
01006     }
01007   }
01008 
01009   lagrangeConstraints_[constraintID] = newconstraint;
01010 
01011   return(0);
01012 }
01013 
01014 //------------------------------------------------------------------------------
01015 int fei::MatrixGraph_Impl2::initPenaltyConstraint(int constraintID,
01016                                             int constraintIDType,
01017                                             int numIDs,
01018                                             const int* idTypes,
01019                                             const int* IDs,
01020                                             const int* fieldIDs)
01021 {
01022   if (rowSpace_.get() == NULL) ERReturn(-1);
01023 
01024   ConstraintType* constraint = getPenaltyConstraint(constraintID);
01025 
01026   ConstraintType* newconstraint =
01027     new ConstraintType(constraintID, constraintIDType,
01028                                               false, //isSlave
01029                                               true,  //isPenalty
01030                                               numIDs, idTypes, IDs, fieldIDs,
01031                                               0, //slaveOffset
01032                                               0, //offsetIntoSlaveField
01033                                               NULL, //weights
01034                                               0.0, //rhsValue
01035                                               rowSpace_.get());
01036 
01037   if (constraint != NULL) {
01038     if (*constraint != *newconstraint) {
01039       delete constraint;
01040     }
01041     else {
01042       delete newconstraint;
01043       return(0);
01044     }
01045   }
01046 
01047   penaltyConstraints_[constraintID] = newconstraint;
01048 
01049   return(0);
01050 }
01051 
01052 //------------------------------------------------------------------------------
01053 bool fei::MatrixGraph_Impl2::hasSlaveDof(int ID, int idType)
01054 {
01055   snl_fei::RecordCollection* collection = NULL;
01056   rowSpace_->getRecordCollection(idType, collection);
01057   if (collection == NULL) {
01058     throw std::runtime_error("fei::MatrixGraph_Impl2::hasSlaveDof: ERROR, unknown idType");
01059   }
01060 
01061   fei::Record<int>* rec = collection->getRecordWithID(ID);
01062 
01063   if (rec == NULL) {
01064     FEI_OSTRINGSTREAM osstr;
01065     osstr << "fei::MatrixGraph_Impl2::hasSlaveDof: ERROR, specified ID ("
01066      << ID << ") not found.";
01067     throw std::runtime_error(osstr.str());
01068   }
01069 
01070   return( rec->hasSlaveDof() );
01071 }
01072 
01073 //------------------------------------------------------------------------------
01074 int fei::MatrixGraph_Impl2::initSlaveConstraint(int numIDs,
01075                                           const int* idTypes,
01076                                           const int* IDs,
01077                                           const int* fieldIDs,
01078                                           int offsetOfSlave,
01079                                           int offsetIntoSlaveField,
01080                                           const double* weights,
01081                                           double rhsValue)
01082 {
01083   if (rowSpace_.get() == NULL) ERReturn(-1);
01084 
01085   FEI_OSTRINGSTREAM idosstr;
01086   idosstr << IDs[offsetOfSlave] << fieldIDs[offsetOfSlave] << offsetIntoSlaveField;
01087   FEI_ISTRINGSTREAM idisstr(idosstr.str());
01088   int crID = IDs[offsetOfSlave];
01089   idisstr >> crID;
01090 
01091   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
01092     FEI_OSTREAM& dbgos = *output_stream_;
01093 
01094     dbgos << dbgprefix_<<"initSlaveConstraint crID=" << crID << ", slaveID="
01095        << IDs[offsetOfSlave];
01096 
01097     if (output_level_ >= fei::FULL_LOGS) {
01098       dbgos << " { ";
01099       for(int ii=0; ii<numIDs; ++ii) {
01100         dbgos << "("<<IDs[ii] << ","<<fieldIDs[ii]<<") ";
01101       }
01102       dbgos << "}"<<FEI_ENDL;
01103     }
01104     else dbgos << FEI_ENDL;
01105   }
01106 
01107   ConstraintType* constraint = getSlaveConstraint(crID);
01108 
01109   ConstraintType* newconstraint =
01110     new ConstraintType(crID, 0,
01111                        true, //isSlave
01112                        false,  //isPenalty
01113                        numIDs, idTypes, IDs, fieldIDs,
01114                        offsetOfSlave,
01115                        offsetIntoSlaveField,
01116                        weights,
01117                        rhsValue,
01118                        rowSpace_.get());
01119  
01120   if (constraint != NULL) {
01121     if (*constraint != *newconstraint) {
01122       if (!constraint->structurallySame(*newconstraint)) {
01123         FEI_OSTRINGSTREAM osstr;
01124         osstr << "fei::MatrixGraph_Impl2::initSlaveConstraint: slave ID "<<IDs[offsetOfSlave]
01125               << " is already constrained, with different connectivity. Changing the"
01126               << " the structure of an existing constraint is not allowed.";
01127         throw std::runtime_error(osstr.str());
01128       }
01129       newSlaveData_ = true;
01130       delete constraint;
01131     }
01132     else {
01133       delete newconstraint;
01134       return(0);
01135     }
01136   }
01137   else {
01138     newSlaveData_ = true;
01139   }
01140 
01141   slaveConstraints_[crID] = newconstraint;
01142 
01143   return(0);
01144 }
01145 
01146 //------------------------------------------------------------------------------
01147 bool
01148 fei::MatrixGraph_Impl2::newSlaveData()
01149 {
01150   bool globalBool = false;
01151   fei::Allreduce(comm_, newSlaveData_, globalBool);
01152   newSlaveData_ = globalBool;
01153 
01154   return(newSlaveData_);
01155 }
01156 
01157 //------------------------------------------------------------------------------
01158 fei::ConstraintType*
01159 fei::MatrixGraph_Impl2::getLagrangeConstraint(int constraintID)
01160 {
01161   std::map<int,ConstraintType*>::iterator
01162     c_iter = lagrangeConstraints_.find(constraintID);
01163   if (c_iter == lagrangeConstraints_.end()) return(NULL);
01164 
01165   return( (*c_iter).second );
01166 }
01167 
01168 //------------------------------------------------------------------------------
01169 fei::ConstraintType*
01170 fei::MatrixGraph_Impl2::getSlaveConstraint(int constraintID)
01171 {
01172   std::map<int,ConstraintType*>::iterator
01173     c_iter = slaveConstraints_.find(constraintID);
01174   if (c_iter == slaveConstraints_.end()) return(NULL);
01175 
01176   return( (*c_iter).second );
01177 }
01178 
01179 //------------------------------------------------------------------------------
01180 fei::ConstraintType*
01181 fei::MatrixGraph_Impl2::getPenaltyConstraint(int constraintID)
01182 {
01183   std::map<int,ConstraintType*>::iterator
01184     c_iter = penaltyConstraints_.find(constraintID);
01185   if (c_iter == penaltyConstraints_.end()) return(NULL);
01186 
01187   return( (*c_iter).second );
01188 }
01189 
01190 //------------------------------------------------------------------------------
01191 int fei::MatrixGraph_Impl2::initComplete()
01192 {
01193   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
01194     (*output_stream_) <<dbgprefix_<< "initComplete" << FEI_ENDL;
01195   }
01196   if (rowSpace_.get() == NULL) {
01197     ERReturn(-1);
01198   }
01199   else {
01200     CHK_ERR( rowSpace_->initComplete() );
01201   }
01202 
01203   if (haveColSpace_ && colSpace_.get() == NULL) ERReturn(-1);
01204   if (haveColSpace_) {
01205     CHK_ERR( colSpace_->initComplete() );
01206   }
01207 
01208   if (rowSpace_->fieldMasks_.size() == 1 &&
01209       rowSpace_->fieldMasks_[0]->getNumFields() == 1) {
01210     simpleProblem_ = true;
01211   }
01212 
01213   std::vector<int>& eqnNums = rowSpace_->getEqnNumbers();
01214   vspcEqnPtr_ = eqnNums.size() > 0 ? &eqnNums[0] : NULL;
01215 
01216   //If there are slave constraints (on any processor), we need to create
01217   //a slave dependency matrix.
01218   localNumSlaves_ = slaveConstraints_.size();
01219   globalNumSlaves_ = 0;
01220   CHK_ERR( fei::GlobalSum(comm_, localNumSlaves_, globalNumSlaves_) );
01221 
01222   if (globalNumSlaves_ > 0) {
01223     //we need to allocate the slave dependency and reduction matrices too.
01224     CHK_ERR( createSlaveMatrices() );
01225   }
01226 
01227   return(0);
01228 }
01229 
01230 //----------------------------------------------------------------------------
01231 int fei::MatrixGraph_Impl2::createAlgebraicGraph(bool blockEntryGraph,
01232                                                  fei::Graph* graph,
01233                                                  bool gatherFromOverlap)
01234 {
01235   //Now run the connectivity blocks, adding vertices (locations of nonzeros) to
01236   //the graph.
01237 
01238   bool wasAlreadyBlockEntry = blockEntryGraph_;
01239 
01240   blockEntryGraph_ = blockEntryGraph;
01241 
01242   for(unsigned i=0; i<sparseBlocks_.size(); ++i) {
01243     fei::ConnectivityBlock& cblock = *(sparseBlocks_[i]);
01244     CHK_ERR( addBlockToGraph_sparse(graph, &cblock) );
01245   }
01246 
01247   std::map<int,fei::ConnectivityBlock*>::const_iterator
01248     cb_iter = connectivityBlocks_.begin(),
01249     cb_end  = connectivityBlocks_.end();
01250 
01251   while(cb_iter != cb_end) {
01252     fei::ConnectivityBlock& cblock = *((*cb_iter).second);
01253 
01254     bool symmetricBlock = cblock.isSymmetric();
01255 
01256     if (symmetricBlock) {
01257       fei::Pattern* pattern = cblock.getRowPattern();
01258       if (pattern == NULL) {
01259         ERReturn(-1);
01260       }
01261 
01262       fei::Pattern::PatternType pType = pattern->getPatternType();
01263       if (pType == fei::Pattern::GENERAL || pType == fei::Pattern::SINGLE_IDTYPE) {
01264         CHK_ERR( addBlockToGraph_multiField_symmetric(graph, &cblock) );
01265       }
01266       else if (pType == fei::Pattern::SIMPLE) {
01267         CHK_ERR( addBlockToGraph_singleField_symmetric(graph, &cblock) );
01268       }
01269       else if (pType == fei::Pattern::NO_FIELD) {
01270         CHK_ERR( addBlockToGraph_noField_symmetric(graph, &cblock) );
01271       }
01272     }
01273     else {
01274       fei::Pattern* pattern = cblock.getRowPattern();
01275       fei::Pattern* colpattern = cblock.getColPattern();
01276       if (pattern == NULL || colpattern == NULL) {
01277         ERReturn(-1);
01278       }
01279 
01280       fei::Pattern::PatternType pType = pattern->getPatternType();
01281       fei::Pattern::PatternType colPType = colpattern->getPatternType();
01282       if (pType == fei::Pattern::SIMPLE && colPType == fei::Pattern::SIMPLE) {
01283         CHK_ERR( addBlockToGraph_singleField_nonsymmetric(graph, &cblock) );
01284       }
01285       else {
01286         CHK_ERR( addBlockToGraph_multiField_nonsymmetric(graph, &cblock) );
01287       }
01288     }
01289     ++cb_iter;
01290   }
01291 
01292   CHK_ERR( addLagrangeConstraintsToGraph(graph) );
01293 
01294   CHK_ERR( addPenaltyConstraintsToGraph(graph) );
01295 
01296   if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
01297     FEI_OSTREAM& os = *output_stream_;
01298     os << dbgprefix_<<"# before graph->gatherFromOverlap()" << FEI_ENDL;
01299     CHK_ERR( graph->writeLocalGraph(os) );
01300     CHK_ERR( graph->writeRemoteGraph(os) );
01301   }
01302 
01303   if (gatherFromOverlap) {
01304     CHK_ERR( graph->gatherFromOverlap() );
01305   }
01306 
01307   if (blockEntryGraph) {
01308     CHK_ERR( exchangeBlkEqnSizes(graph) );
01309   }
01310 
01311   if (output_level_ >= fei::FULL_LOGS && fei::numProcs(comm_)>1 &&
01312       output_stream_ != NULL && gatherFromOverlap) {
01313     FEI_OSTREAM& os = *output_stream_;
01314     os << dbgprefix_<<" after graph->gatherFromOverlap()" << FEI_ENDL;
01315     CHK_ERR( graph->writeLocalGraph(*output_stream_) );
01316   }
01317 
01318   if (!wasAlreadyBlockEntry) {
01319     setIndicesMode(POINT_ENTRY_GRAPH);
01320   }
01321 
01322   return(0);
01323 }
01324 
01325 //----------------------------------------------------------------------------
01326 fei::SharedPtr<fei::SparseRowGraph>
01327 fei::MatrixGraph_Impl2::createGraph(bool blockEntryGraph,
01328                                     bool localRowGraph_includeSharedRows)
01329 {
01330   fei::SharedPtr<fei::SparseRowGraph> localRows;
01331 
01332   std::vector<int> globalOffsets;
01333 
01334   if (blockEntryGraph) {
01335     if (reducer_.get() != NULL) {
01336       throw std::runtime_error("fei::MatrixGraph_Impl2::createGraph ERROR, can't specify both block-entry assembly and slave-constraint reduction.");
01337     }
01338 
01339     rowSpace_->getGlobalBlkIndexOffsets(globalOffsets);
01340   }
01341   else {
01342     rowSpace_->getGlobalIndexOffsets(globalOffsets);
01343   }
01344 
01345   if ((int)globalOffsets.size() < localProc_+2) return localRows;
01346 
01347   int firstOffset = globalOffsets[localProc_];
01348   int lastOffset = globalOffsets[localProc_+1] - 1;
01349 
01350   if (reducer_.get() != NULL) {
01351     std::vector<int>& reduced_eqns = reducer_->getLocalReducedEqns();
01352     firstOffset = reduced_eqns[0];
01353     lastOffset = reduced_eqns[reduced_eqns.size()-1];
01354   }
01355 
01356   fei::SharedPtr<fei::Graph> inner_graph(new fei::Graph_Impl(comm_, firstOffset, lastOffset) );
01357   fei::SharedPtr<fei::Graph> graph;
01358 
01359   if (reducer_.get() == NULL) {
01360     graph = inner_graph;
01361   }
01362   else {
01363     graph.reset( new fei::GraphReducer(reducer_, inner_graph) );
01364   }
01365 
01366   bool gatherFromOverlap = !localRowGraph_includeSharedRows;
01367   int err = createAlgebraicGraph(blockEntryGraph, graph.get(),
01368                                  gatherFromOverlap);
01369   if (err != 0) {
01370     return(localRows);
01371   }
01372 
01373   localRows = fei::createSparseRowGraph(*(inner_graph->getLocalGraph()));
01374 
01375   remotelyOwnedGraphRows_ = fei::createSparseRowGraph(inner_graph->getRemoteGraph());
01376 
01377   remotelyOwnedGraphRows_->blockEntries = blockEntryGraph;
01378 
01379   if (localRowGraph_includeSharedRows &&
01380       remotelyOwnedGraphRows_->rowNumbers.size() > 0) {
01381 
01382     fei::SharedPtr<fei::SparseRowGraph> localPlusShared =
01383       snl_fei::mergeSparseRowGraphs(localRows.get(), remotelyOwnedGraphRows_.get());
01384     localRows = localPlusShared;
01385   }
01386 
01387   return(localRows);
01388 }
01389 
01390 //----------------------------------------------------------------------------
01391 int fei::MatrixGraph_Impl2::exchangeBlkEqnSizes(fei::Graph* graph)
01392 {
01393   //If point-equals-block (meaning block-equations are of size 1) or
01394   //if blockEntryGraph_ is false (meaning we haven't constructed a block-entry-
01395   //graph) then there is nothing to do in this function.
01396   if ( rowSpace_->getPointBlockMap()->ptEqualBlk() ) {
01397     return(0);
01398   }
01399 
01400   snl_fei::BlkSizeMsgHandler blkHandler(rowSpace_.get(), graph, comm_);
01401 
01402   CHK_ERR( blkHandler.do_the_exchange() );
01403 
01404   return(0);
01405 }
01406 
01407 //----------------------------------------------------------------------------
01408 int fei::MatrixGraph_Impl2::createSlaveMatrices()
01409 {
01410   //In this function we extract the dependency matrix of equation numbers 
01411   //from the slave-constraints locally on each processor, then gather those
01412   //slave-equations onto each processor (so that each processor holds ALL
01413   //slave-equations).
01414   //Then we remove any couplings that may exist among the slave-equations and
01415   //finally create a FillableMat object to hold the final dependency matrix D_.
01416   //
01417 
01418   if (!newSlaveData()) {
01419     return(0);
01420   }
01421 
01422   std::vector<int>& eqnNums = rowSpace_->getEqnNumbers();
01423   vspcEqnPtr_ = eqnNums.size() > 0 ? &eqnNums[0] : NULL;
01424 
01425   fei::SharedPtr<fei::FillableMat> local_D(new fei::FillableMat);
01426   fei::SharedPtr<fei::CSVec> local_g(new fei::CSVec);
01427 
01428   std::vector<int> masterEqns;
01429   std::vector<double> masterCoefs;
01430 
01431   std::map<int,ConstraintType*>::const_iterator
01432     cr_iter = slaveConstraints_.begin(),
01433     cr_end  = slaveConstraints_.end();
01434 
01435   for(; cr_iter != cr_end; ++cr_iter) {
01436     ConstraintType* cr = (*cr_iter).second;
01437 
01438     fei::Record<int>* slaveRecord = cr->getSlave();
01439     int slaveID = slaveRecord->getID();
01440     int slaveIDType = cr->getIDType();
01441     int slaveFieldID = cr->getSlaveFieldID();
01442     int offsetIntoSlaveField = cr->getOffsetIntoSlaveField();
01443     int slaveEqn = -1;
01444     CHK_ERR( rowSpace_->getGlobalIndex(slaveIDType, slaveID,
01445                                        slaveFieldID, 0, offsetIntoSlaveField,
01446                                        slaveEqn) );
01447 
01448     std::vector<int>& masterRecords_vec = cr->getMasters();
01449     int* masterRecords = &masterRecords_vec[0];
01450     std::vector<snl_fei::RecordCollection*>& masterRecColls = cr->getMasterRecordCollections();
01451     std::vector<int>& masterIDTypes = cr->getMasterIDTypes();
01452     std::vector<int>& masterFieldIDs = cr->getMasterFieldIDs();
01453     std::vector<double>& masterWeights = cr->getMasterWeights();
01454     double* masterWtPtr = &masterWeights[0];
01455 
01456     masterEqns.resize(masterWeights.size());
01457     masterCoefs.resize(masterWeights.size());
01458 
01459     int* masterEqnsPtr = &(masterEqns[0]);
01460     double* masterCoefsPtr = &(masterCoefs[0]);
01461 
01462     int offset = 0;
01463     for(size_t j=0; j<masterIDTypes.size(); ++j) {
01464       fei::Record<int>* masterRecord = masterRecColls[j]->getRecordWithLocalID(masterRecords[j]);
01465       int* eqnNumbers = vspcEqnPtr_+masterRecord->getOffsetIntoEqnNumbers();
01466       fei::FieldMask* mask = masterRecord->getFieldMask();
01467       int eqnOffset = 0;
01468       if (!simpleProblem_) {
01469         mask->getFieldEqnOffset(masterFieldIDs[j], eqnOffset);
01470       }
01471 
01472       unsigned fieldSize = rowSpace_->getFieldSize(masterFieldIDs[j]);
01473 
01474       int eqn = eqnNumbers[eqnOffset];
01475       for(unsigned k=0; k<fieldSize; ++k) {
01476         masterEqnsPtr[offset++] = eqn+k;
01477       }
01478     }
01479 
01480     double fei_eps = 1.e-49;
01481 
01482     offset = 0;
01483     for(size_t jj=0; jj<masterEqns.size(); ++jj) {
01484       if (std::abs(masterWtPtr[jj]) > fei_eps) {
01485         masterCoefsPtr[offset] = masterWtPtr[jj];
01486         masterEqnsPtr[offset++] = masterEqnsPtr[jj];
01487       }
01488     }
01489 
01490     if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
01491       bool log_eqn = isLogEqn(slaveEqn);
01492       if (!log_eqn) {
01493         for(unsigned ii=0; ii<masterEqns.size(); ++ii) {
01494           if (isLogEqn(masterEqnsPtr[ii])) {
01495             log_eqn = true;
01496             break;
01497           }
01498         }
01499       }
01500 
01501       if (log_eqn) {
01502         FEI_OSTREAM& os = *output_stream_;
01503         os << "createSlaveMatrices: " << slaveEqn << " = ";
01504         for(unsigned ii=0; ii<masterEqns.size(); ++ii) {
01505           if (ii!=0) os << "+ ";
01506           os << masterCoefsPtr[ii]<<"*"<<masterEqnsPtr[ii]<<" ";
01507         }
01508         os << FEI_ENDL;
01509       }
01510     }
01511 
01512     local_D->putRow(slaveEqn, masterEqnsPtr, masterCoefsPtr, offset);
01513 
01514     if (std::abs(cr->getRHSValue()) > fei_eps) {
01515       fei::put_entry(*local_g, slaveEqn, cr->getRHSValue());
01516     }
01517   }
01518 
01519   if (D_.get() == NULL) {
01520     D_.reset(new fei::FillableMat);
01521   }
01522 
01523   if (g_.get() == NULL) {
01524     g_.reset(new fei::CSVec);
01525   }
01526 
01527 #ifndef FEI_SER
01528   if (numProcs_ > 1) {
01529     fei::impl_utils::global_union(comm_, *local_D, *D_);
01530     fei::impl_utils::global_union(comm_, *local_g, *g_);
01531   }
01532   else {
01533     *D_ = *local_D;
01534     *g_ = *local_g;
01535   }
01536 #else
01537   *D_ = *local_D;
01538   *g_ = *local_g;
01539 #endif
01540 
01541   if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
01542     (*output_stream_) << "#  D_ (pre-removeCouplings):"<<FEI_ENDL;
01543     (*output_stream_) << *D_;
01544   }
01545 
01546   int levelsOfCoupling = fei::impl_utils::remove_couplings(*D_);
01547   (void)levelsOfCoupling;
01548 
01549   if (reducer_.get() == NULL) {
01550     reducer_.reset(new fei::Reducer(D_, g_, comm_));
01551 
01552     std::vector<int> indices;
01553     rowSpace_->getIndices_Owned(indices);
01554 
01555     reducer_->setLocalUnreducedEqns(indices);
01556   }
01557   else {
01558     reducer_->initialize();
01559   }
01560 
01561   if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
01562     (*output_stream_) << "#  D_:"<<FEI_ENDL;
01563     (*output_stream_) << *D_;
01564   }
01565 
01566   double fei_eps = 1.e-49;
01567 
01568   g_nonzero_ = false;
01569   for(size_t j=0; j<g_->size(); ++j) {
01570     double coef = g_->coefs()[j];
01571     if (std::abs(coef) > fei_eps) {
01572       g_nonzero_ = true;
01573     }
01574   }
01575 
01576   newSlaveData_ = false;
01577 
01578   return(0);
01579 }
01580 
01581 //----------------------------------------------------------------------------
01582 fei::SharedPtr<fei::Reducer>
01583 fei::MatrixGraph_Impl2::getReducer()
01584 {
01585   return( reducer_ );
01586 }
01587 
01588 //----------------------------------------------------------------------------
01589 fei::SharedPtr<fei::SparseRowGraph>
01590 fei::MatrixGraph_Impl2::getRemotelyOwnedGraphRows()
01591 {
01592   return( remotelyOwnedGraphRows_ );
01593 }
01594 
01595 //----------------------------------------------------------------------------
01596 void fei::MatrixGraph_Impl2::getConstrainedIndices(std::vector<int>& crindices) const
01597 {
01598   if (constrained_indices_.empty()) {
01599     crindices.clear();
01600     return;
01601   }
01602 
01603   std::set<int>::const_iterator
01604     s_iter = constrained_indices_.begin(),
01605     s_end = constrained_indices_.end();
01606 
01607   crindices.resize(constrained_indices_.size());
01608 
01609   int offset = 0;
01610   for(; s_iter != s_end; ++s_iter) {
01611     crindices[offset++] = *s_iter;
01612   }
01613 }
01614 
01615 //----------------------------------------------------------------------------
01616 int fei::MatrixGraph_Impl2::addLagrangeConstraintsToGraph(fei::Graph* graph)
01617 {
01618   std::vector<int> indices;
01619   std::map<int,ConstraintType*>::const_iterator
01620     cr_iter = lagrangeConstraints_.begin(),
01621     cr_end  = lagrangeConstraints_.end();
01622 
01623   constrained_indices_.clear();
01624 
01625   while(cr_iter != cr_end) {
01626     ConstraintType* cr = (*cr_iter).second;
01627     int crID = cr->getConstraintID();
01628 
01629     CHK_ERR( getConstraintConnectivityIndices(cr, indices) );
01630 
01631     int numIndices = indices.size();
01632     int* indicesPtr = &(indices[0]);
01633 
01634     for(int i=0; i<numIndices; ++i) {
01635       constrained_indices_.insert(indicesPtr[i]);
01636     }
01637 
01638     int crEqnRow = -1, tmp;
01639     if (blockEntryGraph_) {
01640       CHK_ERR( rowSpace_->getGlobalBlkIndex(cr->getIDType(),
01641                                                    crID, crEqnRow) );
01642       cr->setBlkEqnNumber(crEqnRow);
01643       CHK_ERR( rowSpace_->getGlobalIndex(cr->getIDType(),
01644                                                 crID, tmp) );
01645       cr->setEqnNumber(tmp);
01646     }
01647     else {
01648       CHK_ERR( rowSpace_->getGlobalIndex(cr->getIDType(),
01649                                                 crID, crEqnRow) );
01650       cr->setEqnNumber(crEqnRow);
01651       CHK_ERR( rowSpace_->getGlobalBlkIndex(cr->getIDType(),
01652                                                    crID, tmp) );
01653       cr->setBlkEqnNumber(tmp);
01654     }
01655 
01656     //now add the row contribution
01657     CHK_ERR( graph->addIndices(crEqnRow, numIndices, &(indices[0])) );
01658 
01659     //Let's add a diagonal entry to the graph for this constraint-equation,
01660     //just in case we need to fiddle with this equation later (e.g. discard the
01661     //constraint equation and replace it with a dirichlet boundary condition...).
01662     CHK_ERR( graph->addIndices(crEqnRow, 1, &crEqnRow) );
01663 
01664     //and finally, add the column contribution (which is simply the transpose
01665     //of the row contribution)
01666     for(int k=0; k<numIndices; ++k) {
01667       CHK_ERR( graph->addIndices(indicesPtr[k], 1, &crEqnRow) );
01668     }
01669 
01670     ++cr_iter;
01671   }
01672 
01673   return(0);
01674 }
01675 
01676 //----------------------------------------------------------------------------
01677 int fei::MatrixGraph_Impl2::
01678 getConstraintConnectivityIndices(ConstraintType* cr,
01679                                  std::vector<int>& globalIndices)
01680 {
01681   std::vector<int>& fieldSizes = tmpIntArray1_;
01682   std::vector<int>& ones = tmpIntArray2_;
01683 
01684   std::vector<int>& constrainedRecords = cr->getMasters();
01685   std::vector<int>& constrainedFieldIDs = cr->getMasterFieldIDs();
01686   std::vector<snl_fei::RecordCollection*>& recordCollections = cr->getMasterRecordCollections();
01687 
01688   int len = constrainedRecords.size();
01689   fieldSizes.resize(len);
01690 
01691   ones.assign(len, 1);
01692 
01693   int numIndices = 0;
01694   if (blockEntryGraph_) {
01695     numIndices = len;
01696   }
01697   else {
01698     for(int j=0; j<len; ++j) {
01699       unsigned fieldSize = rowSpace_->getFieldSize(constrainedFieldIDs[j]);
01700       fieldSizes[j] = fieldSize;
01701       numIndices += fieldSize;
01702     }
01703   }
01704 
01705   globalIndices.resize(numIndices);
01706 
01707   int checkNum;
01708   CHK_ERR( getConnectivityIndices_multiField(&recordCollections[0],
01709                                              &constrainedRecords[0],
01710                                              len, &ones[0],
01711                                              &constrainedFieldIDs[0],
01712                                              &fieldSizes[0],
01713                                              numIndices, &globalIndices[0],
01714                                              checkNum) );
01715   if (numIndices != checkNum) {
01716     ERReturn(-1);
01717   }
01718 
01719   return(0);
01720 }
01721 
01722 //----------------------------------------------------------------------------
01723 int fei::MatrixGraph_Impl2::addPenaltyConstraintsToGraph(fei::Graph* graph)
01724 {
01725   std::vector<int> indices;
01726   std::map<int,ConstraintType*>::const_iterator
01727     cr_iter = penaltyConstraints_.begin(),
01728     cr_end  = penaltyConstraints_.end();
01729 
01730   while(cr_iter != cr_end) {
01731     ConstraintType* cr = (*cr_iter).second;
01732 
01733     CHK_ERR( getConstraintConnectivityIndices(cr, indices) );
01734 
01735     int numIndices = indices.size();
01736 
01737     //now add the symmetric contributions to the graph
01738     CHK_ERR( graph->addSymmetricIndices(numIndices, &(indices[0])) );
01739 
01740     ++cr_iter;
01741   }
01742 
01743   return(0);
01744 }
01745 
01746 //----------------------------------------------------------------------------
01747 int fei::MatrixGraph_Impl2::compareStructure(const fei::MatrixGraph& matrixGraph,
01748                                        bool& equivalent) const
01749 {
01750   equivalent = false;
01751   int localCode = 1;
01752   int globalCode = 0;
01753 
01754   int myNumBlocks = getNumConnectivityBlocks();
01755   int numBlocks = matrixGraph.getNumConnectivityBlocks();
01756 
01757   if (myNumBlocks != numBlocks) {
01758     CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
01759     equivalent = globalCode > 0 ? false : true;
01760     return(0);
01761   }
01762 
01763   if (numBlocks > 0) {
01764     std::vector<int> myBlockIDs;
01765     std::vector<int> blockIDs;
01766 
01767     CHK_ERR( getConnectivityBlockIDs(myBlockIDs) );
01768     CHK_ERR( matrixGraph.getConnectivityBlockIDs(blockIDs) );
01769 
01770     for(int i=0; i<numBlocks; ++i) {
01771       if (myBlockIDs[i] != blockIDs[i]) {
01772         CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
01773         equivalent = globalCode > 0 ? false : true;
01774         return(0);
01775       }
01776 
01777       const fei::ConnectivityBlock* mycblock = getConnectivityBlock(myBlockIDs[i]);
01778       const fei::ConnectivityBlock* cblock = matrixGraph.getConnectivityBlock(blockIDs[i]);
01779 
01780       int myNumLists = mycblock->getNativeConnectivityIDs().size();
01781       int numLists = cblock->getNativeConnectivityIDs().size();
01782 
01783       if (myNumLists != numLists ||
01784           mycblock->isSymmetric() != cblock->isSymmetric()) {
01785         CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
01786         equivalent = globalCode > 0 ? false : true;
01787         return(0);
01788       }
01789 
01790       int myNumIDsPerList = mycblock->getRowPattern()->getNumIDs();
01791       int numIDsPerList = cblock->getRowPattern()->getNumIDs();
01792 
01793       if (myNumIDsPerList != numIDsPerList) {
01794         CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
01795         equivalent = globalCode > 0 ? false : true;
01796         return(0);
01797       }
01798     }
01799   }
01800 
01801   int numMyLagrangeConstraints = getLocalNumLagrangeConstraints();
01802   int numMySlaveConstraints = getGlobalNumSlaveConstraints();
01803   int numLagrangeConstraints = matrixGraph.getLocalNumLagrangeConstraints();
01804   int numSlaveConstraints = matrixGraph.getGlobalNumSlaveConstraints();
01805 
01806   if (numMyLagrangeConstraints != numLagrangeConstraints ||
01807       numMySlaveConstraints != numSlaveConstraints) {
01808     CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
01809     equivalent = globalCode > 0 ? false : true;
01810     return(0);
01811   }
01812 
01813   localCode = 0;
01814   CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
01815   equivalent = globalCode > 0 ? false : true;
01816   return(0);
01817 }
01818 
01819 //----------------------------------------------------------------------------
01820 int fei::MatrixGraph_Impl2::getNumConnectivityBlocks() const
01821 {
01822   return(connectivityBlocks_.size());
01823 }
01824 
01825 //----------------------------------------------------------------------------
01826 int fei::MatrixGraph_Impl2::getConnectivityBlockIDs(std::vector<int>& blockIDs) const
01827 {
01828   blockIDs.resize(connectivityBlocks_.size());
01829 
01830   std::map<int,fei::ConnectivityBlock*>::const_iterator
01831     cdb_iter = connectivityBlocks_.begin();
01832 
01833   for(size_t i=0; i<blockIDs.size(); ++i, ++cdb_iter) {
01834     int blockID = (*cdb_iter).first;
01835     blockIDs[i] = blockID;
01836   }
01837 
01838   return(0);
01839 }
01840 
01841 //----------------------------------------------------------------------------
01842 int fei::MatrixGraph_Impl2::getNumIDsPerConnectivityList(int blockID) const
01843 {
01844   const fei::ConnectivityBlock* cblock = getConnectivityBlock(blockID);
01845   if (cblock == NULL) return(-1);
01846 
01847   const fei::Pattern* pattern = cblock->getRowPattern();
01848   return(pattern->getNumIDs());
01849 }
01850 
01851 //----------------------------------------------------------------------------
01852 int fei::MatrixGraph_Impl2::getConnectivityNumIndices(int blockID) const
01853 {
01854   const fei::ConnectivityBlock* cblock = getConnectivityBlock(blockID);
01855   if (cblock == NULL) return(-1);
01856 
01857   const fei::Pattern* pattern = cblock->getRowPattern();
01858   return( blockEntryGraph_ ?
01859     pattern->getNumIDs() : pattern->getNumIndices());
01860 }
01861 
01862 //----------------------------------------------------------------------------
01863 int fei::MatrixGraph_Impl2::getConnectivityNumIndices(int blockID,
01864                                                   int& numRowIndices,
01865                                                     int& numColIndices)
01866 {
01867   fei::ConnectivityBlock* cblock = getConnectivityBlock(blockID);
01868   if (cblock == NULL) return(-1);
01869 
01870   fei::Pattern* pattern = cblock->getRowPattern();
01871   numRowIndices = pattern->getNumIndices();
01872   fei::Pattern* colpattern = cblock->isSymmetric() ?
01873     pattern : cblock->getColPattern();
01874   numColIndices = colpattern->getNumIndices();
01875 
01876   return(0);
01877 }
01878 
01879 //----------------------------------------------------------------------------
01880 const fei::ConnectivityBlock* fei::MatrixGraph_Impl2::getConnectivityBlock(int blockID) const
01881 {
01882   std::map<int,fei::ConnectivityBlock*>::const_iterator
01883     c_iter = connectivityBlocks_.find(blockID);
01884   if (c_iter == connectivityBlocks_.end()) return(0);
01885 
01886   return( (*c_iter).second );
01887 }
01888 
01889 //----------------------------------------------------------------------------
01890 fei::ConnectivityBlock* fei::MatrixGraph_Impl2::getConnectivityBlock(int blockID)
01891 {
01892   std::map<int,fei::ConnectivityBlock*>::iterator
01893     c_iter = connectivityBlocks_.find(blockID);
01894   if (c_iter == connectivityBlocks_.end()) return(0);
01895 
01896   return( (*c_iter).second );
01897 }
01898 
01899 //----------------------------------------------------------------------------
01900 int fei::MatrixGraph_Impl2::getConnectivityIndices(int blockID,
01901                                              int connectivityID,
01902                                              int indicesAllocLen,
01903                                              int* indices,
01904                                              int& numIndices)
01905 {
01906   fei::ConnectivityBlock* cblock = getConnectivityBlock(blockID);
01907   if (cblock == NULL) return(-1);
01908 
01909   std::vector<int>& eqnNums = rowSpace_->getEqnNumbers();
01910   vspcEqnPtr_ = eqnNums.size() > 0 ? &eqnNums[0] : NULL;
01911 
01912   fei::Pattern* pattern = cblock->getRowPattern();
01913   numIndices = pattern->getNumIndices();
01914 
01915   int len = numIndices > indicesAllocLen ? indicesAllocLen : numIndices;
01916 
01917   int* records = cblock->getRowConnectivity(connectivityID);
01918   if (records == NULL) {
01919     ERReturn(-1);
01920   }
01921 
01922   fei::Pattern::PatternType pType = pattern->getPatternType();
01923 
01924   if (pType == fei::Pattern::GENERAL || pType == fei::Pattern::SINGLE_IDTYPE) {
01925     const int* numFieldsPerID = pattern->getNumFieldsPerID();
01926     const int* fieldIDs = pattern->getFieldIDs();
01927     int totalNumFields = pattern->getTotalNumFields();
01928 
01929     std::vector<int> fieldSizes(totalNumFields);
01930 
01931     for(int j=0; j<totalNumFields; ++j) {
01932       fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
01933                                             colSpace_.get());
01934     }
01935 
01936     CHK_ERR( getConnectivityIndices_multiField(pattern->getRecordCollections(),
01937                                             records, pattern->getNumIDs(),
01938                                                numFieldsPerID,
01939                                                fieldIDs, &fieldSizes[0],
01940                                                len, indices, numIndices) );
01941   }
01942   else if (pType == fei::Pattern::SIMPLE) {
01943     const int* fieldIDs = pattern->getFieldIDs();
01944 
01945     int fieldID = fieldIDs[0];
01946     unsigned fieldSize = snl_fei::getFieldSize(fieldID,
01947                                                rowSpace_.get(),
01948                                                colSpace_.get());
01949 
01950     CHK_ERR( getConnectivityIndices_singleField(pattern->getRecordCollections(),
01951                                             records, pattern->getNumIDs(),
01952                                                 fieldID, fieldSize,
01953                                                 len, indices, numIndices) );
01954   }
01955   else if (pType == fei::Pattern::NO_FIELD) {
01956     CHK_ERR( getConnectivityIndices_noField(pattern->getRecordCollections(),
01957                                             records, pattern->getNumIDs(),
01958                                             len, indices, numIndices) );
01959   }
01960 
01961   return(0);
01962 }
01963 
01964 //----------------------------------------------------------------------------
01965 int fei::MatrixGraph_Impl2::getConnectivityIndices(int blockID,
01966                                              int connectivityID,
01967                                              int rowIndicesAllocLen,
01968                                              int* rowIndices,
01969                                              int& numRowIndices,
01970                                              int colIndicesAllocLen,
01971                                              int* colIndices,
01972                                              int& numColIndices)
01973 {
01974   fei::ConnectivityBlock* cblock = getConnectivityBlock(blockID);
01975   if (cblock == NULL) return(-1);
01976 
01977   std::vector<int>& eqnNums = rowSpace_->getEqnNumbers();
01978   vspcEqnPtr_ = eqnNums.size() > 0 ? &eqnNums[0] : NULL;
01979 
01980   fei::Pattern* pattern = cblock->getRowPattern();
01981   numRowIndices = pattern->getNumIndices();
01982 
01983   int len = numRowIndices > rowIndicesAllocLen ?
01984     rowIndicesAllocLen : numRowIndices;
01985 
01986   int* records = cblock->getRowConnectivity(connectivityID);
01987   if (records == NULL) {
01988     ERReturn(-1);
01989   }
01990 
01991   fei::Pattern::PatternType pType = pattern->getPatternType();
01992 
01993   if (pType == fei::Pattern::GENERAL || pType == fei::Pattern::SINGLE_IDTYPE) {
01994     const int* numFieldsPerID = pattern->getNumFieldsPerID();
01995     const int* fieldIDs = pattern->getFieldIDs();
01996     int totalNumFields = pattern->getTotalNumFields();
01997 
01998     std::vector<int> fieldSizes(totalNumFields);
01999 
02000     for(int j=0; j<totalNumFields; ++j) {
02001       fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
02002                                             colSpace_.get());
02003     }
02004 
02005     CHK_ERR( getConnectivityIndices_multiField(pattern->getRecordCollections(),
02006                                                records, pattern->getNumIDs(),
02007                                                numFieldsPerID,
02008                                                fieldIDs, &fieldSizes[0],
02009                                                len, rowIndices, numRowIndices) );
02010   }
02011   else if (pType == fei::Pattern::SIMPLE) {
02012     const int* fieldIDs = pattern->getFieldIDs();
02013 
02014     int fieldID = fieldIDs[0];
02015     unsigned fieldSize = snl_fei::getFieldSize(fieldID, rowSpace_.get(),
02016                                             colSpace_.get());
02017 
02018     CHK_ERR( getConnectivityIndices_singleField(pattern->getRecordCollections(),
02019                                                 records, pattern->getNumIDs(),
02020                                                 fieldID, fieldSize,
02021                                                 len, rowIndices, numRowIndices) );
02022   }
02023   else if (pType == fei::Pattern::NO_FIELD) {
02024     CHK_ERR( getConnectivityIndices_noField(pattern->getRecordCollections(),
02025                                             records, pattern->getNumIDs(),
02026                                             len, rowIndices, numRowIndices) );
02027   }
02028 
02029   fei::Pattern* colpattern = cblock->isSymmetric() ? 
02030     pattern : cblock->getColPattern();
02031   pType = colpattern->getPatternType();
02032   numColIndices = colpattern->getNumIndices();
02033   len = numColIndices > colIndicesAllocLen ?
02034     colIndicesAllocLen : numColIndices;
02035 
02036   if (!cblock->isSymmetric()) {
02037     records = cblock->getColConnectivity(connectivityID);
02038   }
02039   if (records == NULL) {
02040     return(-1);
02041   }
02042 
02043   if (pType == fei::Pattern::GENERAL || pType == fei::Pattern::SINGLE_IDTYPE) {
02044     const int* numFieldsPerID = colpattern->getNumFieldsPerID();
02045     const int* fieldIDs = colpattern->getFieldIDs();
02046     int totalNumFields = colpattern->getTotalNumFields();
02047 
02048     std::vector<int> fieldSizes(totalNumFields);
02049 
02050     for(int j=0; j<totalNumFields; ++j) {
02051       fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
02052                                             colSpace_.get());
02053     }
02054 
02055     CHK_ERR( getConnectivityIndices_multiField(colpattern->getRecordCollections(),
02056                                                records, colpattern->getNumIDs(),
02057                                                numFieldsPerID,
02058                                                fieldIDs, &fieldSizes[0],
02059                                                len, colIndices, numColIndices) );
02060   }
02061   else if (pType == fei::Pattern::SIMPLE) {
02062     const int* fieldIDs = colpattern->getFieldIDs();
02063 
02064     int fieldID = fieldIDs[0];
02065     unsigned fieldSize = snl_fei::getFieldSize(fieldID, rowSpace_.get(),
02066                                                colSpace_.get());
02067 
02068     CHK_ERR( getConnectivityIndices_singleField(colpattern->getRecordCollections(),
02069                                                 records, colpattern->getNumIDs(),
02070                                                 fieldID, fieldSize,
02071                                                 len, colIndices, numColIndices) );
02072   }
02073   else if (pType == fei::Pattern::NO_FIELD) {
02074     CHK_ERR( getConnectivityIndices_noField(colpattern->getRecordCollections(),
02075                                             records, colpattern->getNumIDs(),
02076                                             len, colIndices, numColIndices) );
02077   }
02078 
02079   return(0);
02080 }
02081 
02082 //----------------------------------------------------------------------------
02083 int fei::MatrixGraph_Impl2::getLocalNumLagrangeConstraints() const
02084 {
02085   return( lagrangeConstraints_.size() );
02086 }
02087 
02088 //----------------------------------------------------------------------------
02089 int fei::MatrixGraph_Impl2::addBlockToGraph_multiField_symmetric(fei::Graph* graph,
02090                                                         fei::ConnectivityBlock* cblock)
02091 {
02092   fei::Pattern* pattern = cblock->getRowPattern();
02093   int j;
02094   int numIDs = pattern->getNumIDs();
02095   int numIndices = blockEntryGraph_ ?
02096     pattern->getNumIDs() : pattern->getNumIndices();
02097   int checkNumIndices = numIndices;
02098   std::vector<int> indices(numIndices);
02099   int* indicesPtr = &indices[0];
02100 
02101   const int* numFieldsPerID = pattern->getNumFieldsPerID();
02102   const int* fieldIDs = pattern->getFieldIDs();
02103   int totalNumFields = pattern->getTotalNumFields();
02104 
02105   std::vector<int> fieldSizes(totalNumFields);
02106 
02107   for(j=0; j<totalNumFields; ++j) {
02108     fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
02109                                           colSpace_.get());
02110   }
02111 
02112   IndexType<int,int>& connIDs = cblock->getNativeConnectivityIDs();
02113   std::vector<int>& values = cblock->getRowConnectivities();
02114 
02115   IndexType<int,int>::iterator
02116     c_iter = connIDs.begin(),
02117     c_end  = connIDs.end();
02118 
02119   for(; c_iter != c_end; ++c_iter) {
02120     int offset = c_iter->second;
02121     int* records = &values[offset*numIDs];
02122 
02123     CHK_ERR( getConnectivityIndices_multiField(pattern->getRecordCollections(),
02124                                                 records, numIDs,
02125                                                  numFieldsPerID,
02126                                                  fieldIDs,
02127                                                  &fieldSizes[0],
02128                                                  numIndices,
02129                                                  indicesPtr,
02130                                                  checkNumIndices) );
02131 
02132     if (checkNumIndices != numIndices) {
02133       ERReturn(-1);
02134     }
02135 
02136     if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
02137       FEI_OSTREAM& os = *output_stream_;
02138 
02139       const snl_fei::RecordCollection*const* recordColls = pattern->getRecordCollections();
02140       unsigned thisoffset = 0;
02141       for(int ii=0; ii<numIDs; ++ii) {
02142         const fei::Record<int>* record = recordColls[ii]->getRecordWithLocalID(records[ii]);
02143         int ID = record->getID();
02144         os << dbgprefix_<<"scatterIndices: ID=" <<ID<<": ";
02145         int num = pattern->getNumIndicesPerID()[ii];
02146         for(int jj=0; jj<num; ++jj) {
02147           os << indicesPtr[thisoffset++] << " ";
02148         }
02149         os << FEI_ENDL;
02150       }
02151 
02152       for(int ii=0; ii<numIndices; ++ii) {
02153         if (isLogEqn(indicesPtr[ii])) {
02154           os << "adding Symm inds: ";
02155           for(int jj=0; jj<numIndices; ++jj) {
02156             os << indicesPtr[jj] << " ";
02157           }
02158           os << FEI_ENDL;
02159           break;
02160         }
02161       }
02162     }
02163 
02164     //now we have the indices array, so we're ready to push it into
02165     //the graph container
02166     if (numIndices == numIDs || !cblock->isDiagonal()) {
02167       CHK_ERR( graph->addSymmetricIndices(numIndices, indicesPtr,
02168                                           cblock->isDiagonal()) );
02169     }
02170     else {
02171       int ioffset = 0;
02172       int* fieldSizesPtr = &fieldSizes[0];
02173       for(int i=0; i<numIDs; ++i) {
02174         CHK_ERR( graph->addSymmetricIndices(fieldSizesPtr[i], &(indicesPtr[ioffset])));
02175         ioffset += fieldSizesPtr[i];
02176       }
02177     }
02178   }
02179 
02180   return(0);
02181 }
02182 
02183 //----------------------------------------------------------------------------
02184 int fei::MatrixGraph_Impl2::addBlockToGraph_multiField_nonsymmetric(fei::Graph* graph,
02185                                                         fei::ConnectivityBlock* cblock)
02186 {
02187   fei::Pattern* pattern = cblock->getRowPattern();
02188   fei::Pattern* colpattern = cblock->getColPattern();
02189   int j;
02190   int numIDs = pattern->getNumIDs();
02191   int numIndices = blockEntryGraph_ ? numIDs : pattern->getNumIndices();
02192   int checkNumIndices = numIndices;
02193   std::vector<int> indices(numIndices);
02194   int* indicesPtr = &indices[0];
02195 
02196   int numColIDs = colpattern->getNumIDs();
02197   int numColIndices = blockEntryGraph_ ? numColIDs : colpattern->getNumIndices();
02198   int checkNumColIndices = numColIndices;
02199   std::vector<int> colindices(numColIndices);
02200   int* colindicesPtr = &colindices[0];
02201 
02202   const int* numFieldsPerID = pattern->getNumFieldsPerID();
02203   const int* fieldIDs = pattern->getFieldIDs();
02204   int totalNumFields = pattern->getTotalNumFields();
02205 
02206   const int* numFieldsPerColID = colpattern->getNumFieldsPerID();
02207   const int* colfieldIDs = colpattern->getFieldIDs();
02208   int totalNumColFields = colpattern->getTotalNumFields();
02209 
02210   std::vector<int> fieldSizes(totalNumFields);
02211   std::vector<int> colfieldSizes(totalNumColFields);
02212 
02213   for(j=0; j<totalNumFields; ++j) {
02214     fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
02215                                           colSpace_.get());
02216   }
02217   for(j=0; j<totalNumColFields; ++j) {
02218     colfieldSizes[j] = snl_fei::getFieldSize(colfieldIDs[j], rowSpace_.get(),
02219                                            colSpace_.get());
02220   }
02221 
02222   IndexType<int,int>& connIDs = cblock->getNativeConnectivityIDs();
02223   std::vector<int>& rowrecords = cblock->getRowConnectivities();
02224   std::vector<int>& colrecords = cblock->getColConnectivities();
02225 
02226   IndexType<int,int>::iterator
02227    c_iter = connIDs.begin(),
02228    c_end  = connIDs.end();
02229 
02230   for(; c_iter != c_end; ++c_iter) {
02231     int offset = c_iter->second;
02232     int* records = &rowrecords[offset*numIDs];
02233 
02234     int* colRecords = &colrecords[offset*numColIDs];
02235 
02236     CHK_ERR( getConnectivityIndices_multiField(pattern->getRecordCollections(),
02237                                                 records, numIDs,
02238                                                  numFieldsPerID,
02239                                                  fieldIDs,
02240                                                  &fieldSizes[0],
02241                                                  numIndices,
02242                                                  indicesPtr,
02243                                                  checkNumIndices) );
02244 
02245     if (checkNumIndices != numIndices) {
02246       ERReturn(-1);
02247     }
02248 
02249       
02250     CHK_ERR( getConnectivityIndices_multiField(colpattern->getRecordCollections(),
02251                                                 colRecords, numColIDs,
02252                                                  numFieldsPerColID,
02253                                                  colfieldIDs,
02254                                                  &colfieldSizes[0],
02255                                                  numColIndices,
02256                                                  colindicesPtr,
02257                                                  checkNumColIndices) );
02258 
02259     if (checkNumColIndices != numColIndices) {
02260       ERReturn(-1);
02261     }
02262 
02263     //now we have the indices arrays, so we're ready to push them into
02264     //the graph container
02265     for(int r=0; r<numIndices; ++r) {
02266       CHK_ERR( graph->addIndices(indicesPtr[r], numColIndices, colindicesPtr));
02267     }
02268   }
02269 
02270   return(0);
02271 }
02272 
02273 //----------------------------------------------------------------------------
02274 int fei::MatrixGraph_Impl2::getConnectivityIndices_multiField(const snl_fei::RecordCollection*const* recordCollections,
02275                                                           int* records,
02276                                                           int numRecords,
02277                                                           const int* numFieldsPerID,
02278                                                           const int* fieldIDs,
02279                                                           const int* fieldSizes,
02280                                                           int indicesAllocLen,
02281                                                           int* indices,
02282                                                           int& numIndices)
02283 {
02284   numIndices = 0;
02285   int fld_offset = 0;
02286 
02287   for(int i=0; i<numRecords; ++i) {
02288     const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
02289     if (record==NULL) continue;
02290 
02291     if (blockEntryGraph_) {
02292       indices[numIndices++] = record->getNumber();
02293       continue;
02294     }
02295 
02296     const fei::FieldMask* fieldMask = record->getFieldMask();
02297     int* eqnNumbers = vspcEqnPtr_ + record->getOffsetIntoEqnNumbers();
02298 
02299     for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
02300       int eqnOffset = 0;
02301       if (!simpleProblem_) {
02302         fieldMask->getFieldEqnOffset(fieldIDs[fld_offset], eqnOffset);
02303       }
02304 
02305       for(int fs=0; fs<fieldSizes[fld_offset]; ++fs) {
02306         indices[numIndices++] = eqnNumbers[eqnOffset+fs];
02307       }
02308 
02309       ++fld_offset;
02310     }
02311   }
02312 
02313   return(0);
02314 }
02315 
02316 //----------------------------------------------------------------------------
02317 int fei::MatrixGraph_Impl2::addBlockToGraph_singleField_symmetric(fei::Graph* graph,
02318                                                         fei::ConnectivityBlock* cblock)
02319 {
02320   fei::Pattern* pattern = cblock->getRowPattern();
02321   int numIDs = pattern->getNumIDs();
02322   int numIndices = blockEntryGraph_ ? numIDs : pattern->getNumIndices();
02323   int checkNumIndices = numIndices;
02324   std::vector<int> indices(numIndices);
02325   int* indicesPtr = &indices[0];
02326 
02327   const int* fieldIDs = pattern->getFieldIDs();
02328 
02329   int fieldID = fieldIDs[0];
02330   unsigned fieldSize = snl_fei::getFieldSize(fieldID, rowSpace_.get(),
02331                                              colSpace_.get());
02332 
02333   IndexType<int,int>& connIDs = cblock->getNativeConnectivityIDs();
02334   std::vector<int>& rowrecords = cblock->getRowConnectivities();
02335 
02336   IndexType<int,int>::iterator
02337     c_iter = connIDs.begin(),
02338     c_end  = connIDs.end();
02339 
02340   for(; c_iter != c_end; ++c_iter) {
02341     int offset = c_iter->second;
02342     int* records = &rowrecords[offset*numIDs];
02343 
02344     CHK_ERR( getConnectivityIndices_singleField(pattern->getRecordCollections(),
02345                                                 records, numIDs,
02346                                                   fieldID, fieldSize,
02347                                                   checkNumIndices,
02348                                                   indicesPtr,
02349                                                   numIndices) );
02350 
02351     if (checkNumIndices != numIndices) {
02352       ERReturn(-1);
02353     }
02354 
02355     if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
02356         for(int ii=0; ii<numIndices; ++ii) {
02357           if (isLogEqn(indicesPtr[ii])) {
02358             FEI_OSTREAM& os = *output_stream_;
02359             os << "adding Symm inds: ";
02360             for(int jj=0; jj<numIndices; ++jj) {
02361               os << indicesPtr[jj] << " ";
02362             }
02363             os << FEI_ENDL;
02364             break;
02365           }
02366         }
02367     }
02368 
02369     //now we have the indices array, so we're ready to push it into
02370     //the graph container
02371     if (numIndices == numIDs || !cblock->isDiagonal()) {
02372       CHK_ERR( graph->addSymmetricIndices(numIndices, indicesPtr,
02373                                             cblock->isDiagonal()));
02374     }
02375     else {
02376       int ioffset = 0;
02377       for(int i=0; i<numIDs; ++i) {
02378         CHK_ERR( graph->addSymmetricIndices(fieldSize, &(indicesPtr[ioffset])));
02379         ioffset += fieldSize;
02380       }
02381     }
02382   }
02383 
02384   return(0);
02385 }
02386 
02387 //----------------------------------------------------------------------------
02388 int fei::MatrixGraph_Impl2::addBlockToGraph_singleField_nonsymmetric(fei::Graph* graph,
02389                                                         fei::ConnectivityBlock* cblock)
02390 {
02391   fei::Pattern* pattern = cblock->getRowPattern();
02392   fei::Pattern* colpattern = cblock->getColPattern();
02393   int numIDs = pattern->getNumIDs();
02394   int numIndices = blockEntryGraph_ ? numIDs : pattern->getNumIndices();
02395   int checkNumIndices = numIndices;
02396   std::vector<int> indices(numIndices);
02397   int* indicesPtr = &indices[0];
02398 
02399   int numColIDs = colpattern->getNumIDs();
02400   int numColIndices = blockEntryGraph_ ?
02401     numColIDs : colpattern->getNumIndices();
02402   int checkNumColIndices = numColIndices;
02403   std::vector<int> colindices(numColIndices);
02404   int* colindicesPtr = &colindices[0];
02405 
02406   const int* fieldIDs = pattern->getFieldIDs();
02407 
02408   int rowFieldID = fieldIDs[0];
02409   int rowFieldSize = snl_fei::getFieldSize(rowFieldID, rowSpace_.get(),
02410                                            colSpace_.get());
02411 
02412   IndexType<int,int>& connIDs = cblock->getNativeConnectivityIDs();
02413   std::vector<int>& rowrecords = cblock->getRowConnectivities();
02414   std::vector<int>& colrecords = cblock->getColConnectivities();
02415 
02416   int colFieldID = colpattern->getFieldIDs()[0];
02417   int colFieldSize = snl_fei::getFieldSize(colFieldID, rowSpace_.get(),
02418                                            colSpace_.get());
02419 
02420   IndexType<int,int>::iterator
02421     c_iter = connIDs.begin(),
02422     c_end  = connIDs.end();
02423 
02424   for(; c_iter != c_end; ++c_iter) {
02425     int offset = c_iter->second;
02426     int* records = &rowrecords[offset*numIDs];
02427 
02428     int* colRecords = &colrecords[offset*numColIDs];
02429 
02430     if (blockEntryGraph_) {
02431       rowSpace_->getGlobalBlkIndicesL(numIDs, pattern->getRecordCollections(),
02432                                       records, checkNumIndices,
02433                                             indicesPtr, numIndices);
02434 
02435       colSpace_->getGlobalBlkIndicesL(numColIDs, colpattern->getRecordCollections(),
02436                                       colRecords, checkNumColIndices,
02437                                             colindicesPtr, numColIndices);
02438     }
02439     else {
02440       rowSpace_->getGlobalIndicesL(numIDs, pattern->getRecordCollections(),
02441                                       records, rowFieldID,
02442                                          rowFieldSize, checkNumIndices,
02443                                          indicesPtr, numIndices);
02444 
02445       colSpace_->getGlobalIndicesL(numColIDs, colpattern->getRecordCollections(),
02446                                      colRecords, colFieldID,
02447                                          colFieldSize, checkNumColIndices,
02448                                          colindicesPtr, numColIndices);
02449     }
02450 
02451     if (checkNumIndices != numIndices || checkNumColIndices != numColIndices) {
02452       ERReturn(-1);
02453     }
02454 
02455     for(int r=0; r<numIndices; ++r) {
02456       CHK_ERR( graph->addIndices(indicesPtr[r], numColIndices, colindicesPtr));
02457     }
02458   }
02459 
02460   return(0);
02461 }
02462 
02463 //----------------------------------------------------------------------------
02464 int fei::MatrixGraph_Impl2::getConnectivityIndices_singleField(const snl_fei::RecordCollection*const* recordCollections,
02465                                                              int* records,
02466                                                              int numRecords,
02467                                                              int fieldID,
02468                                                              int fieldSize,
02469                                                              int indicesAllocLen,
02470                                                              int* indices,
02471                                                              int& numIndices)
02472 {
02473   numIndices = 0;
02474 
02475   for(int i=0; i<numRecords; ++i) {
02476     if (numIndices == indicesAllocLen) break;
02477 
02478     const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
02479 
02480     if (blockEntryGraph_) {
02481       indices[numIndices++] = record->getNumber();
02482       continue;
02483     }
02484 
02485     int* eqnNumbers = vspcEqnPtr_+record->getOffsetIntoEqnNumbers();
02486 
02487     int eqnOffset = 0;
02488     if (!simpleProblem_) {
02489       const fei::FieldMask* fieldMask = record->getFieldMask();
02490       fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
02491     }
02492 
02493     indices[numIndices++] = eqnNumbers[eqnOffset];
02494     if (fieldSize > 1) {
02495       for(int fs=1; fs<fieldSize; ++fs) {
02496         indices[numIndices++] = eqnNumbers[eqnOffset+fs];
02497       }
02498     }
02499   }
02500 
02501   return(0);
02502 }
02503 
02504 //----------------------------------------------------------------------------
02505 int fei::MatrixGraph_Impl2::getConnectivityIndices_noField(const snl_fei::RecordCollection*const* recordCollections,
02506                                                          int* records,
02507                                                          int numRecords,
02508                                                          int indicesAllocLen,
02509                                                          int* indices,
02510                                                          int& numIndices)
02511 {
02512   numIndices = 0;
02513 
02514   for(int i=0; i<numRecords; ++i) {
02515 
02516     const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
02517     int* eqnNumbers = vspcEqnPtr_+record->getOffsetIntoEqnNumbers();
02518 
02519     if (blockEntryGraph_) {
02520       indices[numIndices++] = record->getNumber();
02521     }
02522     else {
02523       indices[numIndices++] = eqnNumbers[0];
02524     }
02525   }
02526 
02527   return(0);
02528 }
02529 
02530 //----------------------------------------------------------------------------
02531 int fei::MatrixGraph_Impl2::addBlockToGraph_noField_symmetric(fei::Graph* graph,
02532                                                         fei::ConnectivityBlock* cblock)
02533 {
02534   fei::Pattern* pattern = cblock->getRowPattern();
02535   int numIDs = pattern->getNumIDs();
02536   int numIndices = pattern->getNumIndices();
02537   std::vector<int> indices(numIndices);
02538   int* indicesPtr = &indices[0];
02539 
02540   IndexType<int,int>& connIDs = cblock->getNativeConnectivityIDs();
02541   std::vector<int>& rowrecords = cblock->getRowConnectivities();
02542 
02543   IndexType<int,int>::iterator
02544     c_iter = connIDs.begin(),
02545     c_end  = connIDs.end();
02546 
02547   for(; c_iter != c_end; ++c_iter) {
02548     int offset = c_iter->second;
02549     int* records = &rowrecords[offset*numIDs];
02550 
02551     int checkNumIndices;
02552     CHK_ERR( getConnectivityIndices_noField(pattern->getRecordCollections(), records, numIDs, numIndices,
02553                                             indicesPtr, checkNumIndices) );
02554 
02555     if (checkNumIndices != numIndices) {
02556       ERReturn(-1);
02557     }
02558 
02559     if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
02560       for(int ii=0; ii<numIndices; ++ii) {
02561         if (isLogEqn(indicesPtr[ii])) {
02562           FEI_OSTREAM& os = *output_stream_;
02563           os << "adding Symm inds: ";
02564           for(int jj=0; jj<numIndices; ++jj) {
02565             os << indicesPtr[jj] << " ";
02566           }
02567           os << FEI_ENDL;
02568           break;
02569         }
02570       }
02571     }
02572 
02573     //now we have the indices array, so we're ready to push it into
02574     //the graph container
02575     CHK_ERR( graph->addSymmetricIndices(numIndices, indicesPtr) );
02576   }
02577 
02578   return(0);
02579 }
02580 
02581 //----------------------------------------------------------------------------
02582 int fei::MatrixGraph_Impl2::addBlockToGraph_sparse(fei::Graph* graph,
02583                                              fei::ConnectivityBlock* cblock)
02584 {
02585   std::vector<int> row_indices;
02586   std::vector<int> indices;
02587 
02588   fei::Pattern* pattern = cblock->getRowPattern();
02589   const snl_fei::RecordCollection*const* recordCollections = pattern->getRecordCollections();
02590 
02591   IndexType<int,int>& connIDs = cblock->getNativeConnectivityIDs();
02592   std::vector<int>& connOffsets = cblock->getConnectivityOffsets();
02593   int* rowrecords = &(cblock->getRowConnectivities()[0]);
02594   int* colrecords = &(cblock->getColConnectivities()[0]);
02595   bool haveField = cblock->haveFieldID();
02596   int fieldID = cblock->fieldID();
02597   int fieldSize = 1;
02598 
02599   if (haveField) {
02600     fieldSize = snl_fei::getFieldSize(fieldID, rowSpace_.get(),
02601                                       colSpace_.get());
02602   }
02603 
02604   IndexType<int,int>::iterator
02605     c_iter = connIDs.begin(),
02606     c_end  = connIDs.end();
02607 
02608   for(; c_iter != c_end; ++c_iter) {
02609     int offset = c_iter->second;
02610     int rowlen = connOffsets[offset+1] - offset;
02611 
02612     int* records = &(rowrecords[offset]);
02613 
02614     int checkNumIndices;
02615     row_indices.resize(fieldSize);
02616 
02617     if (haveField) {
02618       CHK_ERR( getConnectivityIndices_singleField(recordCollections, records, 1,
02619                                                   fieldID, fieldSize,
02620                                                   fieldSize,
02621                                                   &row_indices[0],
02622                                                   checkNumIndices) );
02623     }
02624     else {
02625       CHK_ERR( getConnectivityIndices_noField(recordCollections, records, 1, fieldSize,
02626                                               &row_indices[0],
02627                                               checkNumIndices) );
02628     }
02629 
02630     indices.resize(fieldSize*rowlen);
02631     int* indicesPtr = &indices[0];
02632     int* crecords = &(colrecords[offset]);
02633 
02634     if (haveField) {
02635       CHK_ERR( getConnectivityIndices_singleField(recordCollections, crecords, rowlen,
02636                                                   fieldID, fieldSize,
02637                                                   fieldSize*rowlen,
02638                                                   indicesPtr,
02639                                                   checkNumIndices) );
02640     }
02641     else {
02642       CHK_ERR( getConnectivityIndices_noField(recordCollections, crecords, rowlen, rowlen,
02643                                               indicesPtr, checkNumIndices) );
02644     }
02645     if (checkNumIndices != rowlen) {
02646       ERReturn(-1);
02647     }
02648 
02649     //now we have the indices, so we're ready to push them into
02650     //the graph container
02651     int* row_ind_ptr = &row_indices[0];
02652     for(int i=0; i<fieldSize; ++i) {
02653       CHK_ERR( graph->addIndices(row_ind_ptr[i], fieldSize*rowlen,
02654                                  indicesPtr) );
02655     }
02656   }
02657 
02658   return(0);
02659 }
02660 
02661 //----------------------------------------------------------------------------
02662 void fei::MatrixGraph_Impl2::setName(const char* name)
02663 {
02664   if (name == NULL) return;
02665 
02666   if (name_ == name) return;
02667 
02668   name_ = name;
02669   dbgprefix_ = "MatGraph_"+name_+": ";
02670 }
02671 
02672 //----------------------------------------------------------------------------
02673 void fei::MatrixGraph_Impl2::setIndicesMode(int mode)
02674 {
02675   if (mode == BLOCK_ENTRY_GRAPH) {
02676     blockEntryGraph_ = true;
02677     return;
02678   }
02679 
02680   if (mode == POINT_ENTRY_GRAPH) {
02681     blockEntryGraph_ = false;
02682     return;
02683   }
02684 
02685   voidERReturn;
02686 }
02687 
02688 //----------------------------------------------------------------------------
02689 fei::SharedPtr<fei::FillableMat>
02690 fei::MatrixGraph_Impl2::getSlaveDependencyMatrix()
02691 {
02692   return( D_ );
02693 }
02694 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends