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