44 #define fei_file "fei_LinSysCoreFilter.cpp"
47 #define ASSEMBLE_PUT 0
48 #define ASSEMBLE_SUM 1
58 timesInitializeCalled_(0),
61 newMatrixData_(false),
62 newVectorData_(false),
63 newConstraintData_(false),
65 connectivitiesInitialized_(false),
66 firstRemEqnExchange_(true),
67 needToCallMatrixLoadComplete_(false),
68 resolveConflictRequested_(false),
73 numLocallyOwnedNodes_(0),
75 firstLocalNodeNumber_(-1),
77 tooLateToChooseBlock_(false),
79 localReducedBlkOffset_(0),
80 numLocalReducedEqnBlks_(0),
86 masterRank_(masterRank),
87 problemStructure_(probStruct),
88 matrixAllocated_(false),
94 eqnCommMgr_put_(NULL),
121 char** paramStrings = NULL;
220 std::vector<int> rowLengths;
228 for(
size_t ii=0; ii<rowLengths.size(); ++ii) {
229 numNonzeros += rowLengths[ii];
232 std::vector<int> indices_1D(numNonzeros);
233 std::vector<int*> indices(numReducedEqns);
236 for(
size_t ii=0; ii<rowLengths.size(); ++ii) {
237 indices[ii] = &(indices_1D[offset]);
238 offset += rowLengths[ii];
243 if (maxBlkSize == 1) {
246 debugOutput(
"#LinSysCoreFilter calling point lsc_->setMatrixStructure");
248 &indices[0], &rowLengths[0], &blkSizes[0]) );
251 std::vector<int> blkRowLengths;
252 int* blkIndices_1D = numNonzeros > 0 ?
new int[numNonzeros] : NULL;
258 blkIndices, blkIndices_1D,
259 blkRowLengths, blkSizes) );
263 blkIndices[ii] = &(blkIndices_1D[offset]);
264 offset += blkRowLengths[ii];
267 debugOutput(
"#LinSysCoreFilter calling block lsc_->setMatrixStructure");
269 blkIndices, &blkRowLengths[0], &blkSizes[0]) );
271 if (numLocalReducedEqnBlks_ != 0)
delete [] blkIndices;
272 if (numNonzeros > 0)
delete [] blkIndices_1D;
278 debugOutput(
"#leaving LinSysCoreFilter::initialize");
299 debugOutput(
"#LinSysCoreFilter calling lsc_->setLookup");
308 std::vector<int>& globalBlkEqnOffsets =
320 std::vector<int> reducedEqnOffsets(globalEqnOffsets.size());
321 std::vector<int> reducedBlkEqnOffsets(globalBlkEqnOffsets.size());
322 std::vector<int> reducedNodeOffsets(globalNodeOffsets.size());
329 std::vector<int> tmpSend(2), tmpRecv, tmpRecvLengths;
331 tmpSend[1] = reducedNodeNum;
334 for(
size_t ii=0; ii<globalEqnOffsets.size()-1; ++ii) {
335 reducedEqnOffsets[ii] = tmpRecv[ii*2];
336 reducedNodeOffsets[ii] = tmpRecv[ii*2+1];
341 if (numSlaveEqns > 0) {
342 reducedBlkEqnOffsets[ii] = reducedNodeOffsets[ii];
345 reducedBlkEqnOffsets[ii] = globalBlkEqnOffsets[ii];
352 int blkEqn = globalBlkEqnOffsets[
numProcs_]-1;
353 if (numSlaveEqns > 0) {
354 blkEqn = reducedNodeNum;
359 tmpSend[1] = reducedNodeNum;
367 reducedEqnOffsets[
numProcs_] = tmpSend[0]+1;
368 reducedNodeOffsets[
numProcs_] = tmpSend[1]+1;
369 reducedBlkEqnOffsets[
numProcs_] = tmpSend[2]+1;
371 debugOutput(
"#LinSysCoreFilter calling lsc_->setGlobalOffsets");
373 &reducedNodeOffsets[0],
374 &reducedEqnOffsets[0],
375 &reducedBlkEqnOffsets[0]) );
386 numNodes -= numRemoteNodes;
388 std::vector<int> numElemsPerBlock(numBlocks);
389 std::vector<int> numNodesPerElem(numBlocks);
390 std::vector<int> numDofPerNode(0,1);
392 for(
int blk=0; blk<numBlocks; blk++) {
402 for(
int nn=0; nn<numNodesPerElem[blk]; nn++) {
403 if (fieldsPerNode[nn] <= 0)
ERReturn(-1);
404 numDofPerNode.push_back(0);
405 int indx = numDofPerNode.size()-1;
407 for(
int nf=0; nf<fieldsPerNode[nn]; nf++) {
413 for(
int i=0; i<numBlocks; i++) {
427 numDofPerNode.resize(0);
428 for(
int nn=0; nn<numNodesPerElem[i]; nn++) {
429 if (fieldsPerNode[nn] <= 0)
ERReturn(-1);
430 numDofPerNode.push_back(0);
431 int indx = numDofPerNode.size()-1;
433 for(
int nf=0; nf<fieldsPerNode[nn]; nf++) {
443 for(
int k=0; k<nodesPerElement; k++) {
449 int* nodeNumbers = &cNodeList[0];
450 debugOutput(
"#LinSysCoreFilter calling lsc->setConnectivities");
452 &elemID, &nodeNumbers) );
466 std::vector<int> iwork;
468 std::map<GlobalID,ConstraintType*>::const_iterator
472 while(cr_iter != cr_end) {
474 int numNodesPerCR = constraint.
getMasters().size();
477 std::vector<GlobalID>& nodeID_vec = constraint.
getMasters();
478 GlobalID* nodeIDPtr = &nodeID_vec[0];
480 if ((
int)iwork.size() < 2*numNodesPerCR) {
481 iwork.resize(2*numNodesPerCR);
484 int* nodeList = &(iwork[0]);
486 int* eqnList = nodeList+numNodesPerCR;
489 int* fieldIDs = &fieldIDs_vec[0];
491 for(
int k=0; k<numNodesPerCR; k++) {
512 os <<
"#LinSysCoreFilter calling lsc_->setMultCREqns"<<
FEI_ENDL;
513 os <<
"# multiplier eqn: " << meqn <<
", columns: ";
514 for(
int j=0; j<numNodesPerCR; ++j) os << eqnList[j] <<
" ";
527 debugOutput(
"LinSysCoreFilter calling lscf->setMultCRComplete");
530 fei::console_out() <<
"LinSysCoreFilter::setLinSysCoreCREqns ERROR returned from "
531 <<
"lscf->setMultCRComplete()" <<
FEI_ENDL;
541 while(cr_iter != cr_end) {
543 int numNodesPerCR = crset.
getMasters().size();
545 std::vector<GlobalID>& nodeIDs_vec = crset.
getMasters();
546 GlobalID* nodeIDsPtr = &nodeIDs_vec[0];
548 if ((
int)iwork.size() < 2*numNodesPerCR) {
549 iwork.resize(2*numNodesPerCR);
552 int* nodeList = &(iwork[0]);
554 int* eqnList = nodeList+numNodesPerCR;
557 int* fieldIDs = &fieldIDs_vec[0];
559 for(
int k=0; k<numNodesPerCR; k++) {
580 int fieldID,
int fieldSize,
592 int numParams = fieldSize;
603 for(
int j=0; j<numParams; j++) {
604 cols[j] = eqnNumber + j;
614 int fieldID,
int fieldSize,
615 double* coefs,
int eqn)
626 int numParams = fieldSize;
644 for(
int j=0; j<numParams; j++) {
645 ptRows[j] = eqnNumber + j;
646 values[j] = &(coefs[j]);
656 int fieldID,
int col,
674 for(
int i=0; i<numEqns; i++) {
681 int iField,
int iFieldSize,
684 int jField,
int jFieldSize,
686 double penValue,
double CRValue)
695 int iEqn = -1, jEqn = -1;
699 int iNumParams = iFieldSize;
700 int jNumParams = jFieldSize;
701 if (iNumParams < 1 || jNumParams < 1) {
702 fei::console_out() <<
"FEI ERROR, attempt to store indices for field with non-positive size"
703 <<
" field "<<iField<<
", size "<<iNumParams<<
", field "<<jField<<
", size "
720 for(
int i=0; i<iNumParams; i++) {
721 for(
int j=0; j<jNumParams; j++) {
723 coefs[j] = penValue*iCoefs[i]*jCoefs[j];
729 double rhsValue = penValue*iCoefs[i]*CRValue;
736 int iField,
int iFieldSize,
739 int jField,
int jFieldSize,
741 double penValue,
double CRValue){
748 int iEqn = -1, jEqn = -1;
752 int iNumParams = iFieldSize;
753 int jNumParams = jFieldSize;
754 if (iNumParams < 1 || jNumParams < 1) {
755 fei::console_out() <<
"FEI ERROR, attempt to store indices for field with non-positive size"
756 <<
" field "<<iField<<
", size "<<iNumParams<<
", field "<<jField<<
", size "
767 if ((
int)
dworkSpace_.size() < iNumParams * jNumParams) {
781 for(
int i=0; i<iNumParams; i++) {
782 double* coefPtr = coefList + i*jNumParams;
787 for(
int j=0; j<jNumParams; j++) {
788 if (i==0) cols[j] = jEqn + j;
789 coefPtr[j] = penValue*iCoefs[i]*jCoefs[j];
792 double rhsValue = penValue*iCoefs[i]*CRValue;
815 debugOutput(
"#LinSysCoreFilter leaving resetSystem");
834 debugOutput(
"#LinSysCoreFilter leaving deleteMultCRs");
870 debugOutput(
"#LinSysCoreFilter leaving resetMatrix");
890 debugOutput(
"# LinSysCoreFilter leaving resetRHSVector");
902 os <<
"FEI: resetInitialGuess" <<
FEI_ENDL;
903 os <<
"#value to which initial guess is to be set" <<
FEI_ENDL;
909 if (eqns == NULL || values == NULL)
return(-1);
921 debugOutput(
"# LinSysCoreFilter leaving resetInitialGuess");
930 const int* offsetsIntoField,
931 const double* prescribedValues)
938 fei::console_out() <<
"FEI Warning: loadNodeBCs called for fieldID "<<fieldID
939 <<
", which was defined with size "<<size<<
" (should be positive)."<<
FEI_ENDL;
950 for(
int j=0; j<numNodes; j++) {
952 (*
logStream())<<static_cast<int>(nodeID)<<
" ";
953 (*
logStream())<<offsetsIntoField[j]<<
" "<<prescribedValues[j];
961 offsetsIntoField, prescribedValues);
970 const double *
const *alpha,
971 const double *
const *beta,
972 const double *
const *gamma)
982 for(
int i=0; i<nb; i++) {
1015 const double*
const* elemStiffness,
1016 const double* elemLoad,
1021 << static_cast<int>(elemBlockID) <<
FEI_ENDL
1022 <<
"#elID" << FEI_ENDL << static_cast<int>(elemID) <<
FEI_ENDL;
1025 int numNodes = block->getNumNodesPerElement();
1028 for(
int i=0; i<numNodes; ++i) {
1030 (*
logStream())<<static_cast<int>(nodeID)<<
" ";
1036 elemLoad, elemFormat));
1043 const double*
const* elemStiffness,
1049 <<
"#elID" << FEI_ENDL << static_cast<int>(elemID) <<
FEI_ENDL;
1052 int numNodes = block->getNumNodesPerElement();
1055 for(
int i=0; i<numNodes; ++i) {
1057 (*
logStream())<<static_cast<int>(nodeID)<<
" ";
1070 const double* elemLoad)
1075 <<
"#elID" << FEI_ENDL << static_cast<int>(elemID) <<
FEI_ENDL;
1078 int numNodes = block->getNumNodesPerElement();
1081 for(
int i=0; i<numNodes; ++i) {
1083 (*
logStream())<<static_cast<int>(nodeID)<<
" ";
1096 const double*
const* elemStiffness,
1097 const double* elemLoad,
1108 const double*
const* elemStiffness,
1109 const double* elemLoad,
1131 const double*
const* stiff = NULL;
1132 if (elemStiffness != NULL) stiff = elemStiffness;
1134 const double* load = NULL;
1135 if (elemLoad != NULL) load = elemLoad;
1145 fei::console_out() <<
"LinSysCoreFilter::generalElemInput ERROR, elemFormat="
1146 << elemFormat <<
" not supported."<<
FEI_ENDL;
1160 if (stiff != NULL || load != NULL) {
1164 if (stiff != NULL) {
1165 os <<
"#elem-stiff (after copy into dense-row format)" <<
FEI_ENDL;
1166 for(
int i=0; i<numElemRows; i++) {
1167 const double* stiff_i = stiff[i];
1168 for(
int j=0; j<numElemRows; j++) {
1169 os << stiff_i[j] <<
" ";
1177 for(
int i=0; i<numElemRows; i++) {
1178 os << load[i] <<
" ";
1183 if (stiff != NULL) {
1193 int* blkSizesPtr = blkIndPtr + numBlkElemRows;
1195 bool useBlkEqns =
false;
1196 if (interleave == 0) {
1200 blkIndPtr, blkSizesPtr);
1201 int sumBlkSizes = 0;
1202 for(
int ii=0; ii<numBlkElemRows; ++ii) {
1203 sumBlkSizes += blkSizesPtr[ii];
1205 if (sumBlkSizes == numElemRows) {
1212 interleave, indPtr);
1215 if (stiff != NULL) {
1226 if (interleave == 0) {
1228 numBlkElemRows, blkIndPtr,
1233 numBlkElemRows, blkIndPtr,
1259 const double* rhsEntries)
1269 if (checkNumEqns != numIDs*fieldSize) {
1285 const double* rhsEntries)
1295 if (checkNumEqns != numIDs*fieldSize) {
1312 std::vector<int> essEqns;
1313 std::vector<double> essGamma;
1329 std::vector<double> essAlpha(essEqns.size(), 1);
1333 if (essEqns.size() > 0) {
1336 &essGamma[0], essEqns.size()) );
1339 debugOutput(
"#LinSysCoreFilter leaving implementAllBCs");
1345 const double* alpha,
1346 const double* gamma,
1349 int* cc_eqns =
const_cast<int*
>(eqns);
1350 double* cc_alpha =
const_cast<double*
>(alpha);
1351 double* cc_gamma =
const_cast<double*
>(gamma);
1359 std::vector<int> reducedEqns(numEqns);
1360 for(
int i=0; i<numEqns; i++) {
1372 const int*
const* colIndices,
const int* colIndLens,
1373 const double*
const* BCcoefs)
1380 int* cc_eqns =
const_cast<int*
>(eqns);
1381 int** cc_colIndices =
const_cast<int**
>(colIndices);
1382 int* cc_colIndLens =
const_cast<int*
>(colIndLens);
1383 double** cc_BCcoefs =
const_cast<double**
>(BCcoefs);
1386 cc_colIndLens, cc_BCcoefs) );
1395 if (numMultCRs < 1) {
1399 std::map<GlobalID,ConstraintType*>::const_iterator
1405 std::vector<int>& bcEqnNumbers = bcEqns.
eqnNumbers();
1413 double fei_eps = 1.e-49;
1415 while(cr_iter != cr_end) {
1420 std::vector<GlobalID>& CRNode_vec = multCR.
getMasters();
1421 GlobalID *CRNodePtr = &CRNode_vec[0];
1423 int* CRFieldPtr = &CRField_vec[0];
1425 double* weights = &weights_vec[0];
1428 for(
int j=0; j<lenList; ++j) {
1430 for(
int k=0; k<fieldSize; ++k) {
1431 if (std::abs(weights[offset++] + 1.0) < fei_eps) {
1462 debugOutput(
"#LinSysCoreFilter::exchangeRemoteEquations");
1469 std::vector<int> flags(len), globalFlags(len);
1480 newBCData_ = globalFlags[3] > 0 ?
true :
false;
1491 debugOutput(
"# putting remote contributions into linear system...");
1504 debugOutput(
"#LinSysCoreFilter leaving exchangeRemoteEquations");
1515 std::vector<fei::CSVec*>& recvEqns = eqnCommMgr.
localEqns();
1516 std::vector<std::vector<double>*>& recvRHSs = *(eqnCommMgr.
localRHSsPtr());
1522 double** coefs =
new double*[numRecvEqns];
1524 for(i=0; i<numRecvEqns; i++) {
1525 coefs[i] = &(recvEqns[i]->coefs()[0]);
1528 for(i=0; i<numRecvEqns; i++) {
1530 int eqn = recvEqnNumbers[i];
1532 fei::console_out() <<
"LinSysCoreFilter::unpackRemoteContributions: ERROR, recvEqn "
1539 for(
size_t ii=0; ii<recvEqns[i]->size(); ii++) {
1540 if (coefs[i][ii] > 1.e+200) {
1542 <<
"WARNING, coefs["<<i<<
"]["<<ii<<
"]: " << coefs[i][ii]
1548 if (recvEqns[i]->size() > 0 && newCoefs) {
1551 recvEqns[i]->size(),
1552 &(recvEqns[i]->indices()[0]),
1553 &(coefs[i]), assemblyMode ) );
1560 &eqn, assemblyMode) );
1572 std::vector<double>& essAlpha,
1573 std::vector<double>& essGamma)
1585 std::vector<int>* eqns = &essEqns;
1588 int numEqns = essEqns.size();
1589 eqns =
new std::vector<int>(numEqns);
1590 int* eqnsPtr = &(*eqns)[0];
1592 for(
int ii=0; ii<numEqns; ++ii) {
1603 &essAlpha[0], &essGamma[0],
1607 if (numRemoteEssBCEqns > 0) {
1610 std::vector<int> remEssBCEqnLengths(remEssBCEqnNumbers.size());
1612 int** indices =
new int*[numRemoteEssBCEqns];
1613 double** coefs =
new double*[numRemoteEssBCEqns];
1615 for(
int i=0; i<numRemoteEssBCEqns; i++) {
1616 coefs[i] = &(remEssBCEqns[i]->
coefs()[0]);
1617 indices[i] = &(remEssBCEqns[i]->
indices()[0]);
1618 remEssBCEqnLengths[i] = remEssBCEqns[i]->
size();
1622 &remEssBCEqnNumbers[0], indices,
1623 &remEssBCEqnLengths[0], coefs));
1633 debugOutput(
"#LinSysCoreFilter exchangeRemoteBCs");
1642 const int* CRFields,
1643 const double* CRWeights,
1657 os<<
"#num-nodes"<<FEI_ENDL<<numCRNodes<<
FEI_ENDL;
1660 for(i=0; i<numCRNodes; ++i) {
1662 os << static_cast<int>(nodeID) <<
" ";
1664 os << FEI_ENDL <<
"#fields:"<<
FEI_ENDL;
1665 for(i=0; i<numCRNodes; ++i) os << CRFields[i] <<
" ";
1666 os << FEI_ENDL <<
"#field-sizes:"<<
FEI_ENDL;
1667 for(i=0; i<numCRNodes; ++i) {
1671 os << FEI_ENDL<<
"#weights:"<<
FEI_ENDL;
1673 for(i=0; i<numCRNodes; ++i) {
1675 for(
int j=0; j<size; ++j) {
1676 os << CRWeights[offset++] <<
" ";
1679 os << FEI_ENDL<<
"#CRValue:"<<FEI_ENDL<<CRValue<<
FEI_ENDL;
1689 fei::console_out() <<
"ERROR in FEI, constraint with ID="<<CRID<<
" appears to have"
1690 <<
" a constrained-node list of length "<<lenList<<
", should be > 0."<<
FEI_ENDL;
1697 std::vector<GlobalID>& CRNode_vec = multCR->
getMasters();
1698 GlobalID *CRNodePtr = &CRNode_vec[0];
1700 for(i=0; i<lenList; i++) {
1701 if (CRNodePtr[i] != CRNodes[i]) {
1702 fei::console_out() <<
"ERROR in FEI, constraint with ID="<<CRID<<
" had different node-list"
1703 <<
" in initCRMult than it has in loadCRMult."<<
FEI_ENDL;
1709 int *CRFieldPtr = &CRField_vec[0];
1711 for (i = 0; i < lenList; i++) {
1712 if (CRFieldPtr[i] != CRFields[i]) {
1713 fei::console_out() <<
"ERROR in FEI, constraint with CRID="<<CRID<<
" had different field-list"
1714 <<
" in initCRMult than it has in loadCRMult."<<
FEI_ENDL;
1727 for (i = 0; i < lenList; i++) {
1729 assert(numSolnParams >= 0);
1730 fieldSizes[i] = numSolnParams;
1739 for(i = 0; i < lenList; i++) {
1740 for(
int j = 0; j < fieldSizes[i]; j++) {
1741 CRWeightArray.push_back(CRWeights[offset++]);
1746 catch(std::runtime_error& exc) {
1752 double* CRWeightsPtr = &CRWeightArray[0];
1759 double* zeroPtr = &zero;
1765 for(
int j = 0; j < lenList; j++) {
1766 int myFieldID = CRFields[j];
1773 &(CRWeightsPtr[offset]));
1782 &(CRWeightsPtr[offset]), irow);
1789 offset += fieldSizes[j];
1799 const int* CRFields,
1800 const double* CRWeights,
1816 fei::console_out() <<
"ERROR in FEI, constraint with ID="<<CRID<<
" appears to have"
1817 <<
" a constrained-node list of length "<<lenList<<
", should be > 0."<<
FEI_ENDL;
1824 std::vector<GlobalID>& CRNode_vec = penCR->
getMasters();
1825 GlobalID* CRNodePtr = &CRNode_vec[0];
1827 for(
int j = 0; j < lenList; j++) {
1828 if (CRNodePtr[j] != CRNodes[j]) {
1829 fei::console_out() <<
"ERROR in FEI, constraint with ID="<<CRID<<
" had different node-list"
1830 <<
" in initCRPen than it has in loadCRPen."<<
FEI_ENDL;
1845 for (i = 0; i < lenList; i++) {
1847 assert(numSolnParams >= 0);
1848 fieldSizes[i] = numSolnParams;
1856 for (i = 0; i < lenList; i++) {
1857 for (
int j = 0; j < fieldSizes[i]; j++) {
1858 CRWeightArray.push_back(CRWeights[offset++]);
1863 catch(std::runtime_error& exc) {
1870 double* CRWeightPtr = &CRWeightArray[0];
1872 int ioffset = 0, joffset = 0;
1873 for(i = 0; i < lenList; i++) {
1875 int iField = CRFields[i];
1878 double* iweights = &(CRWeightPtr[ioffset]);
1879 ioffset += fieldSizes[i];
1882 for (
int j = 0; j < lenList; j++) {
1884 int jField = CRFields[j];
1887 double* jweights = &(CRWeightPtr[joffset]);
1888 joffset += fieldSizes[j];
1890 double rhsValue = CRValue;
1891 if (j < lenList-1) {
1898 jNode, jField, fieldSizes[j], jweights,
1899 penValue, rhsValue);
1903 jNode, jField, fieldSizes[j], jweights,
1904 penValue, rhsValue);
1918 if (numParams == 0 || paramStrings == NULL) {
1919 debugOutput(
"#LinSysCoreFilter::parameters--- no parameters.");
1926 if (param1 != NULL) {
1927 if (!strcmp(param1,
"AZ_VBR_MATRIX") ||
1928 !strcmp(param1,
"blockMatrix")) {
1934 numParams, paramStrings);
1935 if (param1 != NULL) {
1936 if (!strcmp(param1,
"AZ_VBR_MATRIX") ||
1937 !strcmp(param1,
"blockMatrix")) {
1944 numParams,paramStrings);
1945 if ( param1 != NULL){
1946 std::string str(param1);
1952 if ( param1 != NULL){
1957 if ( param1 != NULL ){
1958 std::string str(param1);
1966 <<
"# --- numParams: "<< numParams<<
FEI_ENDL;
1967 for(
int i=0; i<numParams; i++){
1968 (*
logStream())<<
"#------ paramStrings["<<i<<
"]: "
1970 if (paramStrings[i][strlen(paramStrings[i])-1] !=
'\n') {
1979 debugOutput(
"#LinSysCoreFilter leaving parameters function");
1988 std::vector<int> flags(len), globalFlags(len);
1999 newBCData_ = globalFlags[3] > 0 ?
true :
false;
2001 bool called_exchange =
false;
2004 called_exchange =
true;
2007 bool called_implbcs =
false;
2010 called_implbcs =
true;
2014 debugOutput(
"#LinSysCoreFilter calling LinSysCore matrixLoadComplete");
2030 int* fieldIDs,
double* norms,
double& residTime)
2038 if (whichNorm < 0 || whichNorm > 2)
return(-1);
2051 fieldIDs, norms, residValues) );
2071 debugOutput(
"#LinSysCoreFilter in solve, calling launchSolver...");
2079 debugOutput(
"#LinSysCoreFilter... back from solver");
2087 if (status != 0)
return(1);
2113 std::vector<int>::iterator iter =
2118 int index = iter -
rhsIDs_.begin();
2128 const int* ptRowNumbers,
2129 const double*
const* coefs,
2132 for(
int i=0; i<numPtRows; i++) {
2133 int row = ptRowNumbers[i];
2134 const double* valptr = coefs[i];
2143 os <<
"# calling sumIntoSystemMatrix, row: " << ptRowNumbers[i]
2145 for(
int j=0; j<numPtRows; ++j) os << ptRowNumbers[j] <<
" ";
2150 numPtRows, ptRowNumbers,
2155 numPtRows, ptRowNumbers,
2165 const int* ptRowNumbers,
2167 const int* blkRowNumbers,
2168 const int* blkRowSizes,
2169 const double*
const* coefs,
2177 for(i=0; i<numPtRows; i++) {
2178 int row = ptRowNumbers[i];
2179 valptr[i] = coefs[i];
2187 numPtRows, ptRowNumbers,
2193 for(i=0; i<numBlkRows; i++) {
2194 int row = ptRowNumbers[offset];
2196 offset += blkRowSizes[i];
2203 os <<
"# calling sumIntoSystemMatrix, row: " << ptRowNumbers[i]
2205 for(
int j=0; j<numPtRows; ++j) os << ptRowNumbers[j] <<
" ";
2210 numPtRows, ptRowNumbers,
2211 1, &(blkRowNumbers[i]),
2212 numBlkRows, blkRowNumbers,
2213 &(valptr[offset])) );
2216 offset += blkRowSizes[i];
2224 int numPtCols,
const int* ptCols,
2225 const double*
const* values,
int mode)
2230 for(
int i=0; i<numPtRows; i++) {
2239 os <<
"# calling sumIntoSystemMatrix, row: " << ptRows[i]
2241 for(
int j=0; j<numPtCols; ++j) os << ptCols[j] <<
" ";
2262 for(
int ii=0; ii<numPtCols; ii++) {
2263 int reducedEqn = -1;
2268 iworkPtr[ii] = reducedEqn;
2271 iworkPtr[ii] = reducedEqn;
2272 iworkPtr2[offset++] = reducedEqn;
2277 for(
int i=0; i<numPtRows; i++) {
2278 int row = ptRows[i];
2282 if (isSlave)
continue;
2287 for(
int j=0; j<numPtCols; j++) {
2290 (*
logStream())<<
"# giveToMatrix remote("<<reducedRow<<
","
2322 for(
int j=0; j<numPtCols; j++) {
2324 int reducedCol = iworkPtr[j];
2325 if (reducedCol<0)
continue;
2327 double* tmpCoef =
const_cast<double*
>(&(values[i][j]));
2330 (*
logStream())<<
"# giveToMatrix local("<<reducedRow
2331 <<
","<<reducedCol<<
","<<*tmpCoef<<
")"<<
FEI_ENDL;
2337 os <<
"# calling sumIntoSystemMatrix, row: " << reducedRow
2338 <<
", columns: " << reducedCol <<
FEI_ENDL;
2353 catch(std::runtime_error& exc) {
2363 int numPtCols,
const int* ptCols,
2364 const double*
const* values,
int mode)
2369 double fei_min = std::numeric_limits<double>::min();
2371 for(
int i=0; i<numPtRows; i++) {
2374 const double* values_i = values[i];
2376 for(
int j=0; j<numPtCols; ++j) {
2377 if (specialCase && std::abs(values_i[j]) < fei_min)
continue;
2379 const double* valPtr = &(values_i[j]);
2401 for(
size_t i=0; i<rowNumbers.size(); ++i) {
2402 int row = rowNumbers[i];
2403 int offset = rowOffsets[i];
2404 int rowlen = rowOffsets[i+1]-offset;
2405 const int* indices = &pckColInds[offset];
2406 const double* coefs = &pckCoefs[offset];
2409 throw std::runtime_error(
"fei::impl_utils::add_to_matrix ERROR in matrix.sumIn.");
2418 const int* rowColOffsets,
const int* ptCols,
2419 int numColsPerRow,
double** values)
2431 for(
int re=0; re<numPtRows; re++) {
2432 int eqn = ptRows[re];
2436 remoteProcEqns.
addEqn(eqn, owner);
2452 std::vector<int>& eqnNumbers = localEqns.eqnNumbers();
2453 fei::CSVec** localEqnsPtr = (localEqns.eqns().size() ? &(localEqns.eqns()[0]) : 0);
2454 std::vector<int> eqnLengths(eqnNumbers.size());
2455 for(
size_t i=0; i<eqnNumbers.size(); ++i) {
2456 eqnLengths[i] = localEqnsPtr[i]->
size();
2469 &remoteProcEqns, &remoteEqns,
false));
2471 std::vector<int>& remEqnNumbers = remoteEqns.eqnNumbers();
2472 fei::CSVec** remEqnsPtr = (remoteEqns.eqns().size() ? &(remoteEqns.eqns()[0]) : 0);
2473 std::vector<fei::CSVec*>& remEqns = remoteEqns.eqns();
2476 for(
int i=0; i<numPtRows; i++) {
2477 int row = ptRows[i];
2482 if (eqnIndex < 0)
continue;
2487 if (ptCols == NULL) {
2488 for(
size_t j=0; j<remEqnsPtr[eqnIndex]->
size(); j++) {
2489 values[i][j] = remEqns[eqnIndex]->coefs()[j];
2494 for(
int j=0; j<numColsPerRow; j++) {
2495 int offset = rowColOffsets[i] + j;
2499 values[i][j] = remEqns[eqnIndex]->coefs()[colIndex];
2504 for(
int i=0; i<numPtRows; i++) {
2505 int row = ptRows[i];
2508 int rowLen = 0, checkRowLen;
2515 std::vector<double> coefs(rowLen);
2516 std::vector<int> indices(rowLen);
2519 if (rowLen != checkRowLen)
ERReturn(-1);
2527 if (ptCols == NULL) {
2528 for(
int j=0; j<rowLen; j++) {
2529 values[i][j] = coefs[j];
2534 for(
int j=0; j<numColsPerRow; j++) {
2535 std::vector<int>::iterator iter =
2536 std::find(indices.begin(), indices.end(), ptCols[rowColOffsets[i]+j]);
2537 if (iter == indices.end()) {
2541 int index = iter - indices.begin();
2542 values[i][j] = coefs[index];
2558 for(
unsigned p=0; p<eqnNumbers.size(); p++) {
2559 for(
unsigned i=0; i<eqnNumbers[p]->size(); i++) {
2560 int eqn = (*(eqnNumbers[p]))[i];
2565 std::vector<int>& eqnDataEqns = eqnData.
eqnNumbers();
2571 if (len <= 0)
continue;
2572 std::vector<double> coefs(len);
2573 std::vector<int> indices(len);
2578 if (outputLen != len)
ERReturn(-1);
2580 CHK_ERR( eqnData.
addEqn(eqn, &coefs[0], &indices[0], len,
false) );
2602 for(
int p=0; p<numSendProcs; p++) {
2603 for(
int i=0; i<eqnsPerProc[p]; i++) {
2610 std::vector<int>& eqnDataEqns = eqnData.
eqnNumbers();
2617 int bogusIndex = 19;
2627 const int* indices,
int mode)
2630 for(
int i=0; i<num; i++) {
2651 for(
int i=0; i<num; i++) {
2654 translateToReducedEqn(indices[i], reducedEqn);
2655 if (isSlave)
continue;
2657 if (reducedEqn < reducedStartRow_ || reducedEqn >
reducedEndRow_) {
2682 const int* indices,
int mode)
2684 for(
int i=0; i<num; i++) {
2700 std::vector<int>& indices = vec.
indices();
2701 std::vector<double>& coefs = vec.
coefs();
2718 for(
int re=0; re<num; re++) {
2719 int eqn = indices[re];
2723 remoteProcEqns.
addEqn(eqn, owner);
2739 std::vector<int>& eqnNumbers = localEqns.eqnNumbers();
2740 fei::CSVec** localEqnsPtr = &(localEqns.eqns()[0]);
2741 std::vector<int> eqnLengths(eqnNumbers.size());
2742 for(
size_t i=0; i<eqnNumbers.size(); ++i) {
2743 eqnLengths[i] = localEqnsPtr[i]->
size();
2756 &remoteProcEqns, &remoteEqns,
false))
2759 std::vector<int>& remEqnNumbers = remoteEqns.eqnNumbers();
2760 std::vector<std::vector<double>*>& remRhsCoefs = *(remoteEqns.rhsCoefsPtr());
2762 for(
int i=0; i<num; i++) {
2763 int row = indices[i];
2766 if (eqnIndex < 0)
continue;
2768 values[i] = (*(remRhsCoefs[eqnIndex]))[0];
2772 for(
int i=0; i<num; i++) {
2818 std::vector<int>* masterEqns = NULL;
2819 std::vector<double>* masterCoefs = NULL;
2823 int len = masterEqns->size();
2828 for(
int i=0; i<len; i++) {
2829 int mEqn = (*masterEqns)[i];
2839 solnValue += coef * (*masterCoefs)[i];
2865 fei::console_out() <<
"LinSysCoreFilter::getSharedRemoteSolnEntry: ERROR, eqn "
2866 << eqnNumber <<
" not found." <<
FEI_ENDL;
2869 solnValue = remoteSoln[index];
2895 (*
logStream())<<
"# entering unpackSolution, outputLevel: "
2908 for(
int i=0; i<numRecvEqns; i++) {
2909 int eqn = recvEqnNumbers[i];
2912 fei::console_out() <<
"LinSysCoreFilter::unpackSolution: ERROR, 'recv' eqn (" << eqn
2913 <<
") out of local range." <<
FEI_ENDL;
2917 double solnValue = 0.0;
2926 debugOutput(
"#LinSysCoreFilter leaving unpackSolution");
2949 if (numActiveNodes <= 0)
return(0);
2951 int numSolnParams = 0;
2960 for(
int i=0; i<numActiveNodes; i++) {
2964 if (offset == numNodes)
break;
2969 offsets[offset++] = numSolnParams;
2993 for(
int j=0; j<numFields; j++) {
3002 for(
int k=0; k<size; k++) {
3004 results[numSolnParams++] = answer;
3010 offsets[numNodes] = numSolnParams;
3026 if (numActiveNodes <= 0)
return(0);
3028 int numSolnParams = 0;
3034 for(
int i=0; i<numActiveNodes; i++) {
3038 if (offset == numNodes)
break;
3043 offsets[offset++] = numSolnParams;
3067 for(
int j=0; j<numFields; j++) {
3075 for(
int k=0; k<size; k++) {
3077 results[numSolnParams++] = answer;
3082 offsets[numNodes] = numSolnParams;
3104 if (numActiveNodes <= 0)
return(0);
3113 fei::console_out() <<
"LinSysCoreFilter::getBlockFieldNodeSolution WARNING: fieldID " << fieldID
3114 <<
" not contained in element-block " << (int)elemBlockID <<
FEI_ENDL;
3121 for(
int i=0; i<numNodes; i++) {
3151 if (!hasField)
continue;
3153 int offset = fieldSize*i;
3154 for(
int j=0; j<fieldSize; j++) {
3155 double answer = 0.0;
3157 results[offset+j] = answer;
3175 if (numActiveNodes <= 0)
return(0);
3183 for(
int i=0; i<numNodes; i++) {
3213 if (!hasField)
continue;
3215 int offset = fieldSize*i;
3216 for(
int j=0; j<fieldSize; j++) {
3217 double answer = 0.0;
3219 results[offset+j] = answer;
3231 const double *estimates) {
3237 if (numActiveNodes <= 0)
return(0);
3249 for(
int i=0; i<numNodes; i++) {
3251 int err = nodeDB.getNodeWithID(nodeIDs[i], node);
3253 if (err != 0)
continue;
3266 int offs = offsets[i];
3268 for(
int j=0; j<numFields; j++) {
3272 for(
int k=0; k<size; k++) {
3275 translateToReducedEqn(fieldEqnNumbers[j]+k, reducedEqn);
3278 &estimates[offs+k], 1) );
3293 const double *estimates)
3299 os <<
"FEI: putBlockFieldNodeSolution" <<
FEI_ENDL;
3300 os <<
"#blkID" << FEI_ENDL << (int)elemBlockID << FEI_ENDL
3301 <<
"#fieldID"<<FEI_ENDL << fieldID << FEI_ENDL
3302 <<
"#fieldSize"<<FEI_ENDL << fieldSize << FEI_ENDL
3303 <<
"#numNodes"<<FEI_ENDL << numNodes << FEI_ENDL
3304 <<
"#nodeIDs" << FEI_ENDL;
3306 for(i=0; i<numNodes; ++i) os << (
int)nodeIDs[i] <<
FEI_ENDL;
3308 for(i=0; i<numNodes*fieldSize; ++i) os << estimates[i] << FEI_ENDL;
3320 std::vector<int> numbers(numNodes);
3325 std::vector<double> data;
3328 assert(fieldSize >= 0);
3329 numbers.resize(numNodes*fieldSize);
3330 data.resize(numNodes*fieldSize);
3335 for(
int i=0; i<numNodes; i++) {
3344 for(
int j=0; j<fieldSize; j++) {
3345 data[count] = estimates[i*fieldSize + j];
3355 &numbers[0], numNodes, estimates));
3368 int& numElemDOFPerElement,
3381 if (numElemDOFPerElement <= 0)
return(0);
3385 std::map<GlobalID,int>& elemIDList = ctable.
elemIDs;
3390 for(
int i=0; i<numElems; i++) {
3391 std::map<GlobalID,int>::const_iterator
3392 iter = elemIDList.find(elemIDs[i]);
3393 if (iter == elemIDList.end())
continue;
3394 int index = iter->second;
3396 int offset = i*numElemDOFPerElement;
3398 for(
int j=0; j<numElemDOFPerElement; j++) {
3399 int eqn = elemDOFEqnNumbers[index] + j;
3403 results[offset+j] = answer;
3415 const double *estimates)
3423 assert(DOFPerElement == dofPerElem);
3424 if (DOFPerElement <= 0)
return(0);
3428 std::map<GlobalID,int>& elemIDList = ctable.
elemIDs;
3433 for(
int i=0; i<numElems; i++) {
3434 std::map<GlobalID,int>::const_iterator
3435 iter = elemIDList.find(elemIDs[i]);
3436 if (iter == elemIDList.end())
continue;
3438 int index = iter->second;
3440 for(
int j=0; j<DOFPerElement; j++) {
3443 translateToReducedEqn(elemDOFEqnNumbers[index] + j, reducedEqn);
3444 double soln = estimates[i*DOFPerElement + j];
3456 double* multipliers)
3459 if (numCRs > multCRsLen) {
3463 std::map<GlobalID, ConstraintType*>::const_iterator
3468 while(cr_iter != cr_end && i < numCRs) {
3471 if (CRID != CRIDs[i]) {
3488 const double *multEstimates)
3492 for(
int j = 0; j < numMultCRs; j++) {
3497 if (eqnNumber < localStartRow_ || eqnNumber >
localEndRow_)
continue;
3508 const double* nodeData)
3519 std::vector<int> nodeNumbers(numNodes);
3521 for(
int i=0; i<numNodes; i++) {
3526 if (nodeNumber < 0) {
3527 fei::console_out() <<
"LinSysCoreFilter::putNodalFieldData ERROR, node with ID "
3528 << (int)nodeIDs[i] <<
" doesn't have an associated nodeNumber "
3529 <<
"assigned. putNodalFieldData shouldn't be called until after the "
3530 <<
"initComplete method has been called." <<
FEI_ENDL;
3534 nodeNumbers[i] = nodeNumber;
3538 &nodeNumbers[0], numNodes, nodeData) );
3547 const double* nodeData)
3558 std::vector<int> eqnNumbers(fieldSize);
3560 for(
int i=0; i<numNodes; i++) {
3566 if (!hasField)
continue;
3568 int reducedEqn = -1;
3570 if (isSlave)
continue;
3574 int localLen = fieldSize;
3575 for(
int j=0; j<fieldSize; j++) {
3576 int thisEqn = reducedEqn+j;
3581 eqnNumbers[j] = reducedEqn+j;
3584 int offset = i*fieldSize;
3586 &nodeData[offset], localLen) );
3595 const int* rowNumbers,
3596 const int* colIndices,
3597 const double*
const* coefs,
3598 bool structurallySymmetric,
3599 int numBlkEqns,
int* blkEqns,
3600 int* blkSizes,
bool useBlkEqns,
3605 bool anySlaves =
false;
3607 if (numSlaveEqns > 0) {
3610 const int* indPtr = colIndices;
3611 for(
int r=0; r<numPtRows; r++) {
3613 if (
rSlave_[r] == 1) anySlaves =
true;
3615 for(
int j=0; j<numPtCols; j++) {
3618 if (isSlave == 1) anySlaves =
true;
3621 if (!structurallySymmetric) indPtr += numPtCols;
3625 if (numSlaveEqns == 0 || !anySlaves) {
3626 if (numSlaveEqns == 0 && structurallySymmetric) {
3629 numBlkEqns, blkEqns, blkSizes,
3641 for(
int i=0; i<numPtRows; i++) {
3642 coefPtr[i] = coefs[i];
3645 if (structurallySymmetric) {
3650 const int* indPtr = colIndices;
3651 for(
int i=0; i<numPtRows; i++) {
3652 int row = rowNumbers[i];
3654 const double* coefPtr1 = coefs[i];
3657 indPtr += numPtCols;
3664 const int* indicesPtr = colIndices;
3665 for(
int i=0; i<numPtRows; i++) {
3666 int row = rowNumbers[i];
3668 const double* coefPtr = coefs[i];
3669 int* colSlave = &
cSlave_[offset];
3670 offset += numPtCols;
3675 for(
int jj=0; jj<numPtCols; jj++) {
3676 int col = indicesPtr[jj];
3686 const int* ii_indicesPtr = colIndices;
3687 for(
int ii=0; ii<numPtRows; ii++) {
3688 int rowi = rowNumbers[ii];
3689 if (
rSlave_[ii] == 1)
continue;
3692 if (index < 0)
continue;
3694 const double* coefs_ii = coefs[ii];
3698 if (!structurallySymmetric) ii_indicesPtr += numPtCols;
3712 if (!structurallySymmetric) indicesPtr += numPtCols;
3772 const double* coefs,
3784 for(
int i = 0; i < numValues; i++) {
3785 int eqn = indices[i];
int storePenNodeData(const NodeDescriptor &iNode, int iField, int iFieldSize, double *iCoefs, const NodeDescriptor &jNode, int jField, int jFieldSize, double *jCoefs, double penValue, double CRValue)
virtual int sumIntoSystemMatrix(int numPtRows, const int *ptRows, int numPtCols, const int *ptCols, int numBlkRows, const int *blkRows, int numBlkCols, const int *blkCols, const double *const *values)=0
virtual int setRHSID(int rhsID)=0
static int exchangeEqnBuffers(MPI_Comm comm, ProcEqns *sendProcEqns, EqnBuffer *sendEqns, ProcEqns *recvProcEqns, EqnBuffer *recvEqns, bool accumulate)
std::vector< fei::CSVec * > & remEssBCEqns()
int giveToLocalReducedMatrix(int numPtRows, const int *ptRows, int numPtCols, const int *ptCols, const double *const *values, int mode)
std::vector< fei::CSVec * > & localEqns()
virtual int loadNodeBCs(int numNodes, const GlobalID *nodeIDs, int fieldID, const int *offsetsIntoField, const double *prescribedValues)
virtual int setMultCREqns(int multCRSetID, int numCRs, int numNodesPerCR, int **nodeNumbers, int **eqnNumbers, int *fieldIDs, int *multiplierEqnNumbers)=0
int getParameters(int &numParams, char **¶mStrings)
int giveToRHS(int num, const double *values, const int *indices, int mode)
int getNumElemDOFPerElement()
std::vector< int > & localEqnNumbers()
int getFieldSize(int fieldID)
int exchangeEqns(std::ostream *dbgOut=NULL)
std::vector< int > iworkSpace_
int assembleEqns(int numPtRows, int numPtCols, const int *rowNumbers, const int *colIndices, const double *const *coefs, bool structurallySymmetric, int numBlkEqns, int *blkEqns, int *blkSizes, bool useBlkEqns, int mode)
bool connectivitiesInitialized_
int giveToLocalReducedRHS(int num, const double *values, const int *indices, int mode)
int Allgatherv(MPI_Comm comm, std::vector< T > &sendbuf, std::vector< int > &recvLengths, std::vector< T > &recvbuf)
virtual int resetMatrix(double s)
int giveToBlkMatrix_symm_noSlaves(int numPtRows, const int *ptRows, int numBlkRows, const int *blkRowNumbers, const int *blkRowSizes, const double *const *coefs, int mode)
std::vector< int > & getMasterFieldIDs()
int getBlockDescriptor_index(int index, BlockDescriptor *&block)
int exchangeRemEssBCs(int *essEqns, int numEssEqns, double *essAlpha, double *essGamma, MPI_Comm comm, std::ostream *dbgOut=NULL)
virtual int getBlockFieldNodeSolution(GlobalID elemBlockID, int fieldID, int numNodes, const GlobalID *nodeIDs, double *results)
bool hasBlockIndex(unsigned blk_idx) const
int assembleReducedEqns()
virtual int sumInElemMatrix(GlobalID elemBlockID, GlobalID elemID, const GlobalID *elemConn, const double *const *elemStiffness, int elemFormat)
virtual int getBlockElemSolution(GlobalID elemBlockID, int numElems, const GlobalID *elemIDs, int &numElemDOFPerElement, double *results)
int addEqn(int eqnNumber, const double *coefs, const int *indices, int len, bool accumulate, bool create_indices_union=false)
std::vector< int > & eqnsPerProcPtr()
void multiply_CSRMat_CSRMat(const CSRMat &A, const CSRMat &B, CSRMat &C, bool storeResultZeros)
int getMasterEqnRHS(int slaveEqn, double &rhsValue)
int getPenConstRecord(int CRID, snl_fei::Constraint< GlobalID > *&penCR)
int getEqnsFromMatrix(ProcEqns &procEqns, EqnBuffer &eqnData)
std::ostream * logStream()
virtual int sumIntoRHSVector(int num, const double *values, const int *indices)=0
virtual int setLookup(Lookup &lookup)=0
const int * getFieldIDList() const
virtual int setLoadVectors(GlobalID elemBlock, int numElems, const GlobalID *elemIDs, const double *const *load, int numEqnsPerElem, const int *const *eqnIndices)=0
virtual int resetSystem(double s)
int getInterleaveStrategy() const
virtual int sumIntoRHS(int IDType, int fieldID, int numIDs, const GlobalID *IDs, const double *rhsEntries)
std::vector< int > & indices()
virtual int putBlockElemSolution(GlobalID elemBlockID, int numElems, const GlobalID *elemIDs, int dofPerElem, const double *estimates)
virtual int setConnectivities(GlobalID elemBlock, int numElements, int numNodesPerElem, const GlobalID *elemIDs, const int *const *connNodes)=0
int giveToMatrix(int numPtRows, const int *ptRows, int numPtCols, const int *ptCols, const double *const *values, int mode)
int firstLocalNodeNumber_
int getMasterEqnCoefs(int slaveEqn, std::vector< double > *&masterCoefs)
virtual int setMultCRComplete()=0
virtual int resetMatrix(double s)=0
void setProcEqnLengths(int *eqnNumbers, int *eqnLengths, int len)
int getNumLocalReducedEqns()
int getEqnNumbers(GlobalID ID, int idType, int fieldID, int &numEqns, int *eqnNumbers)
int mirrorProcEqns(ProcEqns &inProcEqns, ProcEqns &outProcEqns)
virtual int putNodalFieldData(int fieldID, int fieldSize, int *nodeNumbers, int numNodes, const double *data)=0
bool resolveConflictRequested_
void addEqn(int eqnNumber, int proc)
int unpackRemoteContributions(EqnCommMgr &eqnCommMgr, int assemblyMode)
NodeDatabase & getNodeDatabase()
virtual int putNodalFieldData(int fieldID, int numNodes, const GlobalID *nodeIDs, const double *nodeData)
std::vector< int > rowNumbers
std::vector< int > & eqnNumbers()
int resolveConflictingCRs(EqnBuffer &bcEqns)
virtual int implementAllBCs()
virtual int loadCRPen(int CRPenID, int numCRNodes, const GlobalID *CRNodes, const int *CRFields, const double *CRWeights, double CRValue, double penValue)
virtual int putBlockFieldNodeSolution(GlobalID elemBlockID, int fieldID, int numNodes, const GlobalID *nodeIDs, const double *estimates)
const int * getFieldEqnNumbers() const
virtual int loadCRMult(int CRMultID, int numCRNodes, const GlobalID *CRNodes, const int *CRFields, const double *CRWeights, double CRValue)
virtual int getFromRHSVector(int num, double *values, const int *indices)=0
virtual int residualNorm(int whichNorm, int numFields, int *fieldIDs, double *norms, double &residTime)
void multiply_trans_CSRMat_CSVec(const CSRMat &A, const CSVec &x, CSVec &y)
int addRHS(int eqnNumber, int rhsIndex, double value, bool accumulate=true)
fei::FillableMat * getSlaveDependencies()
virtual int parameters(int numParams, const char *const *paramStrings)
virtual int unpackSolution()
void setLinSysCoreCREqns()
std::vector< const double * > dworkSpace2_
#define FEI_BLOCK_DIAGONAL_ROW
bool getFieldEqnNumber(int fieldID, int &eqnNumber) const
int getGlobalMaxBlkSize()
std::vector< double > dworkSpace_
std::vector< int > blkScatterIndices_
EqnCommMgr & getEqnCommMgr()
virtual int getMatrixRowLength(int row, int &length)=0
virtual int parameters(int numParams, const char *const *paramStrings)
virtual int enforceEssentialBC(int *globalEqn, double *alpha, double *gamma, int len)=0
virtual int solve(int &status, double &sTime)
SNL_FEI_Structure * problemStructure_
std::vector< int > elemNumbers
virtual int putBlockNodeSolution(GlobalID elemBlockID, int numNodes, const GlobalID *nodeIDs, const int *offsets, const double *estimates)
int getMasterEqnNumbers(int slaveEqn, std::vector< int > *&masterEqns)
std::vector< int > packedColumnIndices
virtual int resetRHSVector(double s)
std::vector< int > rowOffsets
fei::DirichletBCManager * bcManager_
virtual int setPenCREqns(int penCRSetID, int numCRs, int numNodesPerCR, int **nodeNumbers, int **eqnNumbers, int *fieldIDs)=0
int createEqnCommMgr_put()
virtual int getNodalSolution(int numNodes, const GlobalID *nodeIDs, int *offsets, double *results)
int getFromMatrix(int numPtRows, const int *ptRows, const int *rowColOffsets, const int *ptCols, int numColsPerRow, double **values)
virtual int setGlobalOffsets(int len, int *nodeOffsets, int *eqnOffsets, int *blkEqnOffsets)=0
int ** fieldIDsTablePtr()
void addBCRecords(int numBCs, int IDType, int fieldID, int offsetIntoField, const int *IDs, const double *prescribedValues)
std::vector< GlobalID > & getLocalNodeIDs()
static void copyStiffness(const double *const *elemStiff, int numRows, int elemFormat, double **copy)
int binarySearch(const T &item, const T *list, int len)
std::vector< double > & getMasterWeights()
LinSysCoreFilter(FEI_Implementation *owner, MPI_Comm comm, SNL_FEI_Structure *probStruct, LinearSystemCore *lsc, int masterRank=0)
virtual int putIntoRHSVector(int num, const double *values, const int *indices)=0
int calculateResidualNorms(int whichNorm, int numFields, int *fieldIDs, double *norms, std::vector< double > &residValues)
virtual int resetInitialGuess(double s)
int getMatrixRowLengths(std::vector< int > &rowLengths)
int getConstraintID() const
std::vector< int > & elemDOFEqnNumbers()
#define FEI_ISTRINGSTREAM
int getMatrixStructure(int **colIndices, std::vector< int > &rowLengths)
std::vector< NodeDescriptor * > * elem_conn_ptrs
bool needToCallMatrixLoadComplete_
int finalizeBCEqns(fei::Matrix &matrix, bool throw_if_bc_slave_conflict=false)
void storePenNodeSendData(const NodeDescriptor &iNode, int iField, int iFieldSize, double *iCoefs, const NodeDescriptor &jNode, int jField, int jFieldSize, double *jCoefs, double penValue, double CRValue)
std::vector< std::vector< int > * > & procEqnNumbersPtr()
int assembleRHS(int numValues, const int *indices, const double *coefs, int mode)
void getNodeAtIndex(int i, const NodeDescriptor *&node) const
std::map< GlobalID, snl_fei::Constraint< GlobalID > * > & getPenConstRecords()
int sumIntoMatrix(fei::CSRMat &mat)
int storeNodalColumnCoefs(int eqn, const NodeDescriptor &node, int fieldID, int fieldSize, double *coefs)
std::vector< int > cSlave_
virtual int loadElemBCs(int numElems, const GlobalID *elemIDs, int fieldID, const double *const *alpha, const double *const *beta, const double *const *gamma)
virtual int putIntoSystemMatrix(int numPtRows, const int *ptRows, int numPtCols, const int *ptCols, const double *const *values)=0
GlobalID getGlobalNodeID() const
std::vector< int > & getGlobalEqnOffsets()
void setRHSValue(double rhs)
void getScatterIndices_ID(GlobalID blockID, GlobalID elemID, int interleaveStrategy, int *scatterIndices)
void setEqnCommMgr(EqnCommMgr *eqnCommMgr)
std::map< GlobalID, int > elemIDs
int getEqnSolnEntry(int eqnNumber, double &solnValue)
std::vector< std::vector< double > * > * localRHSsPtr()
const NodeDescriptor & findNodeDescriptor(GlobalID nodeID) const
int addRemoteEqn(int eqnNumber, int destProc, const double *coefs, const int *indices, int num)
int getReducedSolnEntry(int eqnNumber, double &solnValue)
virtual int setMatrixStructure(int **ptColIndices, int *ptRrowLengths, int **blkColIndices, int *blkRowLengths, int *ptRowsPerBlkRow)=0
virtual int setNumRHSVectors(int numRHSs, int *rhsIDs)
virtual int loadComplete()
std::vector< int > & sendEqnNumbersPtr()
std::vector< int > & remEssBCEqnNumbersPtr()
int numLocalReducedEqnBlks_
bool containsField(int fieldID)
virtual int setStiffnessMatrices(GlobalID elemBlock, int numElems, const GlobalID *elemIDs, const double *const *const *stiff, int numEqnsPerElem, const int *const *eqnIndices)=0
int getFromRHS(int num, double *values, const int *indices)
int resetTheMatrix(double s)
int getNumNodesPerElement() const
const char * getParam(const char *key, int numParams, const char *const *paramStrings)
virtual int setCurrentRHS(int rhsID)
std::vector< int > & getGlobalNodeOffsets()
virtual int putNodalFieldSolution(int fieldID, int numNodes, const GlobalID *nodeIDs, const double *nodeData)
int storeNodalRowCoefs(const NodeDescriptor &node, int fieldID, int fieldSize, double *coefs, int eqn)
void storeNodalSendEqn(const NodeDescriptor &node, int fieldID, int col, double *coefs)
std::vector< int > rhsIDs_
virtual int resetRHSVector(double s)=0
int getBlockDescriptor(GlobalID blockID, BlockDescriptor *&block)
SparseRowGraph & getGraph()
std::vector< int > rSlave_
bool translateToReducedEqn(int eqn, int &reducedEqn)
int addIndices(int eqnNumber, const int *indices, int len)
int resetTheRHSVector(double s)
ConnectivityTable & getBlockConnectivity(GlobalID blockID)
std::ostream & console_out()
int timesInitializeCalled_
int GlobalMax(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)
virtual int getNodalFieldSolution(int fieldID, int numNodes, const GlobalID *nodeIDs, double *results)
const NodeDescriptor * findNode(GlobalID nodeID) const
int getSharedRemoteSolnEntry(int eqnNumber, double &solnValue)
virtual ~LinSysCoreFilter()
void separate_BC_eqns(const fei::FillableMat &mat, std::vector< int > &bcEqns, std::vector< double > &bcVals)
virtual int matrixLoadComplete()=0
int generalElemInput(GlobalID elemBlockID, GlobalID elemID, const double *const *elemStiffness, const double *elemLoad, int elemFormat)
int getOwnerProcForEqn(int eqn)
void debugOutput(const char *mesg)
int getAssociatedNodeNumber(int eqnNumber)
virtual int sumInElem(GlobalID elemBlockID, GlobalID elemID, const GlobalID *elemConn, const double *const *elemStiffness, const double *elemLoad, int elemFormat)
void multiply_trans_CSRMat_CSRMat(const CSRMat &A, const CSRMat &B, CSRMat &C, bool storeResultZeros)
int getEqnsFromRHS(ProcEqns &procEqns, EqnBuffer &eqnData)
int localProc(MPI_Comm comm)
void sumInCoef(int row, int col, double coef)
void add_entry(CSVec &vec, int eqn, double coef)
int getNumBlkEqnsPerElement()
int mirrorProcEqnLengths(ProcEqns &inProcEqns, ProcEqns &outProcEqns)
std::vector< int > rowIndices_
virtual int enforceRemoteEssBCs(int numEqns, const int *eqns, const int *const *colIndices, const int *colIndLens, const double *const *coefs)
virtual int getCRMultipliers(int numCRs, const int *CRIDs, double *multipliers)
virtual int deleteMultCRs()
virtual int resetConstraints(double s)=0
GlobalID getGlobalBlockID()
int Bcast(MPI_Comm comm, std::vector< T > &sendbuf, int sourceProc)
void setNumRHSs(int numRHSs)
std::vector< int > scatterIndices_
virtual int putInitialGuess(const int *eqnNumbers, const double *values, int len)=0
int getNumLocalReducedEqnBlks()
virtual int getSolnEntry(int eqnNumber, double &answer)=0
#define FEI_BLOCK_DIAGONAL_COL
int getMultConstRecord(int CRID, snl_fei::Constraint< GlobalID > *&multCR)
virtual int exchangeRemoteEquations()
virtual int enforceRemoteEssBCs(int numEqns, int *globalEqns, int **colIndices, int *colIndLen, double **coefs)=0
int formResidual(double *residValues, int numLocalEqns)
std::vector< int > & getMasters()
int addRemoteRHS(fei::CSVec &vec, int rhsIndex)
int getNumNodeDescriptors() const
virtual int formResidual(double *values, int len)=0
std::map< GlobalID, snl_fei::Constraint< GlobalID > * > & getMultConstRecords()
virtual int getMatrixRow(int row, double *coefs, int *indices, int len, int &rowLength)=0
bool firstRemEqnExchange_
int giveToMatrix_symm_noSlaves(int numPtRows, const int *ptRowNumbers, const double *const *coefs, int mode)
int getNodeWithID(GlobalID nodeID, const NodeDescriptor *&node) const
std::vector< int > & getGlobalBlkEqnOffsets()
int getNumMultConstRecords()
std::vector< double > & getPackedCoefs()
std::vector< int > iworkSpace2_
int getNumEqnsPerElement()
snl_fei::Constraint< GlobalID > ConstraintType
virtual int getBlockNodeSolution(GlobalID elemBlockID, int numNodes, const GlobalID *nodeIDs, int *offsets, double *results)
const char * getParamValue(const char *key, int numParams, const char *const *paramStrings, char separator=' ')
EqnCommMgr * eqnCommMgr_put_
std::vector< GlobalID > & getSharedNodeIDs()
void addSolnValues(int *eqnNumbers, double *values, int num)
int getNodeNumber() const
virtual int enforceEssentialBCs(const int *eqns, const double *alpha, const double *gamma, int numEqns)
double * sendEqnSolnPtr()
virtual int launchSolver(int &solveStatus, int &iterations)=0
virtual int putIntoRHS(int IDType, int fieldID, int numIDs, const GlobalID *IDs, const double *rhsEntries)
virtual int putCRMultipliers(int numMultCRs, const int *CRIDs, const double *multEstimates)
int numProcs(MPI_Comm comm)
std::vector< double > & coefs()
NodeCommMgr & getNodeCommMgr()
virtual int sumInElemRHS(GlobalID elemBlockID, GlobalID elemID, const GlobalID *elemConn, const double *elemLoad)
int getIndexOfBlock(GlobalID blockID) const
int gatherSharedBCs(EqnBuffer &bcEqns)
virtual int exchangeRemoteBCs(std::vector< int > &essEqns, std::vector< double > &essAlpha, std::vector< double > &essGamma)