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