47 #ifndef MUELU_MATLABUTILS_DEF_HPP
48 #define MUELU_MATLABUTILS_DEF_HPP
52 #if !defined(HAVE_MUELU_MATLAB) || !defined(HAVE_MUELU_EPETRA)
53 #error "Muemex types require MATLAB, Epetra and Tpetra."
168 #ifdef HAVE_MUELU_INTREPID2
172 MuemexType getMuemexType<RCP<FieldContainer_ordinal> >() {
return FIELDCONTAINER_ORDINAL; }
204 mxClassID probIDtype = mxGetClassID(mxa);
206 if (probIDtype == mxINT32_CLASS) {
207 rv = *((
int*)mxGetData(mxa));
208 }
else if (probIDtype == mxLOGICAL_CLASS) {
209 rv = (int)*((
bool*)mxGetData(mxa));
210 }
else if (probIDtype == mxDOUBLE_CLASS) {
211 rv = (int)*((
double*)mxGetData(mxa));
212 }
else if (probIDtype == mxUINT32_CLASS) {
213 rv = (int)*((
unsigned int*)mxGetData(mxa));
216 throw std::runtime_error(
"Error: Unrecognized numerical type.");
223 return *((
bool*)mxGetData(mxa));
228 return *((
double*)mxGetPr(mxa));
233 double realpart = real<double>(*((
double*)mxGetPr(mxa)));
234 double imagpart = imag<double>(*((
double*)mxGetPi(mxa)));
241 if (mxGetClassID(mxa) != mxCHAR_CLASS) {
242 throw runtime_error(
"Can't construct string from anything but a char array.");
244 rv = string(mxArrayToString(mxa));
251 int nr = mxGetM(mxa);
252 int nc = mxGetN(mxa);
254 throw std::runtime_error(
"A Xpetra::Map representation from MATLAB must be a single row vector.");
255 double* pr = mxGetPr(mxa);
258 std::vector<mm_GlobalOrd> localGIDs(numGlobalIndices);
259 for (
int i = 0; i < int(numGlobalIndices); i++) {
260 localGIDs[i] = Teuchos::as<mm_GlobalOrd>(pr[i]);
272 throw runtime_error(
"Failed to create Xpetra::Map.");
279 if (mxGetN(mxa) != 1 && mxGetM(mxa) != 1)
280 throw std::runtime_error(
"An OrdinalVector from MATLAB must be a single row or column vector.");
281 mm_GlobalOrd numGlobalIndices = mxGetM(mxa) * mxGetN(mxa);
283 if (mxGetClassID(mxa) != mxINT32_CLASS)
284 throw std::runtime_error(
"Can only construct LOVector with int32 data.");
285 int* array = (
int*)mxGetData(mxa);
287 throw runtime_error(
"Failed to create map for Xpetra ordinal vector.");
290 throw runtime_error(
"Failed to create ordinal vector with Xpetra::VectorFactory.");
291 for (
int i = 0; i < int(numGlobalIndices); i++) {
292 loVec->replaceGlobalValue(i, 0, array[i]);
301 int nr = mxGetM(mxa);
302 int nc = mxGetN(mxa);
303 double* pr = mxGetPr(mxa);
310 }
catch (std::exception& e) {
311 mexPrintf(
"Error constructing Tpetra MultiVector.\n");
312 std::cout << e.what() << std::endl;
321 int nr = mxGetM(mxa);
322 int nc = mxGetN(mxa);
323 double* pr = mxGetPr(mxa);
324 double* pi = mxGetPi(mxa);
330 for (
int n = 0; n < nc; n++) {
331 for (
int m = 0; m < nr; m++) {
332 myArr[n * nr + m] =
complex_t(pr[n * nr + m], pi[n * nr + m]);
337 }
catch (std::exception& e) {
338 mexPrintf(
"Error constructing Tpetra MultiVector.\n");
339 std::cout << e.what() << std::endl;
346 bool success =
false;
355 const size_t numGlobalIndices = mxGetM(mxa);
358 double* valueArray = mxGetPr(mxa);
359 int nc = mxGetN(mxa);
365 rowind = (
int*)mxGetIr(mxa);
366 colptr = (
int*)mxGetJc(mxa);
370 for (
int i = 0; i < nc; i++) {
371 for (
int j = colptr[i]; j < colptr[i + 1]; j++) {
372 rowCounts[rowind[j]]++;
376 for (
int i = 0; i < nc; i++) {
377 for (
int j = colptr[i]; j < colptr[i + 1]; j++) {
382 A->insertGlobalValues(rowind[j], cols, vals);
385 A->fillComplete(domainMap, rowMap);
393 }
catch (std::exception& e) {
395 if (rowind != NULL)
delete[] rowind;
396 if (colptr != NULL)
delete[] colptr;
400 mexPrintf(
"Error while constructing Tpetra matrix:\n");
401 std::cout << e.what() << std::endl;
404 mexErrMsgTxt(
"An error occurred while setting up a Tpetra matrix.\n");
418 double* realArray = mxGetPr(mxa);
419 double* imagArray = mxGetPi(mxa);
422 int nc = mxGetN(mxa);
428 rowind = (
int*)mxGetIr(mxa);
429 colptr = (
int*)mxGetJc(mxa);
433 for (
int i = 0; i < nc; i++) {
434 for (
int j = colptr[i]; j < colptr[i + 1]; j++) {
435 rowCounts[rowind[j]]++;
439 for (
int i = 0; i < nc; i++) {
440 for (
int j = colptr[i]; j < colptr[i + 1]; j++) {
443 complex_t value = std::complex<double>(realArray[j], imagArray[j]);
446 A->insertGlobalValues(rowind[j], cols, vals);
449 A->fillComplete(domainMap, rowMap);
454 }
catch (std::exception& e) {
455 mexPrintf(
"Error while constructing tpetra matrix:\n");
456 std::cout << e.what() << std::endl;
464 return MueLu::TpetraCrs_To_XpetraMatrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tmat);
470 return MueLu::TpetraCrs_To_XpetraMatrix<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tmat);
476 return MueLu::TpetraMultiVector_To_XpetraMultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tpetraMV);
482 return MueLu::TpetraMultiVector_To_XpetraMultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tpetraMV);
491 double* vals = mxGetPr(mxa);
492 int nr = mxGetM(mxa);
493 int nc = mxGetN(mxa);
498 rowind = (
int*)mxGetIr(mxa);
499 colptr = (
int*)mxGetJc(mxa);
506 for (
int i = 0; i < nc; i++) {
507 for (
int j = colptr[i]; j < colptr[i + 1]; j++) {
517 }
catch (std::exception& e) {
518 mexPrintf(
"An error occurred while setting up an Epetra matrix:\n");
519 std::cout << e.what() << std::endl;
526 int nr = mxGetM(mxa);
527 int nc = mxGetN(mxa);
535 if (mxGetNumberOfElements(mxa) != 1)
536 throw runtime_error(
"Aggregates must be individual structs in MATLAB.");
537 if (!mxIsStruct(mxa))
538 throw runtime_error(
"Trying to pull aggregates from non-struct MATLAB object.");
541 const int correctNumFields = 5;
542 if (mxGetNumberOfFields(mxa) != correctNumFields)
543 throw runtime_error(
"Aggregates structure has wrong number of fields.");
545 int nVert = *(
int*)mxGetData(mxGetField(mxa, 0,
"nVertices"));
546 int nAgg = *(
int*)mxGetData(mxGetField(mxa, 0,
"nAggregates"));
550 int myRank = comm->getRank();
554 agg->SetNumAggregates(nAgg);
562 mxArray* vertToAggID_in = mxGetField(mxa, 0,
"vertexToAggID");
563 int* vertToAggID_inArray = (
int*)mxGetData(vertToAggID_in);
564 mxArray* rootNodes_in = mxGetField(mxa, 0,
"rootNodes");
565 int* rootNodes_inArray = (
int*)mxGetData(rootNodes_in);
566 for (
int i = 0; i < nVert; i++) {
567 vertex2AggId[i] = vertToAggID_inArray[i];
568 procWinner[i] = myRank;
569 agg->SetIsRoot(i,
false);
571 for (
int i = 0; i < nAgg; i++)
573 agg->SetIsRoot(rootNodes_inArray[i],
true);
576 agg->ComputeAggregateSizes(
true);
577 agg->AggregatesCrossProcessors(
false);
584 throw runtime_error(
"AmalgamationInfo not supported in Muemex yet.");
589 RCP<MGraph> loadDataFromMatlab<RCP<MGraph> >(
const mxArray* mxa) {
591 mxArray* edges = mxGetField(mxa, 0,
"edges");
592 mxArray* boundaryNodes = mxGetField(mxa, 0,
"boundaryNodes");
594 throw runtime_error(
"Graph structure in MATLAB must have a field called 'edges' (logical sparse matrix)");
595 if (boundaryNodes == NULL)
596 throw runtime_error(
"Graph structure in MATLAB must have a field called 'boundaryNodes' (int32 array containing list of boundary nodes)");
597 int* boundaryList = (
int*)mxGetData(boundaryNodes);
598 if (!mxIsSparse(edges) || mxGetClassID(edges) != mxLOGICAL_CLASS)
599 throw runtime_error(
"Graph edges must be stored as a logical sparse matrix.");
601 mwIndex* matlabColPtrs = mxGetJc(edges);
602 mwIndex* matlabRowIndices = mxGetIr(edges);
609 int nnz = matlabColPtrs[mxGetN(edges)];
610 for (
int i = 0; i < nnz; i++)
611 entriesPerRow[matlabRowIndices[i]]++;
616 for (
int i = 0; i < nRows; i++)
617 rows[i + 1] = rows[i] + entriesPerRow[i];
620 int ncols = mxGetN(edges);
621 for (
int colNum = 0; colNum < ncols; ++colNum) {
622 int ci = matlabColPtrs[colNum];
623 for (
int j = ci; j < Teuchos::as<int>(matlabColPtrs[colNum + 1]); ++j) {
624 int rowNum = matlabRowIndices[j];
625 cols[rows[rowNum] + insertionsPerRow[rowNum]] = colNum;
626 insertionsPerRow[rowNum]++;
631 for (
int i = 0; i < nRows; i++) {
632 if (maxNzPerRow < entriesPerRow[i])
633 maxNzPerRow = entriesPerRow[i];
642 for (
int i = 0; i < nRows; ++i) {
643 tgraph->insertGlobalIndices((
mm_GlobalOrd)i, cols(rows[i], entriesPerRow[i]));
645 tgraph->fillComplete(map, map);
648 int numBoundaryNodes = mxGetNumberOfElements(boundaryNodes);
649 bool* boundaryFlags =
new bool[nRows];
650 for (
int i = 0; i < nRows; i++) {
651 boundaryFlags[i] =
false;
653 for (
int i = 0; i < numBoundaryNodes; i++) {
654 boundaryFlags[boundaryList[i]] =
true;
657 mgraph->SetBoundaryNodeMap(boundaryNodesInput);
661 #ifdef HAVE_MUELU_INTREPID2
664 if (mxGetClassID(mxa) != mxINT32_CLASS)
665 throw runtime_error(
"FieldContainer must have integer storage entries");
667 int* data = (
int*)mxGetData(mxa);
668 int nr = mxGetM(mxa);
669 int nc = mxGetN(mxa);
672 for (
int col = 0; col < nc; col++) {
673 for (
int row = 0; row < nr; row++) {
674 (*fc)(row, col) = data[col * nr + row];
687 mwSize dims[] = {1, 1};
688 mxArray* mxa = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
689 *((
int*)mxGetData(mxa)) = data;
695 mwSize dims[] = {1, 1};
696 mxArray* mxa = mxCreateLogicalArray(2, dims);
697 *((
bool*)mxGetData(mxa)) = data;
703 return mxCreateDoubleScalar(data);
708 mwSize dims[] = {1, 1};
709 mxArray* mxa = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxCOMPLEX);
710 *((
double*)mxGetPr(mxa)) = real<double>(data);
711 *((
double*)mxGetPi(mxa)) = imag<double>(data);
717 return mxCreateString(data.c_str());
723 int nc = data->getGlobalNumElements();
726 double* array = (
double*)malloc(
sizeof(
double) * nr * nc);
727 for (
int col = 0; col < nc; col++) {
729 array[col] = Teuchos::as<double>(gid);
738 mwSize len = data->getGlobalLength();
740 mwSize dimensions[] = {len, 1};
741 mxArray* rv = mxCreateNumericArray(2, dimensions, mxINT32_CLASS, mxREAL);
742 int* dataPtr = (
int*)mxGetData(rv);
744 for (
int i = 0; i < int(data->getGlobalLength()); i++) {
778 Teuchos::rcp_const_cast<
Xpetra_CrsGraph>(data->getCrsGraph())->computeGlobalConstants();
780 int nr = data->getGlobalNumRows();
781 int nc = data->getGlobalNumCols();
782 int nnz = data->getGlobalNumEntries();
784 #ifdef VERBOSE_OUTPUT
788 mxArray* mxa = createMatlabSparse<Scalar>(nr, nc, nnz);
791 for (
int i = 0; i < nc + 1; i++) {
795 size_t maxEntriesPerRow = data->getGlobalMaxNumRowEntries();
798 int* rowProgress =
new int[nc];
800 Scalar* sparseVals =
new Scalar[nnz];
802 if (data->isLocallyIndexed()) {
803 Scalar* rowValArray =
new Scalar[maxEntriesPerRow];
809 data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
810 for (
mm_LocalOrd entry = 0; entry < int(numEntries); entry++)
812 jc[rowIndices[entry] + 1]++;
817 int entriesAccum = 0;
818 for (
int n = 0; n <= nc; n++) {
819 int temp = entriesAccum;
820 entriesAccum += jc[n];
824 for (
int i = 0; i < nc; i++) {
830 data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
835 sparseVals[jc[col] + rowProgress[col]] = rowVals[i];
836 ir[jc[col] + rowProgress[col]] = m;
840 delete[] rowIndicesArray;
845 data->getGlobalRowView(m, rowIndices, rowVals);
847 jc[rowIndices[n] + 1]++;
853 for (
int i = 0; i < nc; i++) {
856 int entriesAccum = 0;
857 for (
int n = 0; n <= nc; n++) {
858 int temp = entriesAccum;
859 entriesAccum += jc[n];
865 data->getGlobalRowView(m, rowIndices, rowVals);
869 sparseVals[jc[col] + rowProgress[col]] = rowVals[i];
870 ir[jc[col] + rowProgress[col]] = m;
876 fillMatlabArray<Scalar>(sparseVals, mxa, nnz);
878 delete[] rowProgress;
887 Teuchos::rcp_const_cast<
Xpetra_CrsGraph>(data->getCrsGraph())->computeGlobalConstants();
889 int nr = data->getGlobalNumRows();
890 int nc = data->getGlobalNumCols();
891 int nnz = data->getGlobalNumEntries();
892 #ifdef VERBOSE_OUTPUT
896 mxArray* mxa = createMatlabSparse<Scalar>(nr, nc, nnz);
899 for (
int i = 0; i < nc + 1; i++) {
902 size_t maxEntriesPerRow = data->getGlobalMaxNumRowEntries();
903 int* rowProgress =
new int[nc];
905 Scalar* sparseVals =
new Scalar[nnz];
907 if (data->isLocallyIndexed()) {
908 Scalar* rowValArray =
new Scalar[maxEntriesPerRow];
914 data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
915 for (
mm_LocalOrd entry = 0; entry < int(numEntries); entry++)
917 jc[rowIndices[entry] + 1]++;
921 int entriesAccum = 0;
922 for (
int n = 0; n <= nc; n++) {
923 int temp = entriesAccum;
924 entriesAccum += jc[n];
928 for (
int i = 0; i < nc; i++) {
934 data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
939 sparseVals[jc[col] + rowProgress[col]] = rowVals[i];
940 ir[jc[col] + rowProgress[col]] = m;
944 delete[] rowIndicesArray;
949 data->getGlobalRowView(m, rowIndices, rowVals);
951 jc[rowIndices[n] + 1]++;
957 for (
int i = 0; i < nc; i++) {
960 int entriesAccum = 0;
961 for (
int n = 0; n <= nc; n++) {
962 int temp = entriesAccum;
963 entriesAccum += jc[n];
969 data->getGlobalRowView(m, rowIndices, rowVals);
973 sparseVals[jc[col] + rowProgress[col]] = rowVals[i];
974 ir[jc[col] + rowProgress[col]] = m;
980 fillMatlabArray<Scalar>(sparseVals, mxa, nnz);
982 delete[] rowProgress;
1012 int nr = data->getGlobalLength();
1013 int nc = data->getNumVectors();
1015 double* array = (
double*)malloc(
sizeof(
double) * nr * nc);
1016 for (
int col = 0; col < nc; col++) {
1018 for (
int row = 0; row < nr; row++) {
1019 array[col * nr + row] = colData[row];
1030 int nr = data->getGlobalLength();
1031 int nc = data->getNumVectors();
1034 for (
int col = 0; col < nc; col++) {
1036 for (
int row = 0; row < nr; row++) {
1037 array[col * nr + row] = colData[row];
1053 mxArray* output = mxCreateDoubleMatrix(data->GlobalLength(), data->NumVectors(), mxREAL);
1054 double* dataPtr = mxGetPr(output);
1055 data->ExtractCopy(dataPtr, data->GlobalLength());
1062 int numNodes = data->GetVertex2AggId()->getData(0).size();
1063 int numAggs = data->GetNumAggregates();
1065 mwSize singleton[] = {1, 1};
1066 dataIn[0] = mxCreateNumericArray(2, singleton, mxINT32_CLASS, mxREAL);
1067 *((
int*)mxGetData(dataIn[0])) = numNodes;
1068 dataIn[1] = mxCreateNumericArray(2, singleton, mxINT32_CLASS, mxREAL);
1069 *((
int*)mxGetData(dataIn[1])) = numAggs;
1070 mwSize nodeArrayDims[] = {(mwSize)numNodes, 1};
1071 dataIn[2] = mxCreateNumericArray(2, nodeArrayDims, mxINT32_CLASS, mxREAL);
1072 int* vtaid = (
int*)mxGetData(dataIn[2]);
1074 for (
int i = 0; i < numNodes; i++) {
1075 vtaid[i] = vertexToAggID[i];
1077 mwSize aggArrayDims[] = {(mwSize)numAggs, 1};
1078 dataIn[3] = mxCreateNumericArray(2, aggArrayDims, mxINT32_CLASS, mxREAL);
1081 for (
int i = 0; i < numNodes; i++) {
1082 if (data->IsRoot(i))
1085 bool reassignRoots =
false;
1086 if (totalRoots != numAggs) {
1088 <<
"Warning: Number of root nodes and number of aggregates do not match." << endl;
1089 cout <<
"Will reassign root nodes when writing aggregates to matlab." << endl
1091 reassignRoots =
true;
1093 int* rn = (
int*)mxGetData(dataIn[3]);
1095 if (reassignRoots) {
1097 int lastFoundNode = 0;
1098 for (
int i = 0; i < numAggs; i++) {
1100 for (
int j = lastFoundNode; j < lastFoundNode + numNodes; j++) {
1101 int index = j % numNodes;
1102 if (vertexToAggID[index] == i) {
1104 lastFoundNode = index;
1107 TEUCHOS_TEST_FOR_EXCEPTION(rn[i] == -1, runtime_error,
"Invalid aggregates: Couldn't find any node in aggregate #" << i <<
".");
1111 for (
int j = 0; j < numNodes; j++) {
1112 if (data->IsRoot(j)) {
1114 throw runtime_error(
"Cannot store invalid aggregates in MATLAB - more root nodes than aggregates.");
1119 if (i + 1 < numAggs)
1120 throw runtime_error(
"Cannot store invalid aggregates in MATLAB - fewer root nodes than aggregates.");
1123 dataIn[4] = mxCreateNumericArray(1, aggArrayDims, mxINT32_CLASS, mxREAL);
1124 int*
as = (
int*)mxGetData(dataIn[4]);
1126 for (
int i = 0; i < numAggs; i++) {
1127 as[i] = aggSizes[i];
1129 mxArray* matlabAggs[1];
1130 int result = mexCallMATLAB(1, matlabAggs, 5, dataIn,
"constructAggregates");
1132 throw runtime_error(
"Matlab encountered an error while constructing aggregates struct.");
1133 return matlabAggs[0];
1138 throw runtime_error(
"AmalgamationInfo not supported in MueMex yet.");
1139 return mxCreateDoubleScalar(0);
1144 int numEntries = (int)data->GetGlobalNumEdges();
1145 int numRows = (int)data->GetDomainMap()->getGlobalNumElements();
1146 mxArray* mat = mxCreateSparseLogicalMatrix(numRows, numRows, numEntries);
1147 mxLogical* outData = (mxLogical*)mxGetData(mat);
1148 mwIndex* rowInds = mxGetIr(mat);
1149 mwIndex* colPtrs = mxGetJc(mat);
1152 int* entriesPerRow =
new int[numRows];
1153 int* entriesPerCol =
new int[numRows];
1154 for (
int i = 0; i < numRows; i++) {
1155 entriesPerRow[i] = 0;
1156 entriesPerCol[i] = 0;
1158 for (
int i = 0; i < numRows; i++) {
1161 entriesPerRow[i] = neighbors.
size();
1162 for (
int j = 0; j < neighbors.
size(); j++) {
1163 entriesPerCol[neighbors[j]]++;
1165 iter += neighbors.
size();
1168 mxLogical** valuesByColumn =
new mxLogical*[numRows];
1169 int* numEnteredPerCol =
new int[numRows];
1171 for (
int i = 0; i < numRows; i++) {
1172 rowIndsByColumn[i] = &rowInds[accum];
1174 valuesByColumn[i] = &outData[accum];
1175 accum += entriesPerCol[i];
1176 if (accum > numEntries)
1177 throw runtime_error(
"potato");
1179 for (
int i = 0; i < numRows; i++) {
1180 numEnteredPerCol[i] = 0;
1184 for (
int row = 0; row < numRows; row++) {
1185 for (
int entryInRow = 0; entryInRow < entriesPerRow[row]; entryInRow++) {
1186 int col = dataCopy[accum];
1188 rowIndsByColumn[col][numEnteredPerCol[col]] = row;
1189 valuesByColumn[col][numEnteredPerCol[col]] = (mxLogical)1;
1190 numEnteredPerCol[col]++;
1194 for (
int col = 0; col < numRows; col++) {
1195 colPtrs[col] = accum;
1196 accum += entriesPerCol[col];
1198 colPtrs[numRows] = accum;
1199 delete[] numEnteredPerCol;
1200 delete[] rowIndsByColumn;
1201 delete[] valuesByColumn;
1203 delete[] entriesPerRow;
1204 delete[] entriesPerCol;
1207 int numBoundaryNodes = 0;
1208 for (
int i = 0; i < boundaryFlags.
size(); i++) {
1209 if (boundaryFlags[i])
1212 cout <<
"Graph has " << numBoundaryNodes <<
" Dirichlet boundary nodes." << endl;
1213 mwSize dims[] = {(mwSize)numBoundaryNodes, 1};
1214 mxArray* boundaryList = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
1215 int* dest = (
int*)mxGetData(boundaryList);
1216 int* destIter = dest;
1217 for (
int i = 0; i < boundaryFlags.
size(); i++) {
1218 if (boundaryFlags[i]) {
1223 mxArray* constructArgs[] = {mat, boundaryList};
1225 mexCallMATLAB(1, out, 2, constructArgs,
"constructGraph");
1229 #ifdef HAVE_MUELU_INTREPID2
1232 int rank = data->rank();
1235 throw std::runtime_error(
"Error: Only rank two FieldContainers are supported.");
1237 int nr = data->extent(0);
1238 int nc = data->extent(1);
1240 mwSize dims[] = {(mwSize)nr, (mwSize)nc};
1241 mxArray* mxa = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
1242 int* array = (
int*)mxGetData(mxa);
1244 for (
int col = 0; col < nc; col++) {
1245 for (
int row = 0; row < nr; row++) {
1246 array[col * nr + row] = (*data)(row, col);
1253 template <
typename T>
1256 data = loadDataFromMatlab<T>(mxa);
1259 template <
typename T>
1261 return saveDataToMatlab<T>(data);
1264 template <
typename T>
1270 template <
typename T>
1274 template <
typename T>
1279 template <
typename T>
1281 this->data = newData;
1288 template <
typename T>
1291 lvl.
Set<T>(name, data, fact);
1294 template <
typename T>
1297 return lvl.
Get<T>(name);
1298 }
catch (std::exception& e) {
1299 throw std::runtime_error(
"Requested custom variable " + name +
" is not in the level.");
1304 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
1306 using namespace std;
1307 using namespace Teuchos;
1314 vector<RCP<MuemexArg> > args;
1315 for (
size_t i = 0; i < needsList.size(); i++) {
1316 if (needsList[i] ==
"A" || needsList[i] ==
"P" || needsList[i] ==
"R" || needsList[i] ==
"Ptent") {
1317 Matrix_t mydata = lvl.
Get<Matrix_t>(needsList[i], factory->
GetFactory(needsList[i]).get());
1319 }
else if (needsList[i] ==
"Nullspace" || needsList[i] ==
"Coordinates") {
1320 MultiVector_t mydata = lvl.
Get<MultiVector_t>(needsList[i], factory->
GetFactory(needsList[i]).get());
1322 }
else if (needsList[i] ==
"Aggregates") {
1323 Aggregates_t mydata = lvl.
Get<Aggregates_t>(needsList[i], factory->
GetFactory(needsList[i]).get());
1325 }
else if (needsList[i] ==
"UnAmalgamationInfo") {
1326 AmalgamationInfo_t mydata = lvl.
Get<AmalgamationInfo_t>(needsList[i], factory->
GetFactory(needsList[i]).get());
1328 }
else if (needsList[i] ==
"Level") {
1331 }
else if (needsList[i] ==
"Graph") {
1332 Graph_t mydata = lvl.
Get<Graph_t>(needsList[i], factory->
GetFactory(needsList[i]).get());
1335 vector<string> words;
1336 string badNameMsg =
"Custom Muemex variables in \"Needs\" list require a type and a name, e.g. \"double myVal\". \n Leading and trailing spaces are OK. \n Don't know how to handle \"" + needsList[i] +
"\".\n";
1338 char* buf = (
char*)malloc(needsList[i].size() + 1);
1339 strcpy(buf, needsList[i].c_str());
1340 for (
char* iter = buf; *iter !=
' '; iter++) {
1343 throw runtime_error(badNameMsg);
1347 const char* wordDelim =
" ";
1348 char* mark = strtok(buf, wordDelim);
1349 while (mark != NULL) {
1350 string wordStr(mark);
1351 words.push_back(wordStr);
1352 mark = strtok(NULL, wordDelim);
1354 if (words.size() != 2) {
1356 throw runtime_error(badNameMsg);
1358 char* typeStr = (
char*)words[0].c_str();
1359 if (strstr(typeStr,
"ordinalvector")) {
1361 LOVector_t mydata = getLevelVariable<LOVector_t>(needsList[i], lvl);
1363 }
else if (strstr(typeStr,
"map")) {
1365 Map_t mydata = getLevelVariable<Map_t>(needsList[i], lvl);
1367 }
else if (strstr(typeStr,
"scalar")) {
1368 Scalar mydata = getLevelVariable<Scalar>(needsList[i], lvl);
1370 }
else if (strstr(typeStr,
"double")) {
1371 double mydata = getLevelVariable<double>(needsList[i], lvl);
1373 }
else if (strstr(typeStr,
"complex")) {
1374 complex_t mydata = getLevelVariable<complex_t>(needsList[i], lvl);
1376 }
else if (strstr(typeStr,
"matrix")) {
1377 Matrix_t mydata = getLevelVariable<Matrix_t>(needsList[i], lvl);
1379 }
else if (strstr(typeStr,
"multivector")) {
1380 MultiVector_t mydata = getLevelVariable<MultiVector_t>(needsList[i], lvl);
1382 }
else if (strstr(typeStr,
"int")) {
1383 int mydata = getLevelVariable<int>(needsList[i], lvl);
1385 }
else if (strstr(typeStr,
"string")) {
1386 string mydata = getLevelVariable<string>(needsList[i], lvl);
1390 throw std::runtime_error(words[0] +
" is not a known variable type.");
1398 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
1400 using namespace std;
1401 using namespace Teuchos;
1408 for (
size_t i = 0; i < size_t(provides.size()); i++) {
1409 if (provides[i] ==
"A" || provides[i] ==
"P" || provides[i] ==
"R" || provides[i] ==
"Ptent") {
1411 lvl.
Set(provides[i], mydata->getData(), factory);
1412 }
else if (provides[i] ==
"Nullspace" || provides[i] ==
"Coordinates") {
1414 lvl.
Set(provides[i], mydata->getData(), factory);
1415 }
else if (provides[i] ==
"Aggregates") {
1417 lvl.
Set(provides[i], mydata->getData(), factory);
1418 }
else if (provides[i] ==
"UnAmalgamationInfo") {
1420 lvl.
Set(provides[i], mydata->getData(), factory);
1421 }
else if (provides[i] ==
"Graph") {
1423 lvl.
Set(provides[i], mydata->getData(), factory);
1425 vector<string> words;
1426 string badNameMsg =
"Custom Muemex variables in \"Provides\" list require a type and a name, e.g. \"double myVal\". \n Leading and trailing spaces are OK. \n Don't know how to handle \"" + provides[i] +
"\".\n";
1428 char* buf = (
char*)malloc(provides[i].size() + 1);
1429 strcpy(buf, provides[i].c_str());
1430 for (
char* iter = buf; *iter !=
' '; iter++) {
1433 throw runtime_error(badNameMsg);
1437 const char* wordDelim =
" ";
1438 char* mark = strtok(buf, wordDelim);
1439 while (mark != NULL) {
1440 string wordStr(mark);
1441 words.push_back(wordStr);
1442 mark = strtok(NULL, wordDelim);
1444 if (words.size() != 2) {
1446 throw runtime_error(badNameMsg);
1448 const char* typeStr = words[0].c_str();
1449 if (strstr(typeStr,
"ordinalvector")) {
1452 addLevelVariable<LOVector_t>(mydata->getData(), words[1], lvl, factory);
1453 }
else if (strstr(typeStr,
"map")) {
1456 addLevelVariable<Map_t>(mydata->getData(), words[1], lvl, factory);
1457 }
else if (strstr(typeStr,
"scalar")) {
1459 addLevelVariable<Scalar>(mydata->getData(), words[1], lvl, factory);
1460 }
else if (strstr(typeStr,
"double")) {
1462 addLevelVariable<double>(mydata->getData(), words[1], lvl, factory);
1463 }
else if (strstr(typeStr,
"complex")) {
1465 addLevelVariable<complex_t>(mydata->getData(), words[1], lvl, factory);
1466 }
else if (strstr(typeStr,
"matrix")) {
1468 addLevelVariable<Matrix_t>(mydata->getData(), words[1], lvl, factory);
1469 }
else if (strstr(typeStr,
"multivector")) {
1471 addLevelVariable<MultiVector_t>(mydata->getData(), words[1], lvl, factory);
1472 }
else if (strstr(typeStr,
"int")) {
1474 addLevelVariable<int>(mydata->getData(), words[1], lvl, factory);
1475 }
else if (strstr(typeStr,
"bool")) {
1477 addLevelVariable<bool>(mydata->getData(), words[1], lvl, factory);
1478 }
else if (strstr(typeStr,
"string")) {
1480 addLevelVariable<string>(mydata->getData(), words[1], lvl, factory);
1483 throw std::runtime_error(words[0] +
" is not a known variable type.");
1494 throw std::runtime_error(
"Muemex does not support long long for global indices");
1499 throw std::runtime_error(
"Muemex does not support long long for global indices");
1504 throw std::runtime_error(
"Muemex does not support long long for global indices");
1509 throw std::runtime_error(
"Muemex does not support long long for global indices");
1513 #endif // HAVE_MUELU_MATLAB error handler
1514 #endif // MUELU_MATLABUTILS_DEF_HPP guard
mxArray * convertToMatlab()
std::vector< Teuchos::RCP< MuemexArg > > processNeeds< double, mm_LocalOrd, long long, mm_node_t >(const Factory *factory, std::string &needsParam, Level &lvl)
template mxArray * saveDataToMatlab(bool &data)
std::vector< std::string > tokenizeList(const std::string ¶ms)
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
mxArray * createMatlabSparse< complex_t >(int numRows, int numCols, int nnz)
MuemexType getMuemexType< string >()
Tpetra::Map::local_ordinal_type mm_LocalOrd
Tpetra::Map::global_ordinal_type mm_GlobalOrd
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
RCP< Xpetra::MultiVector< SC, LO, GO, NO > > TpetraMultiVector_To_XpetraMultiVector(const Teuchos::RCP< Tpetra::MultiVector< SC, LO, GO, NO >> &Vtpetra)
User data are always kept. This flag is set automatically when Level::Set("data", data) is used...
std::string tolower(const std::string &str)
void fillMatlabArray< double >(double *array, const mxArray *mxa, int n)
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
template string loadDataFromMatlab< string >(const mxArray *mxa)
MuemexType getMuemexType(const T &data)
std::vector< Teuchos::RCP< MuemexArg > > processNeeds(const Factory *factory, std::string &needsParam, Level &lvl)
void addLevelVariable(const T &data, std::string &name, Level &lvl, const FactoryBase *fact=NoFactory::get())
int FillComplete(bool OptimizeDataStorage=true)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void fillMatlabArray< complex_t >(complex_t *array, const mxArray *mxa, int n)
mxArray * createMatlabSparse< double >(int numRows, int numCols, int nnz)
MueLu::DefaultScalar Scalar
int * mwIndex_to_int(int N, mwIndex *mwi_array)
Class that holds all level-specific information.
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalGlobal lg=Xpetra::GloballyDistributed)
void processProvides< complex_t, mm_LocalOrd, long long, mm_node_t >(std::vector< Teuchos::RCP< MuemexArg > > &mexOutput, const Factory *factory, std::string &providesParam, Level &lvl)
template complex_t loadDataFromMatlab< complex_t >(const mxArray *mxa)
template int loadDataFromMatlab< int >(const mxArray *mxa)
RCP< Xpetra::Matrix< SC, LO, GO, NO > > TpetraCrs_To_XpetraMatrix(const Teuchos::RCP< Tpetra::CrsMatrix< SC, LO, GO, NO >> &Atpetra)
template bool loadDataFromMatlab< bool >(const mxArray *mxa)
const T & getLevelVariable(std::string &name, Level &lvl)
MuemexType getMuemexType< int >()
mxArray * createMatlabMultiVector< complex_t >(int numRows, int numCols)
std::vector< Teuchos::RCP< MuemexArg > > processNeeds< complex_t, mm_LocalOrd, long long, mm_node_t >(const Factory *factory, std::string &needsParam, Level &lvl)
void AddKeepFlag(const std::string &ename, const FactoryBase *factory=NoFactory::get(), KeepType keep=MueLu::Keep)
MuemexType getMuemexType< bool >()
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
void processProvides< double, mm_LocalOrd, long long, mm_node_t >(std::vector< Teuchos::RCP< MuemexArg > > &mexOutput, const Factory *factory, std::string &providesParam, Level &lvl)
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Teuchos::RCP< const Teuchos::Comm< int > > getDefaultComm()
std::complex< double > complex_t
mxArray * createMatlabMultiVector< double >(int numRows, int numCols)
template double loadDataFromMatlab< double >(const mxArray *mxa)
MueLu::Aggregates< mm_LocalOrd, mm_GlobalOrd, mm_node_t > MAggregates
TypeTo as(const TypeFrom &t)
MuemexType getMuemexType< double >()
Tpetra::Map muemex_map_type
virtual const RCP< const FactoryBase > GetFactory(const std::string &varName) const =0
int GetLevelID() const
Return level number.
MuemexType getMuemexType< complex_t >()
void processProvides(std::vector< Teuchos::RCP< MuemexArg > > &mexOutput, const Factory *factory, std::string &providesParam, Level &lvl)