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 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 
47 #ifndef MUELU_MATLABUTILS_DEF_HPP
48 #define MUELU_MATLABUTILS_DEF_HPP
49 
51 
52 #if !defined(HAVE_MUELU_MATLAB) || !defined(HAVE_MUELU_EPETRA) || !defined(HAVE_MUELU_TPETRA)
53 #error "Muemex types require MATLAB, Epetra and Tpetra."
54 #else
55 
56 using Teuchos::RCP;
57 using Teuchos::rcp;
58 using namespace std;
59 
60 namespace MueLu {
61 
62 extern bool rewrap_ints;
63 
64 /* ******************************* */
65 /* getMuemexType */
66 /* ******************************* */
67 
68 template<typename T> MuemexType getMuemexType(const T & data) {throw std::runtime_error("Unknown Type");}
69 
70 template<> MuemexType getMuemexType(const int & data) {return INT;}
71 template<> MuemexType getMuemexType<int>() {return INT;}
72 template<> MuemexType getMuemexType<bool>() {return BOOL;}
73 
74 template<> MuemexType getMuemexType(const double & data) {return DOUBLE;}
75 template<> MuemexType getMuemexType<double>() {return DOUBLE;}
76 
77 template<> MuemexType getMuemexType(const std::string & data) {return STRING;}
78 template<> MuemexType getMuemexType<string>() {return STRING;}
79 
80 template<> MuemexType getMuemexType(const complex_t& data) {return COMPLEX;}
82 
83 template<> MuemexType getMuemexType(const RCP<Xpetra_map> & data) {return XPETRA_MAP;}
84 template<> MuemexType getMuemexType<RCP<Xpetra_map> >() {return XPETRA_MAP;}
85 
87 template<> MuemexType getMuemexType<RCP<Xpetra_ordinal_vector>>() {return XPETRA_ORDINAL_VECTOR;}
88 
90 template<> MuemexType getMuemexType<RCP<Tpetra_MultiVector_double>>() {return TPETRA_MULTIVECTOR_DOUBLE;}
91 
93 template<> MuemexType getMuemexType<RCP<Tpetra_MultiVector_complex>>() {return TPETRA_MULTIVECTOR_COMPLEX;}
94 
96 template<> MuemexType getMuemexType<RCP<Tpetra_CrsMatrix_double>>() {return TPETRA_MATRIX_DOUBLE;}
97 
99 template<> MuemexType getMuemexType<RCP<Tpetra_CrsMatrix_complex>>() {return TPETRA_MATRIX_COMPLEX;}
100 
102 template<> MuemexType getMuemexType<RCP<Xpetra_MultiVector_double>>() {return XPETRA_MULTIVECTOR_DOUBLE;}
103 
105 template<> MuemexType getMuemexType<RCP<Xpetra_MultiVector_complex>>() {return XPETRA_MULTIVECTOR_COMPLEX;}
106 
108 template<> MuemexType getMuemexType<RCP<Xpetra_Matrix_double>>() {return XPETRA_MATRIX_DOUBLE;}
109 
111 template<> MuemexType getMuemexType<RCP<Xpetra_Matrix_complex>>() {return XPETRA_MATRIX_COMPLEX;}
112 
114 template<> MuemexType getMuemexType<RCP<Epetra_CrsMatrix>>() {return EPETRA_CRSMATRIX;}
115 
117 template<> MuemexType getMuemexType<RCP<Epetra_MultiVector>>() {return EPETRA_MULTIVECTOR;}
118 
119 template<> MuemexType getMuemexType(const RCP<MAggregates>& data) {return AGGREGATES;}
120 template<> MuemexType getMuemexType<RCP<MAggregates>>() {return AGGREGATES;}
121 
122 template<> MuemexType getMuemexType(const RCP<MAmalInfo>& data) {return AMALGAMATION_INFO;}
123 template<> MuemexType getMuemexType<RCP<MAmalInfo>>() {return AMALGAMATION_INFO;}
124 
125 template<> MuemexType getMuemexType(const RCP<MGraph>& data) {return GRAPH;}
126 template<> MuemexType getMuemexType<RCP<MGraph>>() {return GRAPH;}
127 
128 #ifdef HAVE_MUELU_INTREPID2
129 template<> MuemexType getMuemexType(const RCP<FieldContainer_ordinal>& data) {return FIELDCONTAINER_ORDINAL;}
130 template<> MuemexType getMuemexType<RCP<FieldContainer_ordinal>>() {return FIELDCONTAINER_ORDINAL;}
131 #endif
132 
133 /* "prototypes" for specialized functions used in other specialized functions */
134 
135 template<> mxArray* createMatlabSparse<double>(int numRows, int numCols, int nnz);
136 template<> mxArray* createMatlabSparse<complex_t>(int numRows, int numCols, int nnz);
137 template<> mxArray* createMatlabMultiVector<double>(int numRows, int numCols);
138 template<> mxArray* createMatlabMultiVector<complex_t>(int numRows, int numCols);
139 template<> void fillMatlabArray<double>(double* array, const mxArray* mxa, int n);
140 template<> void fillMatlabArray<complex_t>(complex_t* array, const mxArray* mxa, int n);
141 template<> mxArray* saveDataToMatlab(RCP<Xpetra_MultiVector_double>& data);
142 template<> mxArray* saveDataToMatlab(RCP<Xpetra_MultiVector_complex>& data);
143 template<> mxArray* saveDataToMatlab(RCP<Xpetra_Matrix_double>& data);
144 template<> mxArray* saveDataToMatlab(RCP<Xpetra_Matrix_complex>& data);
145 
146 /* ******************************* */
147 /* loadDataFromMatlab */
148 /* ******************************* */
149 
150 template<>
151 int loadDataFromMatlab<int>(const mxArray* mxa)
152 {
153  mxClassID probIDtype = mxGetClassID(mxa);
154  int rv;
155  if(probIDtype == mxINT32_CLASS)
156  {
157  rv = *((int*) mxGetData(mxa));
158  }
159  else if(probIDtype == mxLOGICAL_CLASS)
160  {
161  rv = (int) *((bool*) mxGetData(mxa));
162  }
163  else if(probIDtype == mxDOUBLE_CLASS)
164  {
165  rv = (int) *((double*) mxGetData(mxa));
166  }
167  else if(probIDtype == mxUINT32_CLASS)
168  {
169  rv = (int) *((unsigned int*) mxGetData(mxa));
170  }
171  else
172  {
173  rv = -1;
174  throw std::runtime_error("Error: Unrecognized numerical type.");
175  }
176  return rv;
177 }
178 
179 template<>
180 bool loadDataFromMatlab<bool>(const mxArray* mxa)
181 {
182  return *((bool*) mxGetData(mxa));
183 }
184 
185 template<>
186 double loadDataFromMatlab<double>(const mxArray* mxa)
187 {
188  return *((double*) mxGetPr(mxa));
189 }
190 
191 template<>
193 {
194  double realpart = real<double>(*((double*) mxGetPr(mxa)));
195  double imagpart = imag<double>(*((double*) mxGetPi(mxa)));
196  return complex_t(realpart, imagpart);
197 }
198 
199 template<>
200 string loadDataFromMatlab<string>(const mxArray* mxa)
201 {
202  string rv = "";
203  if (mxGetClassID(mxa) != mxCHAR_CLASS)
204  {
205  throw runtime_error("Can't construct string from anything but a char array.");
206  }
207  rv = string(mxArrayToString(mxa));
208  return rv;
209 }
210 
211 template<>
212 RCP<Xpetra_map> loadDataFromMatlab<RCP<Xpetra_map>>(const mxArray* mxa)
213 {
215  int nr = mxGetM(mxa);
216  int nc = mxGetN(mxa);
217  if(nr != 1)
218  throw std::runtime_error("A Xpetra::Map representation from MATLAB must be a single row vector.");
219  double* pr = mxGetPr(mxa);
220  mm_GlobalOrd numGlobalIndices = nc;
221 
222  std::vector<mm_GlobalOrd> localGIDs(numGlobalIndices);
223  for(int i = 0; i < int(numGlobalIndices); i++) {
224  localGIDs[i] = Teuchos::as<mm_GlobalOrd>(pr[i]);
225  }
226 
227  const Teuchos::ArrayView<const mm_GlobalOrd> localGIDs_view(&localGIDs[0],localGIDs.size());
228  RCP<Xpetra_map> map =
232  localGIDs_view,
233  0, comm);
234 
235  if(map.is_null())
236  throw runtime_error("Failed to create Xpetra::Map.");
237  return map;
238 }
239 
240 template<>
241 RCP<Xpetra_ordinal_vector> loadDataFromMatlab<RCP<Xpetra_ordinal_vector>>(const mxArray* mxa)
242 {
244  if(mxGetN(mxa) != 1 && mxGetM(mxa) != 1)
245  throw std::runtime_error("An OrdinalVector from MATLAB must be a single row or column vector.");
246  mm_GlobalOrd numGlobalIndices = mxGetM(mxa) * mxGetN(mxa);
248  if(mxGetClassID(mxa) != mxINT32_CLASS)
249  throw std::runtime_error("Can only construct LOVector with int32 data.");
250  int* array = (int*) mxGetData(mxa);
251  if(map.is_null())
252  throw runtime_error("Failed to create map for Xpetra ordinal vector.");
254  if(loVec.is_null())
255  throw runtime_error("Failed to create ordinal vector with Xpetra::VectorFactory.");
256  for(int i = 0; i < int(numGlobalIndices); i++)
257  {
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)
265 {
267  try
268  {
269  int nr = mxGetM(mxa);
270  int nc = mxGetN(mxa);
271  double* pr = mxGetPr(mxa);
272  RCP<const Teuchos::Comm<int>> comm = Tpetra::getDefaultComm();
273  //numGlobalIndices for map constructor is the number of rows in matrix/vectors, right?
275  //Allocate a new array of complex values to use with the multivector
276  Teuchos::ArrayView<const double> arrView(pr, nr * nc);
277  mv = rcp(new Tpetra::MultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(map, arrView, size_t(nr), size_t(nc)));
278  }
279  catch(std::exception& e)
280  {
281  mexPrintf("Error constructing Tpetra MultiVector.\n");
282  std::cout << e.what() << std::endl;
283  }
284  return mv;
285 }
286 
287 template<>
288 RCP<Tpetra_MultiVector_complex> loadDataFromMatlab<RCP<Tpetra_MultiVector_complex>>(const mxArray* mxa)
289 {
291  try
292  {
293  int nr = mxGetM(mxa);
294  int nc = mxGetN(mxa);
295  double* pr = mxGetPr(mxa);
296  double* pi = mxGetPi(mxa);
297  RCP<const Teuchos::Comm<int>> comm = Tpetra::getDefaultComm();
298  //numGlobalIndices for map constructor is the number of rows in matrix/vectors, right?
300  //Allocate a new array of complex values to use with the multivector
301  complex_t* myArr = new complex_t[nr * nc];
302  for(int n = 0; n < nc; n++)
303  {
304  for(int m = 0; m < nr; m++)
305  {
306  myArr[n * nr + m] = complex_t(pr[n * nr + m], pi[n * nr + m]);
307  }
308  }
309  Teuchos::ArrayView<complex_t> arrView(myArr, nr * nc);
310  mv = rcp(new Tpetra::MultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(map, arrView, nr, nc));
311  }
312  catch(std::exception& e)
313  {
314  mexPrintf("Error constructing Tpetra MultiVector.\n");
315  std::cout << e.what() << std::endl;
316  }
317  return mv;
318 }
319 
320 template<>
321 RCP<Tpetra_CrsMatrix_double> loadDataFromMatlab<RCP<Tpetra_CrsMatrix_double>>(const mxArray* mxa)
322 {
323  bool success = false;
325 
326  int* colptr = NULL;
327  int* rowind = NULL;
328 
329  try
330  {
332  //numGlobalIndices is just the number of rows in the matrix
333  const Tpetra::global_size_t numGlobalIndices = mxGetM(mxa);
334  const mm_GlobalOrd indexBase = 0;
335  RCP<const muemex_map_type> rowMap = rcp(new muemex_map_type(numGlobalIndices, indexBase, comm));
336  RCP<const muemex_map_type> domainMap = rcp(new muemex_map_type(mxGetN(mxa), indexBase, comm));
337  A = Tpetra::createCrsMatrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(rowMap);
338  double* valueArray = mxGetPr(mxa);
339  int nc = mxGetN(mxa);
340  if(rewrap_ints)
341  {
342  //mwIndex_to_int allocates memory so must delete[] later
343  colptr = mwIndex_to_int(nc + 1, mxGetJc(mxa));
344  rowind = mwIndex_to_int(colptr[nc], mxGetIr(mxa));
345  }
346  else
347  {
348  rowind = (int*) mxGetIr(mxa);
349  colptr = (int*) mxGetJc(mxa);
350  }
351  for(int i = 0; i < nc; i++)
352  {
353  for(int j = colptr[i]; j < colptr[i + 1]; j++)
354  {
355  //'array' of 1 element, containing column (in global matrix).
357  //'array' of 1 element, containing value
358  Teuchos::ArrayView<double> vals = Teuchos::ArrayView<double>(&valueArray[j], 1);
359  A->insertGlobalValues(rowind[j], cols, vals);
360  }
361  }
362  A->fillComplete(domainMap, rowMap);
363  if(rewrap_ints)
364  {
365  delete[] rowind; rowind = NULL;
366  delete[] colptr; colptr = NULL;
367  }
368  success = true;
369  }
370  catch(std::exception& e)
371  {
372  if(rewrap_ints)
373  {
374  if(rowind!=NULL) delete[] rowind;
375  if(colptr!=NULL) delete[] colptr;
376  rowind = NULL;
377  colptr = NULL;
378  }
379  mexPrintf("Error while constructing Tpetra matrix:\n");
380  std::cout << e.what() << std::endl;
381  }
382  if(!success)
383  mexErrMsgTxt("An error occurred while setting up a Tpetra matrix.\n");
384  return A;
385 }
386 
387 template<>
388 RCP<Tpetra_CrsMatrix_complex> loadDataFromMatlab<RCP<Tpetra_CrsMatrix_complex>>(const mxArray* mxa)
389 {
391  //Create a map in order to create the matrix (taken from muelu basic example - complex)
392  try
393  {
394  RCP<const Teuchos::Comm<int>> comm = Tpetra::getDefaultComm();
395  const Tpetra::global_size_t numGlobalIndices = mxGetM(mxa);
396  const mm_GlobalOrd indexBase = 0;
397  RCP<const muemex_map_type> rowMap = rcp(new muemex_map_type(numGlobalIndices, indexBase, comm));
398  RCP<const muemex_map_type> domainMap = rcp(new muemex_map_type(mxGetN(mxa), indexBase, comm));
399  A = Tpetra::createCrsMatrix<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(rowMap);
400  double* realArray = mxGetPr(mxa);
401  double* imagArray = mxGetPi(mxa);
402  int* colptr;
403  int* rowind;
404  int nc = mxGetN(mxa);
405  if(rewrap_ints)
406  {
407  //mwIndex_to_int allocates memory so must delete[] later
408  colptr = mwIndex_to_int(nc + 1, mxGetJc(mxa));
409  rowind = mwIndex_to_int(colptr[nc], mxGetIr(mxa));
410  }
411  else
412  {
413  rowind = (int*) mxGetIr(mxa);
414  colptr = (int*) mxGetJc(mxa);
415  }
416  for(int i = 0; i < nc; i++)
417  {
418  for(int j = colptr[i]; j < colptr[i + 1]; j++)
419  {
420  //here assuming that complex_t will always be defined as std::complex<double>
421  //use 'value' over and over again with Teuchos::ArrayViews to insert into matrix
422  complex_t value = std::complex<double>(realArray[j], imagArray[j]);
425  A->insertGlobalValues(rowind[j], cols, vals);
426  }
427  }
428  A->fillComplete(domainMap, rowMap);
429  if(rewrap_ints)
430  {
431  delete[] rowind;
432  delete[] colptr;
433  }
434  }
435  catch(std::exception& e)
436  {
437  mexPrintf("Error while constructing tpetra matrix:\n");
438  std::cout << e.what() << std::endl;
439  }
440  return A;
441 }
442 
443 template<>
444 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)
445 {
446  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);
447  return MueLu::TpetraCrs_To_XpetraMatrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tmat);
448 }
449 
450 template<>
451 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)
452 {
453  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);
454  return MueLu::TpetraCrs_To_XpetraMatrix<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tmat);
455 }
456 
457 template<>
458 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)
459 {
460  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);
461  return MueLu::TpetraMultiVector_To_XpetraMultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tpetraMV);
462 }
463 
464 template<>
465 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)
466 {
467  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);
468  return MueLu::TpetraMultiVector_To_XpetraMultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tpetraMV);
469 }
470 
471 template<>
472 RCP<Epetra_CrsMatrix> loadDataFromMatlab<RCP<Epetra_CrsMatrix>>(const mxArray* mxa)
473 {
474  RCP<Epetra_CrsMatrix> matrix;
475  try
476  {
477  int* colptr;
478  int* rowind;
479  double* vals = mxGetPr(mxa);
480  int nr = mxGetM(mxa);
481  int nc = mxGetN(mxa);
482  if(rewrap_ints)
483  {
484  colptr = mwIndex_to_int(nc + 1, mxGetJc(mxa));
485  rowind = mwIndex_to_int(colptr[nc], mxGetIr(mxa));
486  }
487  else
488  {
489  rowind = (int*) mxGetIr(mxa);
490  colptr = (int*) mxGetJc(mxa);
491  }
492  Epetra_SerialComm Comm;
493  Epetra_Map RangeMap(nr, 0, Comm);
494  Epetra_Map DomainMap(nc, 0, Comm);
495  matrix = rcp(new Epetra_CrsMatrix(Epetra_DataAccess::Copy, RangeMap, DomainMap, 0));
496  /* Do the matrix assembly */
497  for(int i = 0; i < nc; i++)
498  {
499  for(int j = colptr[i]; j < colptr[i + 1]; j++)
500  {
501  //global row, # of entries, value array, column indices array
502  matrix->InsertGlobalValues(rowind[j], 1, &vals[j], &i);
503  }
504  }
505  matrix->FillComplete(DomainMap, RangeMap);
506  if(rewrap_ints)
507  {
508  delete [] rowind;
509  delete [] colptr;
510  }
511  }
512  catch(std::exception& e)
513  {
514  mexPrintf("An error occurred while setting up an Epetra matrix:\n");
515  std::cout << e.what() << std::endl;
516  }
517  return matrix;
518 }
519 
520 template<>
521 RCP<Epetra_MultiVector> loadDataFromMatlab<RCP<Epetra_MultiVector>>(const mxArray* mxa)
522 {
523  int nr = mxGetM(mxa);
524  int nc = mxGetN(mxa);
525  Epetra_SerialComm Comm;
526  Epetra_BlockMap map(nr * nc, 1, 0, Comm);
527  return rcp(new Epetra_MultiVector(Epetra_DataAccess::Copy, map, mxGetPr(mxa), nr, nc));
528 }
529 
530 template<>
531 RCP<MAggregates> loadDataFromMatlab<RCP<MAggregates>>(const mxArray* mxa)
532 {
533  if(mxGetNumberOfElements(mxa) != 1)
534  throw runtime_error("Aggregates must be individual structs in MATLAB.");
535  if(!mxIsStruct(mxa))
536  throw runtime_error("Trying to pull aggregates from non-struct MATLAB object.");
537  //assume that in matlab aggregate structs will only be stored in a 1x1 array
538  //mxa must have the same fields as the ones declared in constructAggregates function in muelu.m for this to work
539  const int correctNumFields = 5; //change if more fields are added to the aggregates representation in constructAggregates in muelu.m
540  if(mxGetNumberOfFields(mxa) != correctNumFields)
541  throw runtime_error("Aggregates structure has wrong number of fields.");
542  //Pull MuemexData types back out
543  int nVert = *(int*) mxGetData(mxGetField(mxa, 0, "nVertices"));
544  int nAgg = *(int*) mxGetData(mxGetField(mxa, 0, "nAggregates"));
545  //Now have all the data needed to fully reconstruct the aggregate
546  //Use similar approach as UserAggregationFactory (which is written for >1 thread but will just be serial here)
548  int myRank = comm->getRank();
551  RCP<MAggregates> agg = rcp(new MAggregates(map));
552  agg->SetNumAggregates(nAgg);
553  //Get handles for the vertex2AggId and procwinner arrays in reconstituted aggregates object
554  //this is serial so all procwinner values will be same (0)
555  ArrayRCP<mm_LocalOrd> vertex2AggId = agg->GetVertex2AggId()->getDataNonConst(0); //the '0' means first (and only) column of multivector, since is just vector
556  ArrayRCP<mm_LocalOrd> procWinner = agg->GetProcWinner()->getDataNonConst(0);
557  //mm_LocalOrd and int are equivalent, so is ok to talk about aggSize with just 'int'
558  //Deep copy the entire vertex2AggID and isRoot arrays, which are both nVert items long
559  //At the same time, set ProcWinner
560  mxArray* vertToAggID_in = mxGetField(mxa, 0, "vertexToAggID");
561  int* vertToAggID_inArray = (int*) mxGetData(vertToAggID_in);
562  mxArray* rootNodes_in = mxGetField(mxa, 0, "rootNodes");
563  int* rootNodes_inArray = (int*) mxGetData(rootNodes_in);
564  for(int i = 0; i < nVert; i++)
565  {
566  vertex2AggId[i] = vertToAggID_inArray[i];
567  procWinner[i] = myRank; //all nodes are going to be on the same proc
568  agg->SetIsRoot(i, false); //the ones that are root will be set in next loop
569  }
570  for(int i = 0; i < nAgg; i++) //rootNodesToCopy is an array of node IDs which are the roots of their aggs
571  {
572  agg->SetIsRoot(rootNodes_inArray[i], true);
573  }
574  //Now recompute the aggSize array the results in the object
575  agg->ComputeAggregateSizes(true);
576  agg->AggregatesCrossProcessors(false);
577  return agg;
578 }
579 
580 template<>
581 RCP<MAmalInfo> loadDataFromMatlab<RCP<MAmalInfo>>(const mxArray* mxa)
582 {
583  RCP<MAmalInfo> amal;
584  throw runtime_error("AmalgamationInfo not supported in Muemex yet.");
585  return amal;
586 }
587 
588 template<>
589 RCP<MGraph> loadDataFromMatlab<RCP<MGraph>>(const mxArray* mxa)
590 {
591  //mxa must be struct with logical sparse matrix called 'edges' and Nx1 int32 array 'boundaryNodes'
592  mxArray* edges = mxGetField(mxa, 0, "edges");
593  mxArray* boundaryNodes = mxGetField(mxa, 0, "boundaryNodes");
594  if(edges == NULL)
595  throw runtime_error("Graph structure in MATLAB must have a field called 'edges' (logical sparse matrix)");
596  if(boundaryNodes == NULL)
597  throw runtime_error("Graph structure in MATLAB must have a field called 'boundaryNodes' (int32 array containing list of boundary nodes)");
598  int* boundaryList = (int*) mxGetData(boundaryNodes);
599  if(!mxIsSparse(edges) || mxGetClassID(edges) != mxLOGICAL_CLASS)
600  throw runtime_error("Graph edges must be stored as a logical sparse matrix.");
601  // Note that Matlab stores sparse matrices in column major format.
602  mwIndex* matlabColPtrs = mxGetJc(edges);
603  mwIndex* matlabRowIndices = mxGetIr(edges);
604  mm_GlobalOrd nRows = (mm_GlobalOrd) mxGetM(edges);
605 
606  // Create and populate row-major CRS data structures for Xpetra::TpetraCrsGraph.
607 
608  // calculate number of nonzeros in each row
609  Teuchos::Array<int> entriesPerRow(nRows);
610  int nnz = matlabColPtrs[mxGetN(edges)]; //last entry in matlabColPtrs
611  for(int i = 0; i < nnz; i++)
612  entriesPerRow[matlabRowIndices[i]]++;
613  // Populate usual row index array. We don't need this for the Xpetra Graph ctor, but
614  // it's convenient for building up the column index array, which the ctor does need.
615  Teuchos::Array<int> rows(nRows+1);
616  rows[0] = 0;
617  for(int i = 0; i < nRows; i++)
618  rows[i+1] = rows[i] + entriesPerRow[i];
619  Teuchos::Array<int> cols(nnz); //column index array
620  Teuchos::Array<int> insertionsPerRow(nRows,0); //track of #insertions done per row
621  int ncols = mxGetN(edges);
622  for (int colNum=0; colNum<ncols; ++colNum) {
623  int ci = matlabColPtrs[colNum];
624  for (int j=ci; j<Teuchos::as<int>(matlabColPtrs[colNum+1]); ++j) {
625  int rowNum = matlabRowIndices[j];
626  cols[ rows[rowNum] + insertionsPerRow[rowNum] ] = colNum;
627  insertionsPerRow[rowNum]++;
628  }
629  }
630  //Find maximum
631  int maxNzPerRow = 0;
632  for(int i = 0; i < nRows; i++) {
633  if(maxNzPerRow < entriesPerRow[i])
634  maxNzPerRow = entriesPerRow[i];
635  }
636 
639  RCP<MMap> map = rcp(new MMap(nRows, 0, comm));
641  RCP<TpetraGraph> tgraph = rcp(new TpetraGraph(map, (size_t) maxNzPerRow));
642  //Populate tgraph in compressed-row format. Must get each row individually...
643  for(int i = 0; i < nRows; ++i) {
644  tgraph->insertGlobalIndices((mm_GlobalOrd) i, cols(rows[i],entriesPerRow[i]));
645  }
646  tgraph->fillComplete(map, map);
648  //Set boundary nodes
649  int numBoundaryNodes = mxGetNumberOfElements(boundaryNodes);
650  bool* boundaryFlags = new bool[nRows];
651  for(int i = 0; i < nRows; i++)
652  {
653  boundaryFlags[i] = false;
654  }
655  for(int i = 0; i < numBoundaryNodes; i++)
656  {
657  boundaryFlags[boundaryList[i]] = true;
658  }
659  ArrayRCP<bool> boundaryNodesInput(boundaryFlags, 0, nRows, true);
660  mgraph->SetBoundaryNodeMap(boundaryNodesInput);
661  return mgraph;
662 }
663 
664 
665 #ifdef HAVE_MUELU_INTREPID2
666 template<>
667 RCP<FieldContainer_ordinal> loadDataFromMatlab<RCP<FieldContainer_ordinal>>(const mxArray* mxa)
668 {
669  if(mxGetClassID(mxa) != mxINT32_CLASS)
670  throw runtime_error("FieldContainer must have integer storage entries");
671 
672  int *data = (int *) mxGetData(mxa);
673  int nr = mxGetM(mxa);
674  int nc = mxGetN(mxa);
675 
676  RCP<FieldContainer_ordinal> fc = rcp(new FieldContainer_ordinal("FC from Matlab",nr,nc));
677  for(int col = 0; col < nc; col++)
678  {
679  for(int row = 0; row < nr; row++)
680  {
681  (*fc)(row,col) = data[col * nr + row];
682  }
683  }
684  return fc;
685 }
686 #endif
687 
688 /* ******************************* */
689 /* saveDataToMatlab */
690 /* ******************************* */
691 
692 template<>
693 mxArray* saveDataToMatlab(int& data)
694 {
695  mwSize dims[] = {1, 1};
696  mxArray* mxa = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
697  *((int*) mxGetData(mxa)) = data;
698  return mxa;
699 }
700 
701 template<>
702 mxArray* saveDataToMatlab(bool& data)
703 {
704  mwSize dims[] = {1, 1};
705  mxArray* mxa = mxCreateLogicalArray(2, dims);
706  *((bool*) mxGetData(mxa)) = data;
707  return mxa;
708 }
709 
710 template<>
711 mxArray* saveDataToMatlab(double& data)
712 {
713  return mxCreateDoubleScalar(data);
714 }
715 
716 template<>
718 {
719  mwSize dims[] = {1, 1};
720  mxArray* mxa = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxCOMPLEX);
721  *((double*) mxGetPr(mxa)) = real<double>(data);
722  *((double*) mxGetPi(mxa)) = imag<double>(data);
723  return mxa;
724 }
725 
726 template<>
727 mxArray* saveDataToMatlab(string& data)
728 {
729  return mxCreateString(data.c_str());
730 }
731 
732 template<>
734 {
735  //Precondition: Memory has already been allocated by MATLAB for the array.
736  int nc = data->getGlobalNumElements();
737  int nr = 1;
738  mxArray* output = createMatlabMultiVector<double>(nr, nc);
739  double* array = (double*) malloc(sizeof(double) * nr * nc);
740  for(int col = 0; col < nc; col++)
741  {
742  mm_GlobalOrd gid = data->getGlobalElement(col);
743  array[col] = Teuchos::as<double>(gid);
744  }
745  fillMatlabArray<double>(array, output, nc * nr);
746  free(array);
747  return output;
748 }
749 
750 template<>
752 {
753  mwSize len = data->getGlobalLength();
754  //create a single column vector
755  mwSize dimensions[] = {len, 1};
756  mxArray* rv = mxCreateNumericArray(2, dimensions, mxINT32_CLASS, mxREAL);
757  int* dataPtr = (int*) mxGetData(rv);
758  ArrayRCP<const mm_LocalOrd> arr = data->getData(0);
759  for(int i = 0; i < int(data->getGlobalLength()); i++)
760  {
761  dataPtr[i] = arr[i];
762  }
763  return rv;
764 }
765 
766 template<>
767 mxArray* saveDataToMatlab(RCP<Tpetra::MultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>>& data)
768 {
770  return saveDataToMatlab(xmv);
771 }
772 
773 template<>
774 mxArray* saveDataToMatlab(RCP<Tpetra::MultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>>& data)
775 {
777  return saveDataToMatlab(xmv);
778 }
779 
780 template<>
782 {
784  return saveDataToMatlab(xmat);
785 }
786 
787 template<>
789 {
791  return saveDataToMatlab(xmat);
792 }
793 
794 template<>
796 {
797  typedef double Scalar;
798  // Compute global constants, if we need them
799  Teuchos::rcp_const_cast<Xpetra_CrsGraph>(data->getCrsGraph())->computeGlobalConstants();
800 
801  int nr = data->getGlobalNumRows();
802  int nc = data->getGlobalNumCols();
803  int nnz = data->getGlobalNumEntries();
804 
805 #ifdef VERBOSE_OUTPUT
806  RCP<Teuchos::FancyOStream> fancyStream = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
807  mat->describe(*fancyStream, Teuchos::VERB_EXTREME);
808 #endif
809  mxArray* mxa = createMatlabSparse<Scalar>(nr, nc, nnz);
810  mwIndex* ir = mxGetIr(mxa);
811  mwIndex* jc = mxGetJc(mxa);
812  for(int i = 0; i < nc + 1; i++)
813  {
814  jc[i] = 0;
815  }
816 
817  size_t maxEntriesPerRow = data->getGlobalMaxNumRowEntries();
818  if(maxEntriesPerRow == Teuchos::OrdinalTraits<size_t>::invalid() || maxEntriesPerRow == 0) maxEntriesPerRow = data->getNodeMaxNumRowEntries();
819 
820  int* rowProgress = new int[nc];
821  //The array that will be copied to Pr and (if complex) Pi later
822  Scalar* sparseVals = new Scalar[nnz];
823  size_t numEntries;
824  if(data->isLocallyIndexed())
825  {
826  Scalar* rowValArray = new Scalar[maxEntriesPerRow];
827  Teuchos::ArrayView<Scalar> rowVals(rowValArray, maxEntriesPerRow);
828  mm_LocalOrd* rowIndicesArray = new mm_LocalOrd[maxEntriesPerRow];
829  Teuchos::ArrayView<mm_LocalOrd> rowIndices(rowIndicesArray, maxEntriesPerRow);
830  for(mm_LocalOrd m = 0; m < nr; m++) //All rows in the Xpetra matrix
831  {
832  data->getLocalRowCopy(m, rowIndices, rowVals, numEntries); //Get the row
833  for(mm_LocalOrd entry = 0; entry < int(numEntries); entry++) //All entries in row
834  {
835  jc[rowIndices[entry] + 1]++; //for each entry, increase jc for the entry's column
836  }
837  }
838 
839  //now jc holds the number of elements in each column, but needs cumulative sum over all previous columns also
840  int entriesAccum = 0;
841  for(int n = 0; n <= nc; n++)
842  {
843  int temp = entriesAccum;
844  entriesAccum += jc[n];
845  jc[n] += temp;
846  }
847  //Jc now populated with colptrs
848  for(int i = 0; i < nc; i++)
849  {
850  rowProgress[i] = 0;
851  }
852  //Row progress values like jc but keep track as the MATLAB matrix is being filled in
853  for(mm_LocalOrd m = 0; m < nr; m++) //rows
854  {
855  data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
856  for(mm_LocalOrd i = 0; i < int(numEntries); i++) //entries in row m (NOT columns)
857  {
858  //row is m, col is rowIndices[i], val is rowVals[i]
859  mm_LocalOrd col = rowIndices[i];
860  sparseVals[jc[col] + rowProgress[col]] = rowVals[i]; //Set value
861  ir[jc[col] + rowProgress[col]] = m; //Set row at which value occurs
862  rowProgress[col]++;
863  }
864  }
865  delete[] rowIndicesArray;
866  }
867  else
868  {
871  for(mm_GlobalOrd m = 0; m < nr; m++)
872  {
873  data->getGlobalRowView(m, rowIndices, rowVals);
874  for(mm_GlobalOrd n = 0; n < rowIndices.size(); n++)
875  {
876  jc[rowIndices[n] + 1]++;
877  }
878  }
879  //Last element of jc is just nnz
880  jc[nc] = nnz;
881  //Jc now populated with colptrs
882  for(int i = 0; i < nc; i++)
883  {
884  rowProgress[i] = 0;
885  }
886  int entriesAccum = 0;
887  for(int n = 0; n <= nc; n++)
888  {
889  int temp = entriesAccum;
890  entriesAccum += jc[n];
891  jc[n] += temp;
892  }
893  //Row progress values like jc but keep track as the MATLAB matrix is being filled in
894  for(mm_GlobalOrd m = 0; m < nr; m++) //rows
895  {
896  data->getGlobalRowView(m, rowIndices, rowVals);
897  for(mm_LocalOrd i = 0; i < rowIndices.size(); i++) //entries in row m
898  {
899  mm_GlobalOrd col = rowIndices[i]; //row is m, col is rowIndices[i], val is rowVals[i]
900  sparseVals[jc[col] + rowProgress[col]] = rowVals[i]; //Set value
901  ir[jc[col] + rowProgress[col]] = m; //Set row at which value occurs
902  rowProgress[col]++;
903  }
904  }
905  }
906  //finally, copy sparseVals into pr (and pi, if complex)
907  fillMatlabArray<Scalar>(sparseVals, mxa, nnz);
908  delete[] sparseVals;
909  delete[] rowProgress;
910  return mxa;
911 }
912 
913 template<>
915 {
916  typedef complex_t Scalar;
917 
918  // Compute global constants, if we need them
919  Teuchos::rcp_const_cast<Xpetra_CrsGraph>(data->getCrsGraph())->computeGlobalConstants();
920 
921  int nr = data->getGlobalNumRows();
922  int nc = data->getGlobalNumCols();
923  int nnz = data->getGlobalNumEntries();
924 #ifdef VERBOSE_OUTPUT
925  RCP<Teuchos::FancyOStream> fancyStream = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
926  mat->describe(*fancyStream, Teuchos::VERB_EXTREME);
927 #endif
928  mxArray* mxa = createMatlabSparse<Scalar>(nr, nc, nnz);
929  mwIndex* ir = mxGetIr(mxa);
930  mwIndex* jc = mxGetJc(mxa);
931  for(int i = 0; i < nc + 1; i++)
932  {
933  jc[i] = 0;
934  }
935  size_t maxEntriesPerRow = data->getGlobalMaxNumRowEntries();
936  int* rowProgress = new int[nc];
937  //The array that will be copied to Pr and (if complex) Pi later
938  Scalar* sparseVals = new Scalar[nnz];
939  size_t numEntries;
940  if(data->isLocallyIndexed())
941  {
942  Scalar* rowValArray = new Scalar[maxEntriesPerRow];
943  Teuchos::ArrayView<Scalar> rowVals(rowValArray, maxEntriesPerRow);
944  mm_LocalOrd* rowIndicesArray = new mm_LocalOrd[maxEntriesPerRow];
945  Teuchos::ArrayView<mm_LocalOrd> rowIndices(rowIndicesArray, maxEntriesPerRow);
946  for(mm_LocalOrd m = 0; m < nr; m++) //All rows in the Xpetra matrix
947  {
948  data->getLocalRowCopy(m, rowIndices, rowVals, numEntries); //Get the row
949  for(mm_LocalOrd entry = 0; entry < int(numEntries); entry++) //All entries in row
950  {
951  jc[rowIndices[entry] + 1]++; //for each entry, increase jc for the entry's column
952  }
953  }
954  //now jc holds the number of elements in each column, but needs cumulative sum over all previous columns also
955  int entriesAccum = 0;
956  for(int n = 0; n <= nc; n++)
957  {
958  int temp = entriesAccum;
959  entriesAccum += jc[n];
960  jc[n] += temp;
961  }
962  //Jc now populated with colptrs
963  for(int i = 0; i < nc; i++)
964  {
965  rowProgress[i] = 0;
966  }
967  //Row progress values like jc but keep track as the MATLAB matrix is being filled in
968  for(mm_LocalOrd m = 0; m < nr; m++) //rows
969  {
970  data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
971  for(mm_LocalOrd i = 0; i < int(numEntries); i++) //entries in row m (NOT columns)
972  {
973  //row is m, col is rowIndices[i], val is rowVals[i]
974  mm_LocalOrd col = rowIndices[i];
975  sparseVals[jc[col] + rowProgress[col]] = rowVals[i]; //Set value
976  ir[jc[col] + rowProgress[col]] = m; //Set row at which value occurs
977  rowProgress[col]++;
978  }
979  }
980  delete[] rowIndicesArray;
981  }
982  else
983  {
986  for(mm_GlobalOrd m = 0; m < nr; m++)
987  {
988  data->getGlobalRowView(m, rowIndices, rowVals);
989  for(mm_GlobalOrd n = 0; n < rowIndices.size(); n++)
990  {
991  jc[rowIndices[n] + 1]++;
992  }
993  }
994  //Last element of jc is just nnz
995  jc[nc] = nnz;
996  //Jc now populated with colptrs
997  for(int i = 0; i < nc; i++)
998  {
999  rowProgress[i] = 0;
1000  }
1001  int entriesAccum = 0;
1002  for(int n = 0; n <= nc; n++)
1003  {
1004  int temp = entriesAccum;
1005  entriesAccum += jc[n];
1006  jc[n] += temp;
1007  }
1008  //Row progress values like jc but keep track as the MATLAB matrix is being filled in
1009  for(mm_GlobalOrd m = 0; m < nr; m++) //rows
1010  {
1011  data->getGlobalRowView(m, rowIndices, rowVals);
1012  for(mm_LocalOrd i = 0; i < rowIndices.size(); i++) //entries in row m
1013  {
1014  mm_GlobalOrd col = rowIndices[i]; //row is m, col is rowIndices[i], val is rowVals[i]
1015  sparseVals[jc[col] + rowProgress[col]] = rowVals[i]; //Set value
1016  ir[jc[col] + rowProgress[col]] = m; //Set row at which value occurs
1017  rowProgress[col]++;
1018  }
1019  }
1020  }
1021  //finally, copy sparseVals into pr (and pi, if complex)
1022  fillMatlabArray<Scalar>(sparseVals, mxa, nnz);
1023  delete[] sparseVals;
1024  delete[] rowProgress;
1025  return mxa;
1026 }
1027 
1028 /*
1029 template<>
1030 mxArray* saveDataToMatlab(RCP<Xpetra::MultiVector<Scalar, mm_LocalOrd, mm_GlobalOrd, mm_node_t>>& data)
1031 {
1032  //Precondition: Memory has already been allocated by MATLAB for the array.
1033  int nr = data->getGlobalLength();
1034  int nc = data->getNumVectors();
1035  mxArray* output = createMatlabMultiVector<Scalar>(nr, nc);
1036  Scalar* array = (Scalar*) malloc(sizeof(Scalar) * nr * nc);
1037  for(int col = 0; col < nc; col++)
1038  {
1039  Teuchos::ArrayRCP<const Scalar> colData = data->getData(col);
1040  for(int row = 0; row < nr; row++)
1041  {
1042  array[col * nr + row] = colData[row];
1043  }
1044  }
1045  fillMatlabArray<Scalar>(array, output, nc * nr);
1046  free(array);
1047  return output;
1048 }
1049 */
1050 
1051 template<>
1053 {
1054  //Precondition: Memory has already been allocated by MATLAB for the array.
1055  int nr = data->getGlobalLength();
1056  int nc = data->getNumVectors();
1057  mxArray* output = createMatlabMultiVector<double>(nr, nc);
1058  double* array = (double*) malloc(sizeof(double) * nr * nc);
1059  for(int col = 0; col < nc; col++)
1060  {
1061  Teuchos::ArrayRCP<const double> colData = data->getData(col);
1062  for(int row = 0; row < nr; row++)
1063  {
1064  array[col * nr + row] = colData[row];
1065  }
1066  }
1067  fillMatlabArray<double>(array, output, nc * nr);
1068  free(array);
1069  return output;
1070 }
1071 
1072 template<>
1074 {
1075  //Precondition: Memory has already been allocated by MATLAB for the array.
1076  int nr = data->getGlobalLength();
1077  int nc = data->getNumVectors();
1078  mxArray* output = createMatlabMultiVector<complex_t>(nr, nc);
1079  complex_t* array = (complex_t*) malloc(sizeof(complex_t) * nr * nc);
1080  for(int col = 0; col < nc; col++)
1081  {
1082  Teuchos::ArrayRCP<const complex_t> colData = data->getData(col);
1083  for(int row = 0; row < nr; row++)
1084  {
1085  array[col * nr + row] = colData[row];
1086  }
1087  }
1088  fillMatlabArray<complex_t>(array, output, nc * nr);
1089  free(array);
1090  return output;
1091 }
1092 
1093 template<>
1095 {
1096  RCP<Xpetra_Matrix_double> xmat = EpetraCrs_To_XpetraMatrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(data);
1097  return saveDataToMatlab(xmat);
1098 }
1099 
1100 template<>
1102 {
1103  mxArray* output = mxCreateDoubleMatrix(data->GlobalLength(), data->NumVectors(), mxREAL);
1104  double* dataPtr = mxGetPr(output);
1105  data->ExtractCopy(dataPtr, data->GlobalLength());
1106  return output;
1107 }
1108 
1109 template<>
1111 {
1112  //Set up array of inputs for matlab constructAggregates
1113  int numNodes = data->GetVertex2AggId()->getData(0).size();
1114  int numAggs = data->GetNumAggregates();
1115  mxArray* dataIn[5];
1116  mwSize singleton[] = {1, 1};
1117  dataIn[0] = mxCreateNumericArray(2, singleton, mxINT32_CLASS, mxREAL);
1118  *((int*) mxGetData(dataIn[0])) = numNodes;
1119  dataIn[1] = mxCreateNumericArray(2, singleton, mxINT32_CLASS, mxREAL);
1120  *((int*) mxGetData(dataIn[1])) = numAggs;
1121  mwSize nodeArrayDims[] = {(mwSize) numNodes, 1}; //dimensions for Nx1 array, where N is number of nodes (vert2Agg)
1122  dataIn[2] = mxCreateNumericArray(2, nodeArrayDims, mxINT32_CLASS, mxREAL);
1123  int* vtaid = (int*) mxGetData(dataIn[2]);
1124  ArrayRCP<const mm_LocalOrd> vertexToAggID = data->GetVertex2AggId()->getData(0);
1125  for(int i = 0; i < numNodes; i++)
1126  {
1127  vtaid[i] = vertexToAggID[i];
1128  }
1129  mwSize aggArrayDims[] = {(mwSize) numAggs, 1}; //dims for Nx1 array, where N is number of aggregates (rootNodes, aggSizes)
1130  dataIn[3] = mxCreateNumericArray(2, aggArrayDims, mxINT32_CLASS, mxREAL);
1131  //First, find out if the aggregates even have 1 root node per aggregate. If not, assume roots are invalid and assign ourselves
1132  int totalRoots = 0;
1133  for(int i = 0; i < numNodes; i++)
1134  {
1135  if(data->IsRoot(i))
1136  totalRoots++;
1137  }
1138  bool reassignRoots = false;
1139  if(totalRoots != numAggs)
1140  {
1141  cout << endl << "Warning: Number of root nodes and number of aggregates do not match." << endl;
1142  cout << "Will reassign root nodes when writing aggregates to matlab." << endl << endl;
1143  reassignRoots = true;
1144  }
1145  int* rn = (int*) mxGetData(dataIn[3]); //list of root nodes (in no particular order)
1146  {
1147  if(reassignRoots)
1148  {
1149  //For each aggregate, just pick the first node we see in it and set it as root
1150  int lastFoundNode = 0; //heuristic for speed, a node in aggregate N+1 is likely to come very soon after a node in agg N
1151  for(int i = 0; i < numAggs; i++)
1152  {
1153  rn[i] = -1;
1154  for(int j = lastFoundNode; j < lastFoundNode + numNodes; j++)
1155  {
1156  int index = j % numNodes;
1157  if(vertexToAggID[index] == i)
1158  {
1159  rn[i] = index;
1160  lastFoundNode = index;
1161  }
1162  }
1163  TEUCHOS_TEST_FOR_EXCEPTION(rn[i] == -1, runtime_error, "Invalid aggregates: Couldn't find any node in aggregate #" << i << ".");
1164  }
1165  }
1166  else
1167  {
1168  int i = 0; //iterates over aggregate IDs
1169  for(int j = 0; j < numNodes; j++)
1170  {
1171  if(data->IsRoot(j))
1172  {
1173  if(i == numAggs)
1174  throw runtime_error("Cannot store invalid aggregates in MATLAB - more root nodes than aggregates.");
1175  rn[i] = j; //now we know this won't go out of bounds (rn's underlying matlab array is numAggs in length)
1176  i++;
1177  }
1178  }
1179  if(i + 1 < numAggs)
1180  throw runtime_error("Cannot store invalid aggregates in MATLAB - fewer root nodes than aggregates.");
1181  }
1182  }
1183  dataIn[4] = mxCreateNumericArray(1, aggArrayDims, mxINT32_CLASS, mxREAL);
1184  int* as = (int*) mxGetData(dataIn[4]); //list of aggregate sizes
1185  ArrayRCP<mm_LocalOrd> aggSizes = data->ComputeAggregateSizes();
1186  for(int i = 0; i < numAggs; i++)
1187  {
1188  as[i] = aggSizes[i];
1189  }
1190  mxArray* matlabAggs[1];
1191  int result = mexCallMATLAB(1, matlabAggs, 5, dataIn, "constructAggregates");
1192  if(result != 0)
1193  throw runtime_error("Matlab encountered an error while constructing aggregates struct.");
1194  return matlabAggs[0];
1195 }
1196 
1197 template<>
1199 {
1200  throw runtime_error("AmalgamationInfo not supported in MueMex yet.");
1201  return mxCreateDoubleScalar(0);
1202 }
1203 
1204 template<>
1206 {
1207  int numEntries = (int) data->GetGlobalNumEdges();
1208  int numRows = (int) data->GetDomainMap()->getGlobalNumElements(); //assume numRows == numCols
1209  mxArray* mat = mxCreateSparseLogicalMatrix(numRows, numRows, numEntries);
1210  mxLogical* outData = (mxLogical*) mxGetData(mat);
1211  mwIndex* rowInds = mxGetIr(mat);
1212  mwIndex* colPtrs = mxGetJc(mat);
1213  mm_LocalOrd* dataCopy = new mm_LocalOrd[numEntries];
1214  mm_LocalOrd* iter = dataCopy;
1215  int* entriesPerRow = new int[numRows];
1216  int* entriesPerCol = new int[numRows];
1217  for(int i = 0; i < numRows; i++)
1218  {
1219  entriesPerRow[i] = 0;
1220  entriesPerCol[i] = 0;
1221  }
1222  for(int i = 0; i < numRows; i++)
1223  {
1224  ArrayView<const mm_LocalOrd> neighbors = data->getNeighborVertices(i); //neighbors has the column indices for row i
1225  memcpy(iter, neighbors.getRawPtr(), sizeof(mm_LocalOrd) * neighbors.size());
1226  entriesPerRow[i] = neighbors.size();
1227  for(int j = 0; j < neighbors.size(); j++)
1228  {
1229  entriesPerCol[neighbors[j]]++;
1230  }
1231  iter += neighbors.size();
1232  }
1233  mwIndex** rowIndsByColumn = new mwIndex*[numRows]; //rowIndsByColumn[0] points to array of row indices in column 1
1234  mxLogical** valuesByColumn = new mxLogical*[numRows];
1235  int* numEnteredPerCol = new int[numRows];
1236  int accum = 0;
1237  for(int i = 0; i < numRows; i++)
1238  {
1239  rowIndsByColumn[i] = &rowInds[accum];
1240  //cout << "Entries in column " << i << " start at offset " << accum << endl;
1241  valuesByColumn[i] = &outData[accum];
1242  accum += entriesPerCol[i];
1243  if(accum > numEntries)
1244  throw runtime_error("potato");
1245  }
1246  for(int i = 0; i < numRows; i++)
1247  {
1248  numEnteredPerCol[i] = 0; //rowIndsByColumn[n][numEnteredPerCol[n]] gives the next place to put a row index
1249  }
1250  //entriesPerCol now has Jc information (col offsets)
1251  accum = 0; //keep track of cumulative index in dataCopy
1252  for(int row = 0; row < numRows; row++)
1253  {
1254  for(int entryInRow = 0; entryInRow < entriesPerRow[row]; entryInRow++)
1255  {
1256  int col = dataCopy[accum];
1257  accum++;
1258  rowIndsByColumn[col][numEnteredPerCol[col]] = row;
1259  valuesByColumn[col][numEnteredPerCol[col]] = (mxLogical) 1;
1260  numEnteredPerCol[col]++;
1261  }
1262  }
1263  accum = 0; //keep track of total entries over all columns
1264  for(int col = 0; col < numRows; col++)
1265  {
1266  colPtrs[col] = accum;
1267  accum += entriesPerCol[col];
1268  }
1269  colPtrs[numRows] = accum; //the last entry in jc, which is equivalent to numEntries
1270  delete[] numEnteredPerCol;
1271  delete[] rowIndsByColumn;
1272  delete[] valuesByColumn;
1273  delete[] dataCopy;
1274  delete[] entriesPerRow;
1275  delete[] entriesPerCol;
1276  //Construct list of boundary nodes
1277  const ArrayRCP<const bool> boundaryFlags = data->GetBoundaryNodeMap();
1278  int numBoundaryNodes = 0;
1279  for(int i = 0; i < boundaryFlags.size(); i++)
1280  {
1281  if(boundaryFlags[i])
1282  numBoundaryNodes++;
1283  }
1284  cout << "Graph has " << numBoundaryNodes << " Dirichlet boundary nodes." << endl;
1285  mwSize dims[] = {(mwSize) numBoundaryNodes, 1};
1286  mxArray* boundaryList = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
1287  int* dest = (int*) mxGetData(boundaryList);
1288  int* destIter = dest;
1289  for(int i = 0; i < boundaryFlags.size(); i++)
1290  {
1291  if(boundaryFlags[i])
1292  {
1293  *destIter = i;
1294  destIter++;
1295  }
1296  }
1297  mxArray* constructArgs[] = {mat, boundaryList};
1298  mxArray* out[1];
1299  mexCallMATLAB(1, out, 2, constructArgs, "constructGraph");
1300  return out[0];
1301 }
1302 
1303 #ifdef HAVE_MUELU_INTREPID2
1304 template<>
1306 {
1307  int rank = data->rank();
1308  // NOTE: Only supports rank 2 arrays
1309  if(rank!=2)
1310  throw std::runtime_error("Error: Only rank two FieldContainers are supported.");
1311 
1312  int nr = data->dimension(0);
1313  int nc = data->dimension(1);
1314 
1315  mwSize dims[]={(mwSize)nr,(mwSize)nc};
1316  mxArray* mxa = mxCreateNumericArray(2,dims, mxINT32_CLASS, mxREAL);
1317  int *array = (int*) mxGetData(mxa);
1318 
1319  for(int col = 0; col < nc; col++)
1320  {
1321  for(int row = 0; row < nr; row++)
1322  {
1323  array[col * nr + row] = (*data)(row,col);
1324  }
1325  }
1326  return mxa;
1327 }
1328 #endif
1329 
1330 
1331 template<typename T>
1333 {
1334  data = loadDataFromMatlab<T>(mxa);
1335 }
1336 
1337 template<typename T>
1339 {
1340  return saveDataToMatlab<T>(data);
1341 }
1342 
1343 template<typename T>
1344 MuemexData<T>::MuemexData(T& dataToCopy, MuemexType dataType) : MuemexArg(dataType)
1345 {
1346  data = dataToCopy;
1347 }
1348 
1349 template<typename T>
1350 MuemexData<T>::MuemexData(T& dataToCopy) : MuemexData(dataToCopy, getMuemexType(dataToCopy)) {}
1351 
1352 template<typename T>
1354 {
1355  return data;
1356 }
1357 
1358 template<typename T>
1359 void MuemexData<T>::setData(T& newData)
1360 {
1361  this->data = newData;
1362 }
1363 
1364 /* ***************************** */
1365 /* More Template Functions */
1366 /* ***************************** */
1367 
1368 template<typename T>
1369 void addLevelVariable(const T& data, std::string& name, Level& lvl, const Factory * fact)
1370 {
1371  lvl.AddKeepFlag(name, fact, MueLu::UserData);
1372  lvl.Set<T>(name, data, fact);
1373 }
1374 
1375 template<typename T>
1376 const T& getLevelVariable(std::string& name, Level& lvl)
1377 {
1378  try
1379  {
1380  return lvl.Get<T>(name);
1381  }
1382  catch(std::exception& e)
1383  {
1384  throw std::runtime_error("Requested custom variable " + name + " is not in the level.");
1385  }
1386 }
1387 
1388 //Functions used to put data through matlab factories - first arg is "this" pointer of matlab factory
1389 template<typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, typename Node>
1390 std::vector<Teuchos::RCP<MuemexArg>> processNeeds(const Factory* factory, std::string& needsParam, Level& lvl)
1391 {
1392  using namespace std;
1393  using namespace Teuchos;
1397  typedef RCP<AmalgamationInfo<LocalOrdinal, GlobalOrdinal, Node>> AmalgamationInfo_t;
1398  typedef RCP<MGraph> Graph_t;
1399  vector<string> needsList = tokenizeList(needsParam);
1400  vector<RCP<MuemexArg>> args;
1401  for(size_t i = 0; i < needsList.size(); i++)
1402  {
1403  if(needsList[i] == "A" || needsList[i] == "P" || needsList[i] == "R" || needsList[i]=="Ptent")
1404  {
1405  Matrix_t mydata = lvl.Get<Matrix_t>(needsList[i], factory->GetFactory(needsList[i]).get());
1406  args.push_back(rcp(new MuemexData<Matrix_t>(mydata)));
1407  }
1408  else if(needsList[i] == "Nullspace" || needsList[i] == "Coordinates")
1409  {
1410  MultiVector_t mydata = lvl.Get<MultiVector_t>(needsList[i], factory->GetFactory(needsList[i]).get());
1411  args.push_back(rcp(new MuemexData<MultiVector_t>(mydata)));
1412  }
1413  else if(needsList[i] == "Aggregates")
1414  {
1415  Aggregates_t mydata = lvl.Get<Aggregates_t>(needsList[i], factory->GetFactory(needsList[i]).get());
1416  args.push_back(rcp(new MuemexData<Aggregates_t>(mydata)));
1417  }
1418  else if(needsList[i] == "UnAmalgamationInfo")
1419  {
1420  AmalgamationInfo_t mydata = lvl.Get<AmalgamationInfo_t>(needsList[i], factory->GetFactory(needsList[i]).get());
1421  args.push_back(rcp(new MuemexData<AmalgamationInfo_t>(mydata)));
1422  }
1423  else if(needsList[i] == "Level")
1424  {
1425  int levelNum = lvl.GetLevelID();
1426  args.push_back(rcp(new MuemexData<int>(levelNum)));
1427  }
1428  else if(needsList[i] == "Graph")
1429  {
1430  Graph_t mydata = lvl.Get<Graph_t>(needsList[i], factory->GetFactory(needsList[i]).get());
1431  args.push_back(rcp(new MuemexData<Graph_t>(mydata)));
1432  }
1433  else
1434  {
1435  vector<string> words;
1436  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";
1437  //compare type without case sensitivity
1438  char* buf = (char*) malloc(needsList[i].size() + 1);
1439  strcpy(buf, needsList[i].c_str());
1440  for(char* iter = buf; *iter != ' '; iter++)
1441  {
1442  if(*iter == 0)
1443  {
1444  free(buf);
1445  throw runtime_error(badNameMsg);
1446  }
1447  *iter = (char) tolower(*iter);
1448  }
1449  const char* wordDelim = " ";
1450  char* mark = strtok(buf, wordDelim);
1451  while(mark != NULL)
1452  {
1453  string wordStr(mark);
1454  words.push_back(wordStr);
1455  mark = strtok(NULL, wordDelim);
1456  }
1457  if(words.size() != 2)
1458  {
1459  free(buf);
1460  throw runtime_error(badNameMsg);
1461  }
1462  char* typeStr = (char*) words[0].c_str();
1463  if(strstr(typeStr, "ordinalvector"))
1464  {
1466  LOVector_t mydata = getLevelVariable<LOVector_t>(needsList[i], lvl);
1467  args.push_back(rcp(new MuemexData<LOVector_t>(mydata)));
1468  }
1469  else if(strstr(typeStr, "map"))
1470  {
1472  Map_t mydata = getLevelVariable<Map_t>(needsList[i], lvl);
1473  args.push_back(rcp(new MuemexData<Map_t>(mydata)));
1474  }
1475  else if(strstr(typeStr, "scalar"))
1476  {
1477  Scalar mydata = getLevelVariable<Scalar>(needsList[i], lvl);
1478  args.push_back(rcp(new MuemexData<Scalar>(mydata)));
1479  }
1480  else if(strstr(typeStr, "double"))
1481  {
1482  double mydata = getLevelVariable<double>(needsList[i], lvl);
1483  args.push_back(rcp(new MuemexData<double>(mydata)));
1484  }
1485  else if(strstr(typeStr, "complex"))
1486  {
1487  complex_t mydata = getLevelVariable<complex_t>(needsList[i], lvl);
1488  args.push_back(rcp(new MuemexData<complex_t>(mydata)));
1489  }
1490  else if(strstr(typeStr, "matrix"))
1491  {
1492  Matrix_t mydata = getLevelVariable<Matrix_t>(needsList[i], lvl);
1493  args.push_back(rcp(new MuemexData<Matrix_t>(mydata)));
1494  }
1495  else if(strstr(typeStr, "multivector"))
1496  {
1497  MultiVector_t mydata = getLevelVariable<MultiVector_t>(needsList[i], lvl);
1498  args.push_back(rcp(new MuemexData<MultiVector_t>(mydata)));
1499  }
1500  else if(strstr(typeStr, "int"))
1501  {
1502  int mydata = getLevelVariable<int>(needsList[i], lvl);
1503  args.push_back(rcp(new MuemexData<int>(mydata)));
1504  }
1505  else if(strstr(typeStr, "string"))
1506  {
1507  string mydata = getLevelVariable<string>(needsList[i], lvl);
1508  args.push_back(rcp(new MuemexData<string>(mydata)));
1509  }
1510  else
1511  {
1512  free(buf);
1513  throw std::runtime_error(words[0] + " is not a known variable type.");
1514  }
1515  free(buf);
1516  }
1517  }
1518  return args;
1519 }
1520 
1521 template<typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, typename Node>
1522 void processProvides(std::vector<Teuchos::RCP<MuemexArg>>& mexOutput, const Factory* factory, std::string& providesParam, Level& lvl)
1523 {
1524  using namespace std;
1525  using namespace Teuchos;
1529  typedef RCP<AmalgamationInfo<LocalOrdinal, GlobalOrdinal, Node>> AmalgamationInfo_t;
1530  typedef RCP<MGraph> Graph_t;
1531  vector<string> provides = tokenizeList(providesParam);
1532  for(size_t i = 0; i < size_t(provides.size()); i++)
1533  {
1534  if(provides[i] == "A" || provides[i] == "P" || provides[i] == "R" || provides[i]=="Ptent")
1535  {
1536  RCP<MuemexData<Matrix_t>> mydata = Teuchos::rcp_static_cast<MuemexData<Matrix_t>>(mexOutput[i]);
1537  lvl.Set(provides[i], mydata->getData(), factory);
1538  }
1539  else if(provides[i] == "Nullspace" || provides[i] == "Coordinates")
1540  {
1541  RCP<MuemexData<MultiVector_t>> mydata = Teuchos::rcp_static_cast<MuemexData<MultiVector_t>>(mexOutput[i]);
1542  lvl.Set(provides[i], mydata->getData(), factory);
1543  }
1544  else if(provides[i] == "Aggregates")
1545  {
1546  RCP<MuemexData<Aggregates_t>> mydata = Teuchos::rcp_static_cast<MuemexData<Aggregates_t>>(mexOutput[i]);
1547  lvl.Set(provides[i], mydata->getData(), factory);
1548  }
1549  else if(provides[i] == "UnAmalgamationInfo")
1550  {
1551  RCP<MuemexData<AmalgamationInfo_t>> mydata = Teuchos::rcp_static_cast<MuemexData<AmalgamationInfo_t>>(mexOutput[i]);
1552  lvl.Set(provides[i], mydata->getData(), factory);
1553  }
1554  else if(provides[i] == "Graph")
1555  {
1556  RCP<MuemexData<Graph_t>> mydata = Teuchos::rcp_static_cast<MuemexData<Graph_t>>(mexOutput[i]);
1557  lvl.Set(provides[i], mydata->getData(), factory);
1558  }
1559  else
1560  {
1561  vector<string> words;
1562  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";
1563  //compare type without case sensitivity
1564  char* buf = (char*) malloc(provides[i].size() + 1);
1565  strcpy(buf, provides[i].c_str());
1566  for(char* iter = buf; *iter != ' '; iter++)
1567  {
1568  if(*iter == 0)
1569  {
1570  free(buf);
1571  throw runtime_error(badNameMsg);
1572  }
1573  *iter = (char) tolower(*iter);
1574  }
1575  const char* wordDelim = " ";
1576  char* mark = strtok(buf, wordDelim);
1577  while(mark != NULL)
1578  {
1579  string wordStr(mark);
1580  words.push_back(wordStr);
1581  mark = strtok(NULL, wordDelim);
1582  }
1583  if(words.size() != 2)
1584  {
1585  free(buf);
1586  throw runtime_error(badNameMsg);
1587  }
1588  const char* typeStr = words[0].c_str();
1589  if(strstr(typeStr, "ordinalvector"))
1590  {
1592  RCP<MuemexData<LOVector_t>> mydata = Teuchos::rcp_static_cast<MuemexData<LOVector_t>>(mexOutput[i]);
1593  addLevelVariable<LOVector_t>(mydata->getData(), words[1], lvl, factory);
1594  }
1595  else if(strstr(typeStr, "map"))
1596  {
1598  RCP<MuemexData<Map_t>> mydata = Teuchos::rcp_static_cast<MuemexData<Map_t>>(mexOutput[i]);
1599  addLevelVariable<Map_t>(mydata->getData(), words[1], lvl, factory);
1600  }
1601  else if(strstr(typeStr, "scalar"))
1602  {
1603  RCP<MuemexData<Scalar>> mydata = Teuchos::rcp_static_cast<MuemexData<Scalar>>(mexOutput[i]);
1604  addLevelVariable<Scalar>(mydata->getData(), words[1], lvl, factory);
1605  }
1606  else if(strstr(typeStr, "double"))
1607  {
1608  RCP<MuemexData<double>> mydata = Teuchos::rcp_static_cast<MuemexData<double>>(mexOutput[i]);
1609  addLevelVariable<double>(mydata->getData(), words[1], lvl, factory);
1610  }
1611  else if(strstr(typeStr, "complex"))
1612  {
1613  RCP<MuemexData<complex_t>> mydata = Teuchos::rcp_static_cast<MuemexData<complex_t>>(mexOutput[i]);
1614  addLevelVariable<complex_t>(mydata->getData(), words[1], lvl, factory);
1615  }
1616  else if(strstr(typeStr, "matrix"))
1617  {
1618  RCP<MuemexData<Matrix_t>> mydata = Teuchos::rcp_static_cast<MuemexData<Matrix_t>>(mexOutput[i]);
1619  addLevelVariable<Matrix_t>(mydata->getData(), words[1], lvl, factory);
1620  }
1621  else if(strstr(typeStr, "multivector"))
1622  {
1623  RCP<MuemexData<MultiVector_t>> mydata = Teuchos::rcp_static_cast<MuemexData<MultiVector_t>>(mexOutput[i]);
1624  addLevelVariable<MultiVector_t>(mydata->getData(), words[1], lvl, factory);
1625  }
1626  else if(strstr(typeStr, "int"))
1627  {
1628  RCP<MuemexData<int>> mydata = Teuchos::rcp_static_cast<MuemexData<int>>(mexOutput[i]);
1629  addLevelVariable<int>(mydata->getData(), words[1], lvl, factory);
1630  }
1631  else if(strstr(typeStr, "bool"))
1632  {
1633  RCP<MuemexData<bool> > mydata = Teuchos::rcp_static_cast<MuemexData<bool> >(mexOutput[i]);
1634  addLevelVariable<bool>(mydata->getData(), words[1], lvl, factory);
1635  }
1636  else if(strstr(typeStr, "string"))
1637  {
1638  RCP<MuemexData<string>> mydata = Teuchos::rcp_static_cast<MuemexData<string>>(mexOutput[i]);
1639  addLevelVariable<string>(mydata->getData(), words[1], lvl, factory);
1640  }
1641  else
1642  {
1643  free(buf);
1644  throw std::runtime_error(words[0] + " is not a known variable type.");
1645  }
1646  free(buf);
1647  }
1648  }
1649 }
1650 
1651 // Throwable Stubs for long long
1652 
1653 template<>
1654 std::vector<Teuchos::RCP<MuemexArg>> processNeeds<double,mm_LocalOrd,long long,mm_node_t>(const Factory* factory, std::string& needsParam, Level& lvl) {
1655  throw std::runtime_error("Muemex does not support long long for global indices");
1656 }
1657 
1658 template<>
1659 std::vector<Teuchos::RCP<MuemexArg>> processNeeds<complex_t,mm_LocalOrd,long long,mm_node_t>(const Factory* factory, std::string& needsParam, Level& lvl) {
1660  throw std::runtime_error("Muemex does not support long long for global indices");
1661 }
1662 
1663 template<>
1664 void processProvides<double,mm_LocalOrd,long long,mm_node_t>(std::vector<Teuchos::RCP<MuemexArg>>& mexOutput, const Factory* factory, std::string& providesParam, Level& lvl) {
1665  throw std::runtime_error("Muemex does not support long long for global indices");
1666 }
1667 
1668 template<>
1669 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) {
1670  throw std::runtime_error("Muemex does not support long long for global indices");
1671 }
1672 
1673 
1674 }// end namespace
1675 #endif //HAVE_MUELU_MATLAB error handler
1676 #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 >()
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)
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)
size_type size() const
std::vector< Teuchos::RCP< MuemexArg > > processNeeds(const Factory *factory, std::string &needsParam, Level &lvl)
void addLevelVariable(const T &data, std::string &name, Level &lvl, const FactoryBase *fact=NoFactory::get())
int FillComplete(bool OptimizeDataStorage=true)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void fillMatlabArray< complex_t >(complex_t *array, const mxArray *mxa, int n)
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.
Definition: MueLu_Level.hpp:99
template complex_t loadDataFromMatlab< complex_t >(const mxArray *mxa)
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)
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)
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
TransListIter iter
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)
virtual const RCP< const FactoryBase > GetFactory(const std::string &varName) const =0
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
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)
MuemexType getMuemexType< complex_t >()
int mwIndex
RCP< Xpetra::MultiVector< SC, LO, GO, NO > > TpetraMultiVector_To_XpetraMultiVector(const Teuchos::RCP< Tpetra::MultiVector< SC, LO, GO, NO > > &Vtpetra)
bool is_null() const