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