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

Generated on Wed May 12 21:30:41 2010 for FEI by  doxygen 1.4.7