47 #ifndef MUELU_MATLABUTILS_DEF_HPP 48 #define MUELU_MATLABUTILS_DEF_HPP 52 #if !defined(HAVE_MUELU_MATLAB) || !defined(HAVE_MUELU_EPETRA) || !defined(HAVE_MUELU_TPETRA) 53 #error "Muemex types require MATLAB, Epetra and Tpetra." 148 mxClassID probIDtype = mxGetClassID(mxa);
150 if(probIDtype == mxINT32_CLASS)
152 rv = *((
int*) mxGetData(mxa));
154 else if(probIDtype == mxLOGICAL_CLASS)
156 rv = (int) *((
bool*) mxGetData(mxa));
158 else if(probIDtype == mxDOUBLE_CLASS)
160 rv = (int) *((
double*) mxGetData(mxa));
162 else if(probIDtype == mxUINT32_CLASS)
164 rv = (int) *((
unsigned int*) mxGetData(mxa));
169 throw std::runtime_error(
"Error: Unrecognized numerical type.");
177 return *((
bool*) mxGetData(mxa));
183 return *((
double*) mxGetPr(mxa));
189 double realpart = real<double>(*((
double*) mxGetPr(mxa)));
190 double imagpart = imag<double>(*((
double*) mxGetPi(mxa)));
198 if(!mxGetClassID(mxa) != mxCHAR_CLASS)
200 throw runtime_error(
"Can't construct string from anything but a char array.");
202 rv = string(mxArrayToString(mxa));
210 int nr = mxGetM(mxa);
211 int nc = mxGetN(mxa);
213 throw std::runtime_error(
"A Xpetra::Map representation from MATLAB must be a single row vector.");
214 double* pr = mxGetPr(mxa);
217 std::vector<mm_GlobalOrd> localGIDs(numGlobalIndices);
218 for(
int i = 0; i < int(numGlobalIndices); i++) {
219 localGIDs[i] = Teuchos::as<mm_GlobalOrd>(pr[i]);
231 throw runtime_error(
"Failed to create Xpetra::Map.");
239 if(mxGetN(mxa) != 1 && mxGetM(mxa) != 1)
240 throw std::runtime_error(
"An OrdinalVector from MATLAB must be a single row or column vector.");
241 mm_GlobalOrd numGlobalIndices = mxGetM(mxa) * mxGetN(mxa);
243 if(mxGetClassID(mxa) != mxINT32_CLASS)
244 throw std::runtime_error(
"Can only construct LOVector with int32 data.");
245 int* array = (
int*) mxGetData(mxa);
247 throw runtime_error(
"Failed to create map for Xpetra ordinal vector.");
250 throw runtime_error(
"Failed to create ordinal vector with Xpetra::VectorFactory.");
251 for(
int i = 0; i < int(numGlobalIndices); i++)
253 loVec->replaceGlobalValue(i, 0, array[i]);
264 int nr = mxGetM(mxa);
265 int nc = mxGetN(mxa);
266 double* pr = mxGetPr(mxa);
272 mv =
rcp(
new Tpetra::MultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(map, arrView,
size_t(nr),
size_t(nc)));
274 catch(std::exception& e)
276 mexPrintf(
"Error constructing Tpetra MultiVector.\n");
277 std::cout << e.what() << std::endl;
288 int nr = mxGetM(mxa);
289 int nc = mxGetN(mxa);
290 double* pr = mxGetPr(mxa);
291 double* pi = mxGetPi(mxa);
297 for(
int n = 0; n < nc; n++)
299 for(
int m = 0; m < nr; m++)
301 myArr[n * nr + m] =
complex_t(pr[n * nr + m], pi[n * nr + m]);
305 mv =
rcp(
new Tpetra::MultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(map, arrView, nr, nc));
307 catch(std::exception& e)
309 mexPrintf(
"Error constructing Tpetra MultiVector.\n");
310 std::cout << e.what() << std::endl;
318 bool success =
false;
328 const Tpetra::global_size_t numGlobalIndices = mxGetM(mxa);
332 A = Tpetra::createCrsMatrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(rowMap);
333 double* valueArray = mxGetPr(mxa);
334 int nc = mxGetN(mxa);
343 rowind = (
int*) mxGetIr(mxa);
344 colptr = (
int*) mxGetJc(mxa);
346 for(
int i = 0; i < nc; i++)
348 for(
int j = colptr[i]; j < colptr[i + 1]; j++)
354 A->insertGlobalValues(rowind[j], cols, vals);
357 A->fillComplete(domainMap, rowMap);
360 delete[] rowind; rowind = NULL;
361 delete[] colptr; colptr = NULL;
365 catch(std::exception& e)
369 if(rowind!=NULL)
delete[] rowind;
370 if(colptr!=NULL)
delete[] colptr;
374 mexPrintf(
"Error while constructing Tpetra matrix:\n");
375 std::cout << e.what() << std::endl;
378 mexErrMsgTxt(
"An error occurred while setting up a Tpetra matrix.\n");
390 const Tpetra::global_size_t numGlobalIndices = mxGetM(mxa);
394 A = Tpetra::createCrsMatrix<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(rowMap);
395 double* realArray = mxGetPr(mxa);
396 double* imagArray = mxGetPi(mxa);
399 int nc = mxGetN(mxa);
408 rowind = (
int*) mxGetIr(mxa);
409 colptr = (
int*) mxGetJc(mxa);
411 for(
int i = 0; i < nc; i++)
413 for(
int j = colptr[i]; j < colptr[i + 1]; j++)
417 complex_t value = std::complex<double>(realArray[j], imagArray[j]);
420 A->insertGlobalValues(rowind[j], cols, vals);
423 A->fillComplete(domainMap, rowMap);
430 catch(std::exception& e)
432 mexPrintf(
"Error while constructing tpetra matrix:\n");
433 std::cout << e.what() << std::endl;
442 return MueLu::TpetraCrs_To_XpetraMatrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tmat);
449 return MueLu::TpetraCrs_To_XpetraMatrix<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tmat);
456 return MueLu::TpetraMultiVector_To_XpetraMultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tpetraMV);
463 return MueLu::TpetraMultiVector_To_XpetraMultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tpetraMV);
474 double* vals = mxGetPr(mxa);
475 int nr = mxGetM(mxa);
476 int nc = mxGetN(mxa);
484 rowind = (
int*) mxGetIr(mxa);
485 colptr = (
int*) mxGetJc(mxa);
492 for(
int i = 0; i < nc; i++)
494 for(
int j = colptr[i]; j < colptr[i + 1]; j++)
507 catch(std::exception& e)
509 mexPrintf(
"An error occurred while setting up an Epetra matrix:\n");
510 std::cout << e.what() << std::endl;
518 int nr = mxGetM(mxa);
519 int nc = mxGetN(mxa);
528 if(mxGetNumberOfElements(mxa) != 1)
529 throw runtime_error(
"Aggregates must be individual structs in MATLAB.");
531 throw runtime_error(
"Trying to pull aggregates from non-struct MATLAB object.");
534 const int correctNumFields = 5;
535 if(mxGetNumberOfFields(mxa) != correctNumFields)
536 throw runtime_error(
"Aggregates structure has wrong number of fields.");
538 int nVert = *(
int*) mxGetData(mxGetField(mxa, 0,
"nVertices"));
539 int nAgg = *(
int*) mxGetData(mxGetField(mxa, 0,
"nAggregates"));
543 int myRank = comm->getRank();
547 agg->SetNumAggregates(nAgg);
555 mxArray* vertToAggID_in = mxGetField(mxa, 0,
"vertexToAggID");
556 int* vertToAggID_inArray = (
int*) mxGetData(vertToAggID_in);
557 mxArray* rootNodes_in = mxGetField(mxa, 0,
"rootNodes");
558 int* rootNodes_inArray = (
int*) mxGetData(rootNodes_in);
559 for(
int i = 0; i < nVert; i++)
561 vertex2AggId[i] = vertToAggID_inArray[i];
562 procWinner[i] = myRank;
563 agg->SetIsRoot(i,
false);
565 for(
int i = 0; i < nAgg; i++)
567 agg->SetIsRoot(rootNodes_inArray[i],
true);
570 agg->ComputeAggregateSizes(
true,
true);
571 agg->AggregatesCrossProcessors(
false);
579 throw runtime_error(
"AmalgamationInfo not supported in Muemex yet.");
587 mxArray* edges = mxGetField(mxa, 0,
"edges");
588 mxArray* boundaryNodes = mxGetField(mxa, 0,
"boundaryNodes");
590 throw runtime_error(
"Graph structure in MATLAB must have a field called 'edges' (logical sparse matrix)");
591 if(boundaryNodes == NULL)
592 throw runtime_error(
"Graph structure in MATLAB must have a field called 'boundaryNodes' (int32 array containing list of boundary nodes)");
593 int* boundaryList = (
int*) mxGetData(boundaryNodes);
594 if(!mxIsSparse(edges) || mxGetClassID(edges) != mxLOGICAL_CLASS)
595 throw runtime_error(
"Graph edges must be stored as a logical sparse matrix.");
597 mwIndex* matlabColPtrs = mxGetJc(edges);
598 mwIndex* matlabRowIndices = mxGetIr(edges);
605 int nnz = matlabColPtrs[mxGetN(edges)];
606 for(
int i = 0; i < nnz; i++)
607 entriesPerRow[matlabRowIndices[i]]++;
612 for(
int i = 0; i < nRows; i++)
613 rows[i+1] = rows[i] + entriesPerRow[i];
616 int ncols = mxGetN(edges);
617 for (
int colNum=0; colNum<ncols; ++colNum) {
618 int ci = matlabColPtrs[colNum];
619 for (
int j=ci; j<Teuchos::as<int>(matlabColPtrs[colNum+1]); ++j) {
620 int rowNum = matlabRowIndices[j];
621 cols[ rows[rowNum] + insertionsPerRow[rowNum] ] = colNum;
622 insertionsPerRow[rowNum]++;
627 for(
int i = 0; i < nRows; i++) {
628 if(maxNzPerRow < entriesPerRow[i])
629 maxNzPerRow = entriesPerRow[i];
638 for(
int i = 0; i < nRows; ++i) {
639 tgraph->insertGlobalIndices((
mm_GlobalOrd) i, cols(rows[i],entriesPerRow[i]));
641 tgraph->fillComplete(map, map);
644 int numBoundaryNodes = mxGetNumberOfElements(boundaryNodes);
645 bool* boundaryFlags =
new bool[nRows];
646 for(
int i = 0; i < nRows; i++)
648 boundaryFlags[i] =
false;
650 for(
int i = 0; i < numBoundaryNodes; i++)
652 boundaryFlags[boundaryList[i]] =
true;
655 mgraph->SetBoundaryNodeMap(boundaryNodesInput);
666 mwSize dims[] = {1, 1};
667 mxArray* mxa = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
668 *((
int*) mxGetData(mxa)) = data;
675 mwSize dims[] = {1, 1};
676 mxArray* mxa = mxCreateLogicalArray(2, dims);
677 *((
bool*) mxGetData(mxa)) = data;
684 return mxCreateDoubleScalar(data);
690 mwSize dims[] = {1, 1};
691 mxArray* mxa = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxCOMPLEX);
692 *((
double*) mxGetPr(mxa)) = real<double>(data);
693 *((
double*) mxGetPi(mxa)) = imag<double>(data);
700 return mxCreateString(data.c_str());
707 int nc = data->getGlobalNumElements();
710 double* array = (
double*) malloc(
sizeof(
double) * nr * nc);
711 for(
int col = 0; col < nc; col++)
714 array[col] = Teuchos::as<double>(gid);
724 mwSize len = data->getGlobalLength();
726 mwSize dimensions[] = {len, 1};
727 mxArray* rv = mxCreateNumericArray(2, dimensions, mxINT32_CLASS, mxREAL);
728 int* dataPtr = (
int*) mxGetData(rv);
730 for(
int i = 0; i < int(data->getGlobalLength()); i++)
745 mxArray*
saveDataToMatlab(
RCP<Tpetra::MultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>>& data)
768 typedef double Scalar;
769 int nr = data->getGlobalNumRows();
770 int nc = data->getGlobalNumCols();
771 int nnz = data->getGlobalNumEntries();
772 #ifdef VERBOSE_OUTPUT 776 mxArray* mxa = createMatlabSparse<Scalar>(nr, nc, nnz);
779 for(
int i = 0; i < nc + 1; i++)
783 size_t maxEntriesPerRow = data->getGlobalMaxNumRowEntries();
784 int* rowProgress =
new int[nc];
786 Scalar* sparseVals =
new Scalar[nnz];
788 if(data->isLocallyIndexed())
790 Scalar* rowValArray =
new Scalar[maxEntriesPerRow];
796 data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
797 for(
mm_LocalOrd entry = 0; entry < int(numEntries); entry++)
799 jc[rowIndices[entry] + 1]++;
803 int entriesAccum = 0;
804 for(
int n = 0; n <= nc; n++)
806 int temp = entriesAccum;
807 entriesAccum += jc[n];
811 for(
int i = 0; i < nc; i++)
818 data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
823 sparseVals[jc[col] + rowProgress[col]] = rowVals[i];
824 ir[jc[col] + rowProgress[col]] = m;
828 delete[] rowIndicesArray;
836 data->getGlobalRowView(m, rowIndices, rowVals);
839 jc[rowIndices[n] + 1]++;
845 for(
int i = 0; i < nc; i++)
849 int entriesAccum = 0;
850 for(
int n = 0; n <= nc; n++)
852 int temp = entriesAccum;
853 entriesAccum += jc[n];
859 data->getGlobalRowView(m, rowIndices, rowVals);
863 sparseVals[jc[col] + rowProgress[col]] = rowVals[i];
864 ir[jc[col] + rowProgress[col]] = m;
870 fillMatlabArray<Scalar>(sparseVals, mxa, nnz);
872 delete[] rowProgress;
880 int nr = data->getGlobalNumRows();
881 int nc = data->getGlobalNumCols();
882 int nnz = data->getGlobalNumEntries();
883 #ifdef VERBOSE_OUTPUT 887 mxArray* mxa = createMatlabSparse<Scalar>(nr, nc, nnz);
890 for(
int i = 0; i < nc + 1; i++)
894 size_t maxEntriesPerRow = data->getGlobalMaxNumRowEntries();
895 int* rowProgress =
new int[nc];
897 Scalar* sparseVals =
new Scalar[nnz];
899 if(data->isLocallyIndexed())
901 Scalar* rowValArray =
new Scalar[maxEntriesPerRow];
907 data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
908 for(
mm_LocalOrd entry = 0; entry < int(numEntries); entry++)
910 jc[rowIndices[entry] + 1]++;
914 int entriesAccum = 0;
915 for(
int n = 0; n <= nc; n++)
917 int temp = entriesAccum;
918 entriesAccum += jc[n];
922 for(
int i = 0; i < nc; i++)
929 data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
934 sparseVals[jc[col] + rowProgress[col]] = rowVals[i];
935 ir[jc[col] + rowProgress[col]] = m;
939 delete[] rowIndicesArray;
947 data->getGlobalRowView(m, rowIndices, rowVals);
950 jc[rowIndices[n] + 1]++;
956 for(
int i = 0; i < nc; i++)
960 int entriesAccum = 0;
961 for(
int n = 0; n <= nc; n++)
963 int temp = entriesAccum;
964 entriesAccum += jc[n];
970 data->getGlobalRowView(m, rowIndices, rowVals);
974 sparseVals[jc[col] + rowProgress[col]] = rowVals[i];
975 ir[jc[col] + rowProgress[col]] = m;
981 fillMatlabArray<Scalar>(sparseVals, mxa, nnz);
983 delete[] rowProgress;
1014 int nr = data->getGlobalLength();
1015 int nc = data->getNumVectors();
1017 double* array = (
double*) malloc(
sizeof(
double) * nr * nc);
1018 for(
int col = 0; col < nc; col++)
1021 for(
int row = 0; row < nr; row++)
1023 array[col * nr + row] = colData[row];
1035 int nr = data->getGlobalLength();
1036 int nc = data->getNumVectors();
1039 for(
int col = 0; col < nc; col++)
1042 for(
int row = 0; row < nr; row++)
1044 array[col * nr + row] = colData[row];
1063 double* dataPtr = mxGetPr(output);
1072 int numNodes = data->GetVertex2AggId()->getData(0).size();
1073 int numAggs = data->GetNumAggregates();
1075 mwSize singleton[] = {1, 1};
1076 dataIn[0] = mxCreateNumericArray(2, singleton, mxINT32_CLASS, mxREAL);
1077 *((
int*) mxGetData(dataIn[0])) = numNodes;
1078 dataIn[1] = mxCreateNumericArray(2, singleton, mxINT32_CLASS, mxREAL);
1079 *((
int*) mxGetData(dataIn[1])) = numAggs;
1080 mwSize nodeArrayDims[] = {(mwSize) numNodes, 1};
1081 dataIn[2] = mxCreateNumericArray(2, nodeArrayDims, mxINT32_CLASS, mxREAL);
1082 int* vtaid = (
int*) mxGetData(dataIn[2]);
1084 for(
int i = 0; i < numNodes; i++)
1086 vtaid[i] = vertexToAggID[i];
1088 mwSize aggArrayDims[] = {(mwSize) numAggs, 1};
1089 dataIn[3] = mxCreateNumericArray(2, aggArrayDims, mxINT32_CLASS, mxREAL);
1092 for(
int i = 0; i < numNodes; i++)
1097 bool reassignRoots =
false;
1098 if(totalRoots != numAggs)
1100 cout << endl <<
"Warning: Number of root nodes and number of aggregates do not match." << endl;
1101 cout <<
"Will reassign root nodes when writing aggregates to matlab." << endl << endl;
1102 reassignRoots =
true;
1104 int* rn = (
int*) mxGetData(dataIn[3]);
1109 int lastFoundNode = 0;
1110 for(
int i = 0; i < numAggs; i++)
1113 for(
int j = lastFoundNode; j < lastFoundNode + numNodes; j++)
1115 int index = j % numNodes;
1116 if(vertexToAggID[index] == i)
1119 lastFoundNode = index;
1122 TEUCHOS_TEST_FOR_EXCEPTION(rn[i] == -1, runtime_error,
"Invalid aggregates: Couldn't find any node in aggregate #" << i <<
".");
1128 for(
int j = 0; j < numNodes; j++)
1133 throw runtime_error(
"Cannot store invalid aggregates in MATLAB - more root nodes than aggregates.");
1139 throw runtime_error(
"Cannot store invalid aggregates in MATLAB - fewer root nodes than aggregates.");
1142 dataIn[4] = mxCreateNumericArray(1, aggArrayDims, mxINT32_CLASS, mxREAL);
1143 int*
as = (
int*) mxGetData(dataIn[4]);
1145 for(
int i = 0; i < numAggs; i++)
1147 as[i] = aggSizes[i];
1149 mxArray* matlabAggs[1];
1150 int result = mexCallMATLAB(1, matlabAggs, 5, dataIn,
"constructAggregates");
1152 throw runtime_error(
"Matlab encountered an error while constructing aggregates struct.");
1153 return matlabAggs[0];
1159 throw runtime_error(
"AmalgamationInfo not supported in MueMex yet.");
1160 return mxCreateDoubleScalar(0);
1166 int numEntries = (int) data->GetGlobalNumEdges();
1167 int numRows = (int) data->GetDomainMap()->getGlobalNumElements();
1168 mxArray* mat = mxCreateSparseLogicalMatrix(numRows, numRows, numEntries);
1169 mxLogical* outData = (mxLogical*) mxGetData(mat);
1170 mwIndex* rowInds = mxGetIr(mat);
1171 mwIndex* colPtrs = mxGetJc(mat);
1174 int* entriesPerRow =
new int[numRows];
1175 int* entriesPerCol =
new int[numRows];
1176 for(
int i = 0; i < numRows; i++)
1178 entriesPerRow[i] = 0;
1179 entriesPerCol[i] = 0;
1181 for(
int i = 0; i < numRows; i++)
1185 entriesPerRow[i] = neighbors.
size();
1186 for(
int j = 0; j < neighbors.
size(); j++)
1188 entriesPerCol[neighbors[j]]++;
1190 iter += neighbors.
size();
1193 mxLogical** valuesByColumn =
new mxLogical*[numRows];
1194 int* numEnteredPerCol =
new int[numRows];
1196 for(
int i = 0; i < numRows; i++)
1198 rowIndsByColumn[i] = &rowInds[accum];
1200 valuesByColumn[i] = &outData[accum];
1201 accum += entriesPerCol[i];
1202 if(accum > numEntries)
1203 throw runtime_error(
"potato");
1205 for(
int i = 0; i < numRows; i++)
1207 numEnteredPerCol[i] = 0;
1211 for(
int row = 0; row < numRows; row++)
1213 for(
int entryInRow = 0; entryInRow < entriesPerRow[row]; entryInRow++)
1215 int col = dataCopy[accum];
1217 rowIndsByColumn[col][numEnteredPerCol[col]] = row;
1218 valuesByColumn[col][numEnteredPerCol[col]] = (mxLogical) 1;
1219 numEnteredPerCol[col]++;
1223 for(
int col = 0; col < numRows; col++)
1225 colPtrs[col] = accum;
1226 accum += entriesPerCol[col];
1228 colPtrs[numRows] = accum;
1229 delete[] numEnteredPerCol;
1230 delete[] rowIndsByColumn;
1231 delete[] valuesByColumn;
1233 delete[] entriesPerRow;
1234 delete[] entriesPerCol;
1237 int numBoundaryNodes = 0;
1238 for(
int i = 0; i < boundaryFlags.
size(); i++)
1240 if(boundaryFlags[i])
1243 cout <<
"Graph has " << numBoundaryNodes <<
" Dirichlet boundary nodes." << endl;
1244 mwSize dims[] = {(mwSize) numBoundaryNodes, 1};
1245 mxArray* boundaryList = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
1246 int* dest = (
int*) mxGetData(boundaryList);
1247 int* destIter = dest;
1248 for(
int i = 0; i < boundaryFlags.
size(); i++)
1250 if(boundaryFlags[i])
1256 mxArray* constructArgs[] = {mat, boundaryList};
1258 mexCallMATLAB(1, out, 2, constructArgs,
"constructGraph");
1262 template<
typename T>
1265 data = loadDataFromMatlab<T>(mxa);
1268 template<
typename T>
1271 return saveDataToMatlab<T>(data);
1274 template<
typename T>
1280 template<
typename T>
1283 template<
typename T>
1289 template<
typename T>
1299 template<
typename T>
1303 lvl.
Set<T>(name, data, fact);
1306 template<
typename T>
1311 return lvl.
Get<T>(name);
1313 catch(std::exception& e)
1315 throw std::runtime_error(
"Requested custom variable " + name +
" is not in the level.");
1320 template<
typename Scalar =
double,
typename LocalOrdinal = mm_LocalOrd,
typename GlobalOrdinal = mm_GlobalOrd,
typename Node = mm_node_t>
1323 using namespace std;
1331 vector<RCP<MuemexArg>> args;
1332 for(
size_t i = 0; i < needsList.size(); i++)
1334 if(needsList[i] ==
"A" || needsList[i] ==
"P" || needsList[i] ==
"R" || needsList[i]==
"Ptent")
1336 Matrix_t mydata = lvl.
Get<Matrix_t>(needsList[i], factory->
GetFactory(needsList[i]).get());
1339 else if(needsList[i] ==
"Nullspace" || needsList[i] ==
"Coordinates")
1341 MultiVector_t mydata = lvl.
Get<MultiVector_t>(needsList[i], factory->
GetFactory(needsList[i]).get());
1344 else if(needsList[i] ==
"Aggregates")
1346 Aggregates_t mydata = lvl.
Get<Aggregates_t>(needsList[i], factory->
GetFactory(needsList[i]).get());
1349 else if(needsList[i] ==
"UnAmalgamationInfo")
1351 AmalgamationInfo_t mydata = lvl.
Get<AmalgamationInfo_t>(needsList[i], factory->
GetFactory(needsList[i]).get());
1354 else if(needsList[i] ==
"Level")
1359 else if(needsList[i] ==
"Graph")
1361 Graph_t mydata = lvl.
Get<Graph_t>(needsList[i], factory->
GetFactory(needsList[i]).get());
1366 vector<string> words;
1367 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";
1369 char* buf = (
char*) malloc(needsList[i].size() + 1);
1370 strcpy(buf, needsList[i].c_str());
1371 for(
char* iter = buf; *iter !=
' '; iter++)
1376 throw runtime_error(badNameMsg);
1378 *iter = (char) tolower(*iter);
1380 const char* wordDelim =
" ";
1381 char* mark = strtok(buf, wordDelim);
1384 string wordStr(mark);
1385 words.push_back(wordStr);
1386 mark = strtok(NULL, wordDelim);
1388 if(words.size() != 2)
1391 throw runtime_error(badNameMsg);
1393 char* typeStr = (
char*) words[0].c_str();
1394 if(strstr(typeStr,
"ordinalvector"))
1397 LOVector_t mydata = getLevelVariable<LOVector_t>(needsList[i], lvl);
1400 else if(strstr(typeStr,
"map"))
1403 Map_t mydata = getLevelVariable<Map_t>(needsList[i], lvl);
1406 else if(strstr(typeStr,
"scalar"))
1408 Scalar mydata = getLevelVariable<Scalar>(needsList[i], lvl);
1411 else if(strstr(typeStr,
"double"))
1413 double mydata = getLevelVariable<double>(needsList[i], lvl);
1416 else if(strstr(typeStr,
"complex"))
1418 complex_t mydata = getLevelVariable<complex_t>(needsList[i], lvl);
1421 else if(strstr(typeStr,
"matrix"))
1423 Matrix_t mydata = getLevelVariable<Matrix_t>(needsList[i], lvl);
1426 else if(strstr(typeStr,
"multivector"))
1428 MultiVector_t mydata = getLevelVariable<MultiVector_t>(needsList[i], lvl);
1431 else if(strstr(typeStr,
"int"))
1433 int mydata = getLevelVariable<int>(needsList[i], lvl);
1436 else if(strstr(typeStr,
"string"))
1438 string mydata = getLevelVariable<string>(needsList[i], lvl);
1444 throw std::runtime_error(words[0] +
" is not a known variable type.");
1452 template<
typename Scalar =
double,
typename LocalOrdinal = mm_LocalOrd,
typename GlobalOrdinal = mm_GlobalOrd,
typename Node = mm_node_t>
1455 using namespace std;
1463 for(
size_t i = 0; i < size_t(provides.size()); i++)
1465 if(provides[i] ==
"A" || provides[i] ==
"P" || provides[i] ==
"R" || provides[i]==
"Ptent")
1468 lvl.
Set(provides[i], mydata->getData(), factory);
1470 else if(provides[i] ==
"Nullspace" || provides[i] ==
"Coordinates")
1473 lvl.
Set(provides[i], mydata->getData(), factory);
1475 else if(provides[i] ==
"Aggregates")
1478 lvl.
Set(provides[i], mydata->getData(), factory);
1480 else if(provides[i] ==
"UnAmalgamationInfo")
1483 lvl.
Set(provides[i], mydata->getData(), factory);
1485 else if(provides[i] ==
"Graph")
1488 lvl.
Set(provides[i], mydata->getData(), factory);
1492 vector<string> words;
1493 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";
1495 char* buf = (
char*) malloc(provides[i].size() + 1);
1496 strcpy(buf, provides[i].c_str());
1497 for(
char* iter = buf; *iter !=
' '; iter++)
1502 throw runtime_error(badNameMsg);
1504 *iter = (char) tolower(*iter);
1506 const char* wordDelim =
" ";
1507 char* mark = strtok(buf, wordDelim);
1510 string wordStr(mark);
1511 words.push_back(wordStr);
1512 mark = strtok(NULL, wordDelim);
1514 if(words.size() != 2)
1517 throw runtime_error(badNameMsg);
1519 const char* typeStr = words[0].c_str();
1520 if(strstr(typeStr,
"ordinalvector"))
1524 addLevelVariable<LOVector_t>(mydata->getData(), words[1], lvl, factory);
1526 else if(strstr(typeStr,
"map"))
1530 addLevelVariable<Map_t>(mydata->getData(), words[1], lvl, factory);
1532 else if(strstr(typeStr,
"scalar"))
1535 addLevelVariable<Scalar>(mydata->getData(), words[1], lvl, factory);
1537 else if(strstr(typeStr,
"double"))
1540 addLevelVariable<double>(mydata->getData(), words[1], lvl, factory);
1542 else if(strstr(typeStr,
"complex"))
1545 addLevelVariable<complex_t>(mydata->getData(), words[1], lvl, factory);
1547 else if(strstr(typeStr,
"matrix"))
1550 addLevelVariable<Matrix_t>(mydata->getData(), words[1], lvl, factory);
1552 else if(strstr(typeStr,
"multivector"))
1555 addLevelVariable<MultiVector_t>(mydata->getData(), words[1], lvl, factory);
1557 else if(strstr(typeStr,
"int"))
1560 addLevelVariable<int>(mydata->getData(), words[1], lvl, factory);
1562 else if(strstr(typeStr,
"bool"))
1565 addLevelVariable<bool>(mydata->getData(), words[1], lvl, factory);
1567 else if(strstr(typeStr,
"string"))
1570 addLevelVariable<string>(mydata->getData(), words[1], lvl, factory);
1575 throw std::runtime_error(words[0] +
" is not a known variable type.");
1586 throw std::runtime_error(
"Muemex does not support long long for global indices");
1591 throw std::runtime_error(
"Muemex does not support long long for global indices");
1596 throw std::runtime_error(
"Muemex does not support long long for global indices");
1601 throw std::runtime_error(
"Muemex does not support long long for global indices");
1606 #endif //HAVE_MUELU_MATLAB error handler 1607 #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)
std::vector< std::string > tokenizeList(const std::string ¶ms)
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, const Teuchos::RCP< Node > &node=defaultArgNode())
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)
mxArray * saveDataToMatlab(RCP< MGraph > &data)
MuemexType getMuemexType< string >()
MueLu representation of a compressed row storage graph.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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)
User data are always kept. This flag is set automatically when Level::Set("data", data) is used...
Namespace for MueLu classes and methods.
void fillMatlabArray< double >(double *array, const mxArray *mxa, int n)
int ExtractCopy(double *A, int MyLDA) const
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)
int FillComplete(bool OptimizeDataStorage=true)
int GetLevelID() const
Return level number.
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)
void processProvides< double, mm_LocalOrd, long long, mm_node_t >(std::vector< Teuchos::RCP< MuemexArg >> &mexOutput, const Factory *factory, std::string &providesParam, Level &lvl)
mxArray * createMatlabSparse< double >(int numRows, int numCols, int nnz)
int * mwIndex_to_int(int N, mwIndex *mwi_array)
Class that holds all level-specific information.
template complex_t loadDataFromMatlab< complex_t >(const mxArray *mxa)
void addLevelVariable(const T &data, std::string &name, Level &lvl, Factory *fact=NoFactory::get())
template int loadDataFromMatlab< int >(const mxArray *mxa)
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)
MuemexType getMuemexType(const RCP< MGraph > &data)
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)
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > TpetraMultiVector_To_XpetraMultiVector(const Teuchos::RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &Vtpetra)
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
void processProvides(std::vector< Teuchos::RCP< MuemexArg >> &mexOutput, const Factory *factory, std::string &providesParam, Level &lvl)
RCP< Xpetra::Matrix< SC, LO, GO, NO > > TpetraCrs_To_XpetraMatrix(const Teuchos::RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > &Atpetra)
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Default implementation of FactoryAcceptor::GetFactory()
MuemexType getMuemexType< complex_t >()