FEI Version of the Day
fei_LinSysCoreFilter.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 #include <assert.h>
00014 
00015 #include "fei_utils.hpp"
00016 #include <fei_impl_utils.hpp>
00017 
00018 #include "fei_defs.h"
00019 
00020 #include <fei_ostream_ops.hpp>
00021 #include "fei_CommUtils.hpp"
00022 #include "fei_TemplateUtils.hpp"
00023 #include "snl_fei_Constraint.hpp"
00024 typedef snl_fei::Constraint<GlobalID> ConstraintType;
00025 
00026 #include "fei_DirichletBCManager.hpp"
00027 #include "fei_FillableMat.hpp"
00028 #include "fei_CSVec.hpp"
00029 #include "FEI_Implementation.hpp"
00030 #include "fei_EqnCommMgr.hpp"
00031 #include "fei_ConnectivityTable.hpp"
00032 #include "fei_NodeDescriptor.hpp"
00033 #include "fei_NodeDatabase.hpp"
00034 #include "fei_BlockDescriptor.hpp"
00035 #include "SNL_FEI_Structure.hpp"
00036 #include "snl_fei_Utils.hpp"
00037 #include "fei_Data.hpp"
00038 #include "fei_LinearSystemCore.hpp"
00039 #include "fei_LinSysCore_flexible.hpp"
00040 
00041 #include "fei_LinSysCoreFilter.hpp"
00042 
00043 #undef fei_file
00044 #define fei_file "fei_LinSysCoreFilter.cpp"
00045 #include "fei_ErrMacros.hpp"
00046 
00047 #define ASSEMBLE_PUT 0
00048 #define ASSEMBLE_SUM 1
00049 
00050 
00051 //------------------------------------------------------------------------------
00052 LinSysCoreFilter::LinSysCoreFilter(FEI_Implementation* owner,
00053                                    MPI_Comm comm,
00054                                    SNL_FEI_Structure* probStruct,
00055                                    LinearSystemCore* lsc,
00056                                    int masterRank)
00057  : Filter(probStruct),
00058    timesInitializeCalled_(0),
00059    lsc_(lsc),
00060    useLookup_(true),
00061    newMatrixData_(false),
00062    newVectorData_(false),
00063    newConstraintData_(false),
00064    newBCData_(false),
00065    connectivitiesInitialized_(false),
00066    firstRemEqnExchange_(true),
00067    needToCallMatrixLoadComplete_(false),
00068    resolveConflictRequested_(false),
00069    localStartRow_(0),             //
00070    localEndRow_(0),               //Initialize all private variables here,
00071    numLocalEqns_(0),              //in the order that they are declared.
00072    numGlobalEqns_(0),             //(Don't bother initializing those that will
00073    numLocallyOwnedNodes_(0),      //be initialized in the body of the 
00074    numGlobalNodes_(0),            //constructor below.)
00075    firstLocalNodeNumber_(-1),     //
00076    blockMatrix_(false),
00077    tooLateToChooseBlock_(false),
00078    numLocalEqnBlks_(0),
00079    localReducedBlkOffset_(0),
00080    numLocalReducedEqnBlks_(0),
00081    iterations_(0),
00082    currentRHS_(0),
00083    rhsIDs_(),
00084    outputLevel_(0),
00085    comm_(comm),
00086    masterRank_(masterRank),
00087    problemStructure_(probStruct),
00088    matrixAllocated_(false),
00089    rowIndices_(),
00090    rowColOffsets_(0),
00091    colIndices_(0),
00092    nodeIDType_(0),
00093    eqnCommMgr_(NULL),
00094    eqnCommMgr_put_(NULL),
00095    maxElemRows_(0),
00096    scatterIndices_(),
00097    blkScatterIndices_(),
00098    iworkSpace_(),
00099    iworkSpace2_(),
00100    dworkSpace_(),
00101    dworkSpace2_(),
00102    eStiff_(NULL),
00103    eStiff1D_(NULL),
00104    eLoad_(NULL)
00105 {
00106     localRank_ = fei::localProc(comm_);
00107     numProcs_ = fei::numProcs(comm_);
00108 
00109     internalFei_ = 0;
00110 
00111     numRHSs_ = 1;
00112     rhsIDs_.resize(numRHSs_);
00113     rhsIDs_[0] = 0;
00114 
00115     bcManager_ = new fei::DirichletBCManager(probStruct);
00116     eqnCommMgr_ = problemStructure_->getEqnCommMgr().deepCopy();
00117     int err = createEqnCommMgr_put();
00118 
00119    //We need to get the parameters from the owning FEI_Implementation...
00120    int numParams = -1;
00121    char** paramStrings = NULL;
00122    err = owner->getParameters(numParams, paramStrings);
00123 
00124    //Now let's pass them into our own parameter-handling mechanism.
00125    err = parameters(numParams, paramStrings);
00126    if (err != 0) {
00127      FEI_CERR << "LinSysCoreFilter::LinSysCoreFilter ERROR, parameters failed." << FEI_ENDL;
00128      MPI_Abort(comm_, -1);
00129    }
00130 
00131    Kid_ = new fei::FillableMat;
00132    Kdi_ = new fei::FillableMat;
00133    Kdd_ = new fei::FillableMat;
00134    reducedEqnCounter_ = 0;
00135    reducedRHSCounter_ = 0;
00136 
00137    return;
00138 }
00139 
00140 //------------------------------------------------------------------------------
00141 LinSysCoreFilter::~LinSysCoreFilter() {
00142 //
00143 //  Destructor function. Free allocated memory, etc.
00144 //
00145   numRHSs_ = 0;
00146 
00147   delete bcManager_;
00148   delete eqnCommMgr_;
00149   delete eqnCommMgr_put_;
00150 
00151   delete [] eStiff_;
00152   delete [] eStiff1D_;
00153   delete [] eLoad_;
00154 
00155   delete Kid_;
00156   delete Kdi_;
00157   delete Kdd_;
00158 }
00159 
00160 //------------------------------------------------------------------------------
00161 int LinSysCoreFilter::initialize()
00162 {
00163 // Determine final sparsity pattern for setting the structure of the
00164 // underlying sparse matrix.
00165 //
00166   debugOutput("#LinSysCoreFilter::initialize");
00167 
00168   numLocalEqns_ = problemStructure_->getNumLocalEqns();
00169   numLocalEqnBlks_ = problemStructure_->getNumLocalEqnBlks();
00170   numLocalReducedEqnBlks_ = problemStructure_->getNumLocalReducedEqnBlks();
00171 
00172   // now, obtain the global equation info, such as how many equations there
00173   // are globally, and what the local starting and ending row-numbers are.
00174 
00175   // let's also get the number of global nodes, and a first-local-node-number.
00176   // node-number is a globally 0-based number we are assigning to nodes.
00177   // node-numbers are contiguous on a processor -- i.e., a processor owns a
00178   // contiguous block of node-numbers. This provides an easier-to-work-with
00179   // node numbering than the application-supplied nodeIDs which may not be
00180   // assumed to be contiguous or 0-based, or anything else.
00181 
00182   std::vector<int>& eqnOffsets = problemStructure_->getGlobalEqnOffsets();
00183   localStartRow_ = eqnOffsets[localRank_];
00184   localEndRow_ = eqnOffsets[localRank_+1]-1;
00185   numGlobalEqns_ = eqnOffsets[numProcs_];
00186 
00187   std::vector<int>& nodeOffsets = problemStructure_->getGlobalNodeOffsets();
00188   firstLocalNodeNumber_ =  nodeOffsets[localRank_];
00189   numGlobalNodes_ = nodeOffsets[numProcs_];
00190 
00191   //--------------------------------------------------------------------------
00192   //  ----- end active equation calculations -----
00193 
00194   if (eqnCommMgr_ != NULL) delete eqnCommMgr_;
00195   eqnCommMgr_ = NULL;
00196   if (eqnCommMgr_put_ != NULL) delete eqnCommMgr_put_;
00197   eqnCommMgr_put_ = NULL;
00198 
00199   eqnCommMgr_ = problemStructure_->getEqnCommMgr().deepCopy();
00200   if (eqnCommMgr_ == NULL) ERReturn(-1);
00201 
00202   int err = createEqnCommMgr_put();
00203   if (err != 0) ERReturn(err);
00204 
00205   //(we need to set the number of RHSs in the eqn comm manager)
00206   eqnCommMgr_->setNumRHSs(numRHSs_);
00207 
00208   //let's let the underlying linear system know what the global offsets are.
00209   //While we're dealing with global offsets, we'll also calculate the starting
00210   //and ending 'reduced' rows, etc.
00211   CHK_ERR( initLinSysCore() );
00212 
00213   setLinSysCoreCREqns();
00214 
00215   if (timesInitializeCalled_ == 0) {
00216     //
00217     // let's prepare some arrays for handing the matrix structure to
00218     // the linear system.
00219 
00220     std::vector<int> rowLengths;
00221     CHK_ERR( problemStructure_->getMatrixRowLengths(rowLengths) );
00222 
00223     int numReducedEqns = problemStructure_->getNumLocalReducedEqns();
00224     int maxBlkSize = problemStructure_->getGlobalMaxBlkSize();
00225     std::vector<int> blkSizes(numLocalReducedEqnBlks_, 1);
00226 
00227     int numNonzeros = 0;
00228     for(size_t ii=0; ii<rowLengths.size(); ++ii) {
00229       numNonzeros += rowLengths[ii];
00230     }
00231 
00232     std::vector<int> indices_1D(numNonzeros);
00233     std::vector<int*> indices(numReducedEqns);
00234 
00235     int offset = 0;
00236     for(size_t ii=0; ii<rowLengths.size(); ++ii) {
00237       indices[ii] = &(indices_1D[offset]);
00238       offset += rowLengths[ii];
00239     }
00240 
00241     if (maxBlkSize == 0) ERReturn(-1);
00242 
00243     if (maxBlkSize == 1) {
00244       CHK_ERR( problemStructure_->getMatrixStructure(&indices[0], rowLengths) );
00245 
00246       debugOutput("#LinSysCoreFilter calling point lsc_->setMatrixStructure");
00247       CHK_ERR( lsc_->setMatrixStructure(&indices[0], &rowLengths[0],
00248                                         &indices[0], &rowLengths[0], &blkSizes[0]) );
00249     }
00250     else {
00251       std::vector<int> blkRowLengths;
00252       int* blkIndices_1D = numNonzeros > 0 ? new int[numNonzeros] : NULL;
00253       int** blkIndices = numLocalReducedEqnBlks_ ?
00254         new int*[numLocalReducedEqnBlks_] : NULL;
00255       if (blkIndices == NULL && numLocalReducedEqnBlks_ != 0) ERReturn(-1);
00256 
00257       CHK_ERR( problemStructure_->getMatrixStructure(&indices[0], rowLengths,
00258                                                      blkIndices, blkIndices_1D,
00259                                                      blkRowLengths, blkSizes) );
00260 
00261       offset = 0;
00262       for(int ii=0; ii<numLocalReducedEqnBlks_; ++ii) {
00263         blkIndices[ii] = &(blkIndices_1D[offset]);
00264         offset += blkRowLengths[ii];
00265       }
00266 
00267       debugOutput("#LinSysCoreFilter calling block lsc_->setMatrixStructure");
00268       CHK_ERR( lsc_->setMatrixStructure(&indices[0], &rowLengths[0],
00269                                         blkIndices, &blkRowLengths[0], &blkSizes[0]) );
00270 
00271       if (numLocalReducedEqnBlks_ != 0) delete [] blkIndices;
00272       if (numNonzeros > 0) delete [] blkIndices_1D;
00273     }
00274   }
00275 
00276   matrixAllocated_ = true;
00277 
00278   debugOutput("#leaving LinSysCoreFilter::initialize");
00279   ++timesInitializeCalled_;
00280   return(FEI_SUCCESS);
00281 }
00282 
00283 //------------------------------------------------------------------------------
00284 int LinSysCoreFilter::createEqnCommMgr_put()
00285 {
00286   if (eqnCommMgr_put_ != NULL) return(0);
00287 
00288   eqnCommMgr_put_  = eqnCommMgr_->deepCopy();
00289   if (eqnCommMgr_put_ == NULL) ERReturn(-1);
00290 
00291   eqnCommMgr_put_->resetCoefs();
00292   eqnCommMgr_put_->accumulate_ = false;
00293   return(0);
00294 }
00295 
00296 //==============================================================================
00297 int LinSysCoreFilter::initLinSysCore()
00298 {
00299   debugOutput("#LinSysCoreFilter calling lsc_->setLookup");
00300   int err = lsc_->setLookup(*problemStructure_);
00301 
00302   if (err != 0) {
00303     useLookup_ = false;
00304   }
00305 
00306   std::vector<int>& globalNodeOffsets = problemStructure_->getGlobalNodeOffsets();
00307   std::vector<int>& globalEqnOffsets = problemStructure_->getGlobalEqnOffsets();
00308   std::vector<int>& globalBlkEqnOffsets =
00309     problemStructure_->getGlobalBlkEqnOffsets();
00310 
00311   int startRow = localStartRow_;
00312   while(problemStructure_->isSlaveEqn(startRow)) startRow++;
00313   int endRow = localEndRow_;
00314   while(problemStructure_->isSlaveEqn(endRow)) endRow--;
00315 
00316   problemStructure_->translateToReducedEqn(startRow, reducedStartRow_);
00317   problemStructure_->translateToReducedEqn(endRow, reducedEndRow_);
00318   numReducedRows_ = reducedEndRow_ - reducedStartRow_ + 1;
00319 
00320   std::vector<int> reducedEqnOffsets(globalEqnOffsets.size());
00321   std::vector<int> reducedBlkEqnOffsets(globalBlkEqnOffsets.size());
00322   std::vector<int> reducedNodeOffsets(globalNodeOffsets.size());
00323 
00324   int numSlaveEqns = problemStructure_->numSlaveEquations();
00325 
00326   int reducedNodeNum =
00327     problemStructure_->getAssociatedNodeNumber(reducedStartRow_);
00328 
00329   std::vector<int> tmpSend(2), tmpRecv, tmpRecvLengths;
00330   tmpSend[0] = reducedStartRow_;
00331   tmpSend[1] = reducedNodeNum;
00332   CHK_ERR( fei::Allgatherv(comm_, tmpSend, tmpRecvLengths, tmpRecv) );
00333 
00334   for(size_t ii=0; ii<globalEqnOffsets.size()-1; ++ii) {
00335     reducedEqnOffsets[ii] = tmpRecv[ii*2];
00336     reducedNodeOffsets[ii] = tmpRecv[ii*2+1];
00337 
00338     //Major assumption: we're assuming that if there are slave equations, then
00339     //there are no lagrange-multiplier constraints. Because if there are no
00340     //lagrange-multiplier constraints, then blkEqn == nodeNum.
00341     if (numSlaveEqns > 0) {
00342       reducedBlkEqnOffsets[ii] = reducedNodeOffsets[ii];
00343     }
00344     else {
00345       reducedBlkEqnOffsets[ii] = globalBlkEqnOffsets[ii];
00346     }
00347   }
00348 
00349   if (localRank_ == numProcs_-1) {
00350     reducedNodeNum = problemStructure_->
00351       translateToReducedNodeNumber(globalNodeOffsets[numProcs_]-1, localRank_);
00352     int blkEqn = globalBlkEqnOffsets[numProcs_]-1;
00353     if (numSlaveEqns > 0) {
00354       blkEqn = reducedNodeNum;
00355     }
00356 
00357     tmpSend.resize(3);
00358     tmpSend[0] = reducedEndRow_;
00359     tmpSend[1] = reducedNodeNum;
00360     tmpSend[2] = blkEqn;
00361   }
00362   else {
00363     tmpSend.resize(3);
00364   }
00365 
00366   CHK_ERR( fei::Bcast(comm_, tmpSend, numProcs_-1) );
00367   reducedEqnOffsets[numProcs_] = tmpSend[0]+1;
00368   reducedNodeOffsets[numProcs_] = tmpSend[1]+1;
00369   reducedBlkEqnOffsets[numProcs_] = tmpSend[2]+1;
00370 
00371   debugOutput("#LinSysCoreFilter calling lsc_->setGlobalOffsets");
00372   CHK_ERR( lsc_->setGlobalOffsets(numProcs_+1,
00373                                   &reducedNodeOffsets[0],
00374                                   &reducedEqnOffsets[0],
00375                                   &reducedBlkEqnOffsets[0]) );
00376 
00377   if (connectivitiesInitialized_) return(0);
00378 
00379   int numBlocks = problemStructure_->getNumElemBlocks();
00380   NodeDatabase& nodeDB     = problemStructure_->getNodeDatabase();
00381   NodeCommMgr& nodeCommMgr = problemStructure_->getNodeCommMgr();
00382 
00383   int numNodes = nodeDB.getNumNodeDescriptors();
00384   int numRemoteNodes = nodeCommMgr.getSharedNodeIDs().size() -
00385                          nodeCommMgr.getLocalNodeIDs().size();
00386   numNodes -= numRemoteNodes;
00387 
00388   std::vector<int> numElemsPerBlock(numBlocks);
00389   std::vector<int> numNodesPerElem(numBlocks);
00390   std::vector<int> numDofPerNode(0,1);
00391 
00392   for(int blk=0; blk<numBlocks; blk++) {
00393     BlockDescriptor* block = NULL;
00394     CHK_ERR( problemStructure_->getBlockDescriptor_index(blk, block) );
00395 
00396     numElemsPerBlock[blk] = block->getNumElements();
00397     numNodesPerElem[blk]  = block->getNumNodesPerElement();
00398 
00399     int* fieldsPerNode = block->fieldsPerNodePtr();
00400     int** fieldIDsTable = block->fieldIDsTablePtr();
00401 
00402     for(int nn=0; nn<numNodesPerElem[blk]; nn++) {
00403       if (fieldsPerNode[nn] <= 0) ERReturn(-1);
00404       numDofPerNode.push_back(0);
00405       int indx = numDofPerNode.size()-1;
00406       
00407       for(int nf=0; nf<fieldsPerNode[nn]; nf++) {
00408         numDofPerNode[indx] += problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
00409       }
00410     }
00411   }
00412 
00413   for(int i=0; i<numBlocks; i++) {
00414     BlockDescriptor* block = NULL;
00415     CHK_ERR( problemStructure_->getBlockDescriptor_index(i, block) );
00416 
00417     if (block->getNumElements() == 0) continue;
00418 
00419     ConnectivityTable& ctbl =
00420       problemStructure_->getBlockConnectivity(block->getGlobalBlockID());
00421 
00422     std::vector<int> cNodeList(block->getNumNodesPerElement());
00423 
00424     int* fieldsPerNode = block->fieldsPerNodePtr();
00425     int** fieldIDsTable = block->fieldIDsTablePtr();
00426 
00427     numDofPerNode.resize(0);
00428     for(int nn=0; nn<numNodesPerElem[i]; nn++) {
00429       if (fieldsPerNode[nn] <= 0) ERReturn(-1);
00430       numDofPerNode.push_back(0);
00431       int indx = numDofPerNode.size()-1;
00432 
00433       for(int nf=0; nf<fieldsPerNode[nn]; nf++) {
00434         numDofPerNode[indx] += problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
00435       }
00436     }
00437 
00438     int nodesPerElement = block->getNumNodesPerElement();
00439     NodeDescriptor** elemNodePtrs = &((*ctbl.elem_conn_ptrs)[0]);
00440     int offset = 0;
00441     for(int j=0; j<block->getNumElements(); j++) {
00442 
00443       for(int k=0; k<nodesPerElement; k++) {
00444         NodeDescriptor* node = elemNodePtrs[offset++];
00445         cNodeList[k] = node->getNodeNumber();
00446       }
00447 
00448       int elemID = ctbl.elemNumbers[j];
00449       int* nodeNumbers = &cNodeList[0];
00450       debugOutput("#LinSysCoreFilter calling lsc->setConnectivities");
00451       CHK_ERR( lsc_->setConnectivities(i, 1, nodesPerElement,
00452                                        &elemID, &nodeNumbers) );
00453     }
00454   }
00455 
00456   connectivitiesInitialized_ = true;
00457 
00458   return(FEI_SUCCESS);
00459 }
00460 
00461 //==============================================================================
00462 void LinSysCoreFilter::setLinSysCoreCREqns()
00463 {
00464    int err, i=0;
00465 
00466    std::vector<int> iwork;
00467 
00468    std::map<GlobalID,ConstraintType*>::const_iterator
00469      cr_iter = problemStructure_->getMultConstRecords().begin(),
00470      cr_end = problemStructure_->getMultConstRecords().end();
00471 
00472    while(cr_iter != cr_end) {
00473       ConstraintType& constraint = *((*cr_iter).second);
00474       int numNodesPerCR = constraint.getMasters().size();
00475       int meqn = constraint.getEqnNumber();
00476 
00477       std::vector<GlobalID>& nodeID_vec = constraint.getMasters();
00478       GlobalID* nodeIDPtr = &nodeID_vec[0];
00479 
00480       if ((int)iwork.size() < 2*numNodesPerCR) {
00481         iwork.resize(2*numNodesPerCR);
00482       }
00483 
00484       int* nodeList = &(iwork[0]);
00485 
00486       int* eqnList = nodeList+numNodesPerCR;
00487 
00488       std::vector<int>& fieldIDs_vec = constraint.getMasterFieldIDs();
00489       int* fieldIDs = &fieldIDs_vec[0];
00490 
00491       for(int k=0; k<numNodesPerCR; k++) {
00492         const NodeDescriptor& node = Filter::findNodeDescriptor(nodeIDPtr[k]);
00493         nodeList[k] = node.getNodeNumber();
00494 
00495         int eqn = -1;
00496         if (!node.getFieldEqnNumber(fieldIDs[k], eqn)) voidERReturn;
00497 
00498         eqnList[k] = eqn;
00499       }
00500 
00501       int crMultID = constraint.getConstraintID() + i++;
00502       if (Filter::logStream() != NULL) {
00503         FEI_OSTREAM& os = *logStream();
00504         os << "#LinSysCoreFilter calling lsc_->setMultCREqns"<<FEI_ENDL;
00505         os << "#  multiplier eqn: " << meqn << ", columns: ";
00506         for(int j=0; j<numNodesPerCR; ++j) os << eqnList[j] << " ";
00507         os << FEI_ENDL;
00508       }
00509 
00510       err = lsc_->setMultCREqns(crMultID, 1, numNodesPerCR,
00511                                 &nodeList, &eqnList,
00512                                 fieldIDs, &meqn);
00513       if (err) voidERReturn;
00514       ++cr_iter;
00515    }
00516 
00517    LinSysCore_flexible* lscf = dynamic_cast<LinSysCore_flexible*>(lsc_);
00518    if (lscf != NULL) {
00519      debugOutput("LinSysCoreFilter calling lscf->setMultCRComplete");
00520      err = lscf->setMultCRComplete();
00521      if (err != 0) {
00522        FEI_CERR << "LinSysCoreFilter::setLinSysCoreCREqns ERROR returned from "
00523             << "lscf->setMultCRComplete()" << FEI_ENDL;
00524      }
00525    }
00526 
00527    //
00528    //Now the penalty CRs...
00529    //
00530    cr_iter = problemStructure_->getPenConstRecords().begin();
00531    cr_end = problemStructure_->getPenConstRecords().end();
00532 
00533    while(cr_iter != cr_end) {
00534       ConstraintType& crset = *((*cr_iter).second);
00535       int numNodesPerCR = crset.getMasters().size();
00536 
00537       std::vector<GlobalID>& nodeIDs_vec = crset.getMasters();
00538       GlobalID* nodeIDsPtr = &nodeIDs_vec[0];
00539 
00540       if ((int)iwork.size() < 2*numNodesPerCR) {
00541         iwork.resize(2*numNodesPerCR);
00542       }
00543 
00544       int* nodeList = &(iwork[0]);
00545 
00546       int* eqnList = nodeList+numNodesPerCR;
00547 
00548       std::vector<int>& fieldIDs_vec = crset.getMasterFieldIDs();
00549       int* fieldIDs = &fieldIDs_vec[0];
00550 
00551       for(int k=0; k<numNodesPerCR; k++) {
00552         const NodeDescriptor& node = Filter::findNodeDescriptor(nodeIDsPtr[k]);
00553         nodeList[k] = node.getNodeNumber();
00554 
00555         int eqn = -1;
00556         if (!node.getFieldEqnNumber(fieldIDs[k], eqn)) voidERReturn;
00557 
00558         eqnList[k] = eqn;
00559       }
00560 
00561       int crPenID = crset.getConstraintID() + i;
00562       err = lsc_->setPenCREqns(crPenID, 1, numNodesPerCR,
00563                                &nodeList, &eqnList,
00564                                fieldIDs);
00565       if (err) voidERReturn;
00566       ++cr_iter;
00567    }
00568 }
00569 
00570 //==============================================================================
00571 int LinSysCoreFilter::storeNodalColumnCoefs(int eqn, const NodeDescriptor& node,
00572                                             int fieldID, int fieldSize,
00573                                             double* coefs)
00574 {
00575   //
00576   //This function stores the coeficients for 'node' at 'fieldID' at the correct
00577   //column indices in row 'eqn' of the system matrix.
00578   //
00579   if ((localStartRow_ > eqn) || (eqn > localEndRow_)) return(0);
00580 
00581   int eqnNumber = -1;
00582   if (!node.getFieldEqnNumber(fieldID, eqnNumber)) ERReturn(FEI_ID_NOT_FOUND);
00583 
00584   int numParams = fieldSize;
00585   if (numParams < 1) {
00586     return(0);
00587   }
00588 
00589   if ((int)iworkSpace2_.size() < numParams) {
00590     iworkSpace2_.resize(numParams);
00591   }
00592 
00593   int* cols = &iworkSpace2_[0];
00594 
00595   for(int j=0; j<numParams; j++) {
00596     cols[j] = eqnNumber + j;
00597   }
00598 
00599   CHK_ERR( giveToMatrix(1, &eqn, numParams, cols, &coefs, ASSEMBLE_SUM) );
00600 
00601   return(FEI_SUCCESS);
00602 }
00603 
00604 //==============================================================================
00605 int LinSysCoreFilter::storeNodalRowCoefs(const NodeDescriptor& node,
00606                                          int fieldID, int fieldSize,
00607                                          double* coefs, int eqn)
00608 {
00609   //
00610   //This function stores coeficients in the equations for 'node', 'fieldID' at
00611   //column index 'eqn' of the system matrix.
00612   //
00613   int eqnNumber = -1;
00614   if (!node.getFieldEqnNumber(fieldID, eqnNumber)) ERReturn(FEI_ID_NOT_FOUND);
00615 
00616   if ((localStartRow_ > eqnNumber) || (eqnNumber > localEndRow_)) return(0);
00617 
00618   int numParams = fieldSize;
00619 
00620   if (numParams < 1) {
00621     return(0);
00622   }
00623 
00624   if ((int)iworkSpace2_.size() < numParams) {
00625     iworkSpace2_.resize(numParams);
00626   }
00627 
00628   int* ptRows = &iworkSpace2_[0];
00629 
00630   if ((int)dworkSpace2_.size() < numParams) {
00631     dworkSpace2_.resize(numParams);
00632   }
00633 
00634   const double* * values = &dworkSpace2_[0];
00635 
00636   for(int j=0; j<numParams; j++) {
00637     ptRows[j] = eqnNumber + j;
00638     values[j] = &(coefs[j]);
00639   }
00640 
00641   CHK_ERR( giveToMatrix(numParams, ptRows, 1, &eqn, values, ASSEMBLE_SUM) );
00642 
00643   return(FEI_SUCCESS);
00644 }
00645 
00646 //==============================================================================
00647 void LinSysCoreFilter::storeNodalSendEqn(const NodeDescriptor& node,
00648                                          int fieldID, int col,
00649                                          double* coefs)
00650 {
00651   //
00652   //This is a private LinSysCoreFilter function. We can safely assume that
00653   //it will only be called with a node that is not locally owned.
00654   //
00655   //This function tells the eqn comm mgr that we'll be sending contributions
00656   //to column 'col' for the equations associated with 'fieldID', on 'node', on
00657   //node's owning processor.
00658   //
00659   int proc = node.getOwnerProc();
00660 
00661   int eqnNumber = -1;
00662   if (!node.getFieldEqnNumber(fieldID, eqnNumber)) voidERReturn;
00663 
00664   int numEqns = problemStructure_->getFieldSize(fieldID);
00665 
00666   for(int i=0; i<numEqns; i++) {
00667     eqnCommMgr_->addRemoteEqn(eqnNumber+i, proc, &coefs[i], &col, 1);
00668   }
00669 }
00670 
00671 //==============================================================================
00672 void LinSysCoreFilter::storePenNodeSendData(const NodeDescriptor& iNode,
00673                                             int iField, int iFieldSize,
00674                                             double* iCoefs,
00675                                             const NodeDescriptor& jNode,
00676                                             int jField, int jFieldSize,
00677                                             double* jCoefs,
00678                                             double penValue, double CRValue)
00679 {
00680 //
00681 //This function will register with the eqn comm mgr the equations associated
00682 //with iNode, field 'iField' having column indices that are the equations
00683 //associated with jNode, field 'jField', to be sent to the owner of iNode.
00684 //
00685    int proc = iNode.getOwnerProc();
00686 
00687    int iEqn = -1, jEqn = -1;
00688    if (!iNode.getFieldEqnNumber(iField, iEqn)) voidERReturn;
00689    if (!jNode.getFieldEqnNumber(jField, jEqn)) voidERReturn;
00690 
00691    int iNumParams = iFieldSize;
00692    int jNumParams = jFieldSize;
00693    if (iNumParams < 1 || jNumParams < 1) {
00694      FEI_CERR << "FEI ERROR, attempt to store indices for field with non-positive size"
00695           << " field "<<iField<<", size "<<iNumParams<<", field "<<jField<<", size "
00696           << jNumParams<<FEI_ENDL;
00697      voidERReturn;
00698    }
00699 
00700    if ((int)dworkSpace_.size() < jNumParams) {
00701      dworkSpace_.resize(jNumParams);
00702    }
00703 
00704    double* coefs = &dworkSpace_[0];
00705 
00706    if ((int)iworkSpace2_.size() < jNumParams) {
00707      iworkSpace2_.resize(jNumParams);
00708    }
00709 
00710    int* cols = &iworkSpace2_[0];
00711 
00712    for(int i=0; i<iNumParams; i++) {
00713       for(int j=0; j<jNumParams; j++) {
00714          cols[j] = jEqn + j;
00715          coefs[j] = penValue*iCoefs[i]*jCoefs[j];
00716       }
00717 
00718       int row = iEqn + i;
00719       eqnCommMgr_->addRemoteEqn(row, proc, coefs, cols, jNumParams);
00720 
00721       double rhsValue = penValue*iCoefs[i]*CRValue;
00722       eqnCommMgr_->addRemoteRHS(row, proc, currentRHS_, rhsValue);
00723    }
00724 }
00725 
00726 //==============================================================================
00727 int LinSysCoreFilter::storePenNodeData(const NodeDescriptor& iNode,
00728                                        int iField, int iFieldSize,
00729                                        double* iCoefs,
00730                                        const NodeDescriptor& jNode,
00731                                        int jField, int jFieldSize,
00732                                        double* jCoefs,
00733                                        double penValue, double CRValue){
00734 //
00735 //This function will add to the local matrix the penalty constraint equations
00736 //associated with iNode at iField, having column indices that are the
00737 //equations associated with jNode at jField.
00738 //Also, add the penalty contribution to the RHS vector.
00739 //
00740    int iEqn = -1, jEqn = -1;
00741    if (!iNode.getFieldEqnNumber(iField, iEqn)) ERReturn(FEI_ID_NOT_FOUND);
00742    if (!jNode.getFieldEqnNumber(jField, jEqn)) ERReturn(FEI_ID_NOT_FOUND);
00743 
00744    int iNumParams = iFieldSize;
00745    int jNumParams = jFieldSize;
00746    if (iNumParams < 1 || jNumParams < 1) {
00747      FEI_CERR << "FEI ERROR, attempt to store indices for field with non-positive size"
00748           << " field "<<iField<<", size "<<iNumParams<<", field "<<jField<<", size "
00749           << jNumParams<<FEI_ENDL;
00750      ERReturn(-1);
00751    }
00752 
00753    if ((int)dworkSpace2_.size() < iNumParams) {
00754      dworkSpace2_.resize(iNumParams);
00755    }
00756 
00757    const double* * coefs = &dworkSpace2_[0];
00758 
00759    if ((int)dworkSpace_.size() < iNumParams * jNumParams) {
00760      dworkSpace_.resize(iNumParams * jNumParams);
00761    }
00762 
00763    double* coefList = &dworkSpace_[0];
00764 
00765    if ((int)iworkSpace2_.size() < jNumParams+iNumParams) {
00766      iworkSpace2_.resize(jNumParams+iNumParams);
00767    }
00768 
00769    int* cols = &iworkSpace2_[0];
00770    int* rows = &iworkSpace2_[jNumParams];
00771 
00772 
00773    for(int i=0; i<iNumParams; i++) {
00774       double* coefPtr = coefList + i*jNumParams;
00775       coefs[i] = coefPtr;
00776 
00777       rows[i] = iEqn + i;
00778 
00779       for(int j=0; j<jNumParams; j++) {
00780          if (i==0) cols[j] = jEqn + j;
00781          coefPtr[j] = penValue*iCoefs[i]*jCoefs[j];
00782       }
00783 
00784       double rhsValue = penValue*iCoefs[i]*CRValue;
00785       CHK_ERR( giveToRHS(1, &rhsValue, &rows[i], ASSEMBLE_SUM) );
00786    }
00787 
00788    CHK_ERR( giveToMatrix(iNumParams, rows,
00789                          jNumParams, cols, coefs, ASSEMBLE_SUM) );
00790 
00791    return(FEI_SUCCESS);
00792 }
00793 
00794 //------------------------------------------------------------------------------
00795 int LinSysCoreFilter::resetSystem(double s)
00796 {
00797   //
00798   //  This puts the value s throughout both the matrix and the vector.
00799   //
00800   if (Filter::logStream() != NULL) {
00801     (*logStream()) << "FEI: resetSystem" << FEI_ENDL << s << FEI_ENDL;
00802   }
00803 
00804   CHK_ERR( resetTheMatrix(s) );
00805   CHK_ERR( resetTheRHSVector(s) );
00806 
00807   debugOutput("#LinSysCoreFilter  leaving resetSystem");
00808 
00809   return(FEI_SUCCESS);
00810 }
00811 
00812 //------------------------------------------------------------------------------
00813 int LinSysCoreFilter::deleteMultCRs()
00814 {
00815   debugOutput("#LinSysCoreFilter::deleteMultCRs");
00816 
00817   LinSysCore_flexible* lscf = dynamic_cast<LinSysCore_flexible*>(lsc_);
00818   if (lscf == NULL) {
00819 //    FEI_CERR << "FEI::LinSysCoreFilter: ERROR deleteMultCRs requested, but "
00820 //         << "underlying solver doesn't support this operation." << FEI_ENDL;
00821     return(-1);
00822   }
00823 
00824   int err = lscf->resetConstraints(0.0);
00825 
00826   debugOutput("#LinSysCoreFilter leaving deleteMultCRs");
00827 
00828   return(err);
00829 }
00830 
00831 //------------------------------------------------------------------------------
00832 int LinSysCoreFilter::resetTheMatrix(double s)
00833 {
00834   CHK_ERR( lsc_->resetMatrix(s) );
00835 
00836   return(FEI_SUCCESS);
00837 }
00838 
00839 //------------------------------------------------------------------------------
00840 int LinSysCoreFilter::resetTheRHSVector(double s)
00841 {
00842   CHK_ERR( lsc_->resetRHSVector(s) );
00843 
00844   return(FEI_SUCCESS);
00845 }
00846 
00847 //------------------------------------------------------------------------------
00848 int LinSysCoreFilter::resetMatrix(double s) {
00849 //
00850 //  This puts the value s throughout both the matrix and the vector.
00851 //
00852 
00853     debugOutput("FEI: resetMatrix");
00854 
00855     CHK_ERR( resetTheMatrix(s) )
00856 
00857    //clear away any boundary condition data.
00858    bcManager_->clearAllBCs();
00859 
00860     eqnCommMgr_->resetCoefs();
00861 
00862     debugOutput("#LinSysCoreFilter leaving resetMatrix");
00863 
00864     return(FEI_SUCCESS);
00865 }
00866 
00867 //------------------------------------------------------------------------------
00868 int LinSysCoreFilter::resetRHSVector(double s) {
00869 //
00870 //  This puts the value s throughout the rhs vector.
00871 //
00872 
00873     debugOutput("FEI: resetRHSVector");
00874 
00875     CHK_ERR( resetTheRHSVector(s) )
00876 
00877     //clear away any boundary condition data.
00878     bcManager_->clearAllBCs();
00879 
00880     eqnCommMgr_->resetCoefs();
00881 
00882     debugOutput("# LinSysCoreFilter leaving resetRHSVector");
00883 
00884     return(FEI_SUCCESS);
00885 }
00886 
00887 //------------------------------------------------------------------------------
00888 int LinSysCoreFilter::resetInitialGuess(double s) {
00889 //
00890 //  This puts the value s throughout the initial guess (solution) vector.
00891 //
00892   if (Filter::logStream() != NULL) {
00893     FEI_OSTREAM& os = *logStream();
00894     os << "FEI: resetInitialGuess" << FEI_ENDL;
00895     os << "#value to which initial guess is to be set" << FEI_ENDL;
00896     os << s << FEI_ENDL;
00897   }
00898 
00899   int* eqns = new int[numReducedRows_];
00900   double* values = new double[numReducedRows_];
00901   if (eqns == NULL || values == NULL) return(-1);
00902 
00903   for(int i=0; i<numReducedRows_; ++i) {
00904     eqns[i] = reducedStartRow_ + i;
00905     values[i] = s;
00906   }
00907 
00908   CHK_ERR( lsc_->putInitialGuess(eqns, values, numReducedRows_) );
00909 
00910   delete [] eqns;
00911   delete [] values;
00912 
00913   debugOutput("# LinSysCoreFilter leaving resetInitialGuess");
00914 
00915   return(FEI_SUCCESS);
00916 }
00917 
00918 //------------------------------------------------------------------------------
00919 int LinSysCoreFilter::loadNodeBCs(int numNodes,
00920                                   const GlobalID *nodeIDs,
00921                                   int fieldID,
00922                                   const int* offsetsIntoField,
00923                                   const double* prescribedValues)
00924 {
00925   //
00926   //  load boundary condition information for a given set of nodes
00927   //
00928   int size = problemStructure_->getFieldSize(fieldID);
00929   if (size < 1) {
00930     FEI_CERR << "FEI Warning: loadNodeBCs called for fieldID "<<fieldID
00931          <<", which was defined with size "<<size<<" (should be positive)."<<FEI_ENDL;
00932     return(0);
00933   }
00934 
00935   if (Filter::logStream() != NULL) {
00936     (*logStream())<<"FEI: loadNodeBCs"<<FEI_ENDL
00937                      <<"#num-nodes"<<FEI_ENDL<<numNodes<<FEI_ENDL
00938                      <<"#fieldID"<<FEI_ENDL<<fieldID<<FEI_ENDL
00939                      <<"#field-size"<<FEI_ENDL<<size<<FEI_ENDL;
00940     (*logStream())<<"#following lines: nodeID offsetIntoField value "<<FEI_ENDL;
00941 
00942     for(int j=0; j<numNodes; j++) {
00943       GlobalID nodeID = nodeIDs[j];
00944       (*logStream())<<static_cast<int>(nodeID)<<"  ";
00945       (*logStream())<<offsetsIntoField[j]<<"  "<<prescribedValues[j];
00946       (*logStream())<<FEI_ENDL;
00947     }
00948   }
00949 
00950   if (numNodes > 0) newBCData_ = true;
00951 
00952   bcManager_->addBCRecords(numNodes, nodeIDType_, fieldID, nodeIDs,
00953                            offsetsIntoField, prescribedValues);
00954 
00955   return(FEI_SUCCESS);
00956 }
00957 
00958 //------------------------------------------------------------------------------
00959 int LinSysCoreFilter::loadElemBCs(int numElems,
00960                                   const GlobalID *elemIDs,
00961                                   int fieldID,
00962                                   const double *const *alpha,
00963                                   const double *const *beta,
00964                                   const double *const *gamma)
00965 {
00966    return(-1);
00967 }
00968 
00969 //------------------------------------------------------------------------------
00970 void LinSysCoreFilter::allocElemStuff()
00971 {
00972    int nb = problemStructure_->getNumElemBlocks();
00973 
00974    for(int i=0; i<nb; i++) {
00975      BlockDescriptor* block = NULL;
00976      int err = problemStructure_->getBlockDescriptor_index(i, block);
00977      if (err) voidERReturn;
00978 
00979       int numEqns = block->getNumEqnsPerElement();
00980       if (maxElemRows_ < numEqns) maxElemRows_ = numEqns;
00981    }
00982 
00983    scatterIndices_.resize(maxElemRows_);
00984 
00985    if (eStiff_ != NULL) delete [] eStiff_;
00986    if (eStiff1D_ != NULL) delete [] eStiff1D_;
00987    if (eLoad_ != NULL) delete [] eLoad_;
00988 
00989    eStiff_ = new double*[maxElemRows_];
00990    eStiff1D_ = new double[maxElemRows_*maxElemRows_];
00991 
00992    if (eStiff_ == NULL || eStiff1D_ == NULL) voidERReturn
00993 
00994    for(int r=0; r<maxElemRows_; r++) {
00995       eStiff_[r] = eStiff1D_ + r*maxElemRows_;
00996    }
00997 
00998    eLoad_ = new double[maxElemRows_];
00999 
01000    if (eLoad_ == NULL) voidERReturn
01001 }
01002 
01003 //------------------------------------------------------------------------------
01004 int LinSysCoreFilter::sumInElem(GlobalID elemBlockID,
01005                                 GlobalID elemID,
01006                                 const GlobalID* elemConn,
01007                                 const double* const* elemStiffness,
01008                                 const double* elemLoad,
01009                                 int elemFormat)
01010 {
01011   if (Filter::logStream() != NULL && outputLevel_ > 2) {
01012     (*logStream()) << "FEI: sumInElem" << FEI_ENDL <<"#blkID" << FEI_ENDL
01013                       << static_cast<int>(elemBlockID) << FEI_ENDL
01014                       << "#elID" << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
01015     BlockDescriptor* block = NULL;
01016     CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
01017     int numNodes = block->getNumNodesPerElement();
01018     (*logStream()) << "#n-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
01019     (*logStream()) << "#nodes" << FEI_ENDL;
01020     for(int i=0; i<numNodes; ++i) {
01021       GlobalID nodeID = elemConn[i];
01022       (*logStream())<<static_cast<int>(nodeID)<<" ";
01023     }
01024     (*logStream())<<FEI_ENDL;
01025   }
01026 
01027   return(generalElemInput(elemBlockID, elemID, elemConn, elemStiffness,
01028                           elemLoad, elemFormat));
01029 }
01030 
01031 //------------------------------------------------------------------------------
01032 int LinSysCoreFilter::sumInElemMatrix(GlobalID elemBlockID,
01033                                       GlobalID elemID,
01034                                       const GlobalID* elemConn,
01035                                       const double* const* elemStiffness,
01036                                       int elemFormat)
01037 {
01038   if (Filter::logStream() != NULL && outputLevel_ > 2) {
01039     (*logStream()) << "FEI: sumInElemMatrix"<<FEI_ENDL
01040                       << "#blkID" << FEI_ENDL << static_cast<int>(elemBlockID) << FEI_ENDL
01041                       << "#elID" << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
01042     BlockDescriptor* block = NULL;
01043     CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
01044     int numNodes = block->getNumNodesPerElement();
01045     (*logStream()) << "#n-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
01046     (*logStream()) << "#nodes" << FEI_ENDL;
01047     for(int i=0; i<numNodes; ++i) {
01048       GlobalID nodeID = elemConn[i];
01049       (*logStream())<<static_cast<int>(nodeID)<<" ";
01050     }
01051     (*logStream())<<FEI_ENDL;
01052   }
01053 
01054   return(generalElemInput(elemBlockID, elemID, elemConn, elemStiffness,
01055                           NULL, elemFormat));
01056 }
01057 
01058 //------------------------------------------------------------------------------
01059 int LinSysCoreFilter::sumInElemRHS(GlobalID elemBlockID,
01060                            GlobalID elemID,
01061                            const GlobalID* elemConn,
01062                            const double* elemLoad)
01063 {
01064   if (Filter::logStream() != NULL && outputLevel_ > 2) {
01065     (*logStream()) << "FEI: sumInElemRHS"<<FEI_ENDL<<"#blID" << FEI_ENDL
01066                       <<(int)elemBlockID << FEI_ENDL
01067                       << "#elID" << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
01068     BlockDescriptor* block = NULL;
01069     CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
01070     int numNodes = block->getNumNodesPerElement();
01071     (*logStream()) << "#n-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
01072     (*logStream()) << "#nodes" << FEI_ENDL;
01073     for(int i=0; i<numNodes; ++i) {
01074       GlobalID nodeID = elemConn[i];
01075       (*logStream())<<static_cast<int>(nodeID)<<" ";
01076     }
01077     (*logStream())<<FEI_ENDL;
01078   }
01079 
01080   return(generalElemInput(elemBlockID, elemID, elemConn, NULL,
01081                           elemLoad, -1));
01082 }
01083 
01084 //------------------------------------------------------------------------------
01085 int LinSysCoreFilter::generalElemInput(GlobalID elemBlockID,
01086                                        GlobalID elemID,
01087                                        const GlobalID* elemConn,
01088                                        const double* const* elemStiffness,
01089                                        const double* elemLoad,
01090                                        int elemFormat)
01091 {
01092   (void)elemConn;
01093   return(generalElemInput(elemBlockID, elemID, elemStiffness, elemLoad,
01094                           elemFormat) );
01095 }
01096 
01097 //------------------------------------------------------------------------------
01098 int LinSysCoreFilter::generalElemInput(GlobalID elemBlockID,
01099                                        GlobalID elemID,
01100                                        const double* const* elemStiffness,
01101                                        const double* elemLoad,
01102                                        int elemFormat)
01103 {
01104   //first get the block-descriptor for this elemBlockID...
01105 
01106   BlockDescriptor* block = NULL;
01107   CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
01108 
01109   //now allocate our local stiffness/load copy if we haven't already.
01110 
01111   if (maxElemRows_ <= 0) allocElemStuff();
01112 
01113   int numElemRows = block->getNumEqnsPerElement();
01114   int numBlkElemRows = block->getNumBlkEqnsPerElement();
01115   int interleave = block->getInterleaveStrategy();
01116 
01117   //an std::vector.resize operation is free if the size is either shrinking or
01118   //staying the same.
01119 
01120   scatterIndices_.resize(numElemRows);
01121   blkScatterIndices_.resize(numBlkElemRows*2);
01122 
01123   const double* const* stiff = NULL;
01124   if (elemStiffness != NULL) stiff = elemStiffness;
01125 
01126   const double* load = NULL;
01127   if (elemLoad != NULL) load = elemLoad;
01128 
01129   //we'll make a local dense copy of the element stiffness array
01130   //if the stiffness array was passed in using elemFormat != FEI_DENSE_ROW
01131   //AND if the stiffness array is non-null.
01132   //Note that this is not optimal if elemFormat is one of the
01133   //diagonal or block-diagonal formats.
01134   if (elemFormat != FEI_DENSE_ROW && stiff != NULL) {
01135     if (elemFormat == FEI_BLOCK_DIAGONAL_ROW ||
01136         elemFormat == FEI_BLOCK_DIAGONAL_COL) {
01137       FEI_CERR << "LinSysCoreFilter::generalElemInput ERROR, elemFormat="
01138                << elemFormat << " not supported."<<FEI_ENDL;
01139       ERReturn(-1);
01140     }
01141 
01142     Filter::copyStiffness(stiff, numElemRows, elemFormat, eStiff_);
01143     stiff = eStiff_;
01144   }
01145 
01146   if (stiff != NULL) newMatrixData_ = true;
01147   if (load != NULL) newVectorData_ = true;
01148 
01149   if (Filter::logStream() != NULL && outputLevel_ > 2) {
01150     FEI_OSTREAM& os = *logStream();
01151 
01152     if (stiff != NULL || load != NULL) {
01153       os << "#numRows"<< FEI_ENDL << numElemRows << FEI_ENDL;
01154     }
01155 
01156     if (stiff != NULL) {
01157       os << "#elem-stiff (after copy into dense-row format)" << FEI_ENDL;
01158       for(int i=0; i<numElemRows; i++) {
01159         const double* stiff_i = stiff[i];
01160         for(int j=0; j<numElemRows; j++) {
01161           os << stiff_i[j] << " ";
01162         }
01163         os << FEI_ENDL;
01164       }
01165     }
01166 
01167     if (load != NULL) {
01168       os << "#elem-load" << FEI_ENDL;
01169       for(int i=0; i<numElemRows; i++) {
01170         os << load[i] << " ";
01171       }
01172       os<<FEI_ENDL;
01173     }
01174 
01175     if (stiff != NULL) {
01176       os << "#elemformat" << FEI_ENDL << elemFormat << FEI_ENDL;
01177     }
01178   }
01179 
01180   //now let's obtain the scatter indices for assembling the equations
01181   //into their appropriate places in the global matrix and rhs vectors
01182 
01183   int* indPtr = &scatterIndices_[0];
01184   int* blkIndPtr = &blkScatterIndices_[0];
01185   int* blkSizesPtr = blkIndPtr + numBlkElemRows;
01186 
01187   bool useBlkEqns = false;
01188   if (interleave == 0) {
01189     //interleave==0 is node-major, so we'll get the block-indices too.
01190     problemStructure_->getScatterIndices_ID(elemBlockID, elemID,
01191                                             interleave, indPtr,
01192                                             blkIndPtr, blkSizesPtr);
01193     int sumBlkSizes = 0;
01194     for(int ii=0; ii<numBlkElemRows; ++ii) {
01195       sumBlkSizes += blkSizesPtr[ii];
01196     }
01197     if (sumBlkSizes == numElemRows) {
01198       useBlkEqns = true;
01199     }
01200   }
01201   else {
01202     //interleave!=0 is field-major, and we'll only bother with point-indices.
01203     problemStructure_->getScatterIndices_ID(elemBlockID, elemID,
01204                                             interleave, indPtr);
01205   }
01206 
01207   if (stiff != NULL) {
01208     if (problemStructure_->numSlaveEquations() == 0) {
01209       //I'm not checking the return-value (error-code) on this call, because
01210       //I wasn't even calling it until recently, and I'm not sure if all
01211       //LinearSystemCore implementations even have it implemented.
01212       lsc_->setStiffnessMatrices(elemBlockID, 1, &elemID,
01213                                  &stiff, scatterIndices_.size(),
01214                                  &indPtr);
01215     }
01216 
01217     int len = scatterIndices_.size();
01218     if (interleave == 0) {
01219       CHK_ERR( assembleEqns( len, len, indPtr, indPtr, stiff, true,
01220                              numBlkElemRows, blkIndPtr,
01221                              blkSizesPtr, useBlkEqns, ASSEMBLE_SUM ) );
01222     }
01223     else {
01224       CHK_ERR( assembleEqns( len, len, indPtr, indPtr, stiff, true,
01225                              numBlkElemRows, blkIndPtr,
01226                              blkSizesPtr, false, ASSEMBLE_SUM ) );
01227     }
01228   }
01229 
01230   if (load != NULL) {
01231     if (problemStructure_->numSlaveEquations() == 0) {
01232       //I'm not checking the return-value (error-code) on this call, because
01233       //I wasn't even calling it until recently, and I'm not sure if all
01234       //LinearSystemCore implementations even have it implemented.
01235       lsc_->setLoadVectors(elemBlockID, 1, &elemID,
01236                            &load, scatterIndices_.size(),
01237                            &indPtr);
01238     }
01239 
01240     CHK_ERR( assembleRHS(scatterIndices_.size(), indPtr, load, ASSEMBLE_SUM) );
01241   }
01242 
01243   return(FEI_SUCCESS);
01244 }
01245 
01246 //------------------------------------------------------------------------------
01247 int LinSysCoreFilter::putIntoRHS(int IDType,
01248                           int fieldID,
01249                           int numIDs,
01250                           const GlobalID* IDs,
01251                           const double* rhsEntries)
01252 {
01253   int fieldSize = problemStructure_->getFieldSize(fieldID);
01254 
01255   rowIndices_.resize(fieldSize*numIDs);
01256   int checkNumEqns;
01257 
01258   CHK_ERR( problemStructure_->getEqnNumbers(numIDs, IDs, IDType, fieldID,
01259                                             checkNumEqns,
01260                                             &rowIndices_[0]));
01261   if (checkNumEqns != numIDs*fieldSize) {
01262     ERReturn(-1);
01263   }
01264 
01265   CHK_ERR( exchangeRemoteEquations() );
01266 
01267   CHK_ERR(assembleRHS(rowIndices_.size(), &rowIndices_[0], rhsEntries, ASSEMBLE_PUT));
01268 
01269   return(0);
01270 }
01271 
01272 //------------------------------------------------------------------------------
01273 int LinSysCoreFilter::sumIntoRHS(int IDType,
01274                           int fieldID,
01275                           int numIDs,
01276                           const GlobalID* IDs,
01277                           const double* rhsEntries)
01278 {
01279   int fieldSize = problemStructure_->getFieldSize(fieldID);
01280 
01281   rowIndices_.resize(fieldSize*numIDs);
01282   int checkNumEqns;
01283 
01284   CHK_ERR( problemStructure_->getEqnNumbers(numIDs, IDs, IDType, fieldID,
01285                                             checkNumEqns,
01286                                             &rowIndices_[0]));
01287   if (checkNumEqns != numIDs*fieldSize) {
01288     ERReturn(-1);
01289   }
01290 
01291   CHK_ERR(assembleRHS(rowIndices_.size(), &rowIndices_[0], rhsEntries, ASSEMBLE_SUM));
01292 
01293   return(0);
01294 }
01295 
01296 //------------------------------------------------------------------------------
01297 int LinSysCoreFilter::implementAllBCs() {
01298 //
01299 // This function will handle the modifications to the stiffness and load
01300 // necessary to enforce nodal boundary conditions.
01301 //
01302    debugOutput("# implementAllBCs");
01303 
01304    std::vector<int> essEqns;
01305    std::vector<double> essGamma;
01306 
01307    EqnBuffer bcEqns;
01308 
01309    CHK_ERR( bcManager_->finalizeBCEqns(bcEqns) );
01310 
01311    if (resolveConflictRequested_) {
01312      CHK_ERR( resolveConflictingCRs(bcEqns) );
01313    }
01314 
01315    CHK_ERR( eqnCommMgr_->gatherSharedBCs(bcEqns) );
01316 
01317    //now separate the boundary-condition equations into arrays
01318    fei::FillableMat bcEqns_mat(bcEqns);
01319    fei::impl_utils::separate_BC_eqns(bcEqns_mat, essEqns, essGamma);
01320 
01321    std::vector<double> essAlpha(essEqns.size(), 1);
01322 
01323    exchangeRemoteBCs(essEqns, essAlpha, essGamma);
01324 
01325    if (essEqns.size() > 0) {
01326       CHK_ERR( enforceEssentialBCs(&essEqns[0],
01327                                    &essAlpha[0],
01328                                    &essGamma[0], essEqns.size()) );
01329    }
01330 
01331    debugOutput("#LinSysCoreFilter leaving implementAllBCs");
01332    return(FEI_SUCCESS);
01333 }
01334 
01335 //------------------------------------------------------------------------------
01336 int LinSysCoreFilter::enforceEssentialBCs(const int* eqns,
01337                                           const double* alpha,
01338                                           const double* gamma,
01339                                           int numEqns)
01340 {
01341   int* cc_eqns = const_cast<int*>(eqns);
01342   double* cc_alpha = const_cast<double*>(alpha);
01343   double* cc_gamma = const_cast<double*>(gamma);
01344 
01345   if (problemStructure_->numSlaveEquations() == 0) {
01346     CHK_ERR( lsc_->enforceEssentialBC(cc_eqns,
01347                                       cc_alpha, cc_gamma,
01348                                       numEqns) );
01349   }
01350   else {
01351     std::vector<int> reducedEqns(numEqns);
01352     for(int i=0; i<numEqns; i++) {
01353       problemStructure_->translateToReducedEqn(eqns[i], reducedEqns[i]);
01354     }
01355 
01356     CHK_ERR( lsc_->enforceEssentialBC(&reducedEqns[0], cc_alpha, cc_gamma, numEqns) );
01357   }
01358 
01359   return(FEI_SUCCESS);
01360 }
01361 
01362 //------------------------------------------------------------------------------
01363 int LinSysCoreFilter::enforceRemoteEssBCs(int numEqns, const int* eqns, 
01364                            const int* const* colIndices, const int* colIndLens,
01365                            const double* const* BCcoefs)
01366 {
01367   //These eqn-numbers were reduced to the slave-eqn space by
01368   //LinSysCore::exchangeRemoteBCs, which is the only function that calls THIS
01369   //function, so we can simply pass them straight on in to LinearSystemCore.
01370   //
01371 
01372   int* cc_eqns = const_cast<int*>(eqns);
01373   int** cc_colIndices = const_cast<int**>(colIndices);
01374   int* cc_colIndLens = const_cast<int*>(colIndLens);
01375   double** cc_BCcoefs = const_cast<double**>(BCcoefs);
01376 
01377   CHK_ERR( lsc_->enforceRemoteEssBCs(numEqns, cc_eqns, cc_colIndices,
01378                                      cc_colIndLens, cc_BCcoefs) );
01379 
01380   return(FEI_SUCCESS);
01381 }
01382 
01383 //------------------------------------------------------------------------------
01384 int LinSysCoreFilter::resolveConflictingCRs(EqnBuffer& bcEqns)
01385 {
01386   int numMultCRs = problemStructure_->getNumMultConstRecords();
01387   if (numMultCRs < 1) {
01388     return(0);
01389   }
01390 
01391   std::map<GlobalID,ConstraintType*>::const_iterator
01392     cr_iter = problemStructure_->getMultConstRecords().begin(),
01393     cr_end  = problemStructure_->getMultConstRecords().end();
01394 
01395   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
01396 
01397   std::vector<int>& bcEqnNumbers = bcEqns.eqnNumbers();
01398 
01399   double coefs[3];
01400   int indices[3];
01401   indices[0] = 0;
01402   indices[1] = 1;
01403   indices[2] = 2;
01404 
01405   double fei_eps = 1.e-49;
01406 
01407   while(cr_iter != cr_end) {
01408     ConstraintType& multCR = *((*cr_iter).second);
01409 
01410     int lenList = multCR.getMasters().size();
01411 
01412     std::vector<GlobalID>& CRNode_vec = multCR.getMasters();
01413     GlobalID *CRNodePtr = &CRNode_vec[0];
01414     std::vector<int>& CRField_vec = multCR.getMasterFieldIDs();
01415     int* CRFieldPtr = &CRField_vec[0];
01416     std::vector<double>& weights_vec = multCR.getMasterWeights();
01417     double* weights = &weights_vec[0];
01418 
01419     int offset = 0;
01420     for(int j=0; j<lenList; ++j) {
01421       int fieldSize = problemStructure_->getFieldSize(CRFieldPtr[j]);
01422       for(int k=0; k<fieldSize; ++k) {
01423         if (std::abs(weights[offset++] + 1.0) < fei_eps) {
01424           const NodeDescriptor* node = NULL;
01425           CHK_ERR( nodeDB.getNodeWithID(CRNodePtr[j], node) );
01426           int eqn = 0;
01427           node->getFieldEqnNumber(CRFieldPtr[j], eqn);
01428           eqn += k;
01429 
01430           if (fei::binarySearch(eqn, bcEqnNumbers) >= 0) {
01431             coefs[0] = 1.0;
01432             coefs[1] = 0.0;
01433             coefs[2] = 1.0;
01434             int crEqn = multCR.getEqnNumber();
01435             CHK_ERR( bcEqns.addEqn(crEqn, coefs, indices, 3, false) );
01436           }
01437         }
01438       }
01439     }
01440     ++cr_iter;
01441   }
01442 
01443   return(0);
01444 }
01445 
01446 //------------------------------------------------------------------------------
01447 int LinSysCoreFilter::exchangeRemoteEquations()
01448 {
01449   //
01450   // This function is where processors send local contributions to remote
01451   // equations to the owners of those equations, and receive remote
01452   // contributions to local equations.
01453   //
01454   debugOutput("#LinSysCoreFilter::exchangeRemoteEquations");
01455 
01456   if (reducedEqnCounter_ > 0) CHK_ERR( assembleReducedEqns() );
01457 
01458   if (reducedRHSCounter_ > 0) CHK_ERR( assembleReducedRHS() );
01459 
01460   int len = 4;
01461   std::vector<int> flags(len), globalFlags(len);
01462   flags[0] = newMatrixData_ ? 1 : 0;
01463   flags[1] = newVectorData_ ? 1 : 0;
01464   flags[2] = newConstraintData_ ? 1 : 0;
01465   flags[3] = newBCData_ ? 1 : 0;
01466 
01467   CHK_ERR( fei::GlobalMax(comm_, flags, globalFlags) );
01468 
01469   newMatrixData_     = globalFlags[0] > 0 ? true : false;
01470   newVectorData_     = globalFlags[1] > 0 ? true : false;
01471   newConstraintData_ = globalFlags[2] > 0 ? true : false;
01472   newBCData_         = globalFlags[3] > 0 ? true : false;
01473 
01474   if (newMatrixData_ || newVectorData_ || newConstraintData_) { 
01475 
01476     CHK_ERR( eqnCommMgr_->exchangeEqns(logStream()) );
01477 
01478     needToCallMatrixLoadComplete_ = true;
01479 
01480     //so now the remote contributions should be available, let's get them out
01481     //of the eqn comm mgr and put them into our local matrix structure.
01482 
01483     debugOutput("#   putting remote contributions into linear system...");
01484 
01485     CHK_ERR( unpackRemoteContributions(*eqnCommMgr_, ASSEMBLE_SUM) );
01486 
01487     eqnCommMgr_->resetCoefs();
01488 
01489     newMatrixData_ = false;
01490     newVectorData_ = false;
01491     newConstraintData_ = false;
01492   }
01493 
01494   firstRemEqnExchange_ = false;
01495 
01496   debugOutput("#LinSysCoreFilter leaving exchangeRemoteEquations");
01497 
01498   return(FEI_SUCCESS);
01499 }
01500 
01501 //------------------------------------------------------------------------------
01502 int LinSysCoreFilter::unpackRemoteContributions(EqnCommMgr& eqnCommMgr,
01503                                                 int assemblyMode)
01504 {
01505   int numRecvEqns = eqnCommMgr.getNumLocalEqns();
01506   std::vector<int>& recvEqnNumbers = eqnCommMgr.localEqnNumbers();
01507   std::vector<fei::CSVec*>& recvEqns = eqnCommMgr.localEqns();
01508   std::vector<std::vector<double>*>& recvRHSs = *(eqnCommMgr.localRHSsPtr());
01509 
01510   bool newCoefs = eqnCommMgr.newCoefData();
01511   bool newRHSs = eqnCommMgr.newRHSData();
01512 
01513   int i;
01514   double** coefs = new double*[numRecvEqns];
01515 
01516   for(i=0; i<numRecvEqns; i++) {
01517     coefs[i] = &(recvEqns[i]->coefs()[0]);
01518   }
01519 
01520   for(i=0; i<numRecvEqns; i++) {
01521 
01522     int eqn = recvEqnNumbers[i];
01523     if ((reducedStartRow_ > eqn) || (reducedEndRow_ < eqn)) {
01524       FEI_CERR << "LinSysCoreFilter::unpackRemoteContributions: ERROR, recvEqn "
01525            << eqn << " out of range. (localStartRow_: " << reducedStartRow_
01526            << ", localEndRow_: " << reducedEndRow_ << ", localRank_: "
01527            << localRank_ << ")" << FEI_ENDL;
01528       MPI_Abort(comm_, -1);
01529     }
01530 
01531     for(size_t ii=0; ii<recvEqns[i]->size(); ii++) {
01532       if (coefs[i][ii] > 1.e+200) {
01533         FEI_CERR << localRank_ << ": LinSysCoreFilter::unpackRemoteContributions: "
01534              << "WARNING, coefs["<<i<<"]["<<ii<<"]: " << coefs[i][ii]
01535              << FEI_ENDL;
01536         MPI_Abort(comm_, -1);
01537       }
01538     }
01539 
01540     if (recvEqns[i]->size() > 0 && newCoefs) {
01541       //contribute this equation to the matrix,
01542       CHK_ERR( giveToLocalReducedMatrix(1, &(recvEqnNumbers[i]),
01543                                         recvEqns[i]->size(),
01544                                         &(recvEqns[i]->indices()[0]),
01545                                         &(coefs[i]), assemblyMode ) );
01546     }
01547 
01548     //and now the RHS contributions.
01549     if (newRHSs) {
01550       for(int j=0; j<numRHSs_; j++) {
01551         CHK_ERR( giveToLocalReducedRHS(1, &( (*(recvRHSs[i]))[j] ),
01552                                        &eqn, assemblyMode) );
01553       }
01554     }
01555   }
01556 
01557   delete [] coefs;
01558 
01559   return(0);
01560 }
01561 
01562 //------------------------------------------------------------------------------
01563 int LinSysCoreFilter::exchangeRemoteBCs(std::vector<int>& essEqns,
01564                                         std::vector<double>& essAlpha,
01565                                         std::vector<double>& essGamma)
01566 {
01567   //we need to make sure that the right thing happens for essential
01568   //boundary conditions that get applied to nodes on elements that touch
01569   //a processor boundary. (Note that for this case, the BC node itself doesn't
01570   //touch the processor boundary.) For an essential boundary condition, the row
01571   //and column of the corresponding equation must be diagonalized. If there is
01572   //a processor boundary on any side of the element that contains the node,
01573   //then there are column contributions to the matrix on the other processor.
01574   //That other processor must be notified and told to make the adjustments
01575   //necessary to enforce the boundary condition.
01576 
01577   std::vector<int>* eqns = &essEqns;
01578 
01579   if (problemStructure_->numSlaveEquations() > 0) {
01580     int numEqns = essEqns.size();
01581     eqns = new std::vector<int>(numEqns);
01582     int* eqnsPtr = &(*eqns)[0];
01583 
01584     for(int ii=0; ii<numEqns; ++ii) {
01585       problemStructure_->translateToReducedEqn(essEqns[ii], eqnsPtr[ii]);
01586     }
01587   }
01588 
01589   FEI_OSTREAM* dbgOut = NULL;
01590   if (Filter::logStream() != NULL) {
01591     dbgOut = logStream();
01592   }
01593 
01594   eqnCommMgr_->exchangeRemEssBCs(&(*eqns)[0], eqns->size(),
01595                                  &essAlpha[0], &essGamma[0],
01596                                  comm_, dbgOut);
01597 
01598   int numRemoteEssBCEqns = eqnCommMgr_->getNumRemEssBCEqns();
01599   if (numRemoteEssBCEqns > 0) {
01600     std::vector<int>& remEssBCEqnNumbers = eqnCommMgr_->remEssBCEqnNumbersPtr();
01601     fei::CSVec** remEssBCEqns = &(eqnCommMgr_->remEssBCEqns()[0]);
01602     std::vector<int> remEssBCEqnLengths(remEssBCEqnNumbers.size());
01603 
01604     int** indices = new int*[numRemoteEssBCEqns];
01605     double** coefs = new double*[numRemoteEssBCEqns];
01606 
01607     for(int i=0; i<numRemoteEssBCEqns; i++) {
01608       coefs[i] = &(remEssBCEqns[i]->coefs()[0]);
01609       indices[i] = &(remEssBCEqns[i]->indices()[0]);
01610       remEssBCEqnLengths[i] = remEssBCEqns[i]->size();
01611     }
01612 
01613     CHK_ERR( enforceRemoteEssBCs(numRemoteEssBCEqns,
01614                                  &remEssBCEqnNumbers[0], indices,
01615                                  &remEssBCEqnLengths[0], coefs));
01616 
01617     delete [] indices;
01618     delete [] coefs;
01619   }
01620 
01621   if (problemStructure_->numSlaveEquations() > 0) {
01622     delete eqns;
01623   }
01624 
01625   debugOutput("#LinSysCoreFilter exchangeRemoteBCs");
01626 
01627   return(FEI_SUCCESS);
01628 }
01629 
01630 //------------------------------------------------------------------------------
01631 int LinSysCoreFilter::loadCRMult(int CRID,
01632                                  int numCRNodes,
01633                                  const GlobalID* CRNodes,
01634                                  const int* CRFields,
01635                                  const double* CRWeights,
01636                                  double CRValue)
01637 {
01638 //
01639 // Load Lagrange multiplier constraint relation data
01640 //
01641 //   Question: do we really need to pass CRNodes again?  Here, I'm going
01642 //            to ignore it for now (i.e., not store it, but just check it), 
01643 //            as it got passed during the initialization phase, so all we'll 
01644 //            do here is check for errors...
01645 //
01646   if (Filter::logStream() != NULL) {
01647     FEI_OSTREAM& os = *logStream();
01648     os<<"FEI: loadCRMult"<<FEI_ENDL;
01649     os<<"#num-nodes"<<FEI_ENDL<<numCRNodes<<FEI_ENDL;
01650     os<<"#CRNodes:"<<FEI_ENDL;
01651     int i;
01652     for(i=0; i<numCRNodes; ++i) {
01653       GlobalID nodeID = CRNodes[i];
01654       os << static_cast<int>(nodeID) << " ";
01655     }
01656     os << FEI_ENDL << "#fields:"<<FEI_ENDL;
01657     for(i=0; i<numCRNodes; ++i) os << CRFields[i] << " ";
01658     os << FEI_ENDL << "#field-sizes:"<<FEI_ENDL;
01659     for(i=0; i<numCRNodes; ++i) {
01660       int size = problemStructure_->getFieldSize(CRFields[i]);
01661       os << size << " ";
01662     }
01663     os << FEI_ENDL<<"#weights:"<<FEI_ENDL;
01664     int offset = 0;
01665     for(i=0; i<numCRNodes; ++i) {
01666       int size = problemStructure_->getFieldSize(CRFields[i]);
01667       for(int j=0; j<size; ++j) {
01668         os << CRWeights[offset++] << " ";
01669       }
01670     }
01671     os << FEI_ENDL<<"#CRValue:"<<FEI_ENDL<<CRValue<<FEI_ENDL;
01672   }
01673 
01674   ConstraintType* multCR = NULL;
01675   CHK_ERR( problemStructure_->getMultConstRecord(CRID, multCR) );
01676 
01677   int i;
01678 
01679   int lenList = multCR->getMasters().size();
01680   if (lenList < 1) {
01681     FEI_CERR << "ERROR in FEI, constraint with ID="<<CRID<<" appears to have"
01682          <<" a constrained-node list of length "<<lenList<<", should be > 0."<<FEI_ENDL;
01683     ERReturn(-1);
01684   }
01685 
01686   //  recall the data stored earlier and ensure that the passed data (here, the
01687   //  node list) agrees with the initialization data
01688 
01689   std::vector<GlobalID>& CRNode_vec = multCR->getMasters();
01690   GlobalID *CRNodePtr = &CRNode_vec[0];
01691 
01692   for(i=0; i<lenList; i++) {
01693     if (CRNodePtr[i] != CRNodes[i]) {
01694       FEI_CERR << "ERROR in FEI, constraint with ID="<<CRID<<" had different node-list"
01695            << " in initCRMult than it has in loadCRMult."<<FEI_ENDL;
01696       ERReturn(-1);
01697     }
01698   }
01699 
01700   std::vector<int>& CRField_vec = multCR->getMasterFieldIDs();
01701   int *CRFieldPtr = &CRField_vec[0];
01702 
01703   for (i = 0; i < lenList; i++) {
01704     if (CRFieldPtr[i] != CRFields[i]) {
01705       FEI_CERR <<"ERROR in FEI, constraint with CRID="<<CRID<<" had different field-list"
01706            <<" in initCRMult than it has in loadCRMult."<<FEI_ENDL;
01707       ERReturn(-1);
01708     }
01709   }
01710 
01711   newConstraintData_ = true;
01712 
01713   if ((int)iworkSpace_.size() < lenList) {
01714     iworkSpace_.resize(lenList);
01715   }
01716 
01717   int* fieldSizes = &iworkSpace_[0];
01718 
01719   for (i = 0; i < lenList; i++) {
01720     int numSolnParams = problemStructure_->getFieldSize(CRFields[i]);
01721     assert(numSolnParams >= 0);
01722     fieldSizes[i] = numSolnParams;
01723   }
01724 
01725   std::vector<double>& CRWeightArray = multCR->getMasterWeights();
01726 
01727   int offset = 0;
01728 
01729   try {
01730 
01731   for(i = 0; i < lenList; i++) {
01732     for(int j = 0; j < fieldSizes[i]; j++) {
01733       CRWeightArray.push_back(CRWeights[offset++]);
01734     }
01735   }
01736 
01737   }
01738   catch(std::runtime_error& exc) {
01739     FEI_CERR << exc.what() << FEI_ENDL;
01740     ERReturn(-1);
01741   }
01742 
01743   multCR->setRHSValue(CRValue);
01744   double* CRWeightsPtr = &CRWeightArray[0];
01745 
01746 //  next, perform assembly of the various terms into the system arrays
01747 //  (this is a good candidate for a separate function...)
01748 
01749   int irow = multCR->getEqnNumber();
01750   double zero = 0.0;
01751   double* zeroPtr = &zero;
01752   CHK_ERR( giveToMatrix(1, &irow, 1, &irow, &zeroPtr, ASSEMBLE_PUT) );
01753 
01754   CHK_ERR( giveToRHS(1, &(CRValue), &irow, ASSEMBLE_PUT));
01755 
01756   offset = 0;
01757   for(int j = 0; j < lenList; j++) {
01758     int myFieldID = CRFields[j];
01759 
01760     const NodeDescriptor& node = Filter::findNodeDescriptor(CRNodePtr[j]);
01761 
01762     //first, store the column coeficients for equation irow, the
01763     //constraint's equation.
01764     storeNodalColumnCoefs(irow, node, myFieldID, fieldSizes[j],
01765                           &(CRWeightsPtr[offset]));
01766 
01767 
01768     //next, store store the transpose of the above. i.e., column irow,
01769     //in equations associated with 'node' at 'myFieldID'.
01770 
01771     if (node.getOwnerProc() == localRank_) {
01772 
01773       storeNodalRowCoefs(node, myFieldID, fieldSizes[j],
01774                          &(CRWeightsPtr[offset]), irow);
01775     }
01776     else {
01777 
01778       storeNodalSendEqn(node, myFieldID, irow, &(CRWeightsPtr[offset]));
01779     }
01780 
01781     offset += fieldSizes[j];
01782   }
01783     
01784   return(FEI_SUCCESS);
01785 }
01786 
01787 //------------------------------------------------------------------------------
01788 int LinSysCoreFilter::loadCRPen(int CRID, 
01789                                 int numCRNodes,
01790                                 const GlobalID* CRNodes,
01791                                 const int* CRFields,
01792                                 const double* CRWeights,
01793                                 double CRValue,
01794                                 double penValue)
01795 {
01796   //
01797   // Load penalty constraint relation data
01798   //
01799 
01800   debugOutput("FEI: loadCRPen");
01801 
01802   ConstraintType* penCR = NULL;
01803   CHK_ERR( problemStructure_->getPenConstRecord(CRID, penCR) );
01804 
01805   int i;
01806   int lenList = penCR->getMasters().size();
01807   if (lenList < 1) {
01808     FEI_CERR << "ERROR in FEI, constraint with ID="<<CRID<<" appears to have"
01809          <<" a constrained-node list of length "<<lenList<<", should be > 0."<<FEI_ENDL;
01810     ERReturn(-1);
01811   }
01812 
01813   // recall the data stored earlier and ensure that the passed data (here,
01814   // the node list) agrees with the initialization data
01815 
01816   std::vector<GlobalID>& CRNode_vec = penCR->getMasters();
01817   GlobalID* CRNodePtr = &CRNode_vec[0];
01818                                   
01819   for(int j = 0; j < lenList; j++) {
01820     if (CRNodePtr[j] != CRNodes[j]) {
01821       FEI_CERR << "ERROR in FEI, constraint with ID="<<CRID<<" had different node-list"
01822            << " in initCRPen than it has in loadCRPen."<<FEI_ENDL;
01823       ERReturn(-1);
01824     }
01825   }
01826 
01827   newConstraintData_ = true;
01828 
01829   //  store the weights and rhs-value in the constraint records.
01830 
01831   if ((int)iworkSpace_.size() < lenList) {
01832     iworkSpace_.resize(lenList);
01833   }
01834 
01835   int* fieldSizes = &iworkSpace_[0];
01836 
01837   for (i = 0; i < lenList; i++) {
01838     int numSolnParams = problemStructure_->getFieldSize(CRFields[i]);
01839     assert(numSolnParams >= 0);
01840     fieldSizes[i] = numSolnParams;
01841   }
01842 
01843   std::vector<double>& CRWeightArray = penCR->getMasterWeights();
01844 
01845   try {
01846 
01847   int offset = 0;
01848   for (i = 0; i < lenList; i++) {
01849     for (int j = 0; j < fieldSizes[i]; j++) {
01850       CRWeightArray.push_back(CRWeights[offset++]);
01851     }
01852   }
01853 
01854   }
01855   catch(std::runtime_error& exc) {
01856     FEI_CERR << exc.what() << FEI_ENDL;
01857     ERReturn(-1);
01858   }
01859 
01860   penCR->setRHSValue(CRValue);
01861 
01862   double* CRWeightPtr = &CRWeightArray[0];
01863 
01864   int ioffset = 0, joffset = 0;
01865   for(i = 0; i < lenList; i++) {
01866     GlobalID iNodeID = CRNodePtr[i];
01867     int iField = CRFields[i];
01868 
01869     const NodeDescriptor& iNode = Filter::findNodeDescriptor(iNodeID);
01870     double* iweights = &(CRWeightPtr[ioffset]);
01871     ioffset += fieldSizes[i];
01872 
01873     joffset = 0;
01874     for (int j = 0; j < lenList; j++) {
01875       GlobalID jNodeID = CRNodePtr[j];
01876       int jField = CRFields[j];
01877 
01878       const NodeDescriptor& jNode = Filter::findNodeDescriptor(jNodeID);
01879       double* jweights = &(CRWeightPtr[joffset]);
01880       joffset += fieldSizes[j];
01881 
01882       double rhsValue = CRValue;
01883       if (j < lenList-1) {
01884         rhsValue = 0.0;
01885       }
01886 
01887       if (iNode.getOwnerProc() == localRank_) {
01888 
01889         storePenNodeData(iNode, iField, fieldSizes[i], iweights,
01890                          jNode, jField, fieldSizes[j], jweights,
01891                          penValue, rhsValue);
01892       }
01893       else {
01894         storePenNodeSendData(iNode, iField, fieldSizes[i], iweights,
01895                              jNode, jField, fieldSizes[j], jweights,
01896                              penValue, rhsValue);
01897       }
01898     }
01899   }
01900 
01901   return(FEI_SUCCESS);
01902 }
01903 
01904 //------------------------------------------------------------------------------
01905 int LinSysCoreFilter::parameters(int numParams, const char *const* paramStrings) {
01906 //
01907 // this function takes parameters for setting internal things like solver
01908 // and preconditioner choice, etc.
01909 //
01910    if (numParams == 0 || paramStrings == NULL) {
01911       debugOutput("#LinSysCoreFilter::parameters--- no parameters.");
01912    }
01913    else {
01914       const char* param1 = snl_fei::getParamValue("AZ_matrix_type",
01915                                                          numParams,
01916                                                          paramStrings);
01917 
01918       if (param1 != NULL) {
01919         if (!strcmp(param1, "AZ_VBR_MATRIX") ||
01920             !strcmp(param1, "blockMatrix")) {
01921           blockMatrix_ = true;
01922         }        
01923       }
01924       else {
01925         param1 = snl_fei::getParamValue("matrixType",
01926                                                numParams, paramStrings);
01927         if (param1 != NULL) {
01928           if (!strcmp(param1, "AZ_VBR_MATRIX") ||
01929               !strcmp(param1, "blockMatrix")) {
01930             blockMatrix_ = true;
01931           }        
01932         }
01933       }
01934 
01935       param1 = snl_fei::getParamValue("outputLevel",
01936                                              numParams,paramStrings);
01937       if ( param1 != NULL){
01938         std::string str(param1);
01939         FEI_ISTRINGSTREAM isstr(str);
01940         isstr >> outputLevel_;
01941       }
01942 
01943       param1 = snl_fei::getParam("resolveConflict",numParams,paramStrings);
01944       if ( param1 != NULL){
01945          resolveConflictRequested_ = true;
01946       }
01947 
01948       param1 = snl_fei::getParamValue("internalFei", numParams,paramStrings);
01949       if ( param1 != NULL ){
01950         std::string str(param1);
01951         FEI_ISTRINGSTREAM isstr(str);
01952         isstr >> internalFei_;
01953       }
01954 
01955       if (Filter::logStream() != NULL) {
01956 
01957         (*logStream())<<"#LinSysCoreFilter::parameters"<<FEI_ENDL
01958                          <<"# --- numParams: "<< numParams<<FEI_ENDL;
01959          for(int i=0; i<numParams; i++){
01960            (*logStream())<<"#------ paramStrings["<<i<<"]: "
01961                             <<paramStrings[i];
01962            if (paramStrings[i][strlen(paramStrings[i])-1] != '\n') {
01963              (*logStream())<<FEI_ENDL;
01964            }
01965          }
01966       }
01967    }
01968 
01969    CHK_ERR( Filter::parameters(numParams, paramStrings) );
01970 
01971    debugOutput("#LinSysCoreFilter leaving parameters function");
01972 
01973    return(FEI_SUCCESS);
01974 }
01975 
01976 //------------------------------------------------------------------------------
01977 int LinSysCoreFilter::loadComplete()
01978 {
01979   int len = 4;
01980   std::vector<int> flags(len), globalFlags(len);
01981   flags[0] = newMatrixData_ ? 1 : 0;
01982   flags[1] = newVectorData_ ? 1 : 0;
01983   flags[2] = newConstraintData_ ? 1 : 0;
01984   flags[3] = newBCData_ ? 1 : 0;
01985 
01986   CHK_ERR( fei::GlobalMax(comm_, flags, globalFlags) );
01987 
01988   newMatrixData_     = globalFlags[0] > 0 ? true : false;
01989   newVectorData_     = globalFlags[1] > 0 ? true : false;
01990   newConstraintData_ = globalFlags[2] > 0 ? true : false;
01991   newBCData_         = globalFlags[3] > 0 ? true : false;
01992 
01993   bool called_exchange = false;
01994   if (newMatrixData_ || newVectorData_ || newConstraintData_) {
01995     CHK_ERR( exchangeRemoteEquations() );
01996     called_exchange = true;
01997   }
01998 
01999   bool called_implbcs = false;
02000   if (newBCData_) {
02001     CHK_ERR( implementAllBCs() );
02002     called_implbcs = true;
02003   }
02004 
02005   if (called_exchange || called_implbcs ||needToCallMatrixLoadComplete_) {
02006     debugOutput("#LinSysCoreFilter calling LinSysCore matrixLoadComplete");
02007 
02008     CHK_ERR( lsc_->matrixLoadComplete() );
02009     needToCallMatrixLoadComplete_ = false;
02010   }
02011 
02012   newMatrixData_ = false;
02013   newVectorData_ = false;
02014   newConstraintData_ = false;
02015   newBCData_ = false;
02016 
02017   return(0);
02018 }
02019 
02020 //------------------------------------------------------------------------------
02021 int LinSysCoreFilter::residualNorm(int whichNorm, int numFields,
02022                            int* fieldIDs, double* norms, double& residTime)
02023 {
02024   //
02025   //This function can do 3 kinds of norms: infinity-norm (denoted
02026   //by whichNorm==0), 1-norm and 2-norm.
02027   //
02028   debugOutput("FEI: residualNorm");
02029 
02030   if (whichNorm < 0 || whichNorm > 2) return(-1);
02031 
02032   CHK_ERR( loadComplete() );
02033 
02034   std::vector<double> residValues(numReducedRows_, 0.0);
02035 
02036   double start = fei::utils::cpu_time();
02037 
02038   CHK_ERR( formResidual(&(residValues[0]), numReducedRows_) );
02039 
02040   residTime = fei::utils::cpu_time() - start;
02041 
02042   CHK_ERR( Filter::calculateResidualNorms(whichNorm, numFields,
02043                                           fieldIDs, norms, residValues) );
02044 
02045   return(FEI_SUCCESS);
02046 }
02047 
02048 //------------------------------------------------------------------------------
02049 int LinSysCoreFilter::formResidual(double* residValues, int numLocalEqns)
02050 {
02051   CHK_ERR( lsc_->formResidual(residValues, numLocalEqns))
02052 
02053   return(FEI_SUCCESS);
02054 }
02055 
02056 //------------------------------------------------------------------------------
02057 int LinSysCoreFilter::solve(int& status, double& sTime) {
02058 
02059    debugOutput("FEI: solve");
02060 
02061    CHK_ERR( loadComplete() );
02062 
02063    debugOutput("#LinSysCoreFilter in solve, calling launchSolver...");
02064  
02065    double start = fei::utils::cpu_time();
02066 
02067    CHK_ERR( lsc_->launchSolver(status, iterations_) );
02068 
02069    sTime = fei::utils::cpu_time() - start;
02070 
02071    debugOutput("#LinSysCoreFilter... back from solver");
02072  
02073    //now unpack the locally-owned shared entries of the solution vector into
02074    //the eqn-comm-mgr data structures.
02075    CHK_ERR( unpackSolution() );
02076 
02077    debugOutput("#LinSysCoreFilter leaving solve");
02078 
02079    if (status != 0) return(1);
02080    else return(FEI_SUCCESS);
02081 }
02082 
02083 //------------------------------------------------------------------------------
02084 int LinSysCoreFilter::setNumRHSVectors(int numRHSs, int* rhsIDs){
02085 
02086    if (numRHSs < 0) {
02087       FEI_CERR << "LinSysCoreFilter::setNumRHSVectors: ERROR, numRHSs < 0." << FEI_ENDL;
02088       ERReturn(-1);
02089    }
02090 
02091    numRHSs_ = numRHSs;
02092 
02093    rhsIDs_.resize(numRHSs_);
02094    for(int i=0; i<numRHSs_; i++) rhsIDs_[i] = rhsIDs[i];
02095 
02096   //(we need to set the number of RHSs in the eqn comm manager)
02097   eqnCommMgr_->setNumRHSs(numRHSs_);
02098 
02099    return(FEI_SUCCESS);
02100 }
02101 
02102 //------------------------------------------------------------------------------
02103 int LinSysCoreFilter::setCurrentRHS(int rhsID)
02104 {
02105    std::vector<int>::iterator iter =
02106       std::find( rhsIDs_.begin(), rhsIDs_.end(), rhsID);
02107 
02108    if (iter == rhsIDs_.end()) ERReturn(-1)
02109  
02110    int index = iter - rhsIDs_.begin();
02111    currentRHS_ = index;
02112 
02113    lsc_->setRHSID(rhsID);
02114 
02115    return(FEI_SUCCESS);
02116 }
02117 
02118 //------------------------------------------------------------------------------
02119 int LinSysCoreFilter::giveToMatrix_symm_noSlaves(int numPtRows,
02120                                                  const int* ptRowNumbers,
02121                                                  const double* const* coefs,
02122                                                  int mode)
02123 {
02124   for(int i=0; i<numPtRows; i++) {
02125     int row = ptRowNumbers[i];
02126     const double* valptr = coefs[i];
02127     if (row < localStartRow_ || row > localEndRow_) {
02128       eqnCommMgr_->addRemoteEqn(row, valptr, ptRowNumbers, numPtRows);
02129       continue;
02130     }
02131 
02132     if (mode == ASSEMBLE_SUM) {
02133       if (Filter::logStream() != NULL && 0) {
02134         FEI_OSTREAM& os = *logStream();
02135         os << "#  calling sumIntoSystemMatrix, row: " << ptRowNumbers[i]
02136            << ", columns: ";
02137         for(int j=0; j<numPtRows; ++j) os << ptRowNumbers[j] << " ";
02138         os << FEI_ENDL;
02139       }
02140 
02141       CHK_ERR( lsc_->sumIntoSystemMatrix(1, &(ptRowNumbers[i]),
02142                                          numPtRows, ptRowNumbers,
02143                                          &valptr) );
02144     }
02145     else {
02146       CHK_ERR( lsc_->putIntoSystemMatrix(1, &(ptRowNumbers[i]),
02147                                          numPtRows, ptRowNumbers,
02148                                          &valptr) );
02149     }
02150   }
02151 
02152   return(0);
02153 }
02154 
02155 //------------------------------------------------------------------------------
02156 int LinSysCoreFilter::giveToBlkMatrix_symm_noSlaves(int numPtRows,
02157                                                     const int* ptRowNumbers,
02158                                                     int numBlkRows,
02159                                                     const int* blkRowNumbers,
02160                                                     const int* blkRowSizes,
02161                                                     const double* const* coefs,
02162                                                     int mode)
02163 {
02164   int i;
02165   if ((int)dworkSpace2_.size() < numPtRows) {
02166     dworkSpace2_.resize(numPtRows);
02167   }
02168   const double* * valptr = &dworkSpace2_[0];
02169   for(i=0; i<numPtRows; i++) {
02170     int row = ptRowNumbers[i];
02171     valptr[i] = coefs[i];
02172     if (row < localStartRow_ || row > localEndRow_) {
02173       eqnCommMgr_->addRemoteEqn(row, valptr[i], ptRowNumbers, numPtRows);
02174       continue;
02175     }
02176 
02177     if (mode == ASSEMBLE_PUT) {
02178        CHK_ERR( lsc_->putIntoSystemMatrix(1, &(ptRowNumbers[i]),
02179                                          numPtRows, ptRowNumbers,
02180                                          &(valptr[i])) );
02181    }
02182   }
02183 
02184   int offset = 0;
02185   for(i=0; i<numBlkRows; i++) {
02186     int row = ptRowNumbers[offset];
02187     if (row < localStartRow_ || row > localEndRow_) {
02188       offset += blkRowSizes[i];
02189       continue;
02190     }
02191 
02192     if (mode == ASSEMBLE_SUM) {
02193       if (Filter::logStream() != NULL && 0) {
02194         FEI_OSTREAM& os = *logStream();
02195         os << "#  calling sumIntoSystemMatrix, row: " << ptRowNumbers[i]
02196            << ", columns: ";
02197         for(int j=0; j<numPtRows; ++j) os << ptRowNumbers[j] << " ";
02198         os << FEI_ENDL;
02199       }
02200       
02201       CHK_ERR(lsc_->sumIntoSystemMatrix(blkRowSizes[i],&(ptRowNumbers[offset]),
02202                                         numPtRows, ptRowNumbers,
02203                                         1, &(blkRowNumbers[i]),
02204                                         numBlkRows, blkRowNumbers,
02205                                         &(valptr[offset])) );
02206     }
02207     
02208     offset += blkRowSizes[i];
02209   }
02210 
02211   return(0);
02212 }
02213 
02214 //------------------------------------------------------------------------------
02215 int LinSysCoreFilter::giveToMatrix(int numPtRows, const int* ptRows,
02216                            int numPtCols, const int* ptCols,
02217                            const double* const* values, int mode)
02218 {
02219   try {
02220 
02221   if (problemStructure_->numSlaveEquations() == 0) {
02222     for(int i=0; i<numPtRows; i++) {
02223       if (ptRows[i] < localStartRow_ || ptRows[i] > localEndRow_) {
02224         eqnCommMgr_->addRemoteEqn(ptRows[i], values[i], ptCols, numPtCols);
02225         continue;
02226       }
02227 
02228       if (mode == ASSEMBLE_SUM) {
02229         if (Filter::logStream() != NULL && 0) {
02230           FEI_OSTREAM& os = *logStream();
02231           os << "#  calling sumIntoSystemMatrix, row: " << ptRows[i]
02232              << ", columns: ";
02233           for(int j=0; j<numPtCols; ++j) os << ptCols[j] << " ";
02234           os << FEI_ENDL;
02235         }
02236 
02237         CHK_ERR( lsc_->sumIntoSystemMatrix(1, &(ptRows[i]),
02238                                            numPtCols, ptCols,
02239                                            &(values[i])) );
02240       }
02241       else {
02242         CHK_ERR( lsc_->putIntoSystemMatrix(1, &(ptRows[i]),
02243                                            numPtCols, ptCols,
02244                                            &(values[i])) );
02245       }
02246     }
02247   }
02248   else {
02249     iworkSpace_.resize(numPtCols);
02250     iworkSpace2_.resize(numPtCols);
02251     int* iworkPtr = &iworkSpace_[0];
02252     int* iworkPtr2= &iworkSpace2_[0];
02253     int offset = 0;
02254     for(int ii=0; ii<numPtCols; ii++) {
02255       int reducedEqn = -1;
02256       bool isSlave = problemStructure_->translateToReducedEqn(ptCols[ii],
02257                                                               reducedEqn);
02258       if (isSlave) {
02259         reducedEqn = -1;
02260         iworkPtr[ii] = reducedEqn;
02261       }
02262       else {
02263         iworkPtr[ii] = reducedEqn;
02264         iworkPtr2[offset++] = reducedEqn;
02265       }
02266     }
02267     iworkSpace2_.resize(offset);
02268 
02269     for(int i=0; i<numPtRows; i++) {
02270       int row = ptRows[i];
02271 
02272       int reducedRow;
02273       bool isSlave = problemStructure_->translateToReducedEqn(row, reducedRow);
02274       if (isSlave) continue;
02275 
02276       if (reducedStartRow_ > reducedRow || reducedRow > reducedEndRow_) {
02277 
02278         dworkSpace_.resize(0);
02279         for(int j=0; j<numPtCols; j++) {
02280           if (iworkSpace_[j]>=0) {
02281             if (Filter::logStream() != NULL) {
02282               (*logStream())<<"#  giveToMatrix remote("<<reducedRow<<","
02283                             <<iworkSpace_[j]<<","<<values[i][j]<<")"<<FEI_ENDL;
02284             }
02285 
02286             dworkSpace_.push_back(values[i][j]);
02287           }
02288         }
02289 
02290         if (mode == ASSEMBLE_SUM) {
02291           if (Filter::logStream() != NULL) {
02292             (*logStream())<<"sum"<<FEI_ENDL;
02293           }
02294 
02295           eqnCommMgr_->addRemoteEqn(reducedRow,
02296                                     &dworkSpace_[0],
02297                                     &iworkSpace2_[0],
02298                                     iworkSpace2_.size());
02299         }
02300         else {
02301           if (Filter::logStream() != NULL) {
02302             (*logStream())<<"put"<<FEI_ENDL;
02303           }
02304 
02305           eqnCommMgr_put_->addRemoteEqn(reducedRow,
02306                                         &dworkSpace_[0],
02307                                         &iworkSpace2_[0],
02308                                         iworkSpace2_.size());
02309         }
02310 
02311         continue;
02312       }
02313 
02314       for(int j=0; j<numPtCols; j++) {
02315 
02316         int reducedCol = iworkPtr[j];
02317         if (reducedCol<0) continue;
02318 
02319         double* tmpCoef = const_cast<double*>(&(values[i][j]));
02320 
02321         if (Filter::logStream() != NULL) {
02322           (*logStream())<< "#  giveToMatrix local("<<reducedRow
02323                         <<","<<reducedCol<<","<<*tmpCoef<<")"<<FEI_ENDL;
02324         }
02325 
02326         if (mode == ASSEMBLE_SUM) {
02327           if (Filter::logStream() != NULL && 0) {
02328             FEI_OSTREAM& os = *logStream();
02329             os << "#  calling sumIntoSystemMatrix, row: " << reducedRow
02330                << ", columns: " << reducedCol << FEI_ENDL;
02331           }
02332 
02333           CHK_ERR( lsc_->sumIntoSystemMatrix(1, &reducedRow, 1, &reducedCol,
02334                                              &tmpCoef ) );
02335         }
02336         else {
02337           CHK_ERR( lsc_->putIntoSystemMatrix(1, &reducedRow, 1, &reducedCol,
02338                                              &tmpCoef ) );
02339         }
02340       }
02341     }
02342   }
02343 
02344   }
02345   catch(std::runtime_error& exc) {
02346     FEI_CERR << exc.what() << FEI_ENDL;
02347     ERReturn(-1);
02348   }
02349 
02350   return(FEI_SUCCESS);
02351 }
02352 
02353 //------------------------------------------------------------------------------
02354 int LinSysCoreFilter::giveToLocalReducedMatrix(int numPtRows, const int* ptRows,
02355                                        int numPtCols, const int* ptCols,
02356                                        const double* const* values, int mode)
02357 {
02358   bool specialCase = (!firstRemEqnExchange_ && newConstraintData_
02359                       && !newMatrixData_) ? true : false;
02360 
02361   double fei_min = std::numeric_limits<double>::min();
02362 
02363   for(int i=0; i<numPtRows; i++) {
02364 
02365     if (mode == ASSEMBLE_SUM) {
02366       const double* values_i = values[i];
02367 
02368       for(int j=0; j<numPtCols; ++j) {
02369         if (specialCase && std::abs(values_i[j]) < fei_min) continue;
02370 
02371         const double* valPtr = &(values_i[j]);
02372         CHK_ERR( lsc_->sumIntoSystemMatrix(1, &(ptRows[i]), 1, &(ptCols[j]),
02373                                            &valPtr) );
02374       }
02375     }
02376     else {
02377       CHK_ERR( lsc_->putIntoSystemMatrix(1, &(ptRows[i]), numPtCols, ptCols,
02378                                          &(values[i])) );
02379     }
02380   }
02381 
02382   return(FEI_SUCCESS);
02383 }
02384 
02385 //------------------------------------------------------------------------------
02386 int LinSysCoreFilter::sumIntoMatrix(fei::CSRMat& mat)
02387 {
02388   const std::vector<int>& rowNumbers = mat.getGraph().rowNumbers;
02389   const std::vector<int>& rowOffsets = mat.getGraph().rowOffsets;
02390   const std::vector<int>& pckColInds = mat.getGraph().packedColumnIndices;
02391   const std::vector<double>& pckCoefs = mat.getPackedCoefs();
02392 
02393   for(size_t i=0; i<rowNumbers.size(); ++i) {
02394     int row = rowNumbers[i];
02395     int offset = rowOffsets[i];
02396     int rowlen = rowOffsets[i+1]-offset;
02397     const int* indices = &pckColInds[offset];
02398     const double* coefs = &pckCoefs[offset];
02399 
02400     if (giveToMatrix(1, &row, rowlen, indices, &coefs, ASSEMBLE_SUM) != 0) {
02401       throw std::runtime_error("fei::impl_utils::add_to_matrix ERROR in matrix.sumIn.");
02402     }
02403   }
02404 
02405   return(FEI_SUCCESS);
02406 }
02407 
02408 //------------------------------------------------------------------------------
02409 int LinSysCoreFilter::getFromMatrix(int numPtRows, const int* ptRows,
02410                             const int* rowColOffsets, const int* ptCols,
02411                             int numColsPerRow, double** values)
02412 {
02413   //This function may be attempting to retrieve matrix rows that are not
02414   //locally owned. If those rows correspond to finite-element nodes that we
02415   //share, AND if the owning processor is also making this function call, then
02416   //we can communicate with that processor and obtain those matrix rows.
02417   //
02418 
02419   ProcEqns remoteProcEqns;
02420 
02421   //Let's populate this ProcEqns object with the remote equations and the procs
02422   //that we need to receive the remote equations from.
02423   for(int re=0; re<numPtRows; re++) {
02424     int eqn = ptRows[re];
02425     int owner = problemStructure_->getOwnerProcForEqn(eqn);
02426     if (owner == localRank_) continue;
02427 
02428     remoteProcEqns.addEqn(eqn, owner);
02429   }
02430 
02431   //so now we know which of the requested equations are remotely owned, and we
02432   //know which processors own them.
02433   //Next we're going to need to know which locally-owned equations are needed
02434   //by other processors.
02435   ProcEqns localProcEqns;
02436   CHK_ERR( eqnCommMgr_->mirrorProcEqns(remoteProcEqns, localProcEqns) )
02437 
02438   //ok, now we know which local equations we'll need to send, so let's extract
02439   //those from the matrix
02440   EqnBuffer localEqns;
02441   CHK_ERR( getEqnsFromMatrix(localProcEqns, localEqns) )
02442 
02443   //now we can set the lengths in localProcEqns.
02444   std::vector<int>& eqnNumbers = localEqns.eqnNumbers();
02445   fei::CSVec** localEqnsPtr = (localEqns.eqns().size() ? &(localEqns.eqns()[0]) : 0);
02446   std::vector<int> eqnLengths(eqnNumbers.size());
02447   for(size_t i=0; i<eqnNumbers.size(); ++i) {
02448     eqnLengths[i] = localEqnsPtr[i]->size();
02449   }
02450 
02451   localProcEqns.setProcEqnLengths(&eqnNumbers[0], &eqnLengths[0],
02452                                   eqnNumbers.size());
02453 
02454   //now mirror those lengths in the remoteProcEqns objects to get ready for the
02455   //all-to-all exchange of equation data.
02456   CHK_ERR( eqnCommMgr_->mirrorProcEqnLengths(localProcEqns, remoteProcEqns) );
02457 
02458   EqnBuffer remoteEqns;
02459   //we're now ready to do the exchange.
02460   CHK_ERR( EqnCommMgr::exchangeEqnBuffers(comm_, &localProcEqns, &localEqns,
02461                                           &remoteProcEqns, &remoteEqns, false));
02462 
02463   std::vector<int>& remEqnNumbers = remoteEqns.eqnNumbers();
02464   fei::CSVec** remEqnsPtr = (remoteEqns.eqns().size() ? &(remoteEqns.eqns()[0]) : 0);
02465   std::vector<fei::CSVec*>& remEqns   = remoteEqns.eqns();
02466 
02467   //now we're ready to fill the values array with the remote coefficients.
02468   for(int i=0; i<numPtRows; i++) {
02469     int row = ptRows[i];
02470 
02471     int eqnIndex = fei::binarySearch(row, remEqnNumbers);
02472 
02473     //if eqnIndex < 0, this is a local equation, so skip to the next loop iter.
02474     if (eqnIndex < 0) continue;
02475 
02476     //the equation is remote, so stride across it copying out the coefs.
02477     //if ptCols is NULL, then we're going to copy all coefficients (the whole
02478     //row) into 'values'.
02479     if (ptCols == NULL) {
02480       for(size_t j=0; j<remEqnsPtr[eqnIndex]->size(); j++) {
02481         values[i][j] = remEqns[eqnIndex]->coefs()[j];
02482       }
02483       continue;
02484     }
02485 
02486     for(int j=0; j<numColsPerRow; j++) {
02487       int offset = rowColOffsets[i] + j;
02488       int colIndex = fei::binarySearch(ptCols[offset], remEqns[eqnIndex]->indices());
02489       if (colIndex < 0) ERReturn(-1);
02490 
02491       values[i][j] = remEqns[eqnIndex]->coefs()[colIndex];
02492     }
02493   }
02494 
02495   //and now, get the local stuff out of the matrix.
02496   for(int i=0; i<numPtRows; i++) {
02497     int row = ptRows[i];
02498     if (row < localStartRow_ || localEndRow_ < row) continue;
02499 
02500     int rowLen = 0, checkRowLen;
02501     CHK_ERR( lsc_->getMatrixRowLength(row, rowLen) )
02502       if (rowLen <= 0) ERReturn(-1);
02503 
02504     //for each local row, establish some temp arrays and get the row from
02505     //the matrix.
02506 
02507     std::vector<double> coefs(rowLen);
02508     std::vector<int> indices(rowLen);
02509 
02510     CHK_ERR( lsc_->getMatrixRow(row, &coefs[0], &indices[0], rowLen, checkRowLen) );
02511     if (rowLen != checkRowLen) ERReturn(-1);
02512 
02513     //now stride across the list of requested column-indices, and find the
02514     //corresponding location in the matrix row. Copy that location into the
02515     //values array.
02516 
02517     //again, if ptCols is NULL, then we're going to copy all coefficients 
02518     //(the whole row) into 'values'.
02519     if (ptCols == NULL) {
02520       for(int j=0; j<rowLen; j++) {
02521         values[i][j] = coefs[j];
02522       }
02523       continue;
02524     }
02525 
02526     for(int j=0; j<numColsPerRow; j++) {
02527       std::vector<int>::iterator iter =
02528           std::find(indices.begin(), indices.end(), ptCols[rowColOffsets[i]+j]);
02529       if (iter == indices.end()) {
02530         ERReturn(-1);
02531       }
02532 
02533       int index = iter - indices.begin();
02534       values[i][j] = coefs[index];
02535     }
02536   }
02537 
02538   return(FEI_SUCCESS);
02539 }
02540 
02541 //------------------------------------------------------------------------------
02542 int LinSysCoreFilter::getEqnsFromMatrix(ProcEqns& procEqns, EqnBuffer& eqnData)
02543 {
02544   //Given a ProcEqns object containing lists of equation-numbers, get the data
02545   //for those equations from the local portion of the matrix and store that data
02546   //in the eqnData object.
02547 
02548   std::vector<std::vector<int>*>& eqnNumbers = procEqns.procEqnNumbersPtr();
02549 
02550   for(unsigned p=0; p<eqnNumbers.size(); p++) {
02551     for(unsigned i=0; i<eqnNumbers[p]->size(); i++) {
02552       int eqn = (*(eqnNumbers[p]))[i];
02553 
02554       if (localStartRow_ > eqn || localEndRow_ < eqn) continue;
02555 
02556       //if this equation is already in eqnData, then don't put it in again...
02557       std::vector<int>& eqnDataEqns = eqnData.eqnNumbers();
02558       if (fei::binarySearch(eqn, eqnDataEqns) >= 0) continue;
02559 
02560       int len = 0;
02561       CHK_ERR( lsc_->getMatrixRowLength(eqn, len) );
02562 
02563       if (len <= 0) continue;
02564       std::vector<double> coefs(len);
02565       std::vector<int> indices(len);
02566       int outputLen = 0;
02567 
02568       CHK_ERR( lsc_->getMatrixRow(eqn, &coefs[0], &indices[0],
02569                                   len, outputLen) );
02570       if (outputLen != len) ERReturn(-1);
02571 
02572       CHK_ERR( eqnData.addEqn(eqn, &coefs[0], &indices[0], len, false) );
02573     }
02574   }
02575   return(FEI_SUCCESS);
02576 }
02577 
02578 //------------------------------------------------------------------------------
02579 int LinSysCoreFilter::getEqnsFromRHS(ProcEqns& procEqns, EqnBuffer& eqnData)
02580 {
02581   //Given a ProcEqns object containing lists of equation-numbers, get the data
02582   //for those equations from the local portion of the RHS vector and store that
02583   // data in the eqnData object. We're only storing rhs coefs in an EqnBuffer
02584   //that was designed for also storing equations with column-indices. So we'll
02585   //put a bogus column-index in with each equation just to make sure the
02586   //EqnBuffer does the right stuff internally...
02587 
02588   int numSendProcs = procEqns.getNumProcs();
02589   std::vector<int>& eqnsPerProc = procEqns.eqnsPerProcPtr();
02590   std::vector<std::vector<int>*>& eqnNumbers = procEqns.procEqnNumbersPtr();
02591 
02592   eqnData.setNumRHSs(1);
02593 
02594   for(int p=0; p<numSendProcs; p++) {
02595     for(int i=0; i<eqnsPerProc[p]; i++) {
02596       int reducedEqn;
02597       problemStructure_->translateToReducedEqn((*(eqnNumbers[p]))[i], reducedEqn);
02598 
02599       if (reducedStartRow_ > reducedEqn || reducedEndRow_ < reducedEqn) continue;
02600 
02601       //if this equation is already in eqnData, then don't put it in again...
02602       std::vector<int>& eqnDataEqns = eqnData.eqnNumbers();
02603       if (fei::binarySearch(reducedEqn, eqnDataEqns) >= 0) continue;
02604 
02605       double rhsValue;
02606 
02607       CHK_ERR( lsc_->getFromRHSVector(1, &rhsValue, &reducedEqn) );
02608 
02609       int bogusIndex = 19;
02610       CHK_ERR( eqnData.addIndices(reducedEqn, &bogusIndex, 1) );
02611       CHK_ERR( eqnData.addRHS(reducedEqn, 0, rhsValue) );
02612     }
02613   }
02614   return(FEI_SUCCESS);
02615 }
02616 
02617 //------------------------------------------------------------------------------
02618 int LinSysCoreFilter::giveToRHS(int num, const double* values,
02619                         const int* indices, int mode)
02620 {
02621   if (problemStructure_->numSlaveEquations() == 0) {
02622     for(int i=0; i<num; i++) {
02623       if (indices[i] < localStartRow_ || indices[i] > localEndRow_) {
02624         if (mode == ASSEMBLE_SUM) {
02625           eqnCommMgr_->addRemoteRHS(indices[i], currentRHS_, values[i]);
02626         }
02627         else {
02628           eqnCommMgr_put_->addRemoteRHS(indices[i], currentRHS_, values[i]);
02629         }
02630 
02631         continue;
02632       }
02633 
02634       if (mode == ASSEMBLE_SUM) {
02635         CHK_ERR( lsc_->sumIntoRHSVector(1, &(values[i]), &(indices[i])) );
02636       }
02637       else {
02638         CHK_ERR( lsc_->putIntoRHSVector(1, &(values[i]), &(indices[i])) );
02639       }
02640     }
02641   }
02642   else {
02643     for(int i=0; i<num; i++) {
02644       int reducedEqn;
02645       bool isSlave = problemStructure_->
02646         translateToReducedEqn(indices[i], reducedEqn);
02647       if (isSlave) continue;
02648 
02649       if (reducedEqn < reducedStartRow_ || reducedEqn > reducedEndRow_) {
02650         if (mode == ASSEMBLE_SUM) {
02651           eqnCommMgr_->addRemoteRHS(reducedEqn, currentRHS_, values[i]);
02652         }
02653         else {
02654           eqnCommMgr_put_->addRemoteRHS(reducedEqn, currentRHS_, values[i]);
02655         }
02656 
02657         continue;
02658       }
02659 
02660       if (mode == ASSEMBLE_SUM) {
02661         CHK_ERR( lsc_->sumIntoRHSVector(1, &(values[i]), &reducedEqn) );
02662       }
02663       else {
02664         CHK_ERR( lsc_->putIntoRHSVector(1, &(values[i]), &reducedEqn) );
02665       }
02666     }
02667   }
02668 
02669   return(FEI_SUCCESS);
02670 }
02671 
02672 //------------------------------------------------------------------------------
02673 int LinSysCoreFilter::giveToLocalReducedRHS(int num, const double* values,
02674                                     const int* indices, int mode)
02675 {
02676   for(int i=0; i<num; i++) {
02677 
02678     if (mode == ASSEMBLE_SUM) {
02679       CHK_ERR( lsc_->sumIntoRHSVector(1, &(values[i]), &(indices[i])) );
02680     }
02681     else {
02682       CHK_ERR( lsc_->putIntoRHSVector(1, &(values[i]), &(indices[i])) );
02683     }
02684   }
02685 
02686   return(FEI_SUCCESS);
02687 }
02688 
02689 //------------------------------------------------------------------------------
02690 int LinSysCoreFilter::sumIntoRHS(fei::CSVec& vec)
02691 {
02692   std::vector<int>& indices = vec.indices();
02693   std::vector<double>& coefs = vec.coefs();
02694 
02695   CHK_ERR( giveToRHS(indices.size(), &coefs[0], &indices[0], ASSEMBLE_SUM) );
02696 
02697   return(FEI_SUCCESS);
02698 }
02699 
02700 //------------------------------------------------------------------------------
02701 int LinSysCoreFilter::getFromRHS(int num, double* values, const int* indices)
02702 {
02703   //We need to do similar things here as we do in getFromMatrix, with respect to
02704   //communications to obtain values for equations that are remotely owned.
02705 
02706   ProcEqns remoteProcEqns;
02707 
02708   //Let's populate this ProcEqns object with the remote equations and the procs
02709   //that we need to receive the remote equations from.
02710   for(int re=0; re<num; re++) {
02711     int eqn = indices[re];
02712     int owner = problemStructure_->getOwnerProcForEqn(eqn);
02713     if (owner==localRank_) continue;
02714 
02715     remoteProcEqns.addEqn(eqn, owner);
02716   }
02717 
02718   //so now we know which of the requested equations are remotely owned, and we
02719   //know which processors own them.
02720   //Next we're going to need to know which locally-owned equations are needed
02721   //by other processors.
02722   ProcEqns localProcEqns;
02723   CHK_ERR( eqnCommMgr_->mirrorProcEqns(remoteProcEqns, localProcEqns) );
02724 
02725   //ok, now we know which equations we'll need to send, so let's extract
02726   //them from the rhs vector.
02727   EqnBuffer localEqns;
02728   CHK_ERR( getEqnsFromRHS(localProcEqns, localEqns) );
02729 
02730   //now we can set the lengths in localProcEqns.
02731   std::vector<int>& eqnNumbers = localEqns.eqnNumbers();
02732   fei::CSVec** localEqnsPtr = &(localEqns.eqns()[0]);
02733   std::vector<int> eqnLengths(eqnNumbers.size());
02734   for(size_t i=0; i<eqnNumbers.size(); ++i) {
02735     eqnLengths[i] = localEqnsPtr[i]->size();
02736   }
02737 
02738   localProcEqns.setProcEqnLengths(&eqnNumbers[0], &eqnLengths[0],
02739                                   eqnNumbers.size());
02740 
02741   //now mirror those lengths in the remoteProcEqns objects to get ready for the
02742   //all-to-all exchange of equation data.
02743   CHK_ERR( eqnCommMgr_->mirrorProcEqnLengths(localProcEqns, remoteProcEqns) );
02744 
02745   EqnBuffer remoteEqns;
02746   //we're now ready to do the exchange.
02747   CHK_ERR( EqnCommMgr::exchangeEqnBuffers(comm_, &localProcEqns, &localEqns,
02748                                            &remoteProcEqns, &remoteEqns, false))
02749 
02750   //now we're ready to get the rhs data we've received from other processors.
02751   std::vector<int>& remEqnNumbers = remoteEqns.eqnNumbers();
02752   std::vector<std::vector<double>*>& remRhsCoefs = *(remoteEqns.rhsCoefsPtr());
02753 
02754   for(int i=0; i<num; i++) {
02755     int row = indices[i];
02756 
02757     int eqnIndex = fei::binarySearch(row, remEqnNumbers);
02758     if (eqnIndex < 0) continue;
02759 
02760     values[i] = (*(remRhsCoefs[eqnIndex]))[0];
02761   }
02762 
02763   //and now get the local stuff.
02764   for(int i=0; i<num; i++) {
02765     if (indices[i] < localStartRow_ || indices[i] > localEndRow_) continue;
02766 
02767     CHK_ERR( lsc_->getFromRHSVector(num, values, indices) );
02768   }
02769 
02770   return(FEI_SUCCESS);
02771 }
02772 
02773 //------------------------------------------------------------------------------
02774 int LinSysCoreFilter::getEqnSolnEntry(int eqnNumber, double& solnValue)
02775 {
02776   //This function's task is to retrieve the solution-value for a global
02777   //equation-number. eqnNumber may or may not be a slave-equation, and may or
02778   //may not be locally owned. If it is not locally owned, it should at least
02779   //be shared.
02780   //return 0 if the solution is successfully retrieved, otherwise return 1.
02781   //
02782 
02783   //
02784   //First and probably most common case: there are no slave equations.
02785   //
02786   if (problemStructure_->numSlaveEquations() == 0) {
02787     if (localStartRow_ > eqnNumber || eqnNumber > localEndRow_) {
02788       //Dig into the eqn-comm-mgr for the shared-remote solution value.
02789       CHK_ERR( getSharedRemoteSolnEntry(eqnNumber, solnValue) );
02790     }
02791     else {
02792       //It's local, simply get the solution from the assembled linear system.
02793       CHK_ERR( getReducedSolnEntry( eqnNumber, solnValue ) );
02794     }
02795     return(0);
02796   }
02797 
02798   //
02799   //If we reach this point, there are slave equations to account for.
02800   //So first translate this equation into 'assembled-linear-system'
02801   //equation-numbers.
02802   //
02803   int reducedEqn;
02804   bool isSlave = problemStructure_->translateToReducedEqn(eqnNumber,reducedEqn);
02805 
02806   if (isSlave) {
02807     //This is a slave-equation, so construct its solution-value as the linear-
02808     //combination of the master-equations it is defined in terms of.
02809 
02810     std::vector<int>* masterEqns = NULL;
02811     std::vector<double>* masterCoefs = NULL;
02812     CHK_ERR( problemStructure_->getMasterEqnNumbers(eqnNumber, masterEqns) );
02813     CHK_ERR( problemStructure_->getMasterEqnCoefs(eqnNumber, masterCoefs) );
02814 
02815     int len = masterEqns->size();
02816     solnValue = 0.0;
02817     CHK_ERR( problemStructure_->getMasterEqnRHS(eqnNumber, solnValue) );
02818 
02819     double coef = 0.0;
02820     for(int i=0; i<len; i++) {
02821       int mEqn = (*masterEqns)[i];
02822       int mReducedeqn;
02823       problemStructure_->translateToReducedEqn(mEqn, mReducedeqn);
02824 
02825       if (reducedStartRow_ > mReducedeqn || mReducedeqn > reducedEndRow_) {
02826         CHK_ERR( getSharedRemoteSolnEntry(mReducedeqn, coef) );
02827       }
02828       else {
02829         CHK_ERR( getReducedSolnEntry(mReducedeqn, coef) );
02830       }
02831       solnValue += coef * (*masterCoefs)[i];
02832     }
02833   }
02834   else {
02835     //This is not a slave-equation, so retrieve the solution from either the
02836     //assembled linear system or the shared-remote data structures.
02837 
02838     if (reducedStartRow_ > reducedEqn || reducedEqn > reducedEndRow_) {
02839       CHK_ERR( getSharedRemoteSolnEntry(reducedEqn, solnValue) );
02840     }
02841     else {
02842       CHK_ERR( getReducedSolnEntry(reducedEqn, solnValue) );
02843     }
02844   }
02845 
02846   return(0);
02847 }
02848 
02849 //------------------------------------------------------------------------------
02850 int LinSysCoreFilter::getSharedRemoteSolnEntry(int eqnNumber, double& solnValue)
02851 {
02852   std::vector<int>& remoteEqnNumbers = eqnCommMgr_->sendEqnNumbersPtr();
02853   double* remoteSoln = eqnCommMgr_->sendEqnSolnPtr();
02854 
02855   int index = fei::binarySearch(eqnNumber, remoteEqnNumbers);
02856   if (index < 0) {
02857     FEI_CERR << "LinSysCoreFilter::getSharedRemoteSolnEntry: ERROR, eqn "
02858          << eqnNumber << " not found." << FEI_ENDL;
02859     ERReturn(-1);
02860   }
02861   solnValue = remoteSoln[index];
02862   return(0);
02863 }
02864 
02865 //------------------------------------------------------------------------------
02866 int LinSysCoreFilter::getReducedSolnEntry(int eqnNumber, double& solnValue)
02867 {
02868   //We may safely assume that this function is called with 'eqnNumber' that is
02869   //local in the underlying assembled linear system. i.e., it isn't a slave-
02870   //equation, it isn't remotely owned, etc.
02871   //
02872   CHK_ERR( lsc_->getSolnEntry(eqnNumber, solnValue) );
02873 
02874   return(FEI_SUCCESS);
02875 }
02876 
02877 //------------------------------------------------------------------------------
02878 int LinSysCoreFilter::unpackSolution()
02879 {
02880   //
02881   //This function should be called after the solver has returned,
02882   //and we know that there is a solution in the underlying vector.
02883   //This function ensures that any locally-owned shared solution values are
02884   //available on the sharing processors.
02885   //
02886   if (Filter::logStream() != NULL) {
02887     (*logStream())<< "#  entering unpackSolution, outputLevel: "
02888                      <<outputLevel_<<FEI_ENDL;
02889   }
02890 
02891   //what we need to do is as follows.
02892   //The eqn comm mgr has a list of what it calls 'recv eqns'. These are
02893   //equations that we own, for which we received contributions from other
02894   //processors. The solution values corresponding to these equations need
02895   //to be made available to those remote contributing processors.
02896 
02897   int numRecvEqns = eqnCommMgr_->getNumLocalEqns();
02898   std::vector<int>& recvEqnNumbers = eqnCommMgr_->localEqnNumbers();
02899 
02900   for(int i=0; i<numRecvEqns; i++) {
02901     int eqn = recvEqnNumbers[i];
02902 
02903     if ((reducedStartRow_ > eqn) || (reducedEndRow_ < eqn)) {
02904       FEI_CERR << "LinSysCoreFilter::unpackSolution: ERROR, 'recv' eqn (" << eqn
02905            << ") out of local range." << FEI_ENDL;
02906       MPI_Abort(comm_, -1);
02907     }
02908 
02909     double solnValue = 0.0;
02910 
02911     CHK_ERR( getReducedSolnEntry(eqn, solnValue) );
02912 
02913     eqnCommMgr_->addSolnValues(&eqn, &solnValue, 1);
02914   }
02915 
02916   eqnCommMgr_->exchangeSoln();
02917 
02918   debugOutput("#LinSysCoreFilter leaving unpackSolution");
02919   return(FEI_SUCCESS);
02920 }
02921              
02922 //------------------------------------------------------------------------------
02923 void LinSysCoreFilter::setEqnCommMgr(EqnCommMgr* eqnCommMgr)
02924 {
02925   delete eqnCommMgr_;
02926   eqnCommMgr_ = eqnCommMgr;
02927 }
02928 
02929 //------------------------------------------------------------------------------
02930 int LinSysCoreFilter::getBlockNodeSolution(GlobalID elemBlockID,  
02931                                    int numNodes, 
02932                                    const GlobalID *nodeIDs, 
02933                                    int *offsets,
02934                                    double *results) {
02935         
02936    debugOutput("FEI: getBlockNodeSolution");
02937 
02938    int numActiveNodes = problemStructure_->getNumActiveNodes();
02939    NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
02940 
02941    if (numActiveNodes <= 0) return(0);
02942 
02943    int numSolnParams = 0;
02944 
02945    BlockDescriptor* block = NULL;
02946    CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
02947 
02948    //Traverse the node list, checking if nodes are associated with this block.
02949    //If so, put its 'answers' in the results list.
02950 
02951    int offset = 0;
02952    for(int i=0; i<numActiveNodes; i++) {
02953      const NodeDescriptor* node_i = NULL;
02954      CHK_ERR( nodeDB.getNodeAtIndex(i, node_i) );
02955 
02956       if (offset == numNodes) break;
02957 
02958       GlobalID nodeID = nodeIDs[offset];
02959 
02960       //first let's set the offset at which this node's solution coefs start.
02961       offsets[offset++] = numSolnParams;
02962 
02963       const NodeDescriptor* node = NULL;
02964       int err = 0;
02965       //Obtain the NodeDescriptor of nodeID in the activeNodes list...
02966       //Don't call the getActiveNodeDesc_ID function unless we have to.
02967 
02968       if (nodeID == node_i->getGlobalNodeID()) {
02969         node = node_i;
02970       }
02971       else {
02972          err = nodeDB.getNodeWithID(nodeID, node);
02973       }
02974 
02975       //ok. If err is not 0, meaning nodeID is NOT in the
02976       //activeNodes list, then skip to the next loop iteration.
02977 
02978       if (err != 0) {
02979         continue;
02980       }
02981 
02982       int numFields = node->getNumFields();
02983       const int* fieldIDs = node->getFieldIDList();
02984 
02985       for(int j=0; j<numFields; j++) {
02986         if (block->containsField(fieldIDs[j])) {
02987           int size = problemStructure_->getFieldSize(fieldIDs[j]);
02988           assert(size >= 0);
02989 
02990           int thisEqn = -1;
02991           node->getFieldEqnNumber(fieldIDs[j], thisEqn);
02992 
02993           double answer;
02994           for(int k=0; k<size; k++) {
02995             CHK_ERR( getEqnSolnEntry(thisEqn+k, answer) )
02996               results[numSolnParams++] = answer;
02997           }
02998         }
02999       }//for(j<numFields)loop
03000    }
03001 
03002    offsets[numNodes] = numSolnParams;
03003 
03004    return(FEI_SUCCESS);
03005 }
03006 
03007 //------------------------------------------------------------------------------
03008 int LinSysCoreFilter::getNodalSolution(int numNodes, 
03009                                        const GlobalID *nodeIDs, 
03010                                        int *offsets,
03011                                        double *results)
03012 {
03013   debugOutput("FEI: getNodalSolution");
03014 
03015   int numActiveNodes = problemStructure_->getNumActiveNodes();
03016   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
03017 
03018   if (numActiveNodes <= 0) return(0);
03019 
03020   int numSolnParams = 0;
03021 
03022   //Traverse the node list, checking if nodes are local.
03023   //If so, put 'answers' in the results list.
03024 
03025   int offset = 0;
03026   for(int i=0; i<numActiveNodes; i++) {
03027     const NodeDescriptor* node_i = NULL;
03028     CHK_ERR( nodeDB.getNodeAtIndex(i, node_i) );
03029 
03030     if (offset == numNodes) break;
03031 
03032     GlobalID nodeID = nodeIDs[offset];
03033 
03034     //first let's set the offset at which this node's solution coefs start.
03035     offsets[offset++] = numSolnParams;
03036 
03037     const NodeDescriptor* node = NULL;
03038     int err = 0;
03039     //Obtain the NodeDescriptor of nodeID in the activeNodes list...
03040     //Don't call the getNodeWithID function unless we have to.
03041 
03042     if (nodeID == node_i->getGlobalNodeID()) {
03043       node = node_i;
03044     }
03045     else {
03046       err = nodeDB.getNodeWithID(nodeID, node);
03047     }
03048 
03049     //ok. If err is not 0, meaning nodeID is NOT in the
03050     //activeNodes list, then skip to the next loop iteration.
03051 
03052     if (err != 0) {
03053       continue;
03054     }
03055 
03056     int numFields = node->getNumFields();
03057     const int* fieldIDs = node->getFieldIDList();
03058 
03059     for(int j=0; j<numFields; j++) {
03060       int size = problemStructure_->getFieldSize(fieldIDs[j]);
03061       assert(size >= 0);
03062 
03063       int thisEqn = -1;
03064       node->getFieldEqnNumber(fieldIDs[j], thisEqn);
03065 
03066       double answer;
03067       for(int k=0; k<size; k++) {
03068         CHK_ERR( getEqnSolnEntry(thisEqn+k, answer) )
03069           results[numSolnParams++] = answer;
03070       }
03071     }//for(j<numFields)loop
03072   }
03073 
03074   offsets[numNodes] = numSolnParams;
03075 
03076   return(FEI_SUCCESS);
03077 }
03078 
03079 //------------------------------------------------------------------------------
03080 int LinSysCoreFilter::getBlockFieldNodeSolution(GlobalID elemBlockID,
03081                                         int fieldID,
03082                                         int numNodes, 
03083                                         const GlobalID *nodeIDs, 
03084                                         double *results)
03085 {
03086   //Note: if the user-supplied nodeIDs list containts nodes which are not in
03087   //the specified element-block, then the corresponding positions in the
03088   //results array are simply not referenced. This is dangerous behavior that
03089   //hasn't gotten me into trouble yet.
03090 
03091   debugOutput("FEI: getBlockFieldNodeSolution");
03092 
03093   int numActiveNodes = problemStructure_->getNumActiveNodes();
03094   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
03095 
03096   if (numActiveNodes <= 0) return(0);
03097 
03098   BlockDescriptor* block = NULL;
03099   CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
03100 
03101   int fieldSize = problemStructure_->getFieldSize(fieldID);
03102   if (fieldSize <= 0) ERReturn(-1);
03103 
03104   if (!block->containsField(fieldID)) {
03105     FEI_CERR << "LinSysCoreFilter::getBlockFieldNodeSolution WARNING: fieldID " << fieldID
03106          << " not contained in element-block " << (int)elemBlockID << FEI_ENDL;
03107     return(1);
03108   }
03109 
03110    //Traverse the node list, checking if nodes are associated with this block.
03111    //If so, put the answers in the results list.
03112 
03113    for(int i=0; i<numNodes; i++) {
03114      const NodeDescriptor* node_i = NULL;
03115      CHK_ERR( nodeDB.getNodeAtIndex(i, node_i) );
03116 
03117      GlobalID nodeID = nodeIDs[i];
03118 
03119      const NodeDescriptor* node = NULL;
03120      int err = 0;
03121      //Obtain the NodeDescriptor of nodeID in the activeNodes list...
03122      //Don't call the getNodeWithID function unless we have to. (getNodeWithID
03123      //does a binary-search over all local nodeIDs, while getNodeAtIndex is
03124      //a direct lookup.) Often the user supplies a nodeIDs list that is in the
03125      //"natural" order, so we don't need to call getNodeWithID at all.
03126 
03127      if (nodeID == node_i->getGlobalNodeID()) {
03128        node = node_i;
03129      }
03130      else {
03131        err = nodeDB.getNodeWithID(nodeID, node);
03132      }
03133 
03134      //ok. If err is not 0, meaning nodeID is NOT in the
03135      //activeNodes list, then skip to the next loop iteration.
03136 
03137      if (err != 0) {
03138        continue;
03139      }
03140 
03141      int eqnNumber = -1;
03142      bool hasField = node->getFieldEqnNumber(fieldID, eqnNumber);
03143      if (!hasField) continue;
03144 
03145      int offset = fieldSize*i;
03146      for(int j=0; j<fieldSize; j++) {
03147        double answer = 0.0;
03148        CHK_ERR( getEqnSolnEntry(eqnNumber+j, answer) );
03149        results[offset+j] = answer;
03150      }
03151    }
03152 
03153    return(FEI_SUCCESS);
03154 }
03155 
03156 //------------------------------------------------------------------------------
03157 int LinSysCoreFilter::getNodalFieldSolution(int fieldID,
03158                                     int numNodes, 
03159                                     const GlobalID *nodeIDs, 
03160                                     double *results)
03161 {
03162   debugOutput("FEI: getNodalFieldSolution");
03163 
03164   int numActiveNodes = problemStructure_->getNumActiveNodes();
03165   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
03166 
03167   if (numActiveNodes <= 0) return(0);
03168 
03169   int fieldSize = problemStructure_->getFieldSize(fieldID);
03170   if (fieldSize <= 0) ERReturn(-1);
03171 
03172   //Traverse the node list, checking if nodes have the specified field.
03173   //If so, put the answers in the results list.
03174 
03175   for(int i=0; i<numNodes; i++) {
03176     const NodeDescriptor* node_i = NULL;
03177     CHK_ERR( nodeDB.getNodeAtIndex(i, node_i) );
03178 
03179     GlobalID nodeID = nodeIDs[i];
03180 
03181     const NodeDescriptor* node = NULL;
03182     int err = 0;
03183     //Obtain the NodeDescriptor of nodeID in the activeNodes list...
03184     //Don't call the getNodeWithID function unless we have to.
03185 
03186     if (nodeID == node_i->getGlobalNodeID()) {
03187       node = node_i;
03188     }
03189     else {
03190       err = nodeDB.getNodeWithID(nodeID, node);
03191     }
03192 
03193     //ok. If err is not 0, meaning nodeID is NOT in the
03194     //activeNodes list, then skip to the next loop iteration.
03195 
03196     if (err != 0) {
03197       continue;
03198     }
03199 
03200     int eqnNumber = -1;
03201     bool hasField = node->getFieldEqnNumber(fieldID, eqnNumber);
03202 
03203     //If this node doesn't have the specified field, then skip to the
03204     //next loop iteration.
03205     if (!hasField) continue;
03206 
03207     int offset = fieldSize*i;
03208     for(int j=0; j<fieldSize; j++) {
03209       double answer = 0.0;
03210       CHK_ERR( getEqnSolnEntry(eqnNumber+j, answer) );
03211       results[offset+j] = answer;
03212     }
03213   }
03214 
03215   return(FEI_SUCCESS);
03216 }
03217 
03218 //------------------------------------------------------------------------------
03219 int LinSysCoreFilter::putBlockNodeSolution(GlobalID elemBlockID,
03220                                    int numNodes, 
03221                                    const GlobalID *nodeIDs, 
03222                                    const int *offsets,
03223                                    const double *estimates) {
03224         
03225    debugOutput("FEI: putBlockNodeSolution");
03226 
03227    int numActiveNodes = problemStructure_->getNumActiveNodes();
03228 
03229    if (numActiveNodes <= 0) return(0);
03230 
03231    BlockDescriptor* block = NULL;
03232    CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
03233 
03234    NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
03235 
03236    //traverse the node list, checking for nodes associated with this block
03237    //when an associated node is found, put its 'answers' into the linear system.
03238 
03239    for(int i=0; i<numNodes; i++) {
03240      const NodeDescriptor* node = NULL;
03241      int err = nodeDB.getNodeWithID(nodeIDs[i], node);
03242 
03243       if (err != 0) continue;
03244    
03245       if (!node->containedInBlock(elemBlockID)) continue;
03246 
03247       if (node->getOwnerProc() != localRank_) continue;
03248 
03249       int numFields = node->getNumFields();
03250       const int* fieldIDs = node->getFieldIDList();
03251       const int* fieldEqnNumbers = node->getFieldEqnNumbers();
03252 
03253       if (fieldEqnNumbers[0] < localStartRow_ ||
03254           fieldEqnNumbers[0] > localEndRow_) continue;
03255 
03256       int offs = offsets[i];
03257 
03258       for(int j=0; j<numFields; j++) {
03259          int size = problemStructure_->getFieldSize(fieldIDs[j]);
03260 
03261          if (block->containsField(fieldIDs[j])) {
03262             for(int k=0; k<size; k++) {
03263                int reducedEqn;
03264                problemStructure_->
03265                  translateToReducedEqn(fieldEqnNumbers[j]+k, reducedEqn);
03266 
03267                CHK_ERR( lsc_->putInitialGuess(&reducedEqn,
03268                                               &estimates[offs+k], 1) );
03269             }
03270          }
03271          offs += size;
03272       }//for(j<numFields)loop
03273    }
03274 
03275    return(FEI_SUCCESS);
03276 }
03277 
03278 //------------------------------------------------------------------------------
03279 int LinSysCoreFilter::putBlockFieldNodeSolution(GlobalID elemBlockID, 
03280                                         int fieldID, 
03281                                         int numNodes, 
03282                                         const GlobalID *nodeIDs, 
03283                                         const double *estimates)
03284 {
03285    int fieldSize = problemStructure_->getFieldSize(fieldID);
03286 
03287    if (Filter::logStream() != NULL) {
03288      FEI_OSTREAM& os = *logStream();
03289      os << "FEI: putBlockFieldNodeSolution" << FEI_ENDL;
03290      os << "#blkID" << FEI_ENDL << (int)elemBlockID << FEI_ENDL
03291         << "#fieldID"<<FEI_ENDL << fieldID << FEI_ENDL
03292         << "#fieldSize"<<FEI_ENDL << fieldSize << FEI_ENDL
03293         << "#numNodes"<<FEI_ENDL << numNodes << FEI_ENDL
03294         << "#nodeIDs" << FEI_ENDL;
03295      int i;
03296      for(i=0; i<numNodes; ++i) os << (int)nodeIDs[i] << FEI_ENDL;
03297      os << "#estimates" << FEI_ENDL;
03298      for(i=0; i<numNodes*fieldSize; ++i) os << estimates[i] << FEI_ENDL;
03299    }
03300 
03301    BlockDescriptor* block = NULL;
03302    CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
03303    if (!block->containsField(fieldID)) return(1);
03304 
03305    NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
03306 
03307    //if we have a negative fieldID, we'll need a list of length numNodes,
03308    //in which to put nodeNumbers for passing to the solver... 
03309 
03310    std::vector<int> numbers(numNodes);
03311 
03312    //if we have a fieldID >= 0, then our numbers list will hold equation numbers
03313    //and we'll need fieldSize*numNodes of them.
03314 
03315    std::vector<double> data;
03316 
03317    if (fieldID >= 0) {
03318       assert(fieldSize >= 0);
03319       numbers.resize(numNodes*fieldSize);
03320       data.resize(numNodes*fieldSize);
03321    }
03322 
03323    int count = 0;
03324 
03325    for(int i=0; i<numNodes; i++) {
03326      const NodeDescriptor* node = NULL;
03327      CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
03328 
03329       if (fieldID < 0) numbers[count++] = node->getNodeNumber();
03330       else {
03331          int eqn = -1;
03332          if (node->getFieldEqnNumber(fieldID, eqn)) {
03333            if (eqn >= localStartRow_ && eqn <= localEndRow_) {
03334              for(int j=0; j<fieldSize; j++) { 
03335                data[count] = estimates[i*fieldSize + j];
03336                problemStructure_->translateToReducedEqn(eqn+j, numbers[count++]);
03337              }
03338            }
03339          }
03340       }
03341    }
03342 
03343    if (fieldID < 0) {
03344      CHK_ERR( lsc_->putNodalFieldData(fieldID, fieldSize, 
03345                                       &numbers[0], numNodes, estimates));
03346    }
03347    else {
03348      CHK_ERR(lsc_->putInitialGuess(&numbers[0], &data[0], count));
03349    }
03350 
03351    return(FEI_SUCCESS);
03352 }
03353 
03354 //------------------------------------------------------------------------------
03355 int LinSysCoreFilter::getBlockElemSolution(GlobalID elemBlockID,
03356                                    int numElems, 
03357                                    const GlobalID *elemIDs,
03358                                    int& numElemDOFPerElement,
03359                                    double *results)
03360 {
03361 //
03362 //  return the elemental solution parameters associated with a
03363 //  particular block of elements
03364 //
03365    debugOutput("FEI: getBlockElemSolution");
03366 
03367    BlockDescriptor* block = NULL;
03368    CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) )
03369 
03370    numElemDOFPerElement = block->getNumElemDOFPerElement();
03371    if (numElemDOFPerElement <= 0) return(0);
03372 
03373    ConnectivityTable& ctable =
03374      problemStructure_->getBlockConnectivity(elemBlockID);
03375    std::map<GlobalID,int>& elemIDList = ctable.elemIDs;
03376 
03377    std::vector<int>& elemDOFEqnNumbers = block->elemDOFEqnNumbers();
03378    double answer;
03379 
03380    for(int i=0; i<numElems; i++) {
03381       std::map<GlobalID,int>::const_iterator
03382         iter = elemIDList.find(elemIDs[i]);
03383       if (iter == elemIDList.end()) continue;
03384       int index = iter->second;
03385 
03386       int offset = i*numElemDOFPerElement;
03387 
03388       for(int j=0; j<numElemDOFPerElement; j++) {
03389          int eqn = elemDOFEqnNumbers[index] + j;
03390 
03391          CHK_ERR( getEqnSolnEntry(eqn, answer) )
03392 
03393          results[offset+j] = answer;
03394       }
03395    }
03396 
03397    return(FEI_SUCCESS);
03398 }
03399 
03400 //------------------------------------------------------------------------------
03401 int LinSysCoreFilter::putBlockElemSolution(GlobalID elemBlockID,
03402                                    int numElems,
03403                                    const GlobalID *elemIDs,
03404                                    int dofPerElem,
03405                                    const double *estimates)
03406 {
03407    debugOutput("FEI: putBlockElemSolution");
03408 
03409    BlockDescriptor* block = NULL;
03410    CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) )
03411 
03412    int DOFPerElement = block->getNumElemDOFPerElement();
03413    assert(DOFPerElement == dofPerElem);
03414    if (DOFPerElement <= 0) return(0);
03415 
03416    ConnectivityTable& ctable =
03417      problemStructure_->getBlockConnectivity(elemBlockID);
03418    std::map<GlobalID,int>& elemIDList = ctable.elemIDs;
03419 
03420    std::vector<int>& elemDOFEqnNumbers = block->elemDOFEqnNumbers();
03421 
03422 
03423    for(int i=0; i<numElems; i++) {
03424      std::map<GlobalID,int>::const_iterator
03425        iter = elemIDList.find(elemIDs[i]);
03426      if (iter == elemIDList.end()) continue;
03427 
03428      int index = iter->second;
03429 
03430       for(int j=0; j<DOFPerElement; j++) {
03431          int reducedEqn;
03432          problemStructure_->
03433            translateToReducedEqn(elemDOFEqnNumbers[index] + j, reducedEqn);
03434          double soln = estimates[i*DOFPerElement + j];
03435 
03436          CHK_ERR( lsc_->putInitialGuess(&reducedEqn, &soln, 1) );
03437       }
03438    }
03439 
03440    return(FEI_SUCCESS);
03441 }
03442 
03443 //------------------------------------------------------------------------------
03444 int LinSysCoreFilter::getCRMultipliers(int numCRs,
03445                                        const int* CRIDs,
03446                                        double* multipliers)
03447 {
03448   int multCRsLen = problemStructure_->getNumMultConstRecords();
03449   if (numCRs > multCRsLen) {
03450     return(-1);
03451   }
03452 
03453   std::map<GlobalID, ConstraintType*>::const_iterator
03454     cr_iter = problemStructure_->getMultConstRecords().begin(),
03455     cr_end  = problemStructure_->getMultConstRecords().end();
03456 
03457   int i = 0;
03458   while(cr_iter != cr_end && i < numCRs) {
03459     GlobalID CRID = (*cr_iter).first;
03460     ConstraintType* multCR = (*cr_iter).second;
03461     if (CRID != CRIDs[i]) {
03462       CHK_ERR( problemStructure_->getMultConstRecord(CRIDs[i], multCR) );
03463     }
03464 
03465     int eqn = multCR->getEqnNumber();
03466 
03467     CHK_ERR( getEqnSolnEntry(eqn, multipliers[i]) );
03468     ++cr_iter;
03469     ++i;
03470   }
03471 
03472   return(FEI_SUCCESS);
03473 }
03474 
03475 //------------------------------------------------------------------------------
03476 int LinSysCoreFilter::putCRMultipliers(int numMultCRs,
03477                                const int* CRIDs,
03478                                const double *multEstimates)
03479 {
03480   debugOutput("FEI: putCRMultipliers");
03481 
03482   for(int j = 0; j < numMultCRs; j++) {
03483     ConstraintType* multCR = NULL;
03484     CHK_ERR( problemStructure_->getMultConstRecord(CRIDs[j], multCR) );
03485 
03486     int eqnNumber = multCR->getEqnNumber();
03487     if (eqnNumber < localStartRow_ || eqnNumber > localEndRow_) continue;
03488     CHK_ERR( lsc_->putInitialGuess(&eqnNumber, &(multEstimates[j]), 1));
03489   }
03490 
03491   return(FEI_SUCCESS);
03492 }
03493 
03494 //------------------------------------------------------------------------------
03495 int LinSysCoreFilter::putNodalFieldData(int fieldID,
03496                                         int numNodes,
03497                                         const GlobalID* nodeIDs,
03498                                         const double* nodeData)
03499 {
03500   debugOutput("FEI: putNodalFieldData");
03501 
03502   if (fieldID > -1) {
03503     return(putNodalFieldSolution(fieldID, numNodes, nodeIDs, nodeData));
03504   }
03505 
03506   int fieldSize = problemStructure_->getFieldSize(fieldID);
03507   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
03508 
03509   std::vector<int> nodeNumbers(numNodes);
03510 
03511   for(int i=0; i<numNodes; i++) {
03512     const NodeDescriptor* node = NULL;
03513     CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
03514 
03515     int nodeNumber = node->getNodeNumber();
03516     if (nodeNumber < 0) {
03517       FEI_CERR << "LinSysCoreFilter::putNodalFieldData ERROR, node with ID " 
03518            << (int)nodeIDs[i] << " doesn't have an associated nodeNumber "
03519            << "assigned. putNodalFieldData shouldn't be called until after the "
03520            << "initComplete method has been called." << FEI_ENDL;
03521       ERReturn(-1);
03522     }
03523 
03524     nodeNumbers[i] = nodeNumber;
03525   }
03526 
03527   CHK_ERR( lsc_->putNodalFieldData(fieldID, fieldSize,
03528                                    &nodeNumbers[0], numNodes, nodeData) );
03529 
03530   return(0);
03531 }
03532 
03533 //------------------------------------------------------------------------------
03534 int LinSysCoreFilter::putNodalFieldSolution(int fieldID,
03535                                 int numNodes,
03536                                 const GlobalID* nodeIDs,
03537                                 const double* nodeData)
03538 {
03539   debugOutput("FEI: putNodalFieldSolution");
03540 
03541   if (fieldID < 0) {
03542     return(putNodalFieldData(fieldID, numNodes, nodeIDs, nodeData));
03543   }
03544 
03545   int fieldSize = problemStructure_->getFieldSize(fieldID);
03546   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
03547 
03548   std::vector<int> eqnNumbers(fieldSize);
03549 
03550   for(int i=0; i<numNodes; i++) {
03551     const NodeDescriptor* node = NULL;
03552     CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
03553 
03554     int eqn = -1;
03555     bool hasField = node->getFieldEqnNumber(fieldID, eqn);
03556     if (!hasField) continue;
03557 
03558     int reducedEqn = -1;
03559     bool isSlave = problemStructure_->translateToReducedEqn(eqn, reducedEqn);
03560     if (isSlave) continue;
03561 
03562     if (reducedStartRow_ > reducedEqn || reducedEndRow_ < reducedEqn) continue;
03563 
03564     int localLen = fieldSize;
03565     for(int j=0; j<fieldSize; j++) {
03566       int thisEqn = reducedEqn+j;
03567       if (reducedStartRow_ > thisEqn || reducedEndRow_ <thisEqn) {
03568         localLen = j;
03569       }
03570 
03571       eqnNumbers[j] = reducedEqn+j;
03572     }
03573 
03574     int offset = i*fieldSize;
03575     CHK_ERR( lsc_->putInitialGuess(&eqnNumbers[0],
03576                                    &nodeData[offset], localLen) );
03577   }
03578 
03579   return(0);
03580 }
03581 
03582 //------------------------------------------------------------------------------
03583 int LinSysCoreFilter::assembleEqns(int numPtRows, 
03584                                    int numPtCols,
03585                                    const int* rowNumbers,
03586                                    const int* colIndices,
03587                                    const double* const* coefs,
03588                                    bool structurallySymmetric,
03589                                    int numBlkEqns, int* blkEqns,
03590                                    int* blkSizes, bool useBlkEqns,
03591                                    int mode)
03592 {
03593   if (numPtRows == 0) return(FEI_SUCCESS);
03594 
03595   bool anySlaves = false;
03596   int numSlaveEqns = problemStructure_->numSlaveEquations();
03597   if (numSlaveEqns > 0) {
03598     rSlave_.resize(numPtRows);
03599     cSlave_.resize(0);
03600     const int* indPtr = colIndices;
03601     for(int r=0; r<numPtRows; r++) {
03602       rSlave_[r] = problemStructure_->isSlaveEqn(rowNumbers[r]) ? 1 : 0;
03603       if (rSlave_[r] == 1) anySlaves = true;
03604 
03605       for(int j=0; j<numPtCols; j++) {
03606         int isSlave = problemStructure_->isSlaveEqn(indPtr[j]) ? 1 : 0;
03607         cSlave_.push_back(isSlave);
03608         if (isSlave == 1) anySlaves = true;
03609       }
03610 
03611       if (!structurallySymmetric) indPtr += numPtCols;
03612     }
03613   }
03614 
03615   if (numSlaveEqns == 0 || !anySlaves) {
03616     if (numSlaveEqns == 0 && structurallySymmetric) {
03617       if (useBlkEqns) {
03618         CHK_ERR( giveToBlkMatrix_symm_noSlaves(numPtRows, rowNumbers,
03619                                                numBlkEqns, blkEqns, blkSizes,
03620                                                coefs, mode) );
03621       }
03622       else {
03623         CHK_ERR( giveToMatrix_symm_noSlaves(numPtRows, rowNumbers, coefs, mode) );
03624       }
03625     }
03626     else {
03627       if ((int)dworkSpace2_.size() < numPtRows) {
03628         dworkSpace2_.resize(numPtRows);
03629       }
03630       const double* * coefPtr = &dworkSpace2_[0];
03631       for(int i=0; i<numPtRows; i++) {
03632         coefPtr[i] = coefs[i];
03633       }
03634 
03635       if (structurallySymmetric) {
03636         CHK_ERR( giveToMatrix(numPtRows, rowNumbers, numPtRows, rowNumbers,
03637                               coefPtr, mode) );
03638       }
03639       else {
03640         const int* indPtr = colIndices;
03641         for(int i=0; i<numPtRows; i++) {
03642           int row = rowNumbers[i];
03643 
03644           const double* coefPtr1 = coefs[i];
03645 
03646           CHK_ERR(giveToMatrix(1, &row, numPtCols, indPtr, &coefPtr1, mode));
03647           indPtr += numPtCols;
03648         }
03649       }
03650     }
03651   }
03652   else {
03653     int offset = 0;
03654     const int* indicesPtr = colIndices;
03655     for(int i=0; i<numPtRows; i++) {
03656       int row = rowNumbers[i];
03657 
03658       const double* coefPtr = coefs[i];
03659       int* colSlave = &cSlave_[offset];
03660       offset += numPtCols;
03661 
03662       if (rSlave_[i] == 1) {
03663         //Since this is a slave equation, the non-slave columns of this row go
03664         //into 'Kdi_', and the slave columns go into 'Kdd_'.
03665         for(int jj=0; jj<numPtCols; jj++) {
03666           int col = indicesPtr[jj];
03667           if (colSlave[jj]) {
03668             Kdd_->sumInCoef(row, col, coefPtr[jj]);
03669           }
03670           else {
03671             Kdi_->sumInCoef(row, col, coefPtr[jj]);
03672           }
03673         }
03674 
03675         //We also need to put the non-slave rows of column 'row' into 'K_id'.
03676         const int* ii_indicesPtr = colIndices;
03677         for(int ii=0; ii<numPtRows; ii++) {
03678           int rowi = rowNumbers[ii];
03679           if (rSlave_[ii] == 1) continue;
03680 
03681           int index = fei::binarySearch(row, ii_indicesPtr, numPtCols);
03682           if (index < 0) continue;
03683 
03684           const double* coefs_ii = coefs[ii];
03685 
03686           Kid_->sumInCoef(rowi, row, coefs_ii[index]);
03687 
03688           if (!structurallySymmetric) ii_indicesPtr += numPtCols;
03689         }
03690 
03691         reducedEqnCounter_++;
03692 
03693         continue;
03694       }
03695       else {//row is not a slave eqn...
03696 
03697         //put all non-slave columns from this row into the assembled matrix.
03698 
03699         CHK_ERR( giveToMatrix(1, &row, numPtCols, indicesPtr, &coefPtr, mode) );
03700       }
03701 
03702       if (!structurallySymmetric) indicesPtr += numPtCols;
03703     }
03704 
03705     if (reducedEqnCounter_ > 300) CHK_ERR( assembleReducedEqns() );
03706   }
03707 
03708   return(FEI_SUCCESS);
03709 }
03710 
03711 //------------------------------------------------------------------------------
03712 int LinSysCoreFilter::assembleReducedEqns()
03713 {
03714   fei::FillableMat* D = problemStructure_->getSlaveDependencies();
03715 
03716   csrD = *D;
03717   csrKid = *Kid_;
03718   csrKdi = *Kdi_;
03719   csrKdd = *Kdd_;
03720 
03721   //form tmpMat1_ = Kid_*D
03722   fei::multiply_CSRMat_CSRMat(csrKid, csrD, tmpMat1_);
03723 
03724   //form tmpMat2_ = D^T*Kdi_
03725   fei::multiply_trans_CSRMat_CSRMat(csrD, csrKdi, tmpMat2_);
03726 
03727   if (Filter::logStream() != NULL) {
03728     FEI_OSTREAM& os = *Filter::logStream();
03729     os << "#  tmpMat1_"<<FEI_ENDL << tmpMat1_ << FEI_ENDL;
03730     os << "#  tmpMat2_"<<FEI_ENDL << tmpMat2_ << FEI_ENDL;
03731   }
03732 
03733   //accumulate the above two results into the global system matrix.
03734   CHK_ERR( sumIntoMatrix(tmpMat1_) );
03735   CHK_ERR( sumIntoMatrix(tmpMat2_) );
03736 
03737   //form tmpMat1_ = D^T*Kdd_
03738   fei::multiply_trans_CSRMat_CSRMat(csrD, csrKdd, tmpMat1_);
03739 
03740   //form tmpMat2_ = tmpMat1_*D = D^T*Kdd_*D
03741   fei::multiply_CSRMat_CSRMat(tmpMat1_, csrD, tmpMat2_);
03742 
03743   if (Filter::logStream() != NULL) {
03744     FEI_OSTREAM& os = *Filter::logStream();
03745     os << "#  tmpMat2_"<<FEI_ENDL << tmpMat2_ << FEI_ENDL;
03746   }
03747 
03748   //finally, accumulate tmpMat2_ = D^T*Kdd_*D into the global system matrix.
03749   CHK_ERR( sumIntoMatrix(tmpMat2_) );
03750 
03751   Kdi_->clear();
03752   Kid_->clear();
03753   Kdd_->clear();
03754   reducedEqnCounter_ = 0;
03755 
03756   return(0);
03757 }
03758 
03759 //------------------------------------------------------------------------------
03760 int LinSysCoreFilter::assembleRHS(int numValues,
03761                                   const int* indices,
03762                                   const double* coefs,
03763                                   int mode) {
03764 //
03765 //This function hands the data off to the routine that finally
03766 //sticks it into the RHS vector.
03767 //
03768 
03769   if (problemStructure_->numSlaveEquations() == 0) {
03770     CHK_ERR( giveToRHS(numValues, coefs, indices, mode) );
03771     return(FEI_SUCCESS);
03772   }
03773 
03774   for(int i = 0; i < numValues; i++) {
03775     int eqn = indices[i];
03776 
03777     if (problemStructure_->isSlaveEqn(eqn)) {
03778       fei::add_entry( fd_, eqn, coefs[i]);
03779       reducedRHSCounter_++;
03780       continue;
03781     }
03782 
03783     CHK_ERR( giveToRHS(1, &(coefs[i]), &eqn, mode ) );
03784   }
03785 
03786   if (reducedRHSCounter_ > 300) CHK_ERR( assembleReducedRHS() );
03787 
03788   return(FEI_SUCCESS);
03789 }
03790 
03791 //------------------------------------------------------------------------------
03792 int LinSysCoreFilter::assembleReducedRHS()
03793 {
03794   fei::FillableMat* D = problemStructure_->getSlaveDependencies();
03795 
03796   csrD = *D;
03797 
03798   //now form tmpVec1_ = D^T*fd_.
03799   fei::multiply_trans_CSRMat_CSVec(csrD, fd_, tmpVec1_);
03800 
03801   CHK_ERR( sumIntoRHS(tmpVec1_) );
03802 
03803   fd_.clear();
03804   reducedRHSCounter_ = 0;
03805 
03806   return(0);
03807 }
03808 
03809 //==============================================================================
03810 void LinSysCoreFilter::debugOutput(const char* mesg) {
03811    if (Filter::logStream() != NULL) {
03812      (*logStream())<<mesg<<FEI_ENDL;
03813    }
03814 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends