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   IndexType<int,int>& connectivityIDs = connblk->getNativeConnectivityIDs();
00566 
00567   int idOffset = -1;
00568   IndexType<int,int>::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   IndexType<int,int>& connectivityIDs = connblk->getNativeConnectivityIDs();
00818 
00819   int i, idOffset = -1;
00820   IndexType<int,int>::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->getNativeConnectivityIDs().size();
01792       int numLists = cblock->getNativeConnectivityIDs().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   IndexType<int,int>& connIDs = cblock->getNativeConnectivityIDs();
02124   std::vector<int>& values = cblock->getRowConnectivities();
02125 
02126   IndexType<int,int>::iterator
02127     c_iter = connIDs.begin(),
02128     c_end  = connIDs.end();
02129 
02130   for(; c_iter != c_end; ++c_iter) {
02131     int offset = c_iter->second;
02132     int* records = &values[offset*numIDs];
02133 
02134     CHK_ERR( getConnectivityIndices_multiField(pattern->getRecordCollections(),
02135                                                 records, numIDs,
02136                                                  numFieldsPerID,
02137                                                  fieldIDs,
02138                                                  &fieldSizes[0],
02139                                                  numIndices,
02140                                                  indicesPtr,
02141                                                  checkNumIndices) );
02142 
02143     if (checkNumIndices != numIndices) {
02144       ERReturn(-1);
02145     }
02146 
02147     if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
02148       FEI_OSTREAM& os = *output_stream_;
02149 
02150       const snl_fei::RecordCollection*const* recordColls = pattern->getRecordCollections();
02151       unsigned thisoffset = 0;
02152       for(int ii=0; ii<numIDs; ++ii) {
02153         const fei::Record<int>* record = recordColls[ii]->getRecordWithLocalID(records[ii]);
02154         int ID = record->getID();
02155         os << dbgprefix_<<"scatterIndices: ID=" <<ID<<": ";
02156         int num = pattern->getNumIndicesPerID()[ii];
02157         for(int jj=0; jj<num; ++jj) {
02158           os << indicesPtr[thisoffset++] << " ";
02159         }
02160         os << FEI_ENDL;
02161       }
02162 
02163       for(int ii=0; ii<numIndices; ++ii) {
02164         if (isLogEqn(indicesPtr[ii])) {
02165           os << "adding Symm inds: ";
02166           for(int jj=0; jj<numIndices; ++jj) {
02167             os << indicesPtr[jj] << " ";
02168           }
02169           os << FEI_ENDL;
02170           break;
02171         }
02172       }
02173     }
02174 
02175     //now we have the indices array, so we're ready to push it into
02176     //the graph container
02177     if (numIndices == numIDs || !cblock->isDiagonal()) {
02178       CHK_ERR( graph->addSymmetricIndices(numIndices, indicesPtr,
02179                                           cblock->isDiagonal()) );
02180     }
02181     else {
02182       int ioffset = 0;
02183       int* fieldSizesPtr = &fieldSizes[0];
02184       for(int i=0; i<numIDs; ++i) {
02185         CHK_ERR( graph->addSymmetricIndices(fieldSizesPtr[i], &(indicesPtr[ioffset])));
02186         ioffset += fieldSizesPtr[i];
02187       }
02188     }
02189   }
02190 
02191   return(0);
02192 }
02193 
02194 //----------------------------------------------------------------------------
02195 int fei::MatrixGraph_Impl2::addBlockToGraph_multiField_nonsymmetric(fei::Graph* graph,
02196                                                         fei::ConnectivityBlock* cblock)
02197 {
02198   fei::Pattern* pattern = cblock->getRowPattern();
02199   fei::Pattern* colpattern = cblock->getColPattern();
02200   int j;
02201   int numIDs = pattern->getNumIDs();
02202   int numIndices = blockEntryGraph_ ? numIDs : pattern->getNumIndices();
02203   int checkNumIndices = numIndices;
02204   std::vector<int> indices(numIndices);
02205   int* indicesPtr = &indices[0];
02206 
02207   int numColIDs = colpattern->getNumIDs();
02208   int numColIndices = blockEntryGraph_ ? numColIDs : colpattern->getNumIndices();
02209   int checkNumColIndices = numColIndices;
02210   std::vector<int> colindices(numColIndices);
02211   int* colindicesPtr = &colindices[0];
02212 
02213   const int* numFieldsPerID = pattern->getNumFieldsPerID();
02214   const int* fieldIDs = pattern->getFieldIDs();
02215   int totalNumFields = pattern->getTotalNumFields();
02216 
02217   const int* numFieldsPerColID = colpattern->getNumFieldsPerID();
02218   const int* colfieldIDs = colpattern->getFieldIDs();
02219   int totalNumColFields = colpattern->getTotalNumFields();
02220 
02221   std::vector<int> fieldSizes(totalNumFields);
02222   std::vector<int> colfieldSizes(totalNumColFields);
02223 
02224   for(j=0; j<totalNumFields; ++j) {
02225     fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
02226                                           colSpace_.get());
02227   }
02228   for(j=0; j<totalNumColFields; ++j) {
02229     colfieldSizes[j] = snl_fei::getFieldSize(colfieldIDs[j], rowSpace_.get(),
02230                                            colSpace_.get());
02231   }
02232 
02233   IndexType<int,int>& connIDs = cblock->getNativeConnectivityIDs();
02234   std::vector<int>& rowrecords = cblock->getRowConnectivities();
02235   std::vector<int>& colrecords = cblock->getColConnectivities();
02236 
02237   IndexType<int,int>::iterator
02238    c_iter = connIDs.begin(),
02239    c_end  = connIDs.end();
02240 
02241   for(; c_iter != c_end; ++c_iter) {
02242     int offset = c_iter->second;
02243     int* records = &rowrecords[offset*numIDs];
02244 
02245     int* colRecords = &colrecords[offset*numColIDs];
02246 
02247     CHK_ERR( getConnectivityIndices_multiField(pattern->getRecordCollections(),
02248                                                 records, numIDs,
02249                                                  numFieldsPerID,
02250                                                  fieldIDs,
02251                                                  &fieldSizes[0],
02252                                                  numIndices,
02253                                                  indicesPtr,
02254                                                  checkNumIndices) );
02255 
02256     if (checkNumIndices != numIndices) {
02257       ERReturn(-1);
02258     }
02259 
02260       
02261     CHK_ERR( getConnectivityIndices_multiField(colpattern->getRecordCollections(),
02262                                                 colRecords, numColIDs,
02263                                                  numFieldsPerColID,
02264                                                  colfieldIDs,
02265                                                  &colfieldSizes[0],
02266                                                  numColIndices,
02267                                                  colindicesPtr,
02268                                                  checkNumColIndices) );
02269 
02270     if (checkNumColIndices != numColIndices) {
02271       ERReturn(-1);
02272     }
02273 
02274     //now we have the indices arrays, so we're ready to push them into
02275     //the graph container
02276     for(int r=0; r<numIndices; ++r) {
02277       CHK_ERR( graph->addIndices(indicesPtr[r], numColIndices, colindicesPtr));
02278     }
02279   }
02280 
02281   return(0);
02282 }
02283 
02284 //----------------------------------------------------------------------------
02285 int fei::MatrixGraph_Impl2::getConnectivityIndices_multiField(const snl_fei::RecordCollection*const* recordCollections,
02286                                                           int* records,
02287                                                           int numRecords,
02288                                                           const int* numFieldsPerID,
02289                                                           const int* fieldIDs,
02290                                                           const int* fieldSizes,
02291                                                           int indicesAllocLen,
02292                                                           int* indices,
02293                                                           int& numIndices)
02294 {
02295   numIndices = 0;
02296   int fld_offset = 0;
02297 
02298   for(int i=0; i<numRecords; ++i) {
02299     const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
02300     if (record==NULL) continue;
02301 
02302     if (blockEntryGraph_) {
02303       indices[numIndices++] = record->getNumber();
02304       continue;
02305     }
02306 
02307     const fei::FieldMask* fieldMask = record->getFieldMask();
02308     int* eqnNumbers = vspcEqnPtr_ + record->getOffsetIntoEqnNumbers();
02309 
02310     for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
02311       int eqnOffset = 0;
02312       if (!simpleProblem_) {
02313         fieldMask->getFieldEqnOffset(fieldIDs[fld_offset], eqnOffset);
02314       }
02315 
02316       for(int fs=0; fs<fieldSizes[fld_offset]; ++fs) {
02317         indices[numIndices++] = eqnNumbers[eqnOffset+fs];
02318       }
02319 
02320       ++fld_offset;
02321     }
02322   }
02323 
02324   return(0);
02325 }
02326 
02327 //----------------------------------------------------------------------------
02328 int fei::MatrixGraph_Impl2::addBlockToGraph_singleField_symmetric(fei::Graph* graph,
02329                                                         fei::ConnectivityBlock* cblock)
02330 {
02331   fei::Pattern* pattern = cblock->getRowPattern();
02332   int numIDs = pattern->getNumIDs();
02333   int numIndices = blockEntryGraph_ ? numIDs : pattern->getNumIndices();
02334   int checkNumIndices = numIndices;
02335   std::vector<int> indices(numIndices);
02336   int* indicesPtr = &indices[0];
02337 
02338   const int* fieldIDs = pattern->getFieldIDs();
02339 
02340   int fieldID = fieldIDs[0];
02341   unsigned fieldSize = snl_fei::getFieldSize(fieldID, rowSpace_.get(),
02342                                              colSpace_.get());
02343 
02344   IndexType<int,int>& connIDs = cblock->getNativeConnectivityIDs();
02345   std::vector<int>& rowrecords = cblock->getRowConnectivities();
02346 
02347   IndexType<int,int>::iterator
02348     c_iter = connIDs.begin(),
02349     c_end  = connIDs.end();
02350 
02351   for(; c_iter != c_end; ++c_iter) {
02352     int offset = c_iter->second;
02353     int* records = &rowrecords[offset*numIDs];
02354 
02355     CHK_ERR( getConnectivityIndices_singleField(pattern->getRecordCollections(),
02356                                                 records, numIDs,
02357                                                   fieldID, fieldSize,
02358                                                   checkNumIndices,
02359                                                   indicesPtr,
02360                                                   numIndices) );
02361 
02362     if (checkNumIndices != numIndices) {
02363       ERReturn(-1);
02364     }
02365 
02366     if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
02367         for(int ii=0; ii<numIndices; ++ii) {
02368           if (isLogEqn(indicesPtr[ii])) {
02369             FEI_OSTREAM& os = *output_stream_;
02370             os << "adding Symm inds: ";
02371             for(int jj=0; jj<numIndices; ++jj) {
02372               os << indicesPtr[jj] << " ";
02373             }
02374             os << FEI_ENDL;
02375             break;
02376           }
02377         }
02378     }
02379 
02380     //now we have the indices array, so we're ready to push it into
02381     //the graph container
02382     if (numIndices == numIDs || !cblock->isDiagonal()) {
02383       CHK_ERR( graph->addSymmetricIndices(numIndices, indicesPtr,
02384                                             cblock->isDiagonal()));
02385     }
02386     else {
02387       int ioffset = 0;
02388       for(int i=0; i<numIDs; ++i) {
02389         CHK_ERR( graph->addSymmetricIndices(fieldSize, &(indicesPtr[ioffset])));
02390         ioffset += fieldSize;
02391       }
02392     }
02393   }
02394 
02395   return(0);
02396 }
02397 
02398 //----------------------------------------------------------------------------
02399 int fei::MatrixGraph_Impl2::addBlockToGraph_singleField_nonsymmetric(fei::Graph* graph,
02400                                                         fei::ConnectivityBlock* cblock)
02401 {
02402   fei::Pattern* pattern = cblock->getRowPattern();
02403   fei::Pattern* colpattern = cblock->getColPattern();
02404   int numIDs = pattern->getNumIDs();
02405   int numIndices = blockEntryGraph_ ? numIDs : pattern->getNumIndices();
02406   int checkNumIndices = numIndices;
02407   std::vector<int> indices(numIndices);
02408   int* indicesPtr = &indices[0];
02409 
02410   int numColIDs = colpattern->getNumIDs();
02411   int numColIndices = blockEntryGraph_ ?
02412     numColIDs : colpattern->getNumIndices();
02413   int checkNumColIndices = numColIndices;
02414   std::vector<int> colindices(numColIndices);
02415   int* colindicesPtr = &colindices[0];
02416 
02417   const int* fieldIDs = pattern->getFieldIDs();
02418 
02419   int rowFieldID = fieldIDs[0];
02420   int rowFieldSize = snl_fei::getFieldSize(rowFieldID, rowSpace_.get(),
02421                                            colSpace_.get());
02422 
02423   IndexType<int,int>& connIDs = cblock->getNativeConnectivityIDs();
02424   std::vector<int>& rowrecords = cblock->getRowConnectivities();
02425   std::vector<int>& colrecords = cblock->getColConnectivities();
02426 
02427   int colFieldID = colpattern->getFieldIDs()[0];
02428   int colFieldSize = snl_fei::getFieldSize(colFieldID, rowSpace_.get(),
02429                                            colSpace_.get());
02430 
02431   IndexType<int,int>::iterator
02432     c_iter = connIDs.begin(),
02433     c_end  = connIDs.end();
02434 
02435   for(; c_iter != c_end; ++c_iter) {
02436     int offset = c_iter->second;
02437     int* records = &rowrecords[offset*numIDs];
02438 
02439     int* colRecords = &colrecords[offset*numColIDs];
02440 
02441     if (blockEntryGraph_) {
02442       rowSpace_->getGlobalBlkIndicesL(numIDs, pattern->getRecordCollections(),
02443                                       records, checkNumIndices,
02444                                             indicesPtr, numIndices);
02445 
02446       colSpace_->getGlobalBlkIndicesL(numColIDs, colpattern->getRecordCollections(),
02447                                       colRecords, checkNumColIndices,
02448                                             colindicesPtr, numColIndices);
02449     }
02450     else {
02451       rowSpace_->getGlobalIndicesL(numIDs, pattern->getRecordCollections(),
02452                                       records, rowFieldID,
02453                                          rowFieldSize, checkNumIndices,
02454                                          indicesPtr, numIndices);
02455 
02456       colSpace_->getGlobalIndicesL(numColIDs, colpattern->getRecordCollections(),
02457                                      colRecords, colFieldID,
02458                                          colFieldSize, checkNumColIndices,
02459                                          colindicesPtr, numColIndices);
02460     }
02461 
02462     if (checkNumIndices != numIndices || checkNumColIndices != numColIndices) {
02463       ERReturn(-1);
02464     }
02465 
02466     for(int r=0; r<numIndices; ++r) {
02467       CHK_ERR( graph->addIndices(indicesPtr[r], numColIndices, colindicesPtr));
02468     }
02469   }
02470 
02471   return(0);
02472 }
02473 
02474 //----------------------------------------------------------------------------
02475 int fei::MatrixGraph_Impl2::getConnectivityIndices_singleField(const snl_fei::RecordCollection*const* recordCollections,
02476                                                              int* records,
02477                                                              int numRecords,
02478                                                              int fieldID,
02479                                                              int fieldSize,
02480                                                              int indicesAllocLen,
02481                                                              int* indices,
02482                                                              int& numIndices)
02483 {
02484   numIndices = 0;
02485 
02486   for(int i=0; i<numRecords; ++i) {
02487     if (numIndices == indicesAllocLen) break;
02488 
02489     const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
02490 
02491     if (blockEntryGraph_) {
02492       indices[numIndices++] = record->getNumber();
02493       continue;
02494     }
02495 
02496     int* eqnNumbers = vspcEqnPtr_+record->getOffsetIntoEqnNumbers();
02497 
02498     int eqnOffset = 0;
02499     if (!simpleProblem_) {
02500       const fei::FieldMask* fieldMask = record->getFieldMask();
02501       fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
02502     }
02503 
02504     indices[numIndices++] = eqnNumbers[eqnOffset];
02505     if (fieldSize > 1) {
02506       for(int fs=1; fs<fieldSize; ++fs) {
02507         indices[numIndices++] = eqnNumbers[eqnOffset+fs];
02508       }
02509     }
02510   }
02511 
02512   return(0);
02513 }
02514 
02515 //----------------------------------------------------------------------------
02516 int fei::MatrixGraph_Impl2::getConnectivityIndices_noField(const snl_fei::RecordCollection*const* recordCollections,
02517                                                          int* records,
02518                                                          int numRecords,
02519                                                          int indicesAllocLen,
02520                                                          int* indices,
02521                                                          int& numIndices)
02522 {
02523   numIndices = 0;
02524 
02525   for(int i=0; i<numRecords; ++i) {
02526 
02527     const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
02528     int* eqnNumbers = vspcEqnPtr_+record->getOffsetIntoEqnNumbers();
02529 
02530     if (blockEntryGraph_) {
02531       indices[numIndices++] = record->getNumber();
02532     }
02533     else {
02534       indices[numIndices++] = eqnNumbers[0];
02535     }
02536   }
02537 
02538   return(0);
02539 }
02540 
02541 //----------------------------------------------------------------------------
02542 int fei::MatrixGraph_Impl2::addBlockToGraph_noField_symmetric(fei::Graph* graph,
02543                                                         fei::ConnectivityBlock* cblock)
02544 {
02545   fei::Pattern* pattern = cblock->getRowPattern();
02546   int numIDs = pattern->getNumIDs();
02547   int numIndices = pattern->getNumIndices();
02548   std::vector<int> indices(numIndices);
02549   int* indicesPtr = &indices[0];
02550 
02551   IndexType<int,int>& connIDs = cblock->getNativeConnectivityIDs();
02552   std::vector<int>& rowrecords = cblock->getRowConnectivities();
02553 
02554   IndexType<int,int>::iterator
02555     c_iter = connIDs.begin(),
02556     c_end  = connIDs.end();
02557 
02558   for(; c_iter != c_end; ++c_iter) {
02559     int offset = c_iter->second;
02560     int* records = &rowrecords[offset*numIDs];
02561 
02562     int checkNumIndices;
02563     CHK_ERR( getConnectivityIndices_noField(pattern->getRecordCollections(), records, numIDs, numIndices,
02564                                             indicesPtr, checkNumIndices) );
02565 
02566     if (checkNumIndices != numIndices) {
02567       ERReturn(-1);
02568     }
02569 
02570     if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
02571       for(int ii=0; ii<numIndices; ++ii) {
02572         if (isLogEqn(indicesPtr[ii])) {
02573           FEI_OSTREAM& os = *output_stream_;
02574           os << "adding Symm inds: ";
02575           for(int jj=0; jj<numIndices; ++jj) {
02576             os << indicesPtr[jj] << " ";
02577           }
02578           os << FEI_ENDL;
02579           break;
02580         }
02581       }
02582     }
02583 
02584     //now we have the indices array, so we're ready to push it into
02585     //the graph container
02586     CHK_ERR( graph->addSymmetricIndices(numIndices, indicesPtr) );
02587   }
02588 
02589   return(0);
02590 }
02591 
02592 //----------------------------------------------------------------------------
02593 int fei::MatrixGraph_Impl2::addBlockToGraph_sparse(fei::Graph* graph,
02594                                              fei::ConnectivityBlock* cblock)
02595 {
02596   std::vector<int> row_indices;
02597   std::vector<int> indices;
02598 
02599   fei::Pattern* pattern = cblock->getRowPattern();
02600   const snl_fei::RecordCollection*const* recordCollections = pattern->getRecordCollections();
02601 
02602   IndexType<int,int>& connIDs = cblock->getNativeConnectivityIDs();
02603   std::vector<int>& connOffsets = cblock->getConnectivityOffsets();
02604   int* rowrecords = &(cblock->getRowConnectivities()[0]);
02605   int* colrecords = &(cblock->getColConnectivities()[0]);
02606   bool haveField = cblock->haveFieldID();
02607   int fieldID = cblock->fieldID();
02608   int fieldSize = 1;
02609 
02610   if (haveField) {
02611     fieldSize = snl_fei::getFieldSize(fieldID, rowSpace_.get(),
02612                                       colSpace_.get());
02613   }
02614 
02615   IndexType<int,int>::iterator
02616     c_iter = connIDs.begin(),
02617     c_end  = connIDs.end();
02618 
02619   for(; c_iter != c_end; ++c_iter) {
02620     int offset = c_iter->second;
02621     int rowlen = connOffsets[offset+1] - offset;
02622 
02623     int* records = &(rowrecords[offset]);
02624 
02625     int checkNumIndices;
02626     row_indices.resize(fieldSize);
02627 
02628     if (haveField) {
02629       CHK_ERR( getConnectivityIndices_singleField(recordCollections, records, 1,
02630                                                   fieldID, fieldSize,
02631                                                   fieldSize,
02632                                                   &row_indices[0],
02633                                                   checkNumIndices) );
02634     }
02635     else {
02636       CHK_ERR( getConnectivityIndices_noField(recordCollections, records, 1, fieldSize,
02637                                               &row_indices[0],
02638                                               checkNumIndices) );
02639     }
02640 
02641     indices.resize(fieldSize*rowlen);
02642     int* indicesPtr = &indices[0];
02643     int* crecords = &(colrecords[offset]);
02644 
02645     if (haveField) {
02646       CHK_ERR( getConnectivityIndices_singleField(recordCollections, crecords, rowlen,
02647                                                   fieldID, fieldSize,
02648                                                   fieldSize*rowlen,
02649                                                   indicesPtr,
02650                                                   checkNumIndices) );
02651     }
02652     else {
02653       CHK_ERR( getConnectivityIndices_noField(recordCollections, crecords, rowlen, rowlen,
02654                                               indicesPtr, checkNumIndices) );
02655     }
02656     if (checkNumIndices != rowlen) {
02657       ERReturn(-1);
02658     }
02659 
02660     //now we have the indices, so we're ready to push them into
02661     //the graph container
02662     int* row_ind_ptr = &row_indices[0];
02663     for(int i=0; i<fieldSize; ++i) {
02664       CHK_ERR( graph->addIndices(row_ind_ptr[i], fieldSize*rowlen,
02665                                  indicesPtr) );
02666     }
02667   }
02668 
02669   return(0);
02670 }
02671 
02672 //----------------------------------------------------------------------------
02673 void fei::MatrixGraph_Impl2::setName(const char* name)
02674 {
02675   if (name == NULL) return;
02676 
02677   if (name_ == name) return;
02678 
02679   name_ = name;
02680   dbgprefix_ = "MatGraph_"+name_+": ";
02681 }
02682 
02683 //----------------------------------------------------------------------------
02684 void fei::MatrixGraph_Impl2::setIndicesMode(int mode)
02685 {
02686   if (mode == BLOCK_ENTRY_GRAPH) {
02687     blockEntryGraph_ = true;
02688     return;
02689   }
02690 
02691   if (mode == POINT_ENTRY_GRAPH) {
02692     blockEntryGraph_ = false;
02693     return;
02694   }
02695 
02696   voidERReturn;
02697 }
02698 
02699 //----------------------------------------------------------------------------
02700 fei::SharedPtr<fei::FillableMat>
02701 fei::MatrixGraph_Impl2::getSlaveDependencyMatrix()
02702 {
02703   return( D_ );
02704 }
02705 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends