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