10 #ifndef MUELU_MATLABUTILS_DEF_HPP
11 #define MUELU_MATLABUTILS_DEF_HPP
16 #if !defined(HAVE_MUELU_MATLAB) || !defined(HAVE_MUELU_TPETRA)
17 #error "Muemex types require MATLAB and Tpetra."
107 #ifdef HAVE_MUELU_EPETRA
134 #ifdef HAVE_MUELU_INTREPID2
138 MuemexType getMuemexType<RCP<FieldContainer_ordinal> >() {
return FIELDCONTAINER_ORDINAL; }
170 mxClassID probIDtype = mxGetClassID(mxa);
172 if (probIDtype == mxINT32_CLASS) {
173 rv = *((
int*)mxGetData(mxa));
174 }
else if (probIDtype == mxLOGICAL_CLASS) {
175 rv = (int)*((
bool*)mxGetData(mxa));
176 }
else if (probIDtype == mxDOUBLE_CLASS) {
177 rv = (int)*((
double*)mxGetData(mxa));
178 }
else if (probIDtype == mxUINT32_CLASS) {
179 rv = (int)*((
unsigned int*)mxGetData(mxa));
182 throw std::runtime_error(
"Error: Unrecognized numerical type.");
189 return *((
bool*)mxGetData(mxa));
194 return *((
double*)mxGetPr(mxa));
199 double realpart = real<double>(*((
double*)mxGetPr(mxa)));
200 double imagpart = imag<double>(*((
double*)mxGetPi(mxa)));
207 if (mxGetClassID(mxa) != mxCHAR_CLASS) {
208 throw runtime_error(
"Can't construct string from anything but a char array.");
210 rv = string(mxArrayToString(mxa));
217 int nr = mxGetM(mxa);
218 int nc = mxGetN(mxa);
220 throw std::runtime_error(
"A Xpetra::Map representation from MATLAB must be a single row vector.");
221 double* pr = mxGetPr(mxa);
224 std::vector<mm_GlobalOrd> localGIDs(numGlobalIndices);
225 for (
int i = 0; i < int(numGlobalIndices); i++) {
226 localGIDs[i] = Teuchos::as<mm_GlobalOrd>(pr[i]);
238 throw runtime_error(
"Failed to create Xpetra::Map.");
245 if (mxGetN(mxa) != 1 && mxGetM(mxa) != 1)
246 throw std::runtime_error(
"An OrdinalVector from MATLAB must be a single row or column vector.");
247 mm_GlobalOrd numGlobalIndices = mxGetM(mxa) * mxGetN(mxa);
249 if (mxGetClassID(mxa) != mxINT32_CLASS)
250 throw std::runtime_error(
"Can only construct LOVector with int32 data.");
251 int* array = (
int*)mxGetData(mxa);
253 throw runtime_error(
"Failed to create map for Xpetra ordinal vector.");
256 throw runtime_error(
"Failed to create ordinal vector with Xpetra::VectorFactory.");
257 for (
int i = 0; i < int(numGlobalIndices); i++) {
258 loVec->replaceGlobalValue(i, 0, array[i]);
267 int nr = mxGetM(mxa);
268 int nc = mxGetN(mxa);
269 double* pr = mxGetPr(mxa);
276 }
catch (std::exception& e) {
277 mexPrintf(
"Error constructing Tpetra MultiVector.\n");
278 std::cout << e.what() << std::endl;
287 int nr = mxGetM(mxa);
288 int nc = mxGetN(mxa);
289 double* pr = mxGetPr(mxa);
290 double* pi = mxGetPi(mxa);
296 for (
int n = 0; n < nc; n++) {
297 for (
int m = 0; m < nr; m++) {
298 myArr[n * nr + m] =
complex_t(pr[n * nr + m], pi[n * nr + m]);
303 }
catch (std::exception& e) {
304 mexPrintf(
"Error constructing Tpetra MultiVector.\n");
305 std::cout << e.what() << std::endl;
312 bool success =
false;
321 const size_t numGlobalIndices = mxGetM(mxa);
324 double* valueArray = mxGetPr(mxa);
325 int nc = mxGetN(mxa);
331 rowind = (
int*)mxGetIr(mxa);
332 colptr = (
int*)mxGetJc(mxa);
336 for (
int i = 0; i < nc; i++) {
337 for (
int j = colptr[i]; j < colptr[i + 1]; j++) {
338 rowCounts[rowind[j]]++;
342 for (
int i = 0; i < nc; i++) {
343 for (
int j = colptr[i]; j < colptr[i + 1]; j++) {
348 A->insertGlobalValues(rowind[j], cols, vals);
351 A->fillComplete(domainMap, rowMap);
359 }
catch (std::exception& e) {
361 if (rowind != NULL)
delete[] rowind;
362 if (colptr != NULL)
delete[] colptr;
366 mexPrintf(
"Error while constructing Tpetra matrix:\n");
367 std::cout << e.what() << std::endl;
370 mexErrMsgTxt(
"An error occurred while setting up a Tpetra matrix.\n");
384 double* realArray = mxGetPr(mxa);
385 double* imagArray = mxGetPi(mxa);
388 int nc = mxGetN(mxa);
394 rowind = (
int*)mxGetIr(mxa);
395 colptr = (
int*)mxGetJc(mxa);
399 for (
int i = 0; i < nc; i++) {
400 for (
int j = colptr[i]; j < colptr[i + 1]; j++) {
401 rowCounts[rowind[j]]++;
405 for (
int i = 0; i < nc; i++) {
406 for (
int j = colptr[i]; j < colptr[i + 1]; j++) {
409 complex_t value = std::complex<double>(realArray[j], imagArray[j]);
412 A->insertGlobalValues(rowind[j], cols, vals);
415 A->fillComplete(domainMap, rowMap);
420 }
catch (std::exception& e) {
421 mexPrintf(
"Error while constructing tpetra matrix:\n");
422 std::cout << e.what() << std::endl;
430 return MueLu::TpetraCrs_To_XpetraMatrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tmat);
436 return MueLu::TpetraCrs_To_XpetraMatrix<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tmat);
442 return MueLu::TpetraMultiVector_To_XpetraMultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tpetraMV);
448 return MueLu::TpetraMultiVector_To_XpetraMultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tpetraMV);
451 #ifdef HAVE_MUELU_EPETRA
458 double* vals = mxGetPr(mxa);
459 int nr = mxGetM(mxa);
460 int nc = mxGetN(mxa);
465 rowind = (
int*)mxGetIr(mxa);
466 colptr = (
int*)mxGetJc(mxa);
473 for (
int i = 0; i < nc; i++) {
474 for (
int j = colptr[i]; j < colptr[i + 1]; j++) {
484 }
catch (std::exception& e) {
485 mexPrintf(
"An error occurred while setting up an Epetra matrix:\n");
486 std::cout << e.what() << std::endl;
493 int nr = mxGetM(mxa);
494 int nc = mxGetN(mxa);
503 if (mxGetNumberOfElements(mxa) != 1)
504 throw runtime_error(
"Aggregates must be individual structs in MATLAB.");
505 if (!mxIsStruct(mxa))
506 throw runtime_error(
"Trying to pull aggregates from non-struct MATLAB object.");
509 const int correctNumFields = 5;
510 if (mxGetNumberOfFields(mxa) != correctNumFields)
511 throw runtime_error(
"Aggregates structure has wrong number of fields.");
513 int nVert = *(
int*)mxGetData(mxGetField(mxa, 0,
"nVertices"));
514 int nAgg = *(
int*)mxGetData(mxGetField(mxa, 0,
"nAggregates"));
518 int myRank = comm->getRank();
522 agg->SetNumAggregates(nAgg);
530 mxArray* vertToAggID_in = mxGetField(mxa, 0,
"vertexToAggID");
531 int* vertToAggID_inArray = (
int*)mxGetData(vertToAggID_in);
532 mxArray* rootNodes_in = mxGetField(mxa, 0,
"rootNodes");
533 int* rootNodes_inArray = (
int*)mxGetData(rootNodes_in);
534 for (
int i = 0; i < nVert; i++) {
535 vertex2AggId[i] = vertToAggID_inArray[i];
536 procWinner[i] = myRank;
537 agg->SetIsRoot(i,
false);
539 for (
int i = 0; i < nAgg; i++)
541 agg->SetIsRoot(rootNodes_inArray[i],
true);
544 agg->ComputeAggregateSizes(
true);
545 agg->AggregatesCrossProcessors(
false);
552 throw runtime_error(
"AmalgamationInfo not supported in Muemex yet.");
559 mxArray* edges = mxGetField(mxa, 0,
"edges");
560 mxArray* boundaryNodes = mxGetField(mxa, 0,
"boundaryNodes");
562 throw runtime_error(
"Graph structure in MATLAB must have a field called 'edges' (logical sparse matrix)");
563 if (boundaryNodes == NULL)
564 throw runtime_error(
"Graph structure in MATLAB must have a field called 'boundaryNodes' (int32 array containing list of boundary nodes)");
565 int* boundaryList = (
int*)mxGetData(boundaryNodes);
566 if (!mxIsSparse(edges) || mxGetClassID(edges) != mxLOGICAL_CLASS)
567 throw runtime_error(
"Graph edges must be stored as a logical sparse matrix.");
569 mwIndex* matlabColPtrs = mxGetJc(edges);
570 mwIndex* matlabRowIndices = mxGetIr(edges);
577 int nnz = matlabColPtrs[mxGetN(edges)];
578 for (
int i = 0; i < nnz; i++)
579 entriesPerRow[matlabRowIndices[i]]++;
584 for (
int i = 0; i < nRows; i++)
585 rows[i + 1] = rows[i] + entriesPerRow[i];
588 int ncols = mxGetN(edges);
589 for (
int colNum = 0; colNum < ncols; ++colNum) {
590 int ci = matlabColPtrs[colNum];
591 for (
int j = ci; j < Teuchos::as<int>(matlabColPtrs[colNum + 1]); ++j) {
592 int rowNum = matlabRowIndices[j];
593 cols[rows[rowNum] + insertionsPerRow[rowNum]] = colNum;
594 insertionsPerRow[rowNum]++;
599 for (
int i = 0; i < nRows; i++) {
600 if (maxNzPerRow < entriesPerRow[i])
601 maxNzPerRow = entriesPerRow[i];
610 for (
int i = 0; i < nRows; ++i) {
611 tgraph->insertGlobalIndices((
mm_GlobalOrd)i, cols(rows[i], entriesPerRow[i]));
613 tgraph->fillComplete(map, map);
616 int numBoundaryNodes = mxGetNumberOfElements(boundaryNodes);
617 Kokkos::View<bool*, typename mm_node_t::memory_space> boundaryFlags(
"boundaryFlags", nRows);
619 for (
int i = 0; i < nRows; i++) {
620 boundaryFlags[i] =
false;
622 for (
int i = 0; i < numBoundaryNodes; i++) {
623 boundaryFlags[boundaryList[i]] =
true;
625 mgraph->SetBoundaryNodeMap(boundaryFlags);
629 #ifdef HAVE_MUELU_INTREPID2
632 if (mxGetClassID(mxa) != mxINT32_CLASS)
633 throw runtime_error(
"FieldContainer must have integer storage entries");
635 int* data = (
int*)mxGetData(mxa);
636 int nr = mxGetM(mxa);
637 int nc = mxGetN(mxa);
640 for (
int col = 0; col < nc; col++) {
641 for (
int row = 0; row < nr; row++) {
642 (*fc)(row, col) = data[col * nr + row];
655 mwSize dims[] = {1, 1};
656 mxArray* mxa = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
657 *((
int*)mxGetData(mxa)) = data;
663 mwSize dims[] = {1, 1};
664 mxArray* mxa = mxCreateLogicalArray(2, dims);
665 *((
bool*)mxGetData(mxa)) = data;
671 return mxCreateDoubleScalar(data);
676 mwSize dims[] = {1, 1};
677 mxArray* mxa = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxCOMPLEX);
678 *((
double*)mxGetPr(mxa)) = real<double>(data);
679 *((
double*)mxGetPi(mxa)) = imag<double>(data);
685 return mxCreateString(data.c_str());
691 int nc = data->getGlobalNumElements();
694 double* array = (
double*)malloc(
sizeof(
double) * nr * nc);
695 for (
int col = 0; col < nc; col++) {
697 array[col] = Teuchos::as<double>(gid);
706 mwSize len = data->getGlobalLength();
708 mwSize dimensions[] = {len, 1};
709 mxArray* rv = mxCreateNumericArray(2, dimensions, mxINT32_CLASS, mxREAL);
710 int* dataPtr = (
int*)mxGetData(rv);
712 for (
int i = 0; i < int(data->getGlobalLength()); i++) {
746 Teuchos::rcp_const_cast<
Xpetra_CrsGraph>(data->getCrsGraph())->computeGlobalConstants();
748 int nr = data->getGlobalNumRows();
749 int nc = data->getGlobalNumCols();
750 int nnz = data->getGlobalNumEntries();
752 #ifdef VERBOSE_OUTPUT
756 mxArray* mxa = createMatlabSparse<Scalar>(nr, nc, nnz);
759 for (
int i = 0; i < nc + 1; i++) {
763 size_t maxEntriesPerRow = data->getGlobalMaxNumRowEntries();
766 int* rowProgress =
new int[nc];
768 Scalar* sparseVals =
new Scalar[nnz];
770 if (data->isLocallyIndexed()) {
771 Scalar* rowValArray =
new Scalar[maxEntriesPerRow];
777 data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
778 for (
mm_LocalOrd entry = 0; entry < int(numEntries); entry++)
780 jc[rowIndices[entry] + 1]++;
785 int entriesAccum = 0;
786 for (
int n = 0; n <= nc; n++) {
787 int temp = entriesAccum;
788 entriesAccum += jc[n];
792 for (
int i = 0; i < nc; i++) {
798 data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
803 sparseVals[jc[col] + rowProgress[col]] = rowVals[i];
804 ir[jc[col] + rowProgress[col]] = m;
808 delete[] rowIndicesArray;
813 data->getGlobalRowView(m, rowIndices, rowVals);
815 jc[rowIndices[n] + 1]++;
821 for (
int i = 0; i < nc; i++) {
824 int entriesAccum = 0;
825 for (
int n = 0; n <= nc; n++) {
826 int temp = entriesAccum;
827 entriesAccum += jc[n];
833 data->getGlobalRowView(m, rowIndices, rowVals);
837 sparseVals[jc[col] + rowProgress[col]] = rowVals[i];
838 ir[jc[col] + rowProgress[col]] = m;
844 fillMatlabArray<Scalar>(sparseVals, mxa, nnz);
846 delete[] rowProgress;
855 Teuchos::rcp_const_cast<
Xpetra_CrsGraph>(data->getCrsGraph())->computeGlobalConstants();
857 int nr = data->getGlobalNumRows();
858 int nc = data->getGlobalNumCols();
859 int nnz = data->getGlobalNumEntries();
860 #ifdef VERBOSE_OUTPUT
864 mxArray* mxa = createMatlabSparse<Scalar>(nr, nc, nnz);
867 for (
int i = 0; i < nc + 1; i++) {
870 size_t maxEntriesPerRow = data->getGlobalMaxNumRowEntries();
871 int* rowProgress =
new int[nc];
873 Scalar* sparseVals =
new Scalar[nnz];
875 if (data->isLocallyIndexed()) {
876 Scalar* rowValArray =
new Scalar[maxEntriesPerRow];
882 data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
883 for (
mm_LocalOrd entry = 0; entry < int(numEntries); entry++)
885 jc[rowIndices[entry] + 1]++;
889 int entriesAccum = 0;
890 for (
int n = 0; n <= nc; n++) {
891 int temp = entriesAccum;
892 entriesAccum += jc[n];
896 for (
int i = 0; i < nc; i++) {
902 data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
907 sparseVals[jc[col] + rowProgress[col]] = rowVals[i];
908 ir[jc[col] + rowProgress[col]] = m;
912 delete[] rowIndicesArray;
917 data->getGlobalRowView(m, rowIndices, rowVals);
919 jc[rowIndices[n] + 1]++;
925 for (
int i = 0; i < nc; i++) {
928 int entriesAccum = 0;
929 for (
int n = 0; n <= nc; n++) {
930 int temp = entriesAccum;
931 entriesAccum += jc[n];
937 data->getGlobalRowView(m, rowIndices, rowVals);
941 sparseVals[jc[col] + rowProgress[col]] = rowVals[i];
942 ir[jc[col] + rowProgress[col]] = m;
948 fillMatlabArray<Scalar>(sparseVals, mxa, nnz);
950 delete[] rowProgress;
980 int nr = data->getGlobalLength();
981 int nc = data->getNumVectors();
983 double* array = (
double*)malloc(
sizeof(
double) * nr * nc);
984 for (
int col = 0; col < nc; col++) {
986 for (
int row = 0; row < nr; row++) {
987 array[col * nr + row] = colData[row];
998 int nr = data->getGlobalLength();
999 int nc = data->getNumVectors();
1002 for (
int col = 0; col < nc; col++) {
1004 for (
int row = 0; row < nr; row++) {
1005 array[col * nr + row] = colData[row];
1013 #ifdef HAVE_MUELU_EPETRA
1022 mxArray* output = mxCreateDoubleMatrix(data->GlobalLength(), data->NumVectors(), mxREAL);
1023 double* dataPtr = mxGetPr(output);
1024 data->ExtractCopy(dataPtr, data->GlobalLength());
1032 int numNodes = data->GetVertex2AggId()->getData(0).size();
1033 int numAggs = data->GetNumAggregates();
1035 mwSize singleton[] = {1, 1};
1036 dataIn[0] = mxCreateNumericArray(2, singleton, mxINT32_CLASS, mxREAL);
1037 *((
int*)mxGetData(dataIn[0])) = numNodes;
1038 dataIn[1] = mxCreateNumericArray(2, singleton, mxINT32_CLASS, mxREAL);
1039 *((
int*)mxGetData(dataIn[1])) = numAggs;
1040 mwSize nodeArrayDims[] = {(mwSize)numNodes, 1};
1041 dataIn[2] = mxCreateNumericArray(2, nodeArrayDims, mxINT32_CLASS, mxREAL);
1042 int* vtaid = (
int*)mxGetData(dataIn[2]);
1044 for (
int i = 0; i < numNodes; i++) {
1045 vtaid[i] = vertexToAggID[i];
1047 mwSize aggArrayDims[] = {(mwSize)numAggs, 1};
1048 dataIn[3] = mxCreateNumericArray(2, aggArrayDims, mxINT32_CLASS, mxREAL);
1051 for (
int i = 0; i < numNodes; i++) {
1052 if (data->IsRoot(i))
1055 bool reassignRoots =
false;
1056 if (totalRoots != numAggs) {
1058 <<
"Warning: Number of root nodes and number of aggregates do not match." << endl;
1059 cout <<
"Will reassign root nodes when writing aggregates to matlab." << endl
1061 reassignRoots =
true;
1063 int* rn = (
int*)mxGetData(dataIn[3]);
1065 if (reassignRoots) {
1067 int lastFoundNode = 0;
1068 for (
int i = 0; i < numAggs; i++) {
1070 for (
int j = lastFoundNode; j < lastFoundNode + numNodes; j++) {
1071 int index = j % numNodes;
1072 if (vertexToAggID[index] == i) {
1074 lastFoundNode = index;
1077 TEUCHOS_TEST_FOR_EXCEPTION(rn[i] == -1, runtime_error,
"Invalid aggregates: Couldn't find any node in aggregate #" << i <<
".");
1081 for (
int j = 0; j < numNodes; j++) {
1082 if (data->IsRoot(j)) {
1084 throw runtime_error(
"Cannot store invalid aggregates in MATLAB - more root nodes than aggregates.");
1089 if (i + 1 < numAggs)
1090 throw runtime_error(
"Cannot store invalid aggregates in MATLAB - fewer root nodes than aggregates.");
1093 dataIn[4] = mxCreateNumericArray(1, aggArrayDims, mxINT32_CLASS, mxREAL);
1094 int*
as = (
int*)mxGetData(dataIn[4]);
1095 auto aggSizes = data->ComputeAggregateSizes();
1096 for (
int i = 0; i < numAggs; i++) {
1097 as[i] = aggSizes[i];
1100 int result = mexCallMATLAB(1, matlabAggs, 5, dataIn,
"constructAggregates");
1102 throw runtime_error(
"Matlab encountered an error while constructing aggregates struct.");
1103 return matlabAggs[0];
1108 throw runtime_error(
"AmalgamationInfo not supported in MueMex yet.");
1109 return mxCreateDoubleScalar(0);
1114 int numEntries = (int)data->GetGlobalNumEdges();
1115 int numRows = (int)data->GetDomainMap()->getGlobalNumElements();
1116 mxArray* mat = mxCreateSparseLogicalMatrix(numRows, numRows, numEntries);
1117 mxLogical* outData = (mxLogical*)mxGetData(mat);
1118 mwIndex* rowInds = mxGetIr(mat);
1119 mwIndex* colPtrs = mxGetJc(mat);
1122 int* entriesPerRow =
new int[numRows];
1123 int* entriesPerCol =
new int[numRows];
1124 for (
int i = 0; i < numRows; i++) {
1125 entriesPerRow[i] = 0;
1126 entriesPerCol[i] = 0;
1128 for (
int i = 0; i < numRows; i++) {
1131 entriesPerRow[i] = neighbors.
size();
1132 for (
int j = 0; j < neighbors.
size(); j++) {
1133 entriesPerCol[neighbors[j]]++;
1135 iter += neighbors.
size();
1138 mxLogical** valuesByColumn =
new mxLogical*[numRows];
1139 int* numEnteredPerCol =
new int[numRows];
1141 for (
int i = 0; i < numRows; i++) {
1142 rowIndsByColumn[i] = &rowInds[accum];
1144 valuesByColumn[i] = &outData[accum];
1145 accum += entriesPerCol[i];
1146 if (accum > numEntries)
1147 throw runtime_error(
"potato");
1149 for (
int i = 0; i < numRows; i++) {
1150 numEnteredPerCol[i] = 0;
1154 for (
int row = 0; row < numRows; row++) {
1155 for (
int entryInRow = 0; entryInRow < entriesPerRow[row]; entryInRow++) {
1156 int col = dataCopy[accum];
1158 rowIndsByColumn[col][numEnteredPerCol[col]] = row;
1159 valuesByColumn[col][numEnteredPerCol[col]] = (mxLogical)1;
1160 numEnteredPerCol[col]++;
1164 for (
int col = 0; col < numRows; col++) {
1165 colPtrs[col] = accum;
1166 accum += entriesPerCol[col];
1168 colPtrs[numRows] = accum;
1169 delete[] numEnteredPerCol;
1170 delete[] rowIndsByColumn;
1171 delete[] valuesByColumn;
1173 delete[] entriesPerRow;
1174 delete[] entriesPerCol;
1176 auto boundaryFlags = data->GetBoundaryNodeMap();
1177 int numBoundaryNodes = 0;
1178 for (
int i = 0; i < (int)boundaryFlags.size(); i++) {
1179 if (boundaryFlags[i])
1182 cout <<
"Graph has " << numBoundaryNodes <<
" Dirichlet boundary nodes." << endl;
1183 mwSize dims[] = {(mwSize)numBoundaryNodes, 1};
1184 mxArray* boundaryList = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
1185 int* dest = (
int*)mxGetData(boundaryList);
1186 int* destIter = dest;
1187 for (
int i = 0; i < (int)boundaryFlags.size(); i++) {
1188 if (boundaryFlags[i]) {
1193 mxArray* constructArgs[] = {mat, boundaryList};
1195 mexCallMATLAB(1, out, 2, constructArgs,
"constructGraph");
1199 #ifdef HAVE_MUELU_INTREPID2
1202 int rank = data->rank();
1205 throw std::runtime_error(
"Error: Only rank two FieldContainers are supported.");
1207 int nr = data->extent(0);
1208 int nc = data->extent(1);
1210 mwSize dims[] = {(mwSize)nr, (mwSize)nc};
1211 mxArray* mxa = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
1212 int* array = (
int*)mxGetData(mxa);
1214 for (
int col = 0; col < nc; col++) {
1215 for (
int row = 0; row < nr; row++) {
1216 array[col * nr + row] = (*data)(row, col);
1223 template <
typename T>
1226 data = loadDataFromMatlab<T>(mxa);
1229 template <
typename T>
1231 return saveDataToMatlab<T>(data);
1234 template <
typename T>
1240 template <
typename T>
1244 template <
typename T>
1249 template <
typename T>
1251 this->data = newData;
1258 template <
typename T>
1261 lvl.
Set<T>(name, data, fact);
1264 template <
typename T>
1267 return lvl.
Get<T>(name);
1268 }
catch (std::exception& e) {
1269 throw std::runtime_error(
"Requested custom variable " + name +
" is not in the level.");
1274 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
1276 using namespace std;
1277 using namespace Teuchos;
1284 vector<RCP<MuemexArg> > args;
1285 for (
size_t i = 0; i < needsList.size(); i++) {
1286 if (needsList[i] ==
"A" || needsList[i] ==
"P" || needsList[i] ==
"R" || needsList[i] ==
"Ptent") {
1287 Matrix_t mydata = lvl.
Get<Matrix_t>(needsList[i], factory->
GetFactory(needsList[i]).get());
1289 }
else if (needsList[i] ==
"Nullspace" || needsList[i] ==
"Coordinates") {
1290 MultiVector_t mydata = lvl.
Get<MultiVector_t>(needsList[i], factory->
GetFactory(needsList[i]).get());
1292 }
else if (needsList[i] ==
"Aggregates") {
1293 Aggregates_t mydata = lvl.
Get<Aggregates_t>(needsList[i], factory->
GetFactory(needsList[i]).get());
1295 }
else if (needsList[i] ==
"UnAmalgamationInfo") {
1296 AmalgamationInfo_t mydata = lvl.
Get<AmalgamationInfo_t>(needsList[i], factory->
GetFactory(needsList[i]).get());
1298 }
else if (needsList[i] ==
"Level") {
1301 }
else if (needsList[i] ==
"Graph") {
1302 Graph_t mydata = lvl.
Get<Graph_t>(needsList[i], factory->
GetFactory(needsList[i]).get());
1305 vector<string> words;
1306 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";
1308 char* buf = (
char*)malloc(needsList[i].size() + 1);
1309 strcpy(buf, needsList[i].c_str());
1310 for (
char* iter = buf; *iter !=
' '; iter++) {
1313 throw runtime_error(badNameMsg);
1317 const char* wordDelim =
" ";
1318 char* mark = strtok(buf, wordDelim);
1319 while (mark != NULL) {
1320 string wordStr(mark);
1321 words.push_back(wordStr);
1322 mark = strtok(NULL, wordDelim);
1324 if (words.size() != 2) {
1326 throw runtime_error(badNameMsg);
1328 char* typeStr = (
char*)words[0].c_str();
1329 if (strstr(typeStr,
"ordinalvector")) {
1331 LOVector_t mydata = getLevelVariable<LOVector_t>(needsList[i], lvl);
1333 }
else if (strstr(typeStr,
"map")) {
1335 Map_t mydata = getLevelVariable<Map_t>(needsList[i], lvl);
1337 }
else if (strstr(typeStr,
"scalar")) {
1338 Scalar mydata = getLevelVariable<Scalar>(needsList[i], lvl);
1340 }
else if (strstr(typeStr,
"double")) {
1341 double mydata = getLevelVariable<double>(needsList[i], lvl);
1343 }
else if (strstr(typeStr,
"complex")) {
1344 complex_t mydata = getLevelVariable<complex_t>(needsList[i], lvl);
1346 }
else if (strstr(typeStr,
"matrix")) {
1347 Matrix_t mydata = getLevelVariable<Matrix_t>(needsList[i], lvl);
1349 }
else if (strstr(typeStr,
"multivector")) {
1350 MultiVector_t mydata = getLevelVariable<MultiVector_t>(needsList[i], lvl);
1352 }
else if (strstr(typeStr,
"int")) {
1353 int mydata = getLevelVariable<int>(needsList[i], lvl);
1355 }
else if (strstr(typeStr,
"string")) {
1356 string mydata = getLevelVariable<string>(needsList[i], lvl);
1360 throw std::runtime_error(words[0] +
" is not a known variable type.");
1368 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
1370 using namespace std;
1371 using namespace Teuchos;
1378 for (
size_t i = 0; i < size_t(provides.size()); i++) {
1379 if (provides[i] ==
"A" || provides[i] ==
"P" || provides[i] ==
"R" || provides[i] ==
"Ptent") {
1381 lvl.
Set(provides[i], mydata->getData(), factory);
1382 }
else if (provides[i] ==
"Nullspace" || provides[i] ==
"Coordinates") {
1384 lvl.
Set(provides[i], mydata->getData(), factory);
1385 }
else if (provides[i] ==
"Aggregates") {
1387 lvl.
Set(provides[i], mydata->getData(), factory);
1388 }
else if (provides[i] ==
"UnAmalgamationInfo") {
1390 lvl.
Set(provides[i], mydata->getData(), factory);
1391 }
else if (provides[i] ==
"Graph") {
1393 lvl.
Set(provides[i], mydata->getData(), factory);
1395 vector<string> words;
1396 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";
1398 char* buf = (
char*)malloc(provides[i].size() + 1);
1399 strcpy(buf, provides[i].c_str());
1400 for (
char* iter = buf; *iter !=
' '; iter++) {
1403 throw runtime_error(badNameMsg);
1407 const char* wordDelim =
" ";
1408 char* mark = strtok(buf, wordDelim);
1409 while (mark != NULL) {
1410 string wordStr(mark);
1411 words.push_back(wordStr);
1412 mark = strtok(NULL, wordDelim);
1414 if (words.size() != 2) {
1416 throw runtime_error(badNameMsg);
1418 const char* typeStr = words[0].c_str();
1419 if (strstr(typeStr,
"ordinalvector")) {
1422 addLevelVariable<LOVector_t>(mydata->getData(), words[1], lvl, factory);
1423 }
else if (strstr(typeStr,
"map")) {
1426 addLevelVariable<Map_t>(mydata->getData(), words[1], lvl, factory);
1427 }
else if (strstr(typeStr,
"scalar")) {
1429 addLevelVariable<Scalar>(mydata->getData(), words[1], lvl, factory);
1430 }
else if (strstr(typeStr,
"double")) {
1432 addLevelVariable<double>(mydata->getData(), words[1], lvl, factory);
1433 }
else if (strstr(typeStr,
"complex")) {
1435 addLevelVariable<complex_t>(mydata->getData(), words[1], lvl, factory);
1436 }
else if (strstr(typeStr,
"matrix")) {
1438 addLevelVariable<Matrix_t>(mydata->getData(), words[1], lvl, factory);
1439 }
else if (strstr(typeStr,
"multivector")) {
1441 addLevelVariable<MultiVector_t>(mydata->getData(), words[1], lvl, factory);
1442 }
else if (strstr(typeStr,
"int")) {
1444 addLevelVariable<int>(mydata->getData(), words[1], lvl, factory);
1445 }
else if (strstr(typeStr,
"bool")) {
1447 addLevelVariable<bool>(mydata->getData(), words[1], lvl, factory);
1448 }
else if (strstr(typeStr,
"string")) {
1450 addLevelVariable<string>(mydata->getData(), words[1], lvl, factory);
1453 throw std::runtime_error(words[0] +
" is not a known variable type.");
1464 throw std::runtime_error(
"Muemex does not support long long for global indices");
1469 throw std::runtime_error(
"Muemex does not support long long for global indices");
1474 throw std::runtime_error(
"Muemex does not support long long for global indices");
1479 throw std::runtime_error(
"Muemex does not support long long for global indices");
1483 #endif // HAVE_MUELU_MATLAB error handler
1484 #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)
struct mxArray_tag mxArray
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
Lightweight MueLu representation of a compressed row storage graph.
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)
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Default implementation of FactoryAcceptor::GetFactory()