MueLu Version of the Day
Loading...
Searching...
No Matches
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
56using Teuchos::RCP;
57using Teuchos::rcp;
58using namespace std;
59
60namespace MueLu {
61
62extern bool rewrap_ints;
63
64/* ******************************* */
65/* getMuemexType */
66/* ******************************* */
67
68template<typename T> MuemexType getMuemexType(const T & data) {throw std::runtime_error("Unknown Type");}
69
70template<> MuemexType getMuemexType(const int & data) {return INT;}
71template<> MuemexType getMuemexType<int>() {return INT;}
72template<> MuemexType getMuemexType<bool>() {return BOOL;}
73
74template<> MuemexType getMuemexType(const double & data) {return DOUBLE;}
76
77template<> MuemexType getMuemexType(const std::string & data) {return STRING;}
79
80template<> MuemexType getMuemexType(const complex_t& data) {return COMPLEX;}
82
83template<> MuemexType getMuemexType(const RCP<Xpetra_map> & data) {return XPETRA_MAP;}
84template<> MuemexType getMuemexType<RCP<Xpetra_map> >() {return XPETRA_MAP;}
85
86template<> MuemexType getMuemexType(const RCP<Xpetra_ordinal_vector> & data) {return XPETRA_ORDINAL_VECTOR;}
87template<> MuemexType getMuemexType<RCP<Xpetra_ordinal_vector>>() {return XPETRA_ORDINAL_VECTOR;}
88
89template<> MuemexType getMuemexType(const RCP<Tpetra_MultiVector_double> & data) {return TPETRA_MULTIVECTOR_DOUBLE;}
90template<> MuemexType getMuemexType<RCP<Tpetra_MultiVector_double>>() {return TPETRA_MULTIVECTOR_DOUBLE;}
91
92template<> MuemexType getMuemexType(const RCP<Tpetra_MultiVector_complex>& data) {return TPETRA_MULTIVECTOR_COMPLEX;}
93template<> MuemexType getMuemexType<RCP<Tpetra_MultiVector_complex>>() {return TPETRA_MULTIVECTOR_COMPLEX;}
94
95template<> MuemexType getMuemexType(const RCP<Tpetra_CrsMatrix_double> & data) {return TPETRA_MATRIX_DOUBLE;}
96template<> MuemexType getMuemexType<RCP<Tpetra_CrsMatrix_double>>() {return TPETRA_MATRIX_DOUBLE;}
97
98template<> MuemexType getMuemexType(const RCP<Tpetra_CrsMatrix_complex> & data) {return TPETRA_MATRIX_COMPLEX;}
99template<> MuemexType getMuemexType<RCP<Tpetra_CrsMatrix_complex>>() {return TPETRA_MATRIX_COMPLEX;}
100
101template<> MuemexType getMuemexType(const RCP<Xpetra_MultiVector_double> & data) {return XPETRA_MULTIVECTOR_DOUBLE;}
102template<> MuemexType getMuemexType<RCP<Xpetra_MultiVector_double>>() {return XPETRA_MULTIVECTOR_DOUBLE;}
103
104template<> MuemexType getMuemexType(const RCP<Xpetra_MultiVector_complex> & data) {return XPETRA_MULTIVECTOR_COMPLEX;}
105template<> MuemexType getMuemexType<RCP<Xpetra_MultiVector_complex>>() {return XPETRA_MULTIVECTOR_COMPLEX;}
106
107template<> MuemexType getMuemexType(const RCP<Xpetra_Matrix_double> & data) {return XPETRA_MATRIX_DOUBLE;}
108template<> MuemexType getMuemexType<RCP<Xpetra_Matrix_double>>() {return XPETRA_MATRIX_DOUBLE;}
109
110template<> MuemexType getMuemexType(const RCP<Xpetra_Matrix_complex> & data) {return XPETRA_MATRIX_COMPLEX;}
111template<> MuemexType getMuemexType<RCP<Xpetra_Matrix_complex>>() {return XPETRA_MATRIX_COMPLEX;}
112
113template<> MuemexType getMuemexType(const RCP<Epetra_CrsMatrix> & data) {return EPETRA_CRSMATRIX;}
114template<> MuemexType getMuemexType<RCP<Epetra_CrsMatrix>>() {return EPETRA_CRSMATRIX;}
115
116template<> MuemexType getMuemexType(const RCP<Epetra_MultiVector> & data) {return EPETRA_MULTIVECTOR;}
117template<> MuemexType getMuemexType<RCP<Epetra_MultiVector>>() {return EPETRA_MULTIVECTOR;}
118
119template<> MuemexType getMuemexType(const RCP<MAggregates>& data) {return AGGREGATES;}
120template<> MuemexType getMuemexType<RCP<MAggregates>>() {return AGGREGATES;}
121
122template<> MuemexType getMuemexType(const RCP<MAmalInfo>& data) {return AMALGAMATION_INFO;}
123template<> MuemexType getMuemexType<RCP<MAmalInfo>>() {return AMALGAMATION_INFO;}
124
125template<> MuemexType getMuemexType(const RCP<MGraph>& data) {return GRAPH;}
126template<> MuemexType getMuemexType<RCP<MGraph>>() {return GRAPH;}
127
128#ifdef HAVE_MUELU_INTREPID2
129template<> MuemexType getMuemexType(const RCP<FieldContainer_ordinal>& data) {return FIELDCONTAINER_ORDINAL;}
130template<> MuemexType getMuemexType<RCP<FieldContainer_ordinal>>() {return FIELDCONTAINER_ORDINAL;}
131#endif
132
133/* "prototypes" for specialized functions used in other specialized functions */
134
135template<> mxArray* createMatlabSparse<double>(int numRows, int numCols, int nnz);
136template<> mxArray* createMatlabSparse<complex_t>(int numRows, int numCols, int nnz);
137template<> mxArray* createMatlabMultiVector<double>(int numRows, int numCols);
138template<> mxArray* createMatlabMultiVector<complex_t>(int numRows, int numCols);
139template<> void fillMatlabArray<double>(double* array, const mxArray* mxa, int n);
140template<> void fillMatlabArray<complex_t>(complex_t* array, const mxArray* mxa, int n);
141template<> mxArray* saveDataToMatlab(RCP<Xpetra_MultiVector_double>& data);
142template<> mxArray* saveDataToMatlab(RCP<Xpetra_MultiVector_complex>& data);
143template<> mxArray* saveDataToMatlab(RCP<Xpetra_Matrix_double>& data);
144template<> mxArray* saveDataToMatlab(RCP<Xpetra_Matrix_complex>& data);
145
146/* ******************************* */
147/* loadDataFromMatlab */
148/* ******************************* */
149
150template<>
151int 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
179template<>
180bool loadDataFromMatlab<bool>(const mxArray* mxa)
181{
182 return *((bool*) mxGetData(mxa));
183}
184
185template<>
186double loadDataFromMatlab<double>(const mxArray* mxa)
187{
188 return *((double*) mxGetPr(mxa));
189}
190
191template<>
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
199template<>
200string 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
211template<>
212RCP<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
240template<>
241RCP<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
263template<>
264RCP<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
287template<>
288RCP<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
320template<>
321RCP<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
395template<>
396RCP<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
460template<>
461RCP<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
467template<>
468RCP<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
474template<>
475RCP<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
481template<>
482RCP<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
488template<>
489RCP<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 }
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
537template<>
538RCP<Epetra_MultiVector> loadDataFromMatlab<RCP<Epetra_MultiVector>>(const mxArray* mxa)
539{
540 int nr = mxGetM(mxa);
541 int nc = mxGetN(mxa);
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
547template<>
548RCP<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
597template<>
598RCP<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
605template<>
606RCP<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
683template<>
684RCP<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
709template<>
710mxArray* 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
718template<>
719mxArray* 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
727template<>
728mxArray* saveDataToMatlab(double& data)
729{
730 return mxCreateDoubleScalar(data);
731}
732
733template<>
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
743template<>
744mxArray* saveDataToMatlab(string& data)
745{
746 return mxCreateString(data.c_str());
747}
748
749template<>
750mxArray* 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
767template<>
768mxArray* 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
783template<>
784mxArray* 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
790template<>
791mxArray* 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
797template<>
798mxArray* 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
804template<>
805mxArray* 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
811template<>
812mxArray* 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->getLocalMaxNumRowEntries();
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
930template<>
931mxArray* 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/*
1046template<>
1047mxArray* 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
1068template<>
1069mxArray* 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
1089template<>
1090mxArray* 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
1110template<>
1111mxArray* 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
1117template<>
1118mxArray* 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
1126template<>
1127mxArray* 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
1214template<>
1215mxArray* saveDataToMatlab(RCP<MAmalInfo>& data)
1216{
1217 throw runtime_error("AmalgamationInfo not supported in MueMex yet.");
1218 return mxCreateDoubleScalar(0);
1219}
1220
1221template<>
1222mxArray* 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
1321template<>
1322mxArray* 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
1348template<typename T>
1350{
1351 data = loadDataFromMatlab<T>(mxa);
1352}
1353
1354template<typename T>
1356{
1357 return saveDataToMatlab<T>(data);
1358}
1359
1360template<typename T>
1361MuemexData<T>::MuemexData(T& dataToCopy, MuemexType dataType) : MuemexArg(dataType)
1362{
1363 data = dataToCopy;
1364}
1365
1366template<typename T>
1367MuemexData<T>::MuemexData(T& dataToCopy) : MuemexData(dataToCopy, getMuemexType(dataToCopy)) {}
1368
1369template<typename T>
1371{
1372 return data;
1373}
1374
1375template<typename T>
1377{
1378 this->data = newData;
1379}
1380
1381/* ***************************** */
1382/* More Template Functions */
1383/* ***************************** */
1384
1385template<typename T>
1386void 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
1392template<typename T>
1393const 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
1406template<typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, typename Node>
1407std::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
1538template<typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, typename Node>
1539void 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
1670template<>
1671std::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
1675template<>
1676std::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
1680template<>
1681void 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
1685template<>
1686void 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
int mwIndex
MueLu::DefaultScalar Scalar
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Default implementation of FactoryAcceptor::GetFactory()
MueLu representation of a compressed row storage graph.
Class that holds all level-specific information.
int GetLevelID() const
Return level number.
void AddKeepFlag(const std::string &ename, const FactoryBase *factory=NoFactory::get(), KeepType keep=MueLu::Keep)
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Namespace for MueLu classes and methods.
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()
template int loadDataFromMatlab< int >(const mxArray *mxa)
@ UserData
User data are always kept. This flag is set automatically when Level::Set("data", data) is used....
MuemexType getMuemexType< double >()
Tpetra::Map ::local_ordinal_type mm_LocalOrd
RCP< Xpetra::MultiVector< SC, LO, GO, NO > > TpetraMultiVector_To_XpetraMultiVector(const Teuchos::RCP< Tpetra::MultiVector< SC, LO, GO, NO > > &Vtpetra)
MuemexType getMuemexType< complex_t >()
template double loadDataFromMatlab< double >(const mxArray *mxa)
mxArray * createMatlabSparse< double >(int numRows, int numCols, int nnz)
template bool loadDataFromMatlab< bool >(const mxArray *mxa)
MuemexType getMuemexType< int >()
template string loadDataFromMatlab< string >(const mxArray *mxa)
template complex_t loadDataFromMatlab< complex_t >(const mxArray *mxa)
int * mwIndex_to_int(int N, mwIndex *mwi_array)
void processProvides< double, mm_LocalOrd, long long, mm_node_t >(std::vector< Teuchos::RCP< MuemexArg > > &mexOutput, const Factory *factory, std::string &providesParam, Level &lvl)
std::vector< Teuchos::RCP< MuemexArg > > processNeeds(const Factory *factory, std::string &needsParam, Level &lvl)
void fillMatlabArray< double >(double *array, const mxArray *mxa, int n)
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)
void processProvides(std::vector< Teuchos::RCP< MuemexArg > > &mexOutput, const Factory *factory, std::string &providesParam, Level &lvl)
MueLu::Aggregates< mm_LocalOrd, mm_GlobalOrd, mm_node_t > MAggregates
Tpetra::Map muemex_map_type
Tpetra::Map ::global_ordinal_type mm_GlobalOrd
void fillMatlabArray< complex_t >(complex_t *array, const mxArray *mxa, int n)
const T & getLevelVariable(std::string &name, Level &lvl)
RCP< Xpetra::Matrix< SC, LO, GO, NO > > TpetraCrs_To_XpetraMatrix(const Teuchos::RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > &Atpetra)
mxArray * createMatlabMultiVector< complex_t >(int numRows, int numCols)
MuemexType getMuemexType< bool >()
std::complex< double > complex_t
void addLevelVariable(const T &data, std::string &name, Level &lvl, const FactoryBase *fact=NoFactory::get())
MuemexType getMuemexType< string >()
std::vector< std::string > tokenizeList(const std::string &params)
mxArray * createMatlabSparse< complex_t >(int numRows, int numCols, int nnz)
mxArray * createMatlabMultiVector< double >(int numRows, int numCols)
std::vector< Teuchos::RCP< MuemexArg > > processNeeds< double, mm_LocalOrd, long long, mm_node_t >(const Factory *factory, std::string &needsParam, Level &lvl)
template mxArray * saveDataToMatlab(bool &data)