MueLu  Version of the Day
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 /* "prototypes" for specialized functions used in other specialized functions */
129 
130 template<> mxArray* createMatlabSparse<double>(int numRows, int numCols, int nnz);
131 template<> mxArray* createMatlabSparse<complex_t>(int numRows, int numCols, int nnz);
132 template<> mxArray* createMatlabMultiVector<double>(int numRows, int numCols);
133 template<> mxArray* createMatlabMultiVector<complex_t>(int numRows, int numCols);
134 template<> void fillMatlabArray<double>(double* array, const mxArray* mxa, int n);
135 template<> void fillMatlabArray<complex_t>(complex_t* array, const mxArray* mxa, int n);
136 template<> mxArray* saveDataToMatlab(RCP<Xpetra_MultiVector_double>& data);
137 template<> mxArray* saveDataToMatlab(RCP<Xpetra_MultiVector_complex>& data);
138 template<> mxArray* saveDataToMatlab(RCP<Xpetra_Matrix_double>& data);
139 template<> mxArray* saveDataToMatlab(RCP<Xpetra_Matrix_complex>& data);
140 
141 /* ******************************* */
142 /* loadDataFromMatlab */
143 /* ******************************* */
144 
145 template<>
146 int loadDataFromMatlab<int>(const mxArray* mxa)
147 {
148  mxClassID probIDtype = mxGetClassID(mxa);
149  int rv;
150  if(probIDtype == mxINT32_CLASS)
151  {
152  rv = *((int*) mxGetData(mxa));
153  }
154  else if(probIDtype == mxLOGICAL_CLASS)
155  {
156  rv = (int) *((bool*) mxGetData(mxa));
157  }
158  else if(probIDtype == mxDOUBLE_CLASS)
159  {
160  rv = (int) *((double*) mxGetData(mxa));
161  }
162  else if(probIDtype == mxUINT32_CLASS)
163  {
164  rv = (int) *((unsigned int*) mxGetData(mxa));
165  }
166  else
167  {
168  rv = -1;
169  throw std::runtime_error("Error: Unrecognized numerical type.");
170  }
171  return rv;
172 }
173 
174 template<>
175 bool loadDataFromMatlab<bool>(const mxArray* mxa)
176 {
177  return *((bool*) mxGetData(mxa));
178 }
179 
180 template<>
181 double loadDataFromMatlab<double>(const mxArray* mxa)
182 {
183  return *((double*) mxGetPr(mxa));
184 }
185 
186 template<>
188 {
189  double realpart = real<double>(*((double*) mxGetPr(mxa)));
190  double imagpart = imag<double>(*((double*) mxGetPi(mxa)));
191  return complex_t(realpart, imagpart);
192 }
193 
194 template<>
195 string loadDataFromMatlab<string>(const mxArray* mxa)
196 {
197  string rv = "";
198  if(!mxGetClassID(mxa) != mxCHAR_CLASS)
199  {
200  throw runtime_error("Can't construct string from anything but a char array.");
201  }
202  rv = string(mxArrayToString(mxa));
203  return rv;
204 }
205 
206 template<>
207 RCP<Xpetra_map> loadDataFromMatlab<RCP<Xpetra_map>>(const mxArray* mxa)
208 {
210  int nr = mxGetM(mxa);
211  int nc = mxGetN(mxa);
212  if(nr != 1)
213  throw std::runtime_error("A Xpetra::Map representation from MATLAB must be a single row vector.");
214  double* pr = mxGetPr(mxa);
215  mm_GlobalOrd numGlobalIndices = nc;
216 
217  std::vector<mm_GlobalOrd> localGIDs(numGlobalIndices);
218  for(int i = 0; i < int(numGlobalIndices); i++) {
219  localGIDs[i] = Teuchos::as<mm_GlobalOrd>(pr[i]);
220  }
221 
222  const Teuchos::ArrayView<const mm_GlobalOrd> localGIDs_view(&localGIDs[0],localGIDs.size());
223  RCP<Xpetra_map> map =
224  Xpetra::MapFactory<mm_LocalOrd, mm_GlobalOrd, mm_node_t>::Build(
225  Xpetra::UseTpetra,
227  localGIDs_view,
228  0, comm);
229 
230  if(map.is_null())
231  throw runtime_error("Failed to create Xpetra::Map.");
232  return map;
233 }
234 
235 template<>
236 RCP<Xpetra_ordinal_vector> loadDataFromMatlab<RCP<Xpetra_ordinal_vector>>(const mxArray* mxa)
237 {
239  if(mxGetN(mxa) != 1 && mxGetM(mxa) != 1)
240  throw std::runtime_error("An OrdinalVector from MATLAB must be a single row or column vector.");
241  mm_GlobalOrd numGlobalIndices = mxGetM(mxa) * mxGetN(mxa);
242  RCP<Xpetra::Map<mm_LocalOrd, mm_GlobalOrd, mm_node_t>> map = Xpetra::MapFactory<mm_LocalOrd, mm_GlobalOrd, mm_node_t>::Build(Xpetra::UseTpetra, numGlobalIndices, 0, comm);
243  if(mxGetClassID(mxa) != mxINT32_CLASS)
244  throw std::runtime_error("Can only construct LOVector with int32 data.");
245  int* array = (int*) mxGetData(mxa);
246  if(map.is_null())
247  throw runtime_error("Failed to create map for Xpetra ordinal vector.");
248  RCP<Xpetra_ordinal_vector> loVec = Xpetra::VectorFactory<mm_LocalOrd, mm_LocalOrd, mm_GlobalOrd, mm_node_t>::Build(map, false);
249  if(loVec.is_null())
250  throw runtime_error("Failed to create ordinal vector with Xpetra::VectorFactory.");
251  for(int i = 0; i < int(numGlobalIndices); i++)
252  {
253  loVec->replaceGlobalValue(i, 0, array[i]);
254  }
255  return loVec;
256 }
257 
258 template<>
259 RCP<Tpetra_MultiVector_double> loadDataFromMatlab<RCP<Tpetra_MultiVector_double>>(const mxArray* mxa)
260 {
262  try
263  {
264  int nr = mxGetM(mxa);
265  int nc = mxGetN(mxa);
266  double* pr = mxGetPr(mxa);
267  RCP<const Teuchos::Comm<int>> comm = Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
268  //numGlobalIndices for map constructor is the number of rows in matrix/vectors, right?
270  //Allocate a new array of complex values to use with the multivector
271  Teuchos::ArrayView<const double> arrView(pr, nr * nc);
272  mv = rcp(new Tpetra::MultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(map, arrView, size_t(nr), size_t(nc)));
273  }
274  catch(std::exception& e)
275  {
276  mexPrintf("Error constructing Tpetra MultiVector.\n");
277  std::cout << e.what() << std::endl;
278  }
279  return mv;
280 }
281 
282 template<>
283 RCP<Tpetra_MultiVector_complex> loadDataFromMatlab<RCP<Tpetra_MultiVector_complex>>(const mxArray* mxa)
284 {
286  try
287  {
288  int nr = mxGetM(mxa);
289  int nc = mxGetN(mxa);
290  double* pr = mxGetPr(mxa);
291  double* pi = mxGetPi(mxa);
292  RCP<const Teuchos::Comm<int>> comm = Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
293  //numGlobalIndices for map constructor is the number of rows in matrix/vectors, right?
295  //Allocate a new array of complex values to use with the multivector
296  complex_t* myArr = new complex_t[nr * nc];
297  for(int n = 0; n < nc; n++)
298  {
299  for(int m = 0; m < nr; m++)
300  {
301  myArr[n * nr + m] = complex_t(pr[n * nr + m], pi[n * nr + m]);
302  }
303  }
304  Teuchos::ArrayView<complex_t> arrView(myArr, nr * nc);
305  mv = rcp(new Tpetra::MultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(map, arrView, nr, nc));
306  }
307  catch(std::exception& e)
308  {
309  mexPrintf("Error constructing Tpetra MultiVector.\n");
310  std::cout << e.what() << std::endl;
311  }
312  return mv;
313 }
314 
315 template<>
316 RCP<Tpetra_CrsMatrix_double> loadDataFromMatlab<RCP<Tpetra_CrsMatrix_double>>(const mxArray* mxa)
317 {
318  bool success = false;
320 
321  int* colptr = NULL;
322  int* rowind = NULL;
323 
324  try
325  {
327  //numGlobalIndices is just the number of rows in the matrix
328  const Tpetra::global_size_t numGlobalIndices = mxGetM(mxa);
329  const mm_GlobalOrd indexBase = 0;
330  RCP<const muemex_map_type> rowMap = rcp(new muemex_map_type(numGlobalIndices, indexBase, comm));
331  RCP<const muemex_map_type> domainMap = rcp(new muemex_map_type(mxGetN(mxa), indexBase, comm));
332  A = Tpetra::createCrsMatrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(rowMap);
333  double* valueArray = mxGetPr(mxa);
334  int nc = mxGetN(mxa);
335  if(rewrap_ints)
336  {
337  //mwIndex_to_int allocates memory so must delete[] later
338  colptr = mwIndex_to_int(nc + 1, mxGetJc(mxa));
339  rowind = mwIndex_to_int(colptr[nc], mxGetIr(mxa));
340  }
341  else
342  {
343  rowind = (int*) mxGetIr(mxa);
344  colptr = (int*) mxGetJc(mxa);
345  }
346  for(int i = 0; i < nc; i++)
347  {
348  for(int j = colptr[i]; j < colptr[i + 1]; j++)
349  {
350  //'array' of 1 element, containing column (in global matrix).
352  //'array' of 1 element, containing value
353  Teuchos::ArrayView<double> vals = Teuchos::ArrayView<double>(&valueArray[j], 1);
354  A->insertGlobalValues(rowind[j], cols, vals);
355  }
356  }
357  A->fillComplete(domainMap, rowMap);
358  if(rewrap_ints)
359  {
360  delete[] rowind; rowind = NULL;
361  delete[] colptr; colptr = NULL;
362  }
363  success = true;
364  }
365  catch(std::exception& e)
366  {
367  if(rewrap_ints)
368  {
369  if(rowind!=NULL) delete[] rowind;
370  if(colptr!=NULL) delete[] colptr;
371  rowind = NULL;
372  colptr = NULL;
373  }
374  mexPrintf("Error while constructing Tpetra matrix:\n");
375  std::cout << e.what() << std::endl;
376  }
377  if(!success)
378  mexErrMsgTxt("An error occurred while setting up a Tpetra matrix.\n");
379  return A;
380 }
381 
382 template<>
383 RCP<Tpetra_CrsMatrix_complex> loadDataFromMatlab<RCP<Tpetra_CrsMatrix_complex>>(const mxArray* mxa)
384 {
386  //Create a map in order to create the matrix (taken from muelu basic example - complex)
387  try
388  {
389  RCP<const Teuchos::Comm<int>> comm = Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
390  const Tpetra::global_size_t numGlobalIndices = mxGetM(mxa);
391  const mm_GlobalOrd indexBase = 0;
392  RCP<const muemex_map_type> rowMap = rcp(new muemex_map_type(numGlobalIndices, indexBase, comm));
393  RCP<const muemex_map_type> domainMap = rcp(new muemex_map_type(mxGetN(mxa), indexBase, comm));
394  A = Tpetra::createCrsMatrix<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(rowMap);
395  double* realArray = mxGetPr(mxa);
396  double* imagArray = mxGetPi(mxa);
397  int* colptr;
398  int* rowind;
399  int nc = mxGetN(mxa);
400  if(rewrap_ints)
401  {
402  //mwIndex_to_int allocates memory so must delete[] later
403  colptr = mwIndex_to_int(nc + 1, mxGetJc(mxa));
404  rowind = mwIndex_to_int(colptr[nc], mxGetIr(mxa));
405  }
406  else
407  {
408  rowind = (int*) mxGetIr(mxa);
409  colptr = (int*) mxGetJc(mxa);
410  }
411  for(int i = 0; i < nc; i++)
412  {
413  for(int j = colptr[i]; j < colptr[i + 1]; j++)
414  {
415  //here assuming that complex_t will always be defined as std::complex<double>
416  //use 'value' over and over again with Teuchos::ArrayViews to insert into matrix
417  complex_t value = std::complex<double>(realArray[j], imagArray[j]);
420  A->insertGlobalValues(rowind[j], cols, vals);
421  }
422  }
423  A->fillComplete(domainMap, rowMap);
424  if(rewrap_ints)
425  {
426  delete[] rowind;
427  delete[] colptr;
428  }
429  }
430  catch(std::exception& e)
431  {
432  mexPrintf("Error while constructing tpetra matrix:\n");
433  std::cout << e.what() << std::endl;
434  }
435  return A;
436 }
437 
438 template<>
439 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)
440 {
441  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);
442  return MueLu::TpetraCrs_To_XpetraMatrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tmat);
443 }
444 
445 template<>
446 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)
447 {
448  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);
449  return MueLu::TpetraCrs_To_XpetraMatrix<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tmat);
450 }
451 
452 template<>
453 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)
454 {
455  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);
456  return MueLu::TpetraMultiVector_To_XpetraMultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tpetraMV);
457 }
458 
459 template<>
460 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)
461 {
462  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);
463  return MueLu::TpetraMultiVector_To_XpetraMultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(tpetraMV);
464 }
465 
466 template<>
467 RCP<Epetra_CrsMatrix> loadDataFromMatlab<RCP<Epetra_CrsMatrix>>(const mxArray* mxa)
468 {
469  RCP<Epetra_CrsMatrix> matrix;
470  try
471  {
472  int* colptr;
473  int* rowind;
474  double* vals = mxGetPr(mxa);
475  int nr = mxGetM(mxa);
476  int nc = mxGetN(mxa);
477  if(rewrap_ints)
478  {
479  colptr = mwIndex_to_int(nc + 1, mxGetJc(mxa));
480  rowind = mwIndex_to_int(colptr[nc], mxGetIr(mxa));
481  }
482  else
483  {
484  rowind = (int*) mxGetIr(mxa);
485  colptr = (int*) mxGetJc(mxa);
486  }
487  Epetra_SerialComm Comm;
488  Epetra_Map RangeMap(nr, 0, Comm);
489  Epetra_Map DomainMap(nc, 0, Comm);
490  matrix = rcp(new Epetra_CrsMatrix(Epetra_DataAccess::Copy, RangeMap, DomainMap, 0));
491  /* Do the matrix assembly */
492  for(int i = 0; i < nc; i++)
493  {
494  for(int j = colptr[i]; j < colptr[i + 1]; j++)
495  {
496  //global row, # of entries, value array, column indices array
497  matrix->InsertGlobalValues(rowind[j], 1, &vals[j], &i);
498  }
499  }
500  matrix->FillComplete(DomainMap, RangeMap);
501  if(rewrap_ints)
502  {
503  delete [] rowind;
504  delete [] colptr;
505  }
506  }
507  catch(std::exception& e)
508  {
509  mexPrintf("An error occurred while setting up an Epetra matrix:\n");
510  std::cout << e.what() << std::endl;
511  }
512  return matrix;
513 }
514 
515 template<>
516 RCP<Epetra_MultiVector> loadDataFromMatlab<RCP<Epetra_MultiVector>>(const mxArray* mxa)
517 {
518  int nr = mxGetM(mxa);
519  int nc = mxGetN(mxa);
520  Epetra_SerialComm Comm;
521  Epetra_BlockMap map(nr * nc, 1, 0, Comm);
522  return rcp(new Epetra_MultiVector(Epetra_DataAccess::Copy, map, mxGetPr(mxa), nr, nc));
523 }
524 
525 template<>
526 RCP<MAggregates> loadDataFromMatlab<RCP<MAggregates>>(const mxArray* mxa)
527 {
528  if(mxGetNumberOfElements(mxa) != 1)
529  throw runtime_error("Aggregates must be individual structs in MATLAB.");
530  if(!mxIsStruct(mxa))
531  throw runtime_error("Trying to pull aggregates from non-struct MATLAB object.");
532  //assume that in matlab aggregate structs will only be stored in a 1x1 array
533  //mxa must have the same fields as the ones declared in constructAggregates function in muelu.m for this to work
534  const int correctNumFields = 5; //change if more fields are added to the aggregates representation in constructAggregates in muelu.m
535  if(mxGetNumberOfFields(mxa) != correctNumFields)
536  throw runtime_error("Aggregates structure has wrong number of fields.");
537  //Pull MuemexData types back out
538  int nVert = *(int*) mxGetData(mxGetField(mxa, 0, "nVertices"));
539  int nAgg = *(int*) mxGetData(mxGetField(mxa, 0, "nAggregates"));
540  //Now have all the data needed to fully reconstruct the aggregate
541  //Use similar approach as UserAggregationFactory (which is written for >1 thread but will just be serial here)
543  int myRank = comm->getRank();
544  Xpetra::UnderlyingLib lib = Xpetra::UseTpetra;
545  RCP<Xpetra::Map<mm_LocalOrd, mm_GlobalOrd, mm_node_t>> map = Xpetra::MapFactory<mm_LocalOrd, mm_GlobalOrd, mm_node_t>::Build(lib, nVert, 0, comm);
546  RCP<MAggregates> agg = rcp(new MAggregates(map));
547  agg->SetNumAggregates(nAgg);
548  //Get handles for the vertex2AggId and procwinner arrays in reconstituted aggregates object
549  //this is serial so all procwinner values will be same (0)
550  ArrayRCP<mm_LocalOrd> vertex2AggId = agg->GetVertex2AggId()->getDataNonConst(0); //the '0' means first (and only) column of multivector, since is just vector
551  ArrayRCP<mm_LocalOrd> procWinner = agg->GetProcWinner()->getDataNonConst(0);
552  //mm_LocalOrd and int are equivalent, so is ok to talk about aggSize with just 'int'
553  //Deep copy the entire vertex2AggID and isRoot arrays, which are both nVert items long
554  //At the same time, set ProcWinner
555  mxArray* vertToAggID_in = mxGetField(mxa, 0, "vertexToAggID");
556  int* vertToAggID_inArray = (int*) mxGetData(vertToAggID_in);
557  mxArray* rootNodes_in = mxGetField(mxa, 0, "rootNodes");
558  int* rootNodes_inArray = (int*) mxGetData(rootNodes_in);
559  for(int i = 0; i < nVert; i++)
560  {
561  vertex2AggId[i] = vertToAggID_inArray[i];
562  procWinner[i] = myRank; //all nodes are going to be on the same proc
563  agg->SetIsRoot(i, false); //the ones that are root will be set in next loop
564  }
565  for(int i = 0; i < nAgg; i++) //rootNodesToCopy is an array of node IDs which are the roots of their aggs
566  {
567  agg->SetIsRoot(rootNodes_inArray[i], true);
568  }
569  //Now recompute the aggSize array and cache the results in the object
570  agg->ComputeAggregateSizes(true, true);
571  agg->AggregatesCrossProcessors(false);
572  return agg;
573 }
574 
575 template<>
576 RCP<MAmalInfo> loadDataFromMatlab<RCP<MAmalInfo>>(const mxArray* mxa)
577 {
578  RCP<MAmalInfo> amal;
579  throw runtime_error("AmalgamationInfo not supported in Muemex yet.");
580  return amal;
581 }
582 
583 template<>
584 RCP<MGraph> loadDataFromMatlab<RCP<MGraph>>(const mxArray* mxa)
585 {
586  //mxa must be struct with logical sparse matrix called 'edges' and Nx1 int32 array 'boundaryNodes'
587  mxArray* edges = mxGetField(mxa, 0, "edges");
588  mxArray* boundaryNodes = mxGetField(mxa, 0, "boundaryNodes");
589  if(edges == NULL)
590  throw runtime_error("Graph structure in MATLAB must have a field called 'edges' (logical sparse matrix)");
591  if(boundaryNodes == NULL)
592  throw runtime_error("Graph structure in MATLAB must have a field called 'boundaryNodes' (int32 array containing list of boundary nodes)");
593  int* boundaryList = (int*) mxGetData(boundaryNodes);
594  if(!mxIsSparse(edges) || mxGetClassID(edges) != mxLOGICAL_CLASS)
595  throw runtime_error("Graph edges must be stored as a logical sparse matrix.");
596  // Note that Matlab stores sparse matrices in column major format.
597  mwIndex* matlabColPtrs = mxGetJc(edges);
598  mwIndex* matlabRowIndices = mxGetIr(edges);
599  mm_GlobalOrd nRows = (mm_GlobalOrd) mxGetM(edges);
600 
601  // Create and populate row-major CRS data structures for Xpetra::TpetraCrsGraph.
602 
603  // calculate number of nonzeros in each row
604  Teuchos::Array<int> entriesPerRow(nRows);
605  int nnz = matlabColPtrs[mxGetN(edges)]; //last entry in matlabColPtrs
606  for(int i = 0; i < nnz; i++)
607  entriesPerRow[matlabRowIndices[i]]++;
608  // Populate usual row index array. We don't need this for the Xpetra Graph ctor, but
609  // it's convenient for building up the column index array, which the ctor does need.
610  Teuchos::Array<int> rows(nRows+1);
611  rows[0] = 0;
612  for(int i = 0; i < nRows; i++)
613  rows[i+1] = rows[i] + entriesPerRow[i];
614  Teuchos::Array<int> cols(nnz); //column index array
615  Teuchos::Array<int> insertionsPerRow(nRows,0); //track of #insertions done per row
616  int ncols = mxGetN(edges);
617  for (int colNum=0; colNum<ncols; ++colNum) {
618  int ci = matlabColPtrs[colNum];
619  for (int j=ci; j<Teuchos::as<int>(matlabColPtrs[colNum+1]); ++j) {
620  int rowNum = matlabRowIndices[j];
621  cols[ rows[rowNum] + insertionsPerRow[rowNum] ] = colNum;
622  insertionsPerRow[rowNum]++;
623  }
624  }
625  //Find maximum
626  int maxNzPerRow = 0;
627  for(int i = 0; i < nRows; i++) {
628  if(maxNzPerRow < entriesPerRow[i])
629  maxNzPerRow = entriesPerRow[i];
630  }
631 
633  typedef Xpetra::TpetraMap<mm_LocalOrd, mm_GlobalOrd, mm_node_t> MMap;
634  RCP<MMap> map = rcp(new MMap(nRows, 0, comm));
635  typedef Xpetra::TpetraCrsGraph<mm_LocalOrd, mm_GlobalOrd, mm_node_t> TpetraGraph;
636  RCP<TpetraGraph> tgraph = rcp(new TpetraGraph(map, (size_t) maxNzPerRow));
637  //Populate tgraph in compressed-row format. Must get each row individually...
638  for(int i = 0; i < nRows; ++i) {
639  tgraph->insertGlobalIndices((mm_GlobalOrd) i, cols(rows[i],entriesPerRow[i]));
640  }
641  tgraph->fillComplete(map, map);
643  //Set boundary nodes
644  int numBoundaryNodes = mxGetNumberOfElements(boundaryNodes);
645  bool* boundaryFlags = new bool[nRows];
646  for(int i = 0; i < nRows; i++)
647  {
648  boundaryFlags[i] = false;
649  }
650  for(int i = 0; i < numBoundaryNodes; i++)
651  {
652  boundaryFlags[boundaryList[i]] = true;
653  }
654  ArrayRCP<bool> boundaryNodesInput(boundaryFlags, 0, nRows, true);
655  mgraph->SetBoundaryNodeMap(boundaryNodesInput);
656  return mgraph;
657 }
658 
659 /* ******************************* */
660 /* saveDataToMatlab */
661 /* ******************************* */
662 
663 template<>
664 mxArray* saveDataToMatlab(int& data)
665 {
666  mwSize dims[] = {1, 1};
667  mxArray* mxa = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
668  *((int*) mxGetData(mxa)) = data;
669  return mxa;
670 }
671 
672 template<>
673 mxArray* saveDataToMatlab(bool& data)
674 {
675  mwSize dims[] = {1, 1};
676  mxArray* mxa = mxCreateLogicalArray(2, dims);
677  *((bool*) mxGetData(mxa)) = data;
678  return mxa;
679 }
680 
681 template<>
682 mxArray* saveDataToMatlab(double& data)
683 {
684  return mxCreateDoubleScalar(data);
685 }
686 
687 template<>
689 {
690  mwSize dims[] = {1, 1};
691  mxArray* mxa = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxCOMPLEX);
692  *((double*) mxGetPr(mxa)) = real<double>(data);
693  *((double*) mxGetPi(mxa)) = imag<double>(data);
694  return mxa;
695 }
696 
697 template<>
698 mxArray* saveDataToMatlab(string& data)
699 {
700  return mxCreateString(data.c_str());
701 }
702 
703 template<>
705 {
706  //Precondition: Memory has already been allocated by MATLAB for the array.
707  int nc = data->getGlobalNumElements();
708  int nr = 1;
709  mxArray* output = createMatlabMultiVector<double>(nr, nc);
710  double* array = (double*) malloc(sizeof(double) * nr * nc);
711  for(int col = 0; col < nc; col++)
712  {
713  mm_GlobalOrd gid = data->getGlobalElement(col);
714  array[col] = Teuchos::as<double>(gid);
715  }
716  fillMatlabArray<double>(array, output, nc * nr);
717  free(array);
718  return output;
719 }
720 
721 template<>
723 {
724  mwSize len = data->getGlobalLength();
725  //create a single column vector
726  mwSize dimensions[] = {len, 1};
727  mxArray* rv = mxCreateNumericArray(2, dimensions, mxINT32_CLASS, mxREAL);
728  int* dataPtr = (int*) mxGetData(rv);
729  ArrayRCP<const mm_LocalOrd> arr = data->getData(0);
730  for(int i = 0; i < int(data->getGlobalLength()); i++)
731  {
732  dataPtr[i] = arr[i];
733  }
734  return rv;
735 }
736 
737 template<>
738 mxArray* saveDataToMatlab(RCP<Tpetra::MultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>>& data)
739 {
741  return saveDataToMatlab(xmv);
742 }
743 
744 template<>
745 mxArray* saveDataToMatlab(RCP<Tpetra::MultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>>& data)
746 {
748  return saveDataToMatlab(xmv);
749 }
750 
751 template<>
753 {
755  return saveDataToMatlab(xmat);
756 }
757 
758 template<>
760 {
762  return saveDataToMatlab(xmat);
763 }
764 
765 template<>
767 {
768  typedef double Scalar;
769  int nr = data->getGlobalNumRows();
770  int nc = data->getGlobalNumCols();
771  int nnz = data->getGlobalNumEntries();
772 #ifdef VERBOSE_OUTPUT
773  RCP<Teuchos::FancyOStream> fancyStream = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
774  mat->describe(*fancyStream, Teuchos::VERB_EXTREME);
775 #endif
776  mxArray* mxa = createMatlabSparse<Scalar>(nr, nc, nnz);
777  mwIndex* ir = mxGetIr(mxa);
778  mwIndex* jc = mxGetJc(mxa);
779  for(int i = 0; i < nc + 1; i++)
780  {
781  jc[i] = 0;
782  }
783  size_t maxEntriesPerRow = data->getGlobalMaxNumRowEntries();
784  int* rowProgress = new int[nc];
785  //The array that will be copied to Pr and (if complex) Pi later
786  Scalar* sparseVals = new Scalar[nnz];
787  size_t numEntries;
788  if(data->isLocallyIndexed())
789  {
790  Scalar* rowValArray = new Scalar[maxEntriesPerRow];
791  Teuchos::ArrayView<Scalar> rowVals(rowValArray, maxEntriesPerRow);
792  mm_LocalOrd* rowIndicesArray = new mm_LocalOrd[maxEntriesPerRow];
793  Teuchos::ArrayView<mm_LocalOrd> rowIndices(rowIndicesArray, maxEntriesPerRow);
794  for(mm_LocalOrd m = 0; m < nr; m++) //All rows in the Xpetra matrix
795  {
796  data->getLocalRowCopy(m, rowIndices, rowVals, numEntries); //Get the row
797  for(mm_LocalOrd entry = 0; entry < int(numEntries); entry++) //All entries in row
798  {
799  jc[rowIndices[entry] + 1]++; //for each entry, increase jc for the entry's column
800  }
801  }
802  //now jc holds the number of elements in each column, but needs cumulative sum over all previous columns also
803  int entriesAccum = 0;
804  for(int n = 0; n <= nc; n++)
805  {
806  int temp = entriesAccum;
807  entriesAccum += jc[n];
808  jc[n] += temp;
809  }
810  //Jc now populated with colptrs
811  for(int i = 0; i < nc; i++)
812  {
813  rowProgress[i] = 0;
814  }
815  //Row progress values like jc but keep track as the MATLAB matrix is being filled in
816  for(mm_LocalOrd m = 0; m < nr; m++) //rows
817  {
818  data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
819  for(mm_LocalOrd i = 0; i < int(numEntries); i++) //entries in row m (NOT columns)
820  {
821  //row is m, col is rowIndices[i], val is rowVals[i]
822  mm_LocalOrd col = rowIndices[i];
823  sparseVals[jc[col] + rowProgress[col]] = rowVals[i]; //Set value
824  ir[jc[col] + rowProgress[col]] = m; //Set row at which value occurs
825  rowProgress[col]++;
826  }
827  }
828  delete[] rowIndicesArray;
829  }
830  else
831  {
834  for(mm_GlobalOrd m = 0; m < nr; m++)
835  {
836  data->getGlobalRowView(m, rowIndices, rowVals);
837  for(mm_GlobalOrd n = 0; n < rowIndices.size(); n++)
838  {
839  jc[rowIndices[n] + 1]++;
840  }
841  }
842  //Last element of jc is just nnz
843  jc[nc] = nnz;
844  //Jc now populated with colptrs
845  for(int i = 0; i < nc; i++)
846  {
847  rowProgress[i] = 0;
848  }
849  int entriesAccum = 0;
850  for(int n = 0; n <= nc; n++)
851  {
852  int temp = entriesAccum;
853  entriesAccum += jc[n];
854  jc[n] += temp;
855  }
856  //Row progress values like jc but keep track as the MATLAB matrix is being filled in
857  for(mm_GlobalOrd m = 0; m < nr; m++) //rows
858  {
859  data->getGlobalRowView(m, rowIndices, rowVals);
860  for(mm_LocalOrd i = 0; i < rowIndices.size(); i++) //entries in row m
861  {
862  mm_GlobalOrd col = rowIndices[i]; //row is m, col is rowIndices[i], val is rowVals[i]
863  sparseVals[jc[col] + rowProgress[col]] = rowVals[i]; //Set value
864  ir[jc[col] + rowProgress[col]] = m; //Set row at which value occurs
865  rowProgress[col]++;
866  }
867  }
868  }
869  //finally, copy sparseVals into pr (and pi, if complex)
870  fillMatlabArray<Scalar>(sparseVals, mxa, nnz);
871  delete[] sparseVals;
872  delete[] rowProgress;
873  return mxa;
874 }
875 
876 template<>
878 {
879  typedef complex_t Scalar;
880  int nr = data->getGlobalNumRows();
881  int nc = data->getGlobalNumCols();
882  int nnz = data->getGlobalNumEntries();
883 #ifdef VERBOSE_OUTPUT
884  RCP<Teuchos::FancyOStream> fancyStream = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
885  mat->describe(*fancyStream, Teuchos::VERB_EXTREME);
886 #endif
887  mxArray* mxa = createMatlabSparse<Scalar>(nr, nc, nnz);
888  mwIndex* ir = mxGetIr(mxa);
889  mwIndex* jc = mxGetJc(mxa);
890  for(int i = 0; i < nc + 1; i++)
891  {
892  jc[i] = 0;
893  }
894  size_t maxEntriesPerRow = data->getGlobalMaxNumRowEntries();
895  int* rowProgress = new int[nc];
896  //The array that will be copied to Pr and (if complex) Pi later
897  Scalar* sparseVals = new Scalar[nnz];
898  size_t numEntries;
899  if(data->isLocallyIndexed())
900  {
901  Scalar* rowValArray = new Scalar[maxEntriesPerRow];
902  Teuchos::ArrayView<Scalar> rowVals(rowValArray, maxEntriesPerRow);
903  mm_LocalOrd* rowIndicesArray = new mm_LocalOrd[maxEntriesPerRow];
904  Teuchos::ArrayView<mm_LocalOrd> rowIndices(rowIndicesArray, maxEntriesPerRow);
905  for(mm_LocalOrd m = 0; m < nr; m++) //All rows in the Xpetra matrix
906  {
907  data->getLocalRowCopy(m, rowIndices, rowVals, numEntries); //Get the row
908  for(mm_LocalOrd entry = 0; entry < int(numEntries); entry++) //All entries in row
909  {
910  jc[rowIndices[entry] + 1]++; //for each entry, increase jc for the entry's column
911  }
912  }
913  //now jc holds the number of elements in each column, but needs cumulative sum over all previous columns also
914  int entriesAccum = 0;
915  for(int n = 0; n <= nc; n++)
916  {
917  int temp = entriesAccum;
918  entriesAccum += jc[n];
919  jc[n] += temp;
920  }
921  //Jc now populated with colptrs
922  for(int i = 0; i < nc; i++)
923  {
924  rowProgress[i] = 0;
925  }
926  //Row progress values like jc but keep track as the MATLAB matrix is being filled in
927  for(mm_LocalOrd m = 0; m < nr; m++) //rows
928  {
929  data->getLocalRowCopy(m, rowIndices, rowVals, numEntries);
930  for(mm_LocalOrd i = 0; i < int(numEntries); i++) //entries in row m (NOT columns)
931  {
932  //row is m, col is rowIndices[i], val is rowVals[i]
933  mm_LocalOrd col = rowIndices[i];
934  sparseVals[jc[col] + rowProgress[col]] = rowVals[i]; //Set value
935  ir[jc[col] + rowProgress[col]] = m; //Set row at which value occurs
936  rowProgress[col]++;
937  }
938  }
939  delete[] rowIndicesArray;
940  }
941  else
942  {
945  for(mm_GlobalOrd m = 0; m < nr; m++)
946  {
947  data->getGlobalRowView(m, rowIndices, rowVals);
948  for(mm_GlobalOrd n = 0; n < rowIndices.size(); n++)
949  {
950  jc[rowIndices[n] + 1]++;
951  }
952  }
953  //Last element of jc is just nnz
954  jc[nc] = nnz;
955  //Jc now populated with colptrs
956  for(int i = 0; i < nc; i++)
957  {
958  rowProgress[i] = 0;
959  }
960  int entriesAccum = 0;
961  for(int n = 0; n <= nc; n++)
962  {
963  int temp = entriesAccum;
964  entriesAccum += jc[n];
965  jc[n] += temp;
966  }
967  //Row progress values like jc but keep track as the MATLAB matrix is being filled in
968  for(mm_GlobalOrd m = 0; m < nr; m++) //rows
969  {
970  data->getGlobalRowView(m, rowIndices, rowVals);
971  for(mm_LocalOrd i = 0; i < rowIndices.size(); i++) //entries in row m
972  {
973  mm_GlobalOrd col = rowIndices[i]; //row is m, col is rowIndices[i], val is rowVals[i]
974  sparseVals[jc[col] + rowProgress[col]] = rowVals[i]; //Set value
975  ir[jc[col] + rowProgress[col]] = m; //Set row at which value occurs
976  rowProgress[col]++;
977  }
978  }
979  }
980  //finally, copy sparseVals into pr (and pi, if complex)
981  fillMatlabArray<Scalar>(sparseVals, mxa, nnz);
982  delete[] sparseVals;
983  delete[] rowProgress;
984  return mxa;
985 }
986 
987 /*
988 template<>
989 mxArray* saveDataToMatlab(RCP<Xpetra::MultiVector<Scalar, mm_LocalOrd, mm_GlobalOrd, mm_node_t>>& data)
990 {
991  //Precondition: Memory has already been allocated by MATLAB for the array.
992  int nr = data->getGlobalLength();
993  int nc = data->getNumVectors();
994  mxArray* output = createMatlabMultiVector<Scalar>(nr, nc);
995  Scalar* array = (Scalar*) malloc(sizeof(Scalar) * nr * nc);
996  for(int col = 0; col < nc; col++)
997  {
998  Teuchos::ArrayRCP<const Scalar> colData = data->getData(col);
999  for(int row = 0; row < nr; row++)
1000  {
1001  array[col * nr + row] = colData[row];
1002  }
1003  }
1004  fillMatlabArray<Scalar>(array, output, nc * nr);
1005  free(array);
1006  return output;
1007 }
1008 */
1009 
1010 template<>
1011 mxArray* saveDataToMatlab(RCP<Xpetra::MultiVector<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>>& data)
1012 {
1013  //Precondition: Memory has already been allocated by MATLAB for the array.
1014  int nr = data->getGlobalLength();
1015  int nc = data->getNumVectors();
1016  mxArray* output = createMatlabMultiVector<double>(nr, nc);
1017  double* array = (double*) malloc(sizeof(double) * nr * nc);
1018  for(int col = 0; col < nc; col++)
1019  {
1020  Teuchos::ArrayRCP<const double> colData = data->getData(col);
1021  for(int row = 0; row < nr; row++)
1022  {
1023  array[col * nr + row] = colData[row];
1024  }
1025  }
1026  fillMatlabArray<double>(array, output, nc * nr);
1027  free(array);
1028  return output;
1029 }
1030 
1031 template<>
1032 mxArray* saveDataToMatlab(RCP<Xpetra::MultiVector<complex_t, mm_LocalOrd, mm_GlobalOrd, mm_node_t>>& data)
1033 {
1034  //Precondition: Memory has already been allocated by MATLAB for the array.
1035  int nr = data->getGlobalLength();
1036  int nc = data->getNumVectors();
1037  mxArray* output = createMatlabMultiVector<complex_t>(nr, nc);
1038  complex_t* array = (complex_t*) malloc(sizeof(complex_t) * nr * nc);
1039  for(int col = 0; col < nc; col++)
1040  {
1041  Teuchos::ArrayRCP<const complex_t> colData = data->getData(col);
1042  for(int row = 0; row < nr; row++)
1043  {
1044  array[col * nr + row] = colData[row];
1045  }
1046  }
1047  fillMatlabArray<complex_t>(array, output, nc * nr);
1048  free(array);
1049  return output;
1050 }
1051 
1052 template<>
1054 {
1055  RCP<Xpetra_Matrix_double> xmat = EpetraCrs_To_XpetraMatrix<double, mm_LocalOrd, mm_GlobalOrd, mm_node_t>(data);
1056  return saveDataToMatlab(xmat);
1057 }
1058 
1059 template<>
1061 {
1062  mxArray* output = mxCreateDoubleMatrix(data->GlobalLength(), data->NumVectors(), mxREAL);
1063  double* dataPtr = mxGetPr(output);
1064  data->ExtractCopy(dataPtr, data->GlobalLength());
1065  return output;
1066 }
1067 
1068 template<>
1070 {
1071  //Set up array of inputs for matlab constructAggregates
1072  int numNodes = data->GetVertex2AggId()->getData(0).size();
1073  int numAggs = data->GetNumAggregates();
1074  mxArray* dataIn[5];
1075  mwSize singleton[] = {1, 1};
1076  dataIn[0] = mxCreateNumericArray(2, singleton, mxINT32_CLASS, mxREAL);
1077  *((int*) mxGetData(dataIn[0])) = numNodes;
1078  dataIn[1] = mxCreateNumericArray(2, singleton, mxINT32_CLASS, mxREAL);
1079  *((int*) mxGetData(dataIn[1])) = numAggs;
1080  mwSize nodeArrayDims[] = {(mwSize) numNodes, 1}; //dimensions for Nx1 array, where N is number of nodes (vert2Agg)
1081  dataIn[2] = mxCreateNumericArray(2, nodeArrayDims, mxINT32_CLASS, mxREAL);
1082  int* vtaid = (int*) mxGetData(dataIn[2]);
1083  ArrayRCP<const mm_LocalOrd> vertexToAggID = data->GetVertex2AggId()->getData(0);
1084  for(int i = 0; i < numNodes; i++)
1085  {
1086  vtaid[i] = vertexToAggID[i];
1087  }
1088  mwSize aggArrayDims[] = {(mwSize) numAggs, 1}; //dims for Nx1 array, where N is number of aggregates (rootNodes, aggSizes)
1089  dataIn[3] = mxCreateNumericArray(2, aggArrayDims, mxINT32_CLASS, mxREAL);
1090  //First, find out if the aggregates even have 1 root node per aggregate. If not, assume roots are invalid and assign ourselves
1091  int totalRoots = 0;
1092  for(int i = 0; i < numNodes; i++)
1093  {
1094  if(data->IsRoot(i))
1095  totalRoots++;
1096  }
1097  bool reassignRoots = false;
1098  if(totalRoots != numAggs)
1099  {
1100  cout << endl << "Warning: Number of root nodes and number of aggregates do not match." << endl;
1101  cout << "Will reassign root nodes when writing aggregates to matlab." << endl << endl;
1102  reassignRoots = true;
1103  }
1104  int* rn = (int*) mxGetData(dataIn[3]); //list of root nodes (in no particular order)
1105  {
1106  if(reassignRoots)
1107  {
1108  //For each aggregate, just pick the first node we see in it and set it as root
1109  int lastFoundNode = 0; //heuristic for speed, a node in aggregate N+1 is likely to come very soon after a node in agg N
1110  for(int i = 0; i < numAggs; i++)
1111  {
1112  rn[i] = -1;
1113  for(int j = lastFoundNode; j < lastFoundNode + numNodes; j++)
1114  {
1115  int index = j % numNodes;
1116  if(vertexToAggID[index] == i)
1117  {
1118  rn[i] = index;
1119  lastFoundNode = index;
1120  }
1121  }
1122  TEUCHOS_TEST_FOR_EXCEPTION(rn[i] == -1, runtime_error, "Invalid aggregates: Couldn't find any node in aggregate #" << i << ".");
1123  }
1124  }
1125  else
1126  {
1127  int i = 0; //iterates over aggregate IDs
1128  for(int j = 0; j < numNodes; j++)
1129  {
1130  if(data->IsRoot(j))
1131  {
1132  if(i == numAggs)
1133  throw runtime_error("Cannot store invalid aggregates in MATLAB - more root nodes than aggregates.");
1134  rn[i] = j; //now we know this won't go out of bounds (rn's underlying matlab array is numAggs in length)
1135  i++;
1136  }
1137  }
1138  if(i + 1 < numAggs)
1139  throw runtime_error("Cannot store invalid aggregates in MATLAB - fewer root nodes than aggregates.");
1140  }
1141  }
1142  dataIn[4] = mxCreateNumericArray(1, aggArrayDims, mxINT32_CLASS, mxREAL);
1143  int* as = (int*) mxGetData(dataIn[4]); //list of aggregate sizes
1144  ArrayRCP<mm_LocalOrd> aggSizes = data->ComputeAggregateSizes();
1145  for(int i = 0; i < numAggs; i++)
1146  {
1147  as[i] = aggSizes[i];
1148  }
1149  mxArray* matlabAggs[1];
1150  int result = mexCallMATLAB(1, matlabAggs, 5, dataIn, "constructAggregates");
1151  if(result != 0)
1152  throw runtime_error("Matlab encountered an error while constructing aggregates struct.");
1153  return matlabAggs[0];
1154 }
1155 
1156 template<>
1158 {
1159  throw runtime_error("AmalgamationInfo not supported in MueMex yet.");
1160  return mxCreateDoubleScalar(0);
1161 }
1162 
1163 template<>
1165 {
1166  int numEntries = (int) data->GetGlobalNumEdges();
1167  int numRows = (int) data->GetDomainMap()->getGlobalNumElements(); //assume numRows == numCols
1168  mxArray* mat = mxCreateSparseLogicalMatrix(numRows, numRows, numEntries);
1169  mxLogical* outData = (mxLogical*) mxGetData(mat);
1170  mwIndex* rowInds = mxGetIr(mat);
1171  mwIndex* colPtrs = mxGetJc(mat);
1172  mm_LocalOrd* dataCopy = new mm_LocalOrd[numEntries];
1173  mm_LocalOrd* iter = dataCopy;
1174  int* entriesPerRow = new int[numRows];
1175  int* entriesPerCol = new int[numRows];
1176  for(int i = 0; i < numRows; i++)
1177  {
1178  entriesPerRow[i] = 0;
1179  entriesPerCol[i] = 0;
1180  }
1181  for(int i = 0; i < numRows; i++)
1182  {
1183  ArrayView<const mm_LocalOrd> neighbors = data->getNeighborVertices(i); //neighbors has the column indices for row i
1184  memcpy(iter, neighbors.getRawPtr(), sizeof(mm_LocalOrd) * neighbors.size());
1185  entriesPerRow[i] = neighbors.size();
1186  for(int j = 0; j < neighbors.size(); j++)
1187  {
1188  entriesPerCol[neighbors[j]]++;
1189  }
1190  iter += neighbors.size();
1191  }
1192  mwIndex** rowIndsByColumn = new mwIndex*[numRows]; //rowIndsByColumn[0] points to array of row indices in column 1
1193  mxLogical** valuesByColumn = new mxLogical*[numRows];
1194  int* numEnteredPerCol = new int[numRows];
1195  int accum = 0;
1196  for(int i = 0; i < numRows; i++)
1197  {
1198  rowIndsByColumn[i] = &rowInds[accum];
1199  //cout << "Entries in column " << i << " start at offset " << accum << endl;
1200  valuesByColumn[i] = &outData[accum];
1201  accum += entriesPerCol[i];
1202  if(accum > numEntries)
1203  throw runtime_error("potato");
1204  }
1205  for(int i = 0; i < numRows; i++)
1206  {
1207  numEnteredPerCol[i] = 0; //rowIndsByColumn[n][numEnteredPerCol[n]] gives the next place to put a row index
1208  }
1209  //entriesPerCol now has Jc information (col offsets)
1210  accum = 0; //keep track of cumulative index in dataCopy
1211  for(int row = 0; row < numRows; row++)
1212  {
1213  for(int entryInRow = 0; entryInRow < entriesPerRow[row]; entryInRow++)
1214  {
1215  int col = dataCopy[accum];
1216  accum++;
1217  rowIndsByColumn[col][numEnteredPerCol[col]] = row;
1218  valuesByColumn[col][numEnteredPerCol[col]] = (mxLogical) 1;
1219  numEnteredPerCol[col]++;
1220  }
1221  }
1222  accum = 0; //keep track of total entries over all columns
1223  for(int col = 0; col < numRows; col++)
1224  {
1225  colPtrs[col] = accum;
1226  accum += entriesPerCol[col];
1227  }
1228  colPtrs[numRows] = accum; //the last entry in jc, which is equivalent to numEntries
1229  delete[] numEnteredPerCol;
1230  delete[] rowIndsByColumn;
1231  delete[] valuesByColumn;
1232  delete[] dataCopy;
1233  delete[] entriesPerRow;
1234  delete[] entriesPerCol;
1235  //Construct list of boundary nodes
1236  const ArrayRCP<const bool> boundaryFlags = data->GetBoundaryNodeMap();
1237  int numBoundaryNodes = 0;
1238  for(int i = 0; i < boundaryFlags.size(); i++)
1239  {
1240  if(boundaryFlags[i])
1241  numBoundaryNodes++;
1242  }
1243  cout << "Graph has " << numBoundaryNodes << " Dirichlet boundary nodes." << endl;
1244  mwSize dims[] = {(mwSize) numBoundaryNodes, 1};
1245  mxArray* boundaryList = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
1246  int* dest = (int*) mxGetData(boundaryList);
1247  int* destIter = dest;
1248  for(int i = 0; i < boundaryFlags.size(); i++)
1249  {
1250  if(boundaryFlags[i])
1251  {
1252  *destIter = i;
1253  destIter++;
1254  }
1255  }
1256  mxArray* constructArgs[] = {mat, boundaryList};
1257  mxArray* out[1];
1258  mexCallMATLAB(1, out, 2, constructArgs, "constructGraph");
1259  return out[0];
1260 }
1261 
1262 template<typename T>
1264 {
1265  data = loadDataFromMatlab<T>(mxa);
1266 }
1267 
1268 template<typename T>
1270 {
1271  return saveDataToMatlab<T>(data);
1272 }
1273 
1274 template<typename T>
1275 MuemexData<T>::MuemexData(T& dataToCopy, MuemexType dataType) : MuemexArg(dataType)
1276 {
1277  data = dataToCopy;
1278 }
1279 
1280 template<typename T>
1281 MuemexData<T>::MuemexData(T& dataToCopy) : MuemexData(dataToCopy, getMuemexType(dataToCopy)) {}
1282 
1283 template<typename T>
1285 {
1286  return data;
1287 }
1288 
1289 template<typename T>
1290 void MuemexData<T>::setData(T& newData)
1291 {
1292  this->data = data;
1293 }
1294 
1295 /* ***************************** */
1296 /* More Template Functions */
1297 /* ***************************** */
1298 
1299 template<typename T>
1300 void addLevelVariable(const T& data, std::string& name, Level& lvl, const Factory * fact)
1301 {
1302  lvl.AddKeepFlag(name, fact, MueLu::UserData);
1303  lvl.Set<T>(name, data, fact);
1304 }
1305 
1306 template<typename T>
1307 const T& getLevelVariable(std::string& name, Level& lvl)
1308 {
1309  try
1310  {
1311  return lvl.Get<T>(name);
1312  }
1313  catch(std::exception& e)
1314  {
1315  throw std::runtime_error("Requested custom variable " + name + " is not in the level.");
1316  }
1317 }
1318 
1319 //Functions used to put data through matlab factories - first arg is "this" pointer of matlab factory
1320 template<typename Scalar = double, typename LocalOrdinal = mm_LocalOrd, typename GlobalOrdinal = mm_GlobalOrd, typename Node = mm_node_t>
1321 std::vector<Teuchos::RCP<MuemexArg>> processNeeds(const Factory* factory, std::string& needsParam, Level& lvl)
1322 {
1323  using namespace std;
1324  using namespace Teuchos;
1328  typedef RCP<AmalgamationInfo<LocalOrdinal, GlobalOrdinal, Node>> AmalgamationInfo_t;
1329  typedef RCP<MGraph> Graph_t;
1330  vector<string> needsList = tokenizeList(needsParam);
1331  vector<RCP<MuemexArg>> args;
1332  for(size_t i = 0; i < needsList.size(); i++)
1333  {
1334  if(needsList[i] == "A" || needsList[i] == "P" || needsList[i] == "R" || needsList[i]=="Ptent")
1335  {
1336  Matrix_t mydata = lvl.Get<Matrix_t>(needsList[i], factory->GetFactory(needsList[i]).get());
1337  args.push_back(rcp(new MuemexData<Matrix_t>(mydata)));
1338  }
1339  else if(needsList[i] == "Nullspace" || needsList[i] == "Coordinates")
1340  {
1341  MultiVector_t mydata = lvl.Get<MultiVector_t>(needsList[i], factory->GetFactory(needsList[i]).get());
1342  args.push_back(rcp(new MuemexData<MultiVector_t>(mydata)));
1343  }
1344  else if(needsList[i] == "Aggregates")
1345  {
1346  Aggregates_t mydata = lvl.Get<Aggregates_t>(needsList[i], factory->GetFactory(needsList[i]).get());
1347  args.push_back(rcp(new MuemexData<Aggregates_t>(mydata)));
1348  }
1349  else if(needsList[i] == "UnAmalgamationInfo")
1350  {
1351  AmalgamationInfo_t mydata = lvl.Get<AmalgamationInfo_t>(needsList[i], factory->GetFactory(needsList[i]).get());
1352  args.push_back(rcp(new MuemexData<AmalgamationInfo_t>(mydata)));
1353  }
1354  else if(needsList[i] == "Level")
1355  {
1356  int levelNum = lvl.GetLevelID();
1357  args.push_back(rcp(new MuemexData<int>(levelNum)));
1358  }
1359  else if(needsList[i] == "Graph")
1360  {
1361  Graph_t mydata = lvl.Get<Graph_t>(needsList[i], factory->GetFactory(needsList[i]).get());
1362  args.push_back(rcp(new MuemexData<Graph_t>(mydata)));
1363  }
1364  else
1365  {
1366  vector<string> words;
1367  string badNameMsg = "Custom Muemex variables in \"Needs\" list require a type and a name, e.g. \"double myVal\". \n Leading and trailing spaces are OK. \n Don't know how to handle \"" + needsList[i] + "\".\n";
1368  //compare type without case sensitivity
1369  char* buf = (char*) malloc(needsList[i].size() + 1);
1370  strcpy(buf, needsList[i].c_str());
1371  for(char* iter = buf; *iter != ' '; iter++)
1372  {
1373  if(*iter == 0)
1374  {
1375  free(buf);
1376  throw runtime_error(badNameMsg);
1377  }
1378  *iter = (char) tolower(*iter);
1379  }
1380  const char* wordDelim = " ";
1381  char* mark = strtok(buf, wordDelim);
1382  while(mark != NULL)
1383  {
1384  string wordStr(mark);
1385  words.push_back(wordStr);
1386  mark = strtok(NULL, wordDelim);
1387  }
1388  if(words.size() != 2)
1389  {
1390  free(buf);
1391  throw runtime_error(badNameMsg);
1392  }
1393  char* typeStr = (char*) words[0].c_str();
1394  if(strstr(typeStr, "ordinalvector"))
1395  {
1397  LOVector_t mydata = getLevelVariable<LOVector_t>(needsList[i], lvl);
1398  args.push_back(rcp(new MuemexData<LOVector_t>(mydata)));
1399  }
1400  else if(strstr(typeStr, "map"))
1401  {
1403  Map_t mydata = getLevelVariable<Map_t>(needsList[i], lvl);
1404  args.push_back(rcp(new MuemexData<Map_t>(mydata)));
1405  }
1406  else if(strstr(typeStr, "scalar"))
1407  {
1408  Scalar mydata = getLevelVariable<Scalar>(needsList[i], lvl);
1409  args.push_back(rcp(new MuemexData<Scalar>(mydata)));
1410  }
1411  else if(strstr(typeStr, "double"))
1412  {
1413  double mydata = getLevelVariable<double>(needsList[i], lvl);
1414  args.push_back(rcp(new MuemexData<double>(mydata)));
1415  }
1416  else if(strstr(typeStr, "complex"))
1417  {
1418  complex_t mydata = getLevelVariable<complex_t>(needsList[i], lvl);
1419  args.push_back(rcp(new MuemexData<complex_t>(mydata)));
1420  }
1421  else if(strstr(typeStr, "matrix"))
1422  {
1423  Matrix_t mydata = getLevelVariable<Matrix_t>(needsList[i], lvl);
1424  args.push_back(rcp(new MuemexData<Matrix_t>(mydata)));
1425  }
1426  else if(strstr(typeStr, "multivector"))
1427  {
1428  MultiVector_t mydata = getLevelVariable<MultiVector_t>(needsList[i], lvl);
1429  args.push_back(rcp(new MuemexData<MultiVector_t>(mydata)));
1430  }
1431  else if(strstr(typeStr, "int"))
1432  {
1433  int mydata = getLevelVariable<int>(needsList[i], lvl);
1434  args.push_back(rcp(new MuemexData<int>(mydata)));
1435  }
1436  else if(strstr(typeStr, "string"))
1437  {
1438  string mydata = getLevelVariable<string>(needsList[i], lvl);
1439  args.push_back(rcp(new MuemexData<string>(mydata)));
1440  }
1441  else
1442  {
1443  free(buf);
1444  throw std::runtime_error(words[0] + " is not a known variable type.");
1445  }
1446  free(buf);
1447  }
1448  }
1449  return args;
1450 }
1451 
1452 template<typename Scalar = double, typename LocalOrdinal = mm_LocalOrd, typename GlobalOrdinal = mm_GlobalOrd, typename Node = mm_node_t>
1453 void processProvides(std::vector<Teuchos::RCP<MuemexArg>>& mexOutput, const Factory* factory, std::string& providesParam, Level& lvl)
1454 {
1455  using namespace std;
1456  using namespace Teuchos;
1460  typedef RCP<AmalgamationInfo<LocalOrdinal, GlobalOrdinal, Node>> AmalgamationInfo_t;
1461  typedef RCP<MGraph> Graph_t;
1462  vector<string> provides = tokenizeList(providesParam);
1463  for(size_t i = 0; i < size_t(provides.size()); i++)
1464  {
1465  if(provides[i] == "A" || provides[i] == "P" || provides[i] == "R" || provides[i]=="Ptent")
1466  {
1467  RCP<MuemexData<Matrix_t>> mydata = Teuchos::rcp_static_cast<MuemexData<Matrix_t>>(mexOutput[i]);
1468  lvl.Set(provides[i], mydata->getData(), factory);
1469  }
1470  else if(provides[i] == "Nullspace" || provides[i] == "Coordinates")
1471  {
1472  RCP<MuemexData<MultiVector_t>> mydata = Teuchos::rcp_static_cast<MuemexData<MultiVector_t>>(mexOutput[i]);
1473  lvl.Set(provides[i], mydata->getData(), factory);
1474  }
1475  else if(provides[i] == "Aggregates")
1476  {
1477  RCP<MuemexData<Aggregates_t>> mydata = Teuchos::rcp_static_cast<MuemexData<Aggregates_t>>(mexOutput[i]);
1478  lvl.Set(provides[i], mydata->getData(), factory);
1479  }
1480  else if(provides[i] == "UnAmalgamationInfo")
1481  {
1482  RCP<MuemexData<AmalgamationInfo_t>> mydata = Teuchos::rcp_static_cast<MuemexData<AmalgamationInfo_t>>(mexOutput[i]);
1483  lvl.Set(provides[i], mydata->getData(), factory);
1484  }
1485  else if(provides[i] == "Graph")
1486  {
1487  RCP<MuemexData<Graph_t>> mydata = Teuchos::rcp_static_cast<MuemexData<Graph_t>>(mexOutput[i]);
1488  lvl.Set(provides[i], mydata->getData(), factory);
1489  }
1490  else
1491  {
1492  vector<string> words;
1493  string badNameMsg = "Custom Muemex variables in \"Provides\" list require a type and a name, e.g. \"double myVal\". \n Leading and trailing spaces are OK. \n Don't know how to handle \"" + provides[i] + "\".\n";
1494  //compare type without case sensitivity
1495  char* buf = (char*) malloc(provides[i].size() + 1);
1496  strcpy(buf, provides[i].c_str());
1497  for(char* iter = buf; *iter != ' '; iter++)
1498  {
1499  if(*iter == 0)
1500  {
1501  free(buf);
1502  throw runtime_error(badNameMsg);
1503  }
1504  *iter = (char) tolower(*iter);
1505  }
1506  const char* wordDelim = " ";
1507  char* mark = strtok(buf, wordDelim);
1508  while(mark != NULL)
1509  {
1510  string wordStr(mark);
1511  words.push_back(wordStr);
1512  mark = strtok(NULL, wordDelim);
1513  }
1514  if(words.size() != 2)
1515  {
1516  free(buf);
1517  throw runtime_error(badNameMsg);
1518  }
1519  const char* typeStr = words[0].c_str();
1520  if(strstr(typeStr, "ordinalvector"))
1521  {
1523  RCP<MuemexData<LOVector_t>> mydata = Teuchos::rcp_static_cast<MuemexData<LOVector_t>>(mexOutput[i]);
1524  addLevelVariable<LOVector_t>(mydata->getData(), words[1], lvl, factory);
1525  }
1526  else if(strstr(typeStr, "map"))
1527  {
1529  RCP<MuemexData<Map_t>> mydata = Teuchos::rcp_static_cast<MuemexData<Map_t>>(mexOutput[i]);
1530  addLevelVariable<Map_t>(mydata->getData(), words[1], lvl, factory);
1531  }
1532  else if(strstr(typeStr, "scalar"))
1533  {
1534  RCP<MuemexData<Scalar>> mydata = Teuchos::rcp_static_cast<MuemexData<Scalar>>(mexOutput[i]);
1535  addLevelVariable<Scalar>(mydata->getData(), words[1], lvl, factory);
1536  }
1537  else if(strstr(typeStr, "double"))
1538  {
1539  RCP<MuemexData<double>> mydata = Teuchos::rcp_static_cast<MuemexData<double>>(mexOutput[i]);
1540  addLevelVariable<double>(mydata->getData(), words[1], lvl, factory);
1541  }
1542  else if(strstr(typeStr, "complex"))
1543  {
1544  RCP<MuemexData<complex_t>> mydata = Teuchos::rcp_static_cast<MuemexData<complex_t>>(mexOutput[i]);
1545  addLevelVariable<complex_t>(mydata->getData(), words[1], lvl, factory);
1546  }
1547  else if(strstr(typeStr, "matrix"))
1548  {
1549  RCP<MuemexData<Matrix_t>> mydata = Teuchos::rcp_static_cast<MuemexData<Matrix_t>>(mexOutput[i]);
1550  addLevelVariable<Matrix_t>(mydata->getData(), words[1], lvl, factory);
1551  }
1552  else if(strstr(typeStr, "multivector"))
1553  {
1554  RCP<MuemexData<MultiVector_t>> mydata = Teuchos::rcp_static_cast<MuemexData<MultiVector_t>>(mexOutput[i]);
1555  addLevelVariable<MultiVector_t>(mydata->getData(), words[1], lvl, factory);
1556  }
1557  else if(strstr(typeStr, "int"))
1558  {
1559  RCP<MuemexData<int>> mydata = Teuchos::rcp_static_cast<MuemexData<int>>(mexOutput[i]);
1560  addLevelVariable<int>(mydata->getData(), words[1], lvl, factory);
1561  }
1562  else if(strstr(typeStr, "bool"))
1563  {
1564  RCP<MuemexData<bool> > mydata = Teuchos::rcp_static_cast<MuemexData<bool> >(mexOutput[i]);
1565  addLevelVariable<bool>(mydata->getData(), words[1], lvl, factory);
1566  }
1567  else if(strstr(typeStr, "string"))
1568  {
1569  RCP<MuemexData<string>> mydata = Teuchos::rcp_static_cast<MuemexData<string>>(mexOutput[i]);
1570  addLevelVariable<string>(mydata->getData(), words[1], lvl, factory);
1571  }
1572  else
1573  {
1574  free(buf);
1575  throw std::runtime_error(words[0] + " is not a known variable type.");
1576  }
1577  free(buf);
1578  }
1579  }
1580 }
1581 
1582 // Throwable Stubs for long long
1583 
1584 template<>
1585 std::vector<Teuchos::RCP<MuemexArg>> processNeeds<double,mm_LocalOrd,long long,mm_node_t>(const Factory* factory, std::string& needsParam, Level& lvl) {
1586  throw std::runtime_error("Muemex does not support long long for global indices");
1587 }
1588 
1589 template<>
1590 std::vector<Teuchos::RCP<MuemexArg>> processNeeds<complex_t,mm_LocalOrd,long long,mm_node_t>(const Factory* factory, std::string& needsParam, Level& lvl) {
1591  throw std::runtime_error("Muemex does not support long long for global indices");
1592 }
1593 
1594 template<>
1595 void processProvides<double,mm_LocalOrd,long long,mm_node_t>(std::vector<Teuchos::RCP<MuemexArg>>& mexOutput, const Factory* factory, std::string& providesParam, Level& lvl) {
1596  throw std::runtime_error("Muemex does not support long long for global indices");
1597 }
1598 
1599 template<>
1600 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) {
1601  throw std::runtime_error("Muemex does not support long long for global indices");
1602 }
1603 
1604 
1605 }// end namespace
1606 #endif //HAVE_MUELU_MATLAB error handler
1607 #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)
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->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
mxArray * createMatlabSparse< complex_t >(int numRows, int numCols, int nnz)
bool rewrap_ints
mxArray * saveDataToMatlab(RCP< MGraph > &data)
MuemexType getMuemexType< string >()
size_type size() const
MueLu representation of a compressed row storage graph.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void processProvides< complex_t, mm_LocalOrd, long long, mm_node_t >(std::vector< Teuchos::RCP< MuemexArg >> &mexOutput, const Factory *factory, std::string &providesParam, Level &lvl)
User data are always kept. This flag is set automatically when Level::Set("data", data) is used...
size_type size() const
Namespace for MueLu classes and methods.
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)
bool is_null() const
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
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)
void addLevelVariable(const T &data, std::string &name, Level &lvl, Factory *fact=NoFactory::get())
template int loadDataFromMatlab< int >(const mxArray *mxa)
bool is_null() const
template bool loadDataFromMatlab< bool >(const mxArray *mxa)
const T & getLevelVariable(std::string &name, Level &lvl)
MuemexType getMuemexType< int >()
mxArray * createMatlabMultiVector< complex_t >(int numRows, int numCols)
std::vector< Teuchos::RCP< MuemexArg > > processNeeds< complex_t, mm_LocalOrd, long long, mm_node_t >(const Factory *factory, std::string &needsParam, Level &lvl)
MuemexType getMuemexType(const RCP< MGraph > &data)
void AddKeepFlag(const std::string &ename, const FactoryBase *factory=NoFactory::get(), KeepType keep=MueLu::Keep)
MuemexType getMuemexType< bool >()
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > TpetraMultiVector_To_XpetraMultiVector(const Teuchos::RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &Vtpetra)
std::complex< double > complex_t
T * getRawPtr() const
mxArray * createMatlabMultiVector< double >(int numRows, int numCols)
template double loadDataFromMatlab< double >(const mxArray *mxa)
MueLu::Aggregates< mm_LocalOrd, mm_GlobalOrd, mm_node_t > MAggregates
TypeTo as(const TypeFrom &t)
MuemexType getMuemexType< double >()
Tpetra::Map muemex_map_type
void processProvides(std::vector< Teuchos::RCP< MuemexArg >> &mexOutput, const Factory *factory, std::string &providesParam, Level &lvl)
RCP< Xpetra::Matrix< SC, LO, GO, NO > > TpetraCrs_To_XpetraMatrix(const Teuchos::RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > &Atpetra)
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Default implementation of FactoryAcceptor::GetFactory()
MuemexType getMuemexType< complex_t >()
int mwIndex