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