MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_MatlabUtils_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_MATLABUTILS_DEF_HPP
11 #define MUELU_MATLABUTILS_DEF_HPP
12 
14 #include <mex.h>
15 
16 #if !defined(HAVE_MUELU_MATLAB) || !defined(HAVE_MUELU_TPETRA)
17 #error "Muemex types require MATLAB and Tpetra."
18 #else
19 
20 using Teuchos::RCP;
21 using Teuchos::rcp;
22 using namespace std;
23 
24 namespace MueLu {
25 
26 extern bool rewrap_ints;
27 
28 /* ******************************* */
29 /* getMuemexType */
30 /* ******************************* */
31 
32 template <typename T>
33 MuemexType getMuemexType(const T& data) { throw std::runtime_error("Unknown Type"); }
34 
35 template <>
36 MuemexType getMuemexType(const int& data) { return INT; }
37 template <>
39 template <>
41 
42 template <>
43 MuemexType getMuemexType(const double& data) { return DOUBLE; }
44 template <>
46 
47 template <>
48 MuemexType getMuemexType(const std::string& data) { return STRING; }
49 template <>
51 
52 template <>
53 MuemexType getMuemexType(const complex_t& data) { return COMPLEX; }
54 template <>
56 
57 template <>
59 template <>
60 MuemexType getMuemexType<RCP<Xpetra_map> >() { return XPETRA_MAP; }
61 
62 template <>
64 template <>
65 MuemexType getMuemexType<RCP<Xpetra_ordinal_vector> >() { return XPETRA_ORDINAL_VECTOR; }
66 
67 template <>
69 template <>
70 MuemexType getMuemexType<RCP<Tpetra_MultiVector_double> >() { return TPETRA_MULTIVECTOR_DOUBLE; }
71 
72 template <>
74 template <>
75 MuemexType getMuemexType<RCP<Tpetra_MultiVector_complex> >() { return TPETRA_MULTIVECTOR_COMPLEX; }
76 
77 template <>
79 template <>
80 MuemexType getMuemexType<RCP<Tpetra_CrsMatrix_double> >() { return TPETRA_MATRIX_DOUBLE; }
81 
82 template <>
84 template <>
85 MuemexType getMuemexType<RCP<Tpetra_CrsMatrix_complex> >() { return TPETRA_MATRIX_COMPLEX; }
86 
87 template <>
89 template <>
90 MuemexType getMuemexType<RCP<Xpetra_MultiVector_double> >() { return XPETRA_MULTIVECTOR_DOUBLE; }
91 
92 template <>
94 template <>
95 MuemexType getMuemexType<RCP<Xpetra_MultiVector_complex> >() { return XPETRA_MULTIVECTOR_COMPLEX; }
96 
97 template <>
99 template <>
100 MuemexType getMuemexType<RCP<Xpetra_Matrix_double> >() { return XPETRA_MATRIX_DOUBLE; }
101 
102 template <>
104 template <>
105 MuemexType getMuemexType<RCP<Xpetra_Matrix_complex> >() { return XPETRA_MATRIX_COMPLEX; }
106 
107 #ifdef HAVE_MUELU_EPETRA
108 template <>
110 template <>
111 MuemexType getMuemexType<RCP<Epetra_CrsMatrix> >() { return EPETRA_CRSMATRIX; }
112 
113 template <>
115 template <>
116 MuemexType getMuemexType<RCP<Epetra_MultiVector> >() { return EPETRA_MULTIVECTOR; }
117 #endif
118 
119 template <>
121 template <>
122 MuemexType getMuemexType<RCP<MAggregates> >() { return AGGREGATES; }
123 
124 template <>
126 template <>
127 MuemexType getMuemexType<RCP<MAmalInfo> >() { return AMALGAMATION_INFO; }
128 
129 template <>
130 MuemexType getMuemexType(const RCP<MGraph>& data) { return GRAPH; }
131 template <>
132 MuemexType getMuemexType<RCP<MGraph> >() { return GRAPH; }
133 
134 #ifdef HAVE_MUELU_INTREPID2
135 template <>
136 MuemexType getMuemexType(const RCP<FieldContainer_ordinal>& data) { return FIELDCONTAINER_ORDINAL; }
137 template <>
138 MuemexType getMuemexType<RCP<FieldContainer_ordinal> >() { return FIELDCONTAINER_ORDINAL; }
139 #endif
140 
141 /* "prototypes" for specialized functions used in other specialized functions */
142 
143 template <>
144 mxArray* createMatlabSparse<double>(int numRows, int numCols, int nnz);
145 template <>
146 mxArray* createMatlabSparse<complex_t>(int numRows, int numCols, int nnz);
147 template <>
148 mxArray* createMatlabMultiVector<double>(int numRows, int numCols);
149 template <>
150 mxArray* createMatlabMultiVector<complex_t>(int numRows, int numCols);
151 template <>
152 void fillMatlabArray<double>(double* array, const mxArray* mxa, int n);
153 template <>
154 void fillMatlabArray<complex_t>(complex_t* array, const mxArray* mxa, int n);
155 template <>
157 template <>
159 template <>
161 template <>
163 
164 /* ******************************* */
165 /* loadDataFromMatlab */
166 /* ******************************* */
167 
168 template <>
170  mxClassID probIDtype = mxGetClassID(mxa);
171  int rv;
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));
180  } else {
181  rv = -1;
182  throw std::runtime_error("Error: Unrecognized numerical type.");
183  }
184  return rv;
185 }
186 
187 template <>
189  return *((bool*)mxGetData(mxa));
190 }
191 
192 template <>
194  return *((double*)mxGetPr(mxa));
195 }
196 
197 template <>
199  double realpart = real<double>(*((double*)mxGetPr(mxa)));
200  double imagpart = imag<double>(*((double*)mxGetPi(mxa)));
201  return complex_t(realpart, imagpart);
202 }
203 
204 template <>
206  string rv = "";
207  if (mxGetClassID(mxa) != mxCHAR_CLASS) {
208  throw runtime_error("Can't construct string from anything but a char array.");
209  }
210  rv = string(mxArrayToString(mxa));
211  return rv;
212 }
213 
214 template <>
215 RCP<Xpetra_map> loadDataFromMatlab<RCP<Xpetra_map> >(const mxArray* mxa) {
217  int nr = mxGetM(mxa);
218  int nc = mxGetN(mxa);
219  if (nr != 1)
220  throw std::runtime_error("A Xpetra::Map representation from MATLAB must be a single row vector.");
221  double* pr = mxGetPr(mxa);
222  mm_GlobalOrd numGlobalIndices = nc;
223 
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]);
227  }
228 
229  const Teuchos::ArrayView<const mm_GlobalOrd> localGIDs_view(&localGIDs[0], localGIDs.size());
230  RCP<Xpetra_map> map =
234  localGIDs_view,
235  0, comm);
236 
237  if (map.is_null())
238  throw runtime_error("Failed to create Xpetra::Map.");
239  return map;
240 }
241 
242 template <>
243 RCP<Xpetra_ordinal_vector> loadDataFromMatlab<RCP<Xpetra_ordinal_vector> >(const mxArray* mxa) {
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);
252  if (map.is_null())
253  throw runtime_error("Failed to create map for Xpetra ordinal vector.");
255  if (loVec.is_null())
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]);
259  }
260  return loVec;
261 }
262 
263 template <>
264 RCP<Tpetra_MultiVector_double> loadDataFromMatlab<RCP<Tpetra_MultiVector_double> >(const mxArray* mxa) {
266  try {
267  int nr = mxGetM(mxa);
268  int nc = mxGetN(mxa);
269  double* pr = mxGetPr(mxa);
271  // numGlobalIndices for map constructor is the number of rows in matrix/vectors, right?
273  // Allocate a new array of complex values to use with the multivector
274  Teuchos::ArrayView<const double> arrView(pr, nr * nc);
275  mv = rcp(new Tpetra::MultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(map, arrView, size_t(nr), size_t(nc)));
276  } catch (std::exception& e) {
277  mexPrintf("Error constructing Tpetra MultiVector.\n");
278  std::cout << e.what() << std::endl;
279  }
280  return mv;
281 }
282 
283 template <>
284 RCP<Tpetra_MultiVector_complex> loadDataFromMatlab<RCP<Tpetra_MultiVector_complex> >(const mxArray* mxa) {
286  try {
287  int nr = mxGetM(mxa);
288  int nc = mxGetN(mxa);
289  double* pr = mxGetPr(mxa);
290  double* pi = mxGetPi(mxa);
292  // numGlobalIndices for map constructor is the number of rows in matrix/vectors, right?
294  // Allocate a new array of complex values to use with the multivector
295  complex_t* myArr = new complex_t[nr * nc];
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]);
299  }
300  }
301  Teuchos::ArrayView<complex_t> arrView(myArr, nr * nc);
303  } catch (std::exception& e) {
304  mexPrintf("Error constructing Tpetra MultiVector.\n");
305  std::cout << e.what() << std::endl;
306  }
307  return mv;
308 }
309 
310 template <>
311 RCP<Tpetra_CrsMatrix_double> loadDataFromMatlab<RCP<Tpetra_CrsMatrix_double> >(const mxArray* mxa) {
312  bool success = false;
314 
315  int* colptr = NULL;
316  int* rowind = NULL;
317 
318  try {
320  // numGlobalIndices is just the number of rows in the matrix
321  const size_t numGlobalIndices = mxGetM(mxa);
322  RCP<const muemex_map_type> rowMap = rcp(new muemex_map_type(numGlobalIndices, 0, comm));
323  RCP<const muemex_map_type> domainMap = rcp(new muemex_map_type(mxGetN(mxa), 0, comm));
324  double* valueArray = mxGetPr(mxa);
325  int nc = mxGetN(mxa);
326  if (rewrap_ints) {
327  // mwIndex_to_int allocates memory so must delete[] later
328  colptr = mwIndex_to_int(nc + 1, mxGetJc(mxa));
329  rowind = mwIndex_to_int(colptr[nc], mxGetIr(mxa));
330  } else {
331  rowind = (int*)mxGetIr(mxa);
332  colptr = (int*)mxGetJc(mxa);
333  }
334  // Need this to convert CSC colptrs to CRS row counts
335  Teuchos::Array<size_t> rowCounts(numGlobalIndices);
336  for (int i = 0; i < nc; i++) {
337  for (int j = colptr[i]; j < colptr[i + 1]; j++) {
338  rowCounts[rowind[j]]++;
339  }
340  }
342  for (int i = 0; i < nc; i++) {
343  for (int j = colptr[i]; j < colptr[i + 1]; j++) {
344  //'array' of 1 element, containing column (in global matrix).
346  //'array' of 1 element, containing value
347  Teuchos::ArrayView<double> vals = Teuchos::ArrayView<double>(&valueArray[j], 1);
348  A->insertGlobalValues(rowind[j], cols, vals);
349  }
350  }
351  A->fillComplete(domainMap, rowMap);
352  if (rewrap_ints) {
353  delete[] rowind;
354  rowind = NULL;
355  delete[] colptr;
356  colptr = NULL;
357  }
358  success = true;
359  } catch (std::exception& e) {
360  if (rewrap_ints) {
361  if (rowind != NULL) delete[] rowind;
362  if (colptr != NULL) delete[] colptr;
363  rowind = NULL;
364  colptr = NULL;
365  }
366  mexPrintf("Error while constructing Tpetra matrix:\n");
367  std::cout << e.what() << std::endl;
368  }
369  if (!success)
370  mexErrMsgTxt("An error occurred while setting up a Tpetra matrix.\n");
371  return A;
372 }
373 
374 template <>
375 RCP<Tpetra_CrsMatrix_complex> loadDataFromMatlab<RCP<Tpetra_CrsMatrix_complex> >(const mxArray* mxa) {
377  // Create a map in order to create the matrix (taken from muelu basic example - complex)
378  try {
380  const Tpetra::global_size_t numGlobalIndices = mxGetM(mxa);
381  const mm_GlobalOrd indexBase = 0;
382  RCP<const muemex_map_type> rowMap = rcp(new muemex_map_type(numGlobalIndices, indexBase, comm));
383  RCP<const muemex_map_type> domainMap = rcp(new muemex_map_type(mxGetN(mxa), indexBase, comm));
384  double* realArray = mxGetPr(mxa);
385  double* imagArray = mxGetPi(mxa);
386  int* colptr;
387  int* rowind;
388  int nc = mxGetN(mxa);
389  if (rewrap_ints) {
390  // mwIndex_to_int allocates memory so must delete[] later
391  colptr = mwIndex_to_int(nc + 1, mxGetJc(mxa));
392  rowind = mwIndex_to_int(colptr[nc], mxGetIr(mxa));
393  } else {
394  rowind = (int*)mxGetIr(mxa);
395  colptr = (int*)mxGetJc(mxa);
396  }
397  // Need this to convert CSC colptrs to CRS row counts
398  Teuchos::Array<size_t> rowCounts(numGlobalIndices);
399  for (int i = 0; i < nc; i++) {
400  for (int j = colptr[i]; j < colptr[i + 1]; j++) {
401  rowCounts[rowind[j]]++;
402  }
403  }
405  for (int i = 0; i < nc; i++) {
406  for (int j = colptr[i]; j < colptr[i + 1]; j++) {
407  // here assuming that complex_t will always be defined as std::complex<double>
408  // use 'value' over and over again with Teuchos::ArrayViews to insert into matrix
409  complex_t value = std::complex<double>(realArray[j], imagArray[j]);
412  A->insertGlobalValues(rowind[j], cols, vals);
413  }
414  }
415  A->fillComplete(domainMap, rowMap);
416  if (rewrap_ints) {
417  delete[] rowind;
418  delete[] colptr;
419  }
420  } catch (std::exception& e) {
421  mexPrintf("Error while constructing tpetra matrix:\n");
422  std::cout << e.what() << std::endl;
423  }
424  return A;
425 }
426 
427 template <>
428 RCP<Xpetra::Matrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > loadDataFromMatlab<RCP<Xpetra::Matrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > >(const mxArray* mxa) {
429  RCP<Tpetra::CrsMatrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > tmat = loadDataFromMatlab<RCP<Tpetra::CrsMatrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > >(mxa);
430  return MueLu::TpetraCrs_To_XpetraMatrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tmat);
431 }
432 
433 template <>
434 RCP<Xpetra::Matrix<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > loadDataFromMatlab<RCP<Xpetra::Matrix<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > >(const mxArray* mxa) {
435  RCP<Tpetra::CrsMatrix<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > tmat = loadDataFromMatlab<RCP<Tpetra::CrsMatrix<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > >(mxa);
436  return MueLu::TpetraCrs_To_XpetraMatrix<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tmat);
437 }
438 
439 template <>
440 RCP<Xpetra::MultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > loadDataFromMatlab<RCP<Xpetra::MultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > >(const mxArray* mxa) {
441  RCP<Tpetra::MultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > tpetraMV = loadDataFromMatlab<RCP<Tpetra::MultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > >(mxa);
442  return MueLu::TpetraMultiVector_To_XpetraMultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tpetraMV);
443 }
444 
445 template <>
446 RCP<Xpetra::MultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > loadDataFromMatlab<RCP<Xpetra::MultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > >(const mxArray* mxa) {
447  RCP<Tpetra::MultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > tpetraMV = loadDataFromMatlab<RCP<Tpetra::MultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t> > >(mxa);
448  return MueLu::TpetraMultiVector_To_XpetraMultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tpetraMV);
449 }
450 
451 #ifdef HAVE_MUELU_EPETRA
452 template <>
453 RCP<Epetra_CrsMatrix> loadDataFromMatlab<RCP<Epetra_CrsMatrix> >(const mxArray* mxa) {
454  RCP<Epetra_CrsMatrix> matrix;
455  try {
456  int* colptr;
457  int* rowind;
458  double* vals = mxGetPr(mxa);
459  int nr = mxGetM(mxa);
460  int nc = mxGetN(mxa);
461  if (rewrap_ints) {
462  colptr = mwIndex_to_int(nc + 1, mxGetJc(mxa));
463  rowind = mwIndex_to_int(colptr[nc], mxGetIr(mxa));
464  } else {
465  rowind = (int*)mxGetIr(mxa);
466  colptr = (int*)mxGetJc(mxa);
467  }
468  Epetra_SerialComm Comm;
469  Epetra_Map RangeMap(nr, 0, Comm);
470  Epetra_Map DomainMap(nc, 0, Comm);
471  matrix = rcp(new Epetra_CrsMatrix(Epetra_DataAccess::Copy, RangeMap, DomainMap, 0));
472  /* Do the matrix assembly */
473  for (int i = 0; i < nc; i++) {
474  for (int j = colptr[i]; j < colptr[i + 1]; j++) {
475  // global row, # of entries, value array, column indices array
476  matrix->InsertGlobalValues(rowind[j], 1, &vals[j], &i);
477  }
478  }
479  matrix->FillComplete(DomainMap, RangeMap);
480  if (rewrap_ints) {
481  delete[] rowind;
482  delete[] colptr;
483  }
484  } catch (std::exception& e) {
485  mexPrintf("An error occurred while setting up an Epetra matrix:\n");
486  std::cout << e.what() << std::endl;
487  }
488  return matrix;
489 }
490 
491 template <>
492 RCP<Epetra_MultiVector> loadDataFromMatlab<RCP<Epetra_MultiVector> >(const mxArray* mxa) {
493  int nr = mxGetM(mxa);
494  int nc = mxGetN(mxa);
495  Epetra_SerialComm Comm;
496  Epetra_BlockMap map(nr * nc, 1, 0, Comm);
497  return rcp(new Epetra_MultiVector(Epetra_DataAccess::Copy, map, mxGetPr(mxa), nr, nc));
498 }
499 #endif
500 
501 template <>
502 RCP<MAggregates> loadDataFromMatlab<RCP<MAggregates> >(const mxArray* 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.");
507  // assume that in matlab aggregate structs will only be stored in a 1x1 array
508  // mxa must have the same fields as the ones declared in constructAggregates function in muelu.m for this to work
509  const int correctNumFields = 5; // change if more fields are added to the aggregates representation in constructAggregates in muelu.m
510  if (mxGetNumberOfFields(mxa) != correctNumFields)
511  throw runtime_error("Aggregates structure has wrong number of fields.");
512  // Pull MuemexData types back out
513  int nVert = *(int*)mxGetData(mxGetField(mxa, 0, "nVertices"));
514  int nAgg = *(int*)mxGetData(mxGetField(mxa, 0, "nAggregates"));
515  // Now have all the data needed to fully reconstruct the aggregate
516  // Use similar approach as UserAggregationFactory (which is written for >1 thread but will just be serial here)
518  int myRank = comm->getRank();
521  RCP<MAggregates> agg = rcp(new MAggregates(map));
522  agg->SetNumAggregates(nAgg);
523  // Get handles for the vertex2AggId and procwinner arrays in reconstituted aggregates object
524  // this is serial so all procwinner values will be same (0)
525  ArrayRCP<mm_LocalOrd> vertex2AggId = agg->GetVertex2AggId()->getDataNonConst(0); // the '0' means first (and only) column of multivector, since is just vector
526  ArrayRCP<mm_LocalOrd> procWinner = agg->GetProcWinner()->getDataNonConst(0);
527  // mm_LocalOrd and int are equivalent, so is ok to talk about aggSize with just 'int'
528  // Deep copy the entire vertex2AggID and isRoot arrays, which are both nVert items long
529  // At the same time, set ProcWinner
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; // all nodes are going to be on the same proc
537  agg->SetIsRoot(i, false); // the ones that are root will be set in next loop
538  }
539  for (int i = 0; i < nAgg; i++) // rootNodesToCopy is an array of node IDs which are the roots of their aggs
540  {
541  agg->SetIsRoot(rootNodes_inArray[i], true);
542  }
543  // Now recompute the aggSize array the results in the object
544  agg->ComputeAggregateSizes(true);
545  agg->AggregatesCrossProcessors(false);
546  return agg;
547 }
548 
549 template <>
550 RCP<MAmalInfo> loadDataFromMatlab<RCP<MAmalInfo> >(const mxArray* mxa) {
551  RCP<MAmalInfo> amal;
552  throw runtime_error("AmalgamationInfo not supported in Muemex yet.");
553  return amal;
554 }
555 
556 template <>
557 RCP<MGraph> loadDataFromMatlab<RCP<MGraph> >(const mxArray* mxa) {
558  // mxa must be struct with logical sparse matrix called 'edges' and Nx1 int32 array 'boundaryNodes'
559  mxArray* edges = mxGetField(mxa, 0, "edges");
560  mxArray* boundaryNodes = mxGetField(mxa, 0, "boundaryNodes");
561  if (edges == NULL)
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.");
568  // Note that Matlab stores sparse matrices in column major format.
569  mwIndex* matlabColPtrs = mxGetJc(edges);
570  mwIndex* matlabRowIndices = mxGetIr(edges);
571  mm_GlobalOrd nRows = (mm_GlobalOrd)mxGetM(edges);
572 
573  // Create and populate row-major CRS data structures for Xpetra::TpetraCrsGraph.
574 
575  // calculate number of nonzeros in each row
576  Teuchos::Array<int> entriesPerRow(nRows);
577  int nnz = matlabColPtrs[mxGetN(edges)]; // last entry in matlabColPtrs
578  for (int i = 0; i < nnz; i++)
579  entriesPerRow[matlabRowIndices[i]]++;
580  // Populate usual row index array. We don't need this for the Xpetra Graph ctor, but
581  // it's convenient for building up the column index array, which the ctor does need.
582  Teuchos::Array<int> rows(nRows + 1);
583  rows[0] = 0;
584  for (int i = 0; i < nRows; i++)
585  rows[i + 1] = rows[i] + entriesPerRow[i];
586  Teuchos::Array<int> cols(nnz); // column index array
587  Teuchos::Array<int> insertionsPerRow(nRows, 0); // track of #insertions done per row
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]++;
595  }
596  }
597  // Find maximum
598  int maxNzPerRow = 0;
599  for (int i = 0; i < nRows; i++) {
600  if (maxNzPerRow < entriesPerRow[i])
601  maxNzPerRow = entriesPerRow[i];
602  }
603 
606  RCP<MMap> map = rcp(new MMap(nRows, 0, comm));
608  RCP<TpetraGraph> tgraph = rcp(new TpetraGraph(map, (size_t)maxNzPerRow));
609  // Populate tgraph in compressed-row format. Must get each row individually...
610  for (int i = 0; i < nRows; ++i) {
611  tgraph->insertGlobalIndices((mm_GlobalOrd)i, cols(rows[i], entriesPerRow[i]));
612  }
613  tgraph->fillComplete(map, map);
615  // Set boundary nodes
616  int numBoundaryNodes = mxGetNumberOfElements(boundaryNodes);
617  Kokkos::View<bool*, typename mm_node_t::memory_space> boundaryFlags("boundaryFlags", nRows);
618  // NOTE: This will not work correctly for non-CPU Node types
619  for (int i = 0; i < nRows; i++) {
620  boundaryFlags[i] = false;
621  }
622  for (int i = 0; i < numBoundaryNodes; i++) {
623  boundaryFlags[boundaryList[i]] = true;
624  }
625  mgraph->SetBoundaryNodeMap(boundaryFlags);
626  return mgraph;
627 }
628 
629 #ifdef HAVE_MUELU_INTREPID2
630 template <>
631 RCP<FieldContainer_ordinal> loadDataFromMatlab<RCP<FieldContainer_ordinal> >(const mxArray* mxa) {
632  if (mxGetClassID(mxa) != mxINT32_CLASS)
633  throw runtime_error("FieldContainer must have integer storage entries");
634 
635  int* data = (int*)mxGetData(mxa);
636  int nr = mxGetM(mxa);
637  int nc = mxGetN(mxa);
638 
639  RCP<FieldContainer_ordinal> fc = rcp(new FieldContainer_ordinal("FC from Matlab", nr, nc));
640  for (int col = 0; col < nc; col++) {
641  for (int row = 0; row < nr; row++) {
642  (*fc)(row, col) = data[col * nr + row];
643  }
644  }
645  return fc;
646 }
647 #endif
648 
649 /* ******************************* */
650 /* saveDataToMatlab */
651 /* ******************************* */
652 
653 template <>
655  mwSize dims[] = {1, 1};
656  mxArray* mxa = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
657  *((int*)mxGetData(mxa)) = data;
658  return mxa;
659 }
660 
661 template <>
662 mxArray* saveDataToMatlab(bool& data) {
663  mwSize dims[] = {1, 1};
664  mxArray* mxa = mxCreateLogicalArray(2, dims);
665  *((bool*)mxGetData(mxa)) = data;
666  return mxa;
667 }
668 
669 template <>
670 mxArray* saveDataToMatlab(double& data) {
671  return mxCreateDoubleScalar(data);
672 }
673 
674 template <>
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);
680  return mxa;
681 }
682 
683 template <>
684 mxArray* saveDataToMatlab(string& data) {
685  return mxCreateString(data.c_str());
686 }
687 
688 template <>
690  // Precondition: Memory has already been allocated by MATLAB for the array.
691  int nc = data->getGlobalNumElements();
692  int nr = 1;
693  mxArray* output = createMatlabMultiVector<double>(nr, nc);
694  double* array = (double*)malloc(sizeof(double) * nr * nc);
695  for (int col = 0; col < nc; col++) {
696  mm_GlobalOrd gid = data->getGlobalElement(col);
697  array[col] = Teuchos::as<double>(gid);
698  }
699  fillMatlabArray<double>(array, output, nc * nr);
700  free(array);
701  return output;
702 }
703 
704 template <>
706  mwSize len = data->getGlobalLength();
707  // create a single column vector
708  mwSize dimensions[] = {len, 1};
709  mxArray* rv = mxCreateNumericArray(2, dimensions, mxINT32_CLASS, mxREAL);
710  int* dataPtr = (int*)mxGetData(rv);
711  ArrayRCP<const mm_LocalOrd> arr = data->getData(0);
712  for (int i = 0; i < int(data->getGlobalLength()); i++) {
713  dataPtr[i] = arr[i];
714  }
715  return rv;
716 }
717 
718 template <>
721  return saveDataToMatlab(xmv);
722 }
723 
724 template <>
727  return saveDataToMatlab(xmv);
728 }
729 
730 template <>
733  return saveDataToMatlab(xmat);
734 }
735 
736 template <>
739  return saveDataToMatlab(xmat);
740 }
741 
742 template <>
744  typedef double Scalar;
745  // Compute global constants, if we need them
746  Teuchos::rcp_const_cast<Xpetra_CrsGraph>(data->getCrsGraph())->computeGlobalConstants();
747 
748  int nr = data->getGlobalNumRows();
749  int nc = data->getGlobalNumCols();
750  int nnz = data->getGlobalNumEntries();
751 
752 #ifdef VERBOSE_OUTPUT
753  RCP<Teuchos::FancyOStream> fancyStream = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
754  mat->describe(*fancyStream, Teuchos::VERB_EXTREME);
755 #endif
756  mxArray* mxa = createMatlabSparse<Scalar>(nr, nc, nnz);
757  mwIndex* ir = mxGetIr(mxa);
758  mwIndex* jc = mxGetJc(mxa);
759  for (int i = 0; i < nc + 1; i++) {
760  jc[i] = 0;
761  }
762 
763  size_t maxEntriesPerRow = data->getGlobalMaxNumRowEntries();
764  if (maxEntriesPerRow == Teuchos::OrdinalTraits<size_t>::invalid() || maxEntriesPerRow == 0) maxEntriesPerRow = data->getLocalMaxNumRowEntries();
765 
766  int* rowProgress = new int[nc];
767  // The array that will be copied to Pr and (if complex) Pi later
768  Scalar* sparseVals = new Scalar[nnz];
769  size_t numEntries;
770  if (data->isLocallyIndexed()) {
771  Scalar* rowValArray = new Scalar[maxEntriesPerRow];
772  Teuchos::ArrayView<Scalar> rowVals(rowValArray, maxEntriesPerRow);
773  mm_LocalOrd* rowIndicesArray = new mm_LocalOrd[maxEntriesPerRow];
774  Teuchos::ArrayView<mm_LocalOrd> rowIndices(rowIndicesArray, maxEntriesPerRow);
775  for (mm_LocalOrd m = 0; m < nr; m++) // All rows in the Xpetra matrix
776  {
777  data->getLocalRowCopy(m, rowIndices, rowVals, numEntries); // Get the row
778  for (mm_LocalOrd entry = 0; entry < int(numEntries); entry++) // All entries in row
779  {
780  jc[rowIndices[entry] + 1]++; // for each entry, increase jc for the entry's column
781  }
782  }
783 
784  // now jc holds the number of elements in each column, but needs cumulative sum over all previous columns also
785  int entriesAccum = 0;
786  for (int n = 0; n <= nc; n++) {
787  int temp = entriesAccum;
788  entriesAccum += jc[n];
789  jc[n] += temp;
790  }
791  // Jc now populated with colptrs
792  for (int i = 0; i < nc; i++) {
793  rowProgress[i] = 0;
794  }
795  // Row progress values like jc but keep track as the MATLAB matrix is being filled in
796  for (mm_LocalOrd m = 0; m < nr; m++) // rows
797  {
798  data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
799  for (mm_LocalOrd i = 0; i < int(numEntries); i++) // entries in row m (NOT columns)
800  {
801  // row is m, col is rowIndices[i], val is rowVals[i]
802  mm_LocalOrd col = rowIndices[i];
803  sparseVals[jc[col] + rowProgress[col]] = rowVals[i]; // Set value
804  ir[jc[col] + rowProgress[col]] = m; // Set row at which value occurs
805  rowProgress[col]++;
806  }
807  }
808  delete[] rowIndicesArray;
809  } else {
812  for (mm_GlobalOrd m = 0; m < nr; m++) {
813  data->getGlobalRowView(m, rowIndices, rowVals);
814  for (mm_GlobalOrd n = 0; n < rowIndices.size(); n++) {
815  jc[rowIndices[n] + 1]++;
816  }
817  }
818  // Last element of jc is just nnz
819  jc[nc] = nnz;
820  // Jc now populated with colptrs
821  for (int i = 0; i < nc; i++) {
822  rowProgress[i] = 0;
823  }
824  int entriesAccum = 0;
825  for (int n = 0; n <= nc; n++) {
826  int temp = entriesAccum;
827  entriesAccum += jc[n];
828  jc[n] += temp;
829  }
830  // Row progress values like jc but keep track as the MATLAB matrix is being filled in
831  for (mm_GlobalOrd m = 0; m < nr; m++) // rows
832  {
833  data->getGlobalRowView(m, rowIndices, rowVals);
834  for (mm_LocalOrd i = 0; i < rowIndices.size(); i++) // entries in row m
835  {
836  mm_GlobalOrd col = rowIndices[i]; // row is m, col is rowIndices[i], val is rowVals[i]
837  sparseVals[jc[col] + rowProgress[col]] = rowVals[i]; // Set value
838  ir[jc[col] + rowProgress[col]] = m; // Set row at which value occurs
839  rowProgress[col]++;
840  }
841  }
842  }
843  // finally, copy sparseVals into pr (and pi, if complex)
844  fillMatlabArray<Scalar>(sparseVals, mxa, nnz);
845  delete[] sparseVals;
846  delete[] rowProgress;
847  return mxa;
848 }
849 
850 template <>
852  typedef complex_t Scalar;
853 
854  // Compute global constants, if we need them
855  Teuchos::rcp_const_cast<Xpetra_CrsGraph>(data->getCrsGraph())->computeGlobalConstants();
856 
857  int nr = data->getGlobalNumRows();
858  int nc = data->getGlobalNumCols();
859  int nnz = data->getGlobalNumEntries();
860 #ifdef VERBOSE_OUTPUT
861  RCP<Teuchos::FancyOStream> fancyStream = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
862  mat->describe(*fancyStream, Teuchos::VERB_EXTREME);
863 #endif
864  mxArray* mxa = createMatlabSparse<Scalar>(nr, nc, nnz);
865  mwIndex* ir = mxGetIr(mxa);
866  mwIndex* jc = mxGetJc(mxa);
867  for (int i = 0; i < nc + 1; i++) {
868  jc[i] = 0;
869  }
870  size_t maxEntriesPerRow = data->getGlobalMaxNumRowEntries();
871  int* rowProgress = new int[nc];
872  // The array that will be copied to Pr and (if complex) Pi later
873  Scalar* sparseVals = new Scalar[nnz];
874  size_t numEntries;
875  if (data->isLocallyIndexed()) {
876  Scalar* rowValArray = new Scalar[maxEntriesPerRow];
877  Teuchos::ArrayView<Scalar> rowVals(rowValArray, maxEntriesPerRow);
878  mm_LocalOrd* rowIndicesArray = new mm_LocalOrd[maxEntriesPerRow];
879  Teuchos::ArrayView<mm_LocalOrd> rowIndices(rowIndicesArray, maxEntriesPerRow);
880  for (mm_LocalOrd m = 0; m < nr; m++) // All rows in the Xpetra matrix
881  {
882  data->getLocalRowCopy(m, rowIndices, rowVals, numEntries); // Get the row
883  for (mm_LocalOrd entry = 0; entry < int(numEntries); entry++) // All entries in row
884  {
885  jc[rowIndices[entry] + 1]++; // for each entry, increase jc for the entry's column
886  }
887  }
888  // now jc holds the number of elements in each column, but needs cumulative sum over all previous columns also
889  int entriesAccum = 0;
890  for (int n = 0; n <= nc; n++) {
891  int temp = entriesAccum;
892  entriesAccum += jc[n];
893  jc[n] += temp;
894  }
895  // Jc now populated with colptrs
896  for (int i = 0; i < nc; i++) {
897  rowProgress[i] = 0;
898  }
899  // Row progress values like jc but keep track as the MATLAB matrix is being filled in
900  for (mm_LocalOrd m = 0; m < nr; m++) // rows
901  {
902  data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
903  for (mm_LocalOrd i = 0; i < int(numEntries); i++) // entries in row m (NOT columns)
904  {
905  // row is m, col is rowIndices[i], val is rowVals[i]
906  mm_LocalOrd col = rowIndices[i];
907  sparseVals[jc[col] + rowProgress[col]] = rowVals[i]; // Set value
908  ir[jc[col] + rowProgress[col]] = m; // Set row at which value occurs
909  rowProgress[col]++;
910  }
911  }
912  delete[] rowIndicesArray;
913  } else {
916  for (mm_GlobalOrd m = 0; m < nr; m++) {
917  data->getGlobalRowView(m, rowIndices, rowVals);
918  for (mm_GlobalOrd n = 0; n < rowIndices.size(); n++) {
919  jc[rowIndices[n] + 1]++;
920  }
921  }
922  // Last element of jc is just nnz
923  jc[nc] = nnz;
924  // Jc now populated with colptrs
925  for (int i = 0; i < nc; i++) {
926  rowProgress[i] = 0;
927  }
928  int entriesAccum = 0;
929  for (int n = 0; n <= nc; n++) {
930  int temp = entriesAccum;
931  entriesAccum += jc[n];
932  jc[n] += temp;
933  }
934  // Row progress values like jc but keep track as the MATLAB matrix is being filled in
935  for (mm_GlobalOrd m = 0; m < nr; m++) // rows
936  {
937  data->getGlobalRowView(m, rowIndices, rowVals);
938  for (mm_LocalOrd i = 0; i < rowIndices.size(); i++) // entries in row m
939  {
940  mm_GlobalOrd col = rowIndices[i]; // row is m, col is rowIndices[i], val is rowVals[i]
941  sparseVals[jc[col] + rowProgress[col]] = rowVals[i]; // Set value
942  ir[jc[col] + rowProgress[col]] = m; // Set row at which value occurs
943  rowProgress[col]++;
944  }
945  }
946  }
947  // finally, copy sparseVals into pr (and pi, if complex)
948  fillMatlabArray<Scalar>(sparseVals, mxa, nnz);
949  delete[] sparseVals;
950  delete[] rowProgress;
951  return mxa;
952 }
953 
954 /*
955 template<>
956 mxArray* saveDataToMatlab(RCP<Xpetra::MultiVector<Scalar, mm_LocalOrd, mm_GlobalOrd, mm_node_t>>& data)
957 {
958  //Precondition: Memory has already been allocated by MATLAB for the array.
959  int nr = data->getGlobalLength();
960  int nc = data->getNumVectors();
961  mxArray* output = createMatlabMultiVector<Scalar>(nr, nc);
962  Scalar* array = (Scalar*) malloc(sizeof(Scalar) * nr * nc);
963  for(int col = 0; col < nc; col++)
964  {
965  Teuchos::ArrayRCP<const Scalar> colData = data->getData(col);
966  for(int row = 0; row < nr; row++)
967  {
968  array[col * nr + row] = colData[row];
969  }
970  }
971  fillMatlabArray<Scalar>(array, output, nc * nr);
972  free(array);
973  return output;
974 }
975 */
976 
977 template <>
979  // Precondition: Memory has already been allocated by MATLAB for the array.
980  int nr = data->getGlobalLength();
981  int nc = data->getNumVectors();
982  mxArray* output = createMatlabMultiVector<double>(nr, nc);
983  double* array = (double*)malloc(sizeof(double) * nr * nc);
984  for (int col = 0; col < nc; col++) {
985  Teuchos::ArrayRCP<const double> colData = data->getData(col);
986  for (int row = 0; row < nr; row++) {
987  array[col * nr + row] = colData[row];
988  }
989  }
990  fillMatlabArray<double>(array, output, nc * nr);
991  free(array);
992  return output;
993 }
994 
995 template <>
997  // Precondition: Memory has already been allocated by MATLAB for the array.
998  int nr = data->getGlobalLength();
999  int nc = data->getNumVectors();
1000  mxArray* output = createMatlabMultiVector<complex_t>(nr, nc);
1001  complex_t* array = (complex_t*)malloc(sizeof(complex_t) * nr * nc);
1002  for (int col = 0; col < nc; col++) {
1003  Teuchos::ArrayRCP<const complex_t> colData = data->getData(col);
1004  for (int row = 0; row < nr; row++) {
1005  array[col * nr + row] = colData[row];
1006  }
1007  }
1008  fillMatlabArray<complex_t>(array, output, nc * nr);
1009  free(array);
1010  return output;
1011 }
1012 
1013 #ifdef HAVE_MUELU_EPETRA
1014 template <>
1016  RCP<Xpetra_Matrix_double> xmat = EpetraCrs_To_XpetraMatrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(data);
1017  return saveDataToMatlab(xmat);
1018 }
1019 
1020 template <>
1022  mxArray* output = mxCreateDoubleMatrix(data->GlobalLength(), data->NumVectors(), mxREAL);
1023  double* dataPtr = mxGetPr(output);
1024  data->ExtractCopy(dataPtr, data->GlobalLength());
1025  return output;
1026 }
1027 #endif
1028 
1029 template <>
1031  // Set up array of inputs for matlab constructAggregates
1032  int numNodes = data->GetVertex2AggId()->getData(0).size();
1033  int numAggs = data->GetNumAggregates();
1034  mxArray* dataIn[5];
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}; // dimensions for Nx1 array, where N is number of nodes (vert2Agg)
1041  dataIn[2] = mxCreateNumericArray(2, nodeArrayDims, mxINT32_CLASS, mxREAL);
1042  int* vtaid = (int*)mxGetData(dataIn[2]);
1043  ArrayRCP<const mm_LocalOrd> vertexToAggID = data->GetVertex2AggId()->getData(0);
1044  for (int i = 0; i < numNodes; i++) {
1045  vtaid[i] = vertexToAggID[i];
1046  }
1047  mwSize aggArrayDims[] = {(mwSize)numAggs, 1}; // dims for Nx1 array, where N is number of aggregates (rootNodes, aggSizes)
1048  dataIn[3] = mxCreateNumericArray(2, aggArrayDims, mxINT32_CLASS, mxREAL);
1049  // First, find out if the aggregates even have 1 root node per aggregate. If not, assume roots are invalid and assign ourselves
1050  int totalRoots = 0;
1051  for (int i = 0; i < numNodes; i++) {
1052  if (data->IsRoot(i))
1053  totalRoots++;
1054  }
1055  bool reassignRoots = false;
1056  if (totalRoots != numAggs) {
1057  cout << endl
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
1060  << endl;
1061  reassignRoots = true;
1062  }
1063  int* rn = (int*)mxGetData(dataIn[3]); // list of root nodes (in no particular order)
1064  {
1065  if (reassignRoots) {
1066  // For each aggregate, just pick the first node we see in it and set it as root
1067  int lastFoundNode = 0; // heuristic for speed, a node in aggregate N+1 is likely to come very soon after a node in agg N
1068  for (int i = 0; i < numAggs; i++) {
1069  rn[i] = -1;
1070  for (int j = lastFoundNode; j < lastFoundNode + numNodes; j++) {
1071  int index = j % numNodes;
1072  if (vertexToAggID[index] == i) {
1073  rn[i] = index;
1074  lastFoundNode = index;
1075  }
1076  }
1077  TEUCHOS_TEST_FOR_EXCEPTION(rn[i] == -1, runtime_error, "Invalid aggregates: Couldn't find any node in aggregate #" << i << ".");
1078  }
1079  } else {
1080  int i = 0; // iterates over aggregate IDs
1081  for (int j = 0; j < numNodes; j++) {
1082  if (data->IsRoot(j)) {
1083  if (i == numAggs)
1084  throw runtime_error("Cannot store invalid aggregates in MATLAB - more root nodes than aggregates.");
1085  rn[i] = j; // now we know this won't go out of bounds (rn's underlying matlab array is numAggs in length)
1086  i++;
1087  }
1088  }
1089  if (i + 1 < numAggs)
1090  throw runtime_error("Cannot store invalid aggregates in MATLAB - fewer root nodes than aggregates.");
1091  }
1092  }
1093  dataIn[4] = mxCreateNumericArray(1, aggArrayDims, mxINT32_CLASS, mxREAL);
1094  int* as = (int*)mxGetData(dataIn[4]); // list of aggregate sizes
1095  auto aggSizes = data->ComputeAggregateSizes();
1096  for (int i = 0; i < numAggs; i++) {
1097  as[i] = aggSizes[i];
1098  }
1099  mxArray* matlabAggs[1];
1100  int result = mexCallMATLAB(1, matlabAggs, 5, dataIn, "constructAggregates");
1101  if (result != 0)
1102  throw runtime_error("Matlab encountered an error while constructing aggregates struct.");
1103  return matlabAggs[0];
1104 }
1105 
1106 template <>
1108  throw runtime_error("AmalgamationInfo not supported in MueMex yet.");
1109  return mxCreateDoubleScalar(0);
1110 }
1111 
1112 template <>
1114  int numEntries = (int)data->GetGlobalNumEdges();
1115  int numRows = (int)data->GetDomainMap()->getGlobalNumElements(); // assume numRows == numCols
1116  mxArray* mat = mxCreateSparseLogicalMatrix(numRows, numRows, numEntries);
1117  mxLogical* outData = (mxLogical*)mxGetData(mat);
1118  mwIndex* rowInds = mxGetIr(mat);
1119  mwIndex* colPtrs = mxGetJc(mat);
1120  mm_LocalOrd* dataCopy = new mm_LocalOrd[numEntries];
1121  mm_LocalOrd* iter = dataCopy;
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;
1127  }
1128  for (int i = 0; i < numRows; i++) {
1129  ArrayView<typename MGraph::local_ordinal_type> neighbors = data->getNeighborVertices_av(i); // neighbors has the column indices for row i
1130  memcpy(iter, neighbors.getRawPtr(), sizeof(mm_LocalOrd) * neighbors.size());
1131  entriesPerRow[i] = neighbors.size();
1132  for (int j = 0; j < neighbors.size(); j++) {
1133  entriesPerCol[neighbors[j]]++;
1134  }
1135  iter += neighbors.size();
1136  }
1137  mwIndex** rowIndsByColumn = new mwIndex*[numRows]; // rowIndsByColumn[0] points to array of row indices in column 1
1138  mxLogical** valuesByColumn = new mxLogical*[numRows];
1139  int* numEnteredPerCol = new int[numRows];
1140  int accum = 0;
1141  for (int i = 0; i < numRows; i++) {
1142  rowIndsByColumn[i] = &rowInds[accum];
1143  // cout << "Entries in column " << i << " start at offset " << accum << endl;
1144  valuesByColumn[i] = &outData[accum];
1145  accum += entriesPerCol[i];
1146  if (accum > numEntries)
1147  throw runtime_error("potato");
1148  }
1149  for (int i = 0; i < numRows; i++) {
1150  numEnteredPerCol[i] = 0; // rowIndsByColumn[n][numEnteredPerCol[n]] gives the next place to put a row index
1151  }
1152  // entriesPerCol now has Jc information (col offsets)
1153  accum = 0; // keep track of cumulative index in dataCopy
1154  for (int row = 0; row < numRows; row++) {
1155  for (int entryInRow = 0; entryInRow < entriesPerRow[row]; entryInRow++) {
1156  int col = dataCopy[accum];
1157  accum++;
1158  rowIndsByColumn[col][numEnteredPerCol[col]] = row;
1159  valuesByColumn[col][numEnteredPerCol[col]] = (mxLogical)1;
1160  numEnteredPerCol[col]++;
1161  }
1162  }
1163  accum = 0; // keep track of total entries over all columns
1164  for (int col = 0; col < numRows; col++) {
1165  colPtrs[col] = accum;
1166  accum += entriesPerCol[col];
1167  }
1168  colPtrs[numRows] = accum; // the last entry in jc, which is equivalent to numEntries
1169  delete[] numEnteredPerCol;
1170  delete[] rowIndsByColumn;
1171  delete[] valuesByColumn;
1172  delete[] dataCopy;
1173  delete[] entriesPerRow;
1174  delete[] entriesPerCol;
1175  // Construct list of boundary nodes
1176  auto boundaryFlags = data->GetBoundaryNodeMap();
1177  int numBoundaryNodes = 0;
1178  for (int i = 0; i < (int)boundaryFlags.size(); i++) {
1179  if (boundaryFlags[i])
1180  numBoundaryNodes++;
1181  }
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]) {
1189  *destIter = i;
1190  destIter++;
1191  }
1192  }
1193  mxArray* constructArgs[] = {mat, boundaryList};
1194  mxArray* out[1];
1195  mexCallMATLAB(1, out, 2, constructArgs, "constructGraph");
1196  return out[0];
1197 }
1198 
1199 #ifdef HAVE_MUELU_INTREPID2
1200 template <>
1202  int rank = data->rank();
1203  // NOTE: Only supports rank 2 arrays
1204  if (rank != 2)
1205  throw std::runtime_error("Error: Only rank two FieldContainers are supported.");
1206 
1207  int nr = data->extent(0);
1208  int nc = data->extent(1);
1209 
1210  mwSize dims[] = {(mwSize)nr, (mwSize)nc};
1211  mxArray* mxa = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
1212  int* array = (int*)mxGetData(mxa);
1213 
1214  for (int col = 0; col < nc; col++) {
1215  for (int row = 0; row < nr; row++) {
1216  array[col * nr + row] = (*data)(row, col);
1217  }
1218  }
1219  return mxa;
1220 }
1221 #endif
1222 
1223 template <typename T>
1225  : MuemexArg(getMuemexType<T>()) {
1226  data = loadDataFromMatlab<T>(mxa);
1227 }
1228 
1229 template <typename T>
1231  return saveDataToMatlab<T>(data);
1232 }
1233 
1234 template <typename T>
1235 MuemexData<T>::MuemexData(T& dataToCopy, MuemexType dataType)
1236  : MuemexArg(dataType) {
1237  data = dataToCopy;
1238 }
1239 
1240 template <typename T>
1242  : MuemexData(dataToCopy, getMuemexType(dataToCopy)) {}
1243 
1244 template <typename T>
1246  return data;
1247 }
1248 
1249 template <typename T>
1250 void MuemexData<T>::setData(T& newData) {
1251  this->data = newData;
1252 }
1253 
1254 /* ***************************** */
1255 /* More Template Functions */
1256 /* ***************************** */
1257 
1258 template <typename T>
1259 void addLevelVariable(const T& data, std::string& name, Level& lvl, const Factory* fact) {
1260  lvl.AddKeepFlag(name, fact, MueLu::UserData);
1261  lvl.Set<T>(name, data, fact);
1262 }
1263 
1264 template <typename T>
1265 const T& getLevelVariable(std::string& name, Level& lvl) {
1266  try {
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.");
1270  }
1271 }
1272 
1273 // Functions used to put data through matlab factories - first arg is "this" pointer of matlab factory
1274 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, typename Node>
1275 std::vector<Teuchos::RCP<MuemexArg> > processNeeds(const Factory* factory, std::string& needsParam, Level& lvl) {
1276  using namespace std;
1277  using namespace Teuchos;
1281  typedef RCP<AmalgamationInfo<LocalOrdinal, GlobalOrdinal, Node> > AmalgamationInfo_t;
1282  typedef RCP<MGraph> Graph_t;
1283  vector<string> needsList = tokenizeList(needsParam);
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());
1288  args.push_back(rcp(new MuemexData<Matrix_t>(mydata)));
1289  } else if (needsList[i] == "Nullspace" || needsList[i] == "Coordinates") {
1290  MultiVector_t mydata = lvl.Get<MultiVector_t>(needsList[i], factory->GetFactory(needsList[i]).get());
1291  args.push_back(rcp(new MuemexData<MultiVector_t>(mydata)));
1292  } else if (needsList[i] == "Aggregates") {
1293  Aggregates_t mydata = lvl.Get<Aggregates_t>(needsList[i], factory->GetFactory(needsList[i]).get());
1294  args.push_back(rcp(new MuemexData<Aggregates_t>(mydata)));
1295  } else if (needsList[i] == "UnAmalgamationInfo") {
1296  AmalgamationInfo_t mydata = lvl.Get<AmalgamationInfo_t>(needsList[i], factory->GetFactory(needsList[i]).get());
1297  args.push_back(rcp(new MuemexData<AmalgamationInfo_t>(mydata)));
1298  } else if (needsList[i] == "Level") {
1299  int levelNum = lvl.GetLevelID();
1300  args.push_back(rcp(new MuemexData<int>(levelNum)));
1301  } else if (needsList[i] == "Graph") {
1302  Graph_t mydata = lvl.Get<Graph_t>(needsList[i], factory->GetFactory(needsList[i]).get());
1303  args.push_back(rcp(new MuemexData<Graph_t>(mydata)));
1304  } else {
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";
1307  // compare type without case sensitivity
1308  char* buf = (char*)malloc(needsList[i].size() + 1);
1309  strcpy(buf, needsList[i].c_str());
1310  for (char* iter = buf; *iter != ' '; iter++) {
1311  if (*iter == 0) {
1312  free(buf);
1313  throw runtime_error(badNameMsg);
1314  }
1315  *iter = (char)tolower(*iter);
1316  }
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);
1323  }
1324  if (words.size() != 2) {
1325  free(buf);
1326  throw runtime_error(badNameMsg);
1327  }
1328  char* typeStr = (char*)words[0].c_str();
1329  if (strstr(typeStr, "ordinalvector")) {
1331  LOVector_t mydata = getLevelVariable<LOVector_t>(needsList[i], lvl);
1332  args.push_back(rcp(new MuemexData<LOVector_t>(mydata)));
1333  } else if (strstr(typeStr, "map")) {
1335  Map_t mydata = getLevelVariable<Map_t>(needsList[i], lvl);
1336  args.push_back(rcp(new MuemexData<Map_t>(mydata)));
1337  } else if (strstr(typeStr, "scalar")) {
1338  Scalar mydata = getLevelVariable<Scalar>(needsList[i], lvl);
1339  args.push_back(rcp(new MuemexData<Scalar>(mydata)));
1340  } else if (strstr(typeStr, "double")) {
1341  double mydata = getLevelVariable<double>(needsList[i], lvl);
1342  args.push_back(rcp(new MuemexData<double>(mydata)));
1343  } else if (strstr(typeStr, "complex")) {
1344  complex_t mydata = getLevelVariable<complex_t>(needsList[i], lvl);
1345  args.push_back(rcp(new MuemexData<complex_t>(mydata)));
1346  } else if (strstr(typeStr, "matrix")) {
1347  Matrix_t mydata = getLevelVariable<Matrix_t>(needsList[i], lvl);
1348  args.push_back(rcp(new MuemexData<Matrix_t>(mydata)));
1349  } else if (strstr(typeStr, "multivector")) {
1350  MultiVector_t mydata = getLevelVariable<MultiVector_t>(needsList[i], lvl);
1351  args.push_back(rcp(new MuemexData<MultiVector_t>(mydata)));
1352  } else if (strstr(typeStr, "int")) {
1353  int mydata = getLevelVariable<int>(needsList[i], lvl);
1354  args.push_back(rcp(new MuemexData<int>(mydata)));
1355  } else if (strstr(typeStr, "string")) {
1356  string mydata = getLevelVariable<string>(needsList[i], lvl);
1357  args.push_back(rcp(new MuemexData<string>(mydata)));
1358  } else {
1359  free(buf);
1360  throw std::runtime_error(words[0] + " is not a known variable type.");
1361  }
1362  free(buf);
1363  }
1364  }
1365  return args;
1366 }
1367 
1368 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, typename Node>
1369 void processProvides(std::vector<Teuchos::RCP<MuemexArg> >& mexOutput, const Factory* factory, std::string& providesParam, Level& lvl) {
1370  using namespace std;
1371  using namespace Teuchos;
1375  typedef RCP<AmalgamationInfo<LocalOrdinal, GlobalOrdinal, Node> > AmalgamationInfo_t;
1376  typedef RCP<MGraph> Graph_t;
1377  vector<string> provides = tokenizeList(providesParam);
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") {
1380  RCP<MuemexData<Matrix_t> > mydata = Teuchos::rcp_static_cast<MuemexData<Matrix_t> >(mexOutput[i]);
1381  lvl.Set(provides[i], mydata->getData(), factory);
1382  } else if (provides[i] == "Nullspace" || provides[i] == "Coordinates") {
1383  RCP<MuemexData<MultiVector_t> > mydata = Teuchos::rcp_static_cast<MuemexData<MultiVector_t> >(mexOutput[i]);
1384  lvl.Set(provides[i], mydata->getData(), factory);
1385  } else if (provides[i] == "Aggregates") {
1386  RCP<MuemexData<Aggregates_t> > mydata = Teuchos::rcp_static_cast<MuemexData<Aggregates_t> >(mexOutput[i]);
1387  lvl.Set(provides[i], mydata->getData(), factory);
1388  } else if (provides[i] == "UnAmalgamationInfo") {
1389  RCP<MuemexData<AmalgamationInfo_t> > mydata = Teuchos::rcp_static_cast<MuemexData<AmalgamationInfo_t> >(mexOutput[i]);
1390  lvl.Set(provides[i], mydata->getData(), factory);
1391  } else if (provides[i] == "Graph") {
1392  RCP<MuemexData<Graph_t> > mydata = Teuchos::rcp_static_cast<MuemexData<Graph_t> >(mexOutput[i]);
1393  lvl.Set(provides[i], mydata->getData(), factory);
1394  } else {
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";
1397  // compare type without case sensitivity
1398  char* buf = (char*)malloc(provides[i].size() + 1);
1399  strcpy(buf, provides[i].c_str());
1400  for (char* iter = buf; *iter != ' '; iter++) {
1401  if (*iter == 0) {
1402  free(buf);
1403  throw runtime_error(badNameMsg);
1404  }
1405  *iter = (char)tolower(*iter);
1406  }
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);
1413  }
1414  if (words.size() != 2) {
1415  free(buf);
1416  throw runtime_error(badNameMsg);
1417  }
1418  const char* typeStr = words[0].c_str();
1419  if (strstr(typeStr, "ordinalvector")) {
1421  RCP<MuemexData<LOVector_t> > mydata = Teuchos::rcp_static_cast<MuemexData<LOVector_t> >(mexOutput[i]);
1422  addLevelVariable<LOVector_t>(mydata->getData(), words[1], lvl, factory);
1423  } else if (strstr(typeStr, "map")) {
1425  RCP<MuemexData<Map_t> > mydata = Teuchos::rcp_static_cast<MuemexData<Map_t> >(mexOutput[i]);
1426  addLevelVariable<Map_t>(mydata->getData(), words[1], lvl, factory);
1427  } else if (strstr(typeStr, "scalar")) {
1428  RCP<MuemexData<Scalar> > mydata = Teuchos::rcp_static_cast<MuemexData<Scalar> >(mexOutput[i]);
1429  addLevelVariable<Scalar>(mydata->getData(), words[1], lvl, factory);
1430  } else if (strstr(typeStr, "double")) {
1431  RCP<MuemexData<double> > mydata = Teuchos::rcp_static_cast<MuemexData<double> >(mexOutput[i]);
1432  addLevelVariable<double>(mydata->getData(), words[1], lvl, factory);
1433  } else if (strstr(typeStr, "complex")) {
1434  RCP<MuemexData<complex_t> > mydata = Teuchos::rcp_static_cast<MuemexData<complex_t> >(mexOutput[i]);
1435  addLevelVariable<complex_t>(mydata->getData(), words[1], lvl, factory);
1436  } else if (strstr(typeStr, "matrix")) {
1437  RCP<MuemexData<Matrix_t> > mydata = Teuchos::rcp_static_cast<MuemexData<Matrix_t> >(mexOutput[i]);
1438  addLevelVariable<Matrix_t>(mydata->getData(), words[1], lvl, factory);
1439  } else if (strstr(typeStr, "multivector")) {
1440  RCP<MuemexData<MultiVector_t> > mydata = Teuchos::rcp_static_cast<MuemexData<MultiVector_t> >(mexOutput[i]);
1441  addLevelVariable<MultiVector_t>(mydata->getData(), words[1], lvl, factory);
1442  } else if (strstr(typeStr, "int")) {
1443  RCP<MuemexData<int> > mydata = Teuchos::rcp_static_cast<MuemexData<int> >(mexOutput[i]);
1444  addLevelVariable<int>(mydata->getData(), words[1], lvl, factory);
1445  } else if (strstr(typeStr, "bool")) {
1446  RCP<MuemexData<bool> > mydata = Teuchos::rcp_static_cast<MuemexData<bool> >(mexOutput[i]);
1447  addLevelVariable<bool>(mydata->getData(), words[1], lvl, factory);
1448  } else if (strstr(typeStr, "string")) {
1449  RCP<MuemexData<string> > mydata = Teuchos::rcp_static_cast<MuemexData<string> >(mexOutput[i]);
1450  addLevelVariable<string>(mydata->getData(), words[1], lvl, factory);
1451  } else {
1452  free(buf);
1453  throw std::runtime_error(words[0] + " is not a known variable type.");
1454  }
1455  free(buf);
1456  }
1457  }
1458 }
1459 
1460 // Throwable Stubs for long long
1461 
1462 template <>
1463 std::vector<Teuchos::RCP<MuemexArg> > processNeeds<double, mm_LocalOrd, long long, mm_node_t>(const Factory* factory, std::string& needsParam, Level& lvl) {
1464  throw std::runtime_error("Muemex does not support long long for global indices");
1465 }
1466 
1467 template <>
1468 std::vector<Teuchos::RCP<MuemexArg> > processNeeds<complex_t, mm_LocalOrd, long long, mm_node_t>(const Factory* factory, std::string& needsParam, Level& lvl) {
1469  throw std::runtime_error("Muemex does not support long long for global indices");
1470 }
1471 
1472 template <>
1473 void processProvides<double, mm_LocalOrd, long long, mm_node_t>(std::vector<Teuchos::RCP<MuemexArg> >& mexOutput, const Factory* factory, std::string& providesParam, Level& lvl) {
1474  throw std::runtime_error("Muemex does not support long long for global indices");
1475 }
1476 
1477 template <>
1478 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) {
1479  throw std::runtime_error("Muemex does not support long long for global indices");
1480 }
1481 
1482 } // namespace MueLu
1483 #endif // HAVE_MUELU_MATLAB error handler
1484 #endif // MUELU_MATLABUTILS_DEF_HPP guard
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 &params)
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-&gt;Get&lt; RCP&lt;Matrix&gt; &gt;(&quot;A&quot;, factory) if factory == NULL =&gt; use default factory.
mxArray * createMatlabSparse< complex_t >(int numRows, int numCols, int nnz)
bool rewrap_ints
bool is_null() const
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)
size_type size() const
User data are always kept. This flag is set automatically when Level::Set(&quot;data&quot;, 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)
size_t global_size_t
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.
Definition: MueLu_Level.hpp:63
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)
T * getRawPtr() const
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)
TransListIter iter
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.
Definition: MueLu_Level.cpp:51
MuemexType getMuemexType< complex_t >()
void processProvides(std::vector< Teuchos::RCP< MuemexArg > > &mexOutput, const Factory *factory, std::string &providesParam, Level &lvl)
int mwIndex
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Default implementation of FactoryAcceptor::GetFactory()
bool is_null() const