xref: /petsc/src/mat/impls/aij/seq/matlab/aijmatlab.c (revision a96a251defc07734f733d1abf27b5765e87fdc8b)
1 #define PETSCMAT_DLL
2 
3 /*
4         Provides an interface for the Matlab engine sparse solver
5 
6 */
7 #include "src/mat/impls/aij/seq/aij.h"
8 
9 #include "engine.h"   /* Matlab include file */
10 #include "mex.h"      /* Matlab include file */
11 
12 typedef struct {
13   PetscErrorCode (*MatDuplicate)(Mat,MatDuplicateOption,Mat*);
14   PetscErrorCode (*MatView)(Mat,PetscViewer);
15   PetscErrorCode (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*);
16   PetscErrorCode (*MatILUDTFactor)(Mat,IS,IS,MatFactorInfo*,Mat*);
17   PetscErrorCode (*MatDestroy)(Mat);
18 } Mat_Matlab;
19 
20 
21 EXTERN_C_BEGIN
22 #undef __FUNCT__
23 #define __FUNCT__ "MatMatlabEnginePut_Matlab"
24 PetscErrorCode PETSCMAT_DLLEXPORT MatMatlabEnginePut_Matlab(PetscObject obj,void *mengine)
25 {
26   PetscErrorCode ierr;
27   Mat         B = (Mat)obj;
28   mxArray     *mat;
29   Mat_SeqAIJ  *aij = (Mat_SeqAIJ*)B->data;
30 
31   PetscFunctionBegin;
32   mat  = mxCreateSparse(B->n,B->m,aij->nz,mxREAL);
33   ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));CHKERRQ(ierr);
34   /* Matlab stores by column, not row so we pass in the transpose of the matrix */
35   ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr);
36   ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr);
37 
38   /* Matlab indices start at 0 for sparse (what a surprise) */
39 
40   ierr = PetscObjectName(obj);CHKERRQ(ierr);
41   engPutVariable((Engine *)mengine,obj->name,mat);
42   PetscFunctionReturn(0);
43 }
44 EXTERN_C_END
45 
46 EXTERN_C_BEGIN
47 #undef __FUNCT__
48 #define __FUNCT__ "MatMatlabEngineGet_Matlab"
49 PetscErrorCode PETSCMAT_DLLEXPORT MatMatlabEngineGet_Matlab(PetscObject obj,void *mengine)
50 {
51   PetscErrorCode ierr;
52   int            ii;
53   Mat            mat = (Mat)obj;
54   Mat_SeqAIJ     *aij = (Mat_SeqAIJ*)mat->data;
55   mxArray        *mmat;
56 
57   PetscFunctionBegin;
58   if (!aij->singlemalloc) {
59     ierr = PetscFree(aij->a);CHKERRQ(ierr);
60     ierr = PetscFree(aij->j);CHKERRQ(ierr);
61     ierr = PetscFree(aij->i);CHKERRQ(ierr);
62   } else {
63     ierr = PetscFree3(aij->a,aij->j,aij->i);CHKERRQ(ierr);
64   }
65 
66   mmat = engGetVariable((Engine *)mengine,obj->name);
67 
68   aij->nz           = (mxGetJc(mmat))[mat->m];
69   ierr  = PetscMalloc3(aij->nz,PetscScalar,&aij->a,aij->nz,PetscInt,&aij->j,mat->m+1,PetscInt,&aij->i);CHKERRQ(ierr);
70   aij->singlemalloc = PETSC_TRUE;
71   aij->freedata     = PETSC_TRUE;
72 
73   ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr);
74   /* Matlab stores by column, not row so we pass in the transpose of the matrix */
75   ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr);
76   ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr);
77 
78   for (ii=0; ii<mat->m; ii++) {
79     aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii];
80   }
81 
82   ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
83   ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
84 
85   PetscFunctionReturn(0);
86 }
87 EXTERN_C_END
88 
89 EXTERN_C_BEGIN
90 #undef __FUNCT__
91 #define __FUNCT__ "MatConvert_Matlab_SeqAIJ"
92 PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_Matlab_SeqAIJ(Mat A,const MatType type,MatReuse reuse,Mat *newmat)
93 {
94   PetscErrorCode ierr;
95   Mat            B=*newmat;
96   Mat_Matlab    *lu=(Mat_Matlab*)A->spptr;
97 
98   PetscFunctionBegin;
99   if (reuse == MAT_INITIAL_MATRIX) {
100     ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr);
101   }
102   B->ops->duplicate        = lu->MatDuplicate;
103   B->ops->view             = lu->MatView;
104   B->ops->lufactorsymbolic = lu->MatLUFactorSymbolic;
105   B->ops->iludtfactor      = lu->MatILUDTFactor;
106   B->ops->destroy          = lu->MatDestroy;
107 
108   ierr = PetscFree(lu);CHKERRQ(ierr);
109 
110   ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_matlab_C","",PETSC_NULL);CHKERRQ(ierr);
111   ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_matlab_seqaij_C","",PETSC_NULL);CHKERRQ(ierr);
112   ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C","",PETSC_NULL);CHKERRQ(ierr);
113   ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C","",PETSC_NULL);CHKERRQ(ierr);
114 
115   ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr);
116   *newmat = B;
117   PetscFunctionReturn(0);
118 }
119 EXTERN_C_END
120 
121 #undef __FUNCT__
122 #define __FUNCT__ "MatDestroy_Matlab"
123 PetscErrorCode MatDestroy_Matlab(Mat A)
124 {
125   PetscErrorCode ierr;
126 
127   PetscFunctionBegin;
128   ierr = MatConvert_Matlab_SeqAIJ(A,MATSEQAIJ,&A);CHKERRQ(ierr);
129   ierr = (*A->ops->destroy)(A);CHKERRQ(ierr);
130   PetscFunctionReturn(0);
131 }
132 
133 #undef __FUNCT__
134 #define __FUNCT__ "MatSolve_Matlab"
135 PetscErrorCode MatSolve_Matlab(Mat A,Vec b,Vec x)
136 {
137   PetscErrorCode ierr;
138   char            *_A,*_b,*_x;
139 
140   PetscFunctionBegin;
141   /* make sure objects have names; use default if not */
142   ierr = PetscObjectName((PetscObject)b);CHKERRQ(ierr);
143   ierr = PetscObjectName((PetscObject)x);CHKERRQ(ierr);
144 
145   ierr = PetscObjectGetName((PetscObject)A,&_A);CHKERRQ(ierr);
146   ierr = PetscObjectGetName((PetscObject)b,&_b);CHKERRQ(ierr);
147   ierr = PetscObjectGetName((PetscObject)x,&_x);CHKERRQ(ierr);
148   ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)b);CHKERRQ(ierr);
149   ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = u%s\\(l%s\\(p%s*%s));",_x,_A,_A,_A,_b);CHKERRQ(ierr);
150   ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_b);CHKERRQ(ierr);
151   /* ierr = PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(A->comm),stdout);CHKERRQ(ierr);  */
152   ierr = PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)x);CHKERRQ(ierr);
153   PetscFunctionReturn(0);
154 }
155 
156 #undef __FUNCT__
157 #define __FUNCT__ "MatLUFactorNumeric_Matlab"
158 PetscErrorCode MatLUFactorNumeric_Matlab(Mat A,MatFactorInfo *info,Mat *F)
159 {
160   Mat_SeqAIJ      *f = (Mat_SeqAIJ*)(*F)->data;
161   PetscErrorCode ierr;
162   size_t          len;
163   char            *_A,*name;
164 
165   PetscFunctionBegin;
166   ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);CHKERRQ(ierr);
167   _A   = A->name;
168   ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"[l_%s,u_%s,p_%s] = lu(%s',%g);",_A,_A,_A,_A,info->dtcol);CHKERRQ(ierr);
169   ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_A);CHKERRQ(ierr);
170   ierr = PetscStrlen(_A,&len);CHKERRQ(ierr);
171   ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr);
172   sprintf(name,"_%s",_A);
173   ierr = PetscObjectSetName((PetscObject)*F,name);CHKERRQ(ierr);
174   ierr = PetscFree(name);CHKERRQ(ierr);
175   PetscFunctionReturn(0);
176 }
177 
178 #undef __FUNCT__
179 #define __FUNCT__ "MatLUFactorSymbolic_Matlab"
180 PetscErrorCode MatLUFactorSymbolic_Matlab(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F)
181 {
182   PetscErrorCode ierr;
183   Mat_SeqAIJ *f;
184 
185   PetscFunctionBegin;
186   if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square");
187   ierr                       = MatCreate(A->comm,A->m,A->n,A->m,A->n,F);CHKERRQ(ierr);
188   ierr                       = MatSetType(*F,A->type_name);CHKERRQ(ierr);
189   ierr                       = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr);
190   (*F)->ops->solve           = MatSolve_Matlab;
191   (*F)->ops->lufactornumeric = MatLUFactorNumeric_Matlab;
192   (*F)->factor               = FACTOR_LU;
193   f                          = (Mat_SeqAIJ*)(*F)->data;
194   f->lu_dtcol = info->dtcol;
195   PetscFunctionReturn(0);
196 }
197 
198 /* ---------------------------------------------------------------------------------*/
199 #undef __FUNCT__
200 #define __FUNCT__ "MatSolve_Matlab_QR"
201 PetscErrorCode MatSolve_Matlab_QR(Mat A,Vec b,Vec x)
202 {
203   PetscErrorCode ierr;
204   char            *_A,*_b,*_x;
205 
206   PetscFunctionBegin;
207   /* make sure objects have names; use default if not */
208   ierr = PetscObjectName((PetscObject)b);CHKERRQ(ierr);
209   ierr = PetscObjectName((PetscObject)x);CHKERRQ(ierr);
210 
211   ierr = PetscObjectGetName((PetscObject)A,&_A);CHKERRQ(ierr);
212   ierr = PetscObjectGetName((PetscObject)b,&_b);CHKERRQ(ierr);
213   ierr = PetscObjectGetName((PetscObject)x,&_x);CHKERRQ(ierr);
214   ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)b);CHKERRQ(ierr);
215   ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = r%s\\(r%s'\\(%s*%s));",_x,_A,_A,_A+1,_b);CHKERRQ(ierr);
216   ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_b);CHKERRQ(ierr);
217   /* ierr = PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(A->comm),stdout);CHKERRQ(ierr);  */
218   ierr = PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)x);CHKERRQ(ierr);
219   PetscFunctionReturn(0);
220 }
221 
222 #undef __FUNCT__
223 #define __FUNCT__ "MatLUFactorNumeric_Matlab_QR"
224 PetscErrorCode MatLUFactorNumeric_Matlab_QR(Mat A,MatFactorInfo *info,Mat *F)
225 {
226   PetscErrorCode ierr;
227   size_t          len;
228   char            *_A,*name;
229 
230   PetscFunctionBegin;
231   ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);CHKERRQ(ierr);
232   _A   = A->name;
233   ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"r_%s = qr(%s');",_A,_A);CHKERRQ(ierr);
234   ierr = PetscStrlen(_A,&len);CHKERRQ(ierr);
235   ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr);
236   sprintf(name,"_%s",_A);
237   ierr = PetscObjectSetName((PetscObject)*F,name);CHKERRQ(ierr);
238   ierr = PetscFree(name);CHKERRQ(ierr);
239   PetscFunctionReturn(0);
240 }
241 
242 #undef __FUNCT__
243 #define __FUNCT__ "MatLUFactorSymbolic_Matlab_QR"
244 PetscErrorCode MatLUFactorSymbolic_Matlab_QR(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F)
245 {
246   PetscErrorCode ierr;
247 
248   PetscFunctionBegin;
249   if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square");
250   ierr                       = MatCreate(A->comm,A->m,A->n,A->m,A->n,F);CHKERRQ(ierr);
251   ierr                       = MatSetType(*F,A->type_name);CHKERRQ(ierr);
252   ierr                       = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr);
253   (*F)->ops->solve           = MatSolve_Matlab_QR;
254   (*F)->ops->lufactornumeric = MatLUFactorNumeric_Matlab_QR;
255   (*F)->factor               = FACTOR_LU;
256   (*F)->assembled            = PETSC_TRUE;  /* required by -ksp_view */
257 
258   PetscFunctionReturn(0);
259 }
260 
261 /* --------------------------------------------------------------------------------*/
262 #undef __FUNCT__
263 #define __FUNCT__ "MatILUDTFactor_Matlab"
264 PetscErrorCode MatILUDTFactor_Matlab(Mat A,IS isrow,IS iscol,MatFactorInfo *info,Mat *F)
265 {
266   PetscErrorCode ierr;
267   size_t     len;
268   char       *_A,*name;
269 
270   PetscFunctionBegin;
271   if (info->dt == PETSC_DEFAULT)      info->dt      = .005;
272   if (info->dtcol == PETSC_DEFAULT)   info->dtcol   = .01;
273   if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square");
274   ierr                       = MatCreate(A->comm,A->m,A->n,A->m,A->n,F);CHKERRQ(ierr);
275   ierr                       = MatSetType(*F,A->type_name);CHKERRQ(ierr);
276   ierr                       = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr);
277   (*F)->ops->solve           = MatSolve_Matlab;
278   (*F)->factor               = FACTOR_LU;
279   ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);CHKERRQ(ierr);
280   _A   = A->name;
281   ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,info->dtcol);CHKERRQ(ierr);
282   ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"[l_%s,u_%s,p_%s] = luinc(%s',info_%s);",_A,_A,_A,_A,_A);CHKERRQ(ierr);
283   ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_A);CHKERRQ(ierr);
284 
285   ierr = PetscStrlen(_A,&len);CHKERRQ(ierr);
286   ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr);
287   sprintf(name,"_%s",_A);
288   ierr = PetscObjectSetName((PetscObject)*F,name);CHKERRQ(ierr);
289   ierr = PetscFree(name);CHKERRQ(ierr);
290   PetscFunctionReturn(0);
291 }
292 
293 #undef __FUNCT__
294 #define __FUNCT__ "MatFactorInfo_Matlab"
295 PetscErrorCode MatFactorInfo_Matlab(Mat A,PetscViewer viewer)
296 {
297   PetscErrorCode ierr;
298 
299   PetscFunctionBegin;
300   ierr = PetscViewerASCIIPrintf(viewer,"Matlab run parameters:  -- not written yet!\n");CHKERRQ(ierr);
301   PetscFunctionReturn(0);
302 }
303 
304 #undef __FUNCT__
305 #define __FUNCT__ "MatView_Matlab"
306 PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer) {
307   PetscErrorCode    ierr;
308   PetscTruth        iascii;
309   PetscViewerFormat format;
310   Mat_Matlab        *lu=(Mat_Matlab*)(A->spptr);
311 
312   PetscFunctionBegin;
313   ierr = (*lu->MatView)(A,viewer);CHKERRQ(ierr);
314   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr);
315   if (iascii) {
316     ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
317     if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
318       ierr = MatFactorInfo_Matlab(A,viewer);
319     }
320   }
321   PetscFunctionReturn(0);
322 }
323 
324 #undef __FUNCT__
325 #define __FUNCT__ "MatDuplicate_Matlab"
326 PetscErrorCode MatDuplicate_Matlab(Mat A, MatDuplicateOption op, Mat *M) {
327   PetscErrorCode ierr;
328   Mat_Matlab *lu=(Mat_Matlab*)A->spptr;
329 
330   PetscFunctionBegin;
331   ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr);
332   ierr = PetscMemcpy((*M)->spptr,lu,sizeof(Mat_Matlab));CHKERRQ(ierr);
333   PetscFunctionReturn(0);
334 }
335 
336 EXTERN_C_BEGIN
337 #undef __FUNCT__
338 #define __FUNCT__ "MatConvert_SeqAIJ_Matlab"
339 PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_Matlab(Mat A,const MatType type,MatReuse reuse,Mat *newmat)
340 {
341   /* This routine is only called to convert to MATMATLAB */
342   /* from MATSEQAIJ, so we will ignore 'MatType type'. */
343   PetscErrorCode ierr;
344   Mat            B=*newmat;
345   Mat_Matlab     *lu;
346   PetscTruth     qr;
347 
348   PetscFunctionBegin;
349   if (reuse == MAT_INITIAL_MATRIX) {
350     ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr);
351   }
352 
353   ierr = PetscNew(Mat_Matlab,&lu);CHKERRQ(ierr);
354   lu->MatDuplicate         = A->ops->duplicate;
355   lu->MatView              = A->ops->view;
356   lu->MatLUFactorSymbolic  = A->ops->lufactorsymbolic;
357   lu->MatILUDTFactor       = A->ops->iludtfactor;
358   lu->MatDestroy           = A->ops->destroy;
359 
360   B->spptr                 = (void*)lu;
361   B->ops->duplicate        = MatDuplicate_Matlab;
362   B->ops->view             = MatView_Matlab;
363   B->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab;
364   B->ops->iludtfactor      = MatILUDTFactor_Matlab;
365   B->ops->destroy          = MatDestroy_Matlab;
366 
367   ierr = PetscOptionsHasName(A->prefix,"-mat_matlab_qr",&qr);CHKERRQ(ierr);
368   if (qr) {
369     B->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab_QR;
370     ierr = PetscLogInfo((0,"MatConvert_SeqAIJ_Matlab:Using Matlab QR with iterative refinement for LU factorization and solves\n"));CHKERRQ(ierr);
371   } else {
372     ierr = PetscLogInfo((0,"MatConvert_SeqAIJ_Matlab:Using Matlab for LU factorizations and solves.\n"));CHKERRQ(ierr);
373   }
374   ierr = PetscLogInfo((0,"MatConvert_SeqAIJ_Matlab:Using Matlab for ILUDT factorizations and solves.\n"));CHKERRQ(ierr);
375 
376   ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_matlab_C",
377                                            "MatConvert_SeqAIJ_Matlab",MatConvert_SeqAIJ_Matlab);CHKERRQ(ierr);
378   ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_matlab_seqaij_C",
379                                            "MatConvert_Matlab_SeqAIJ",MatConvert_Matlab_SeqAIJ);CHKERRQ(ierr);
380   ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C",
381                                            "MatMatlabEnginePut_Matlab",MatMatlabEnginePut_Matlab);CHKERRQ(ierr);
382   ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C",
383                                            "MatMatlabEngineGet_Matlab",MatMatlabEngineGet_Matlab);CHKERRQ(ierr);
384   ierr = PetscObjectChangeTypeName((PetscObject)B,MATMATLAB);CHKERRQ(ierr);
385   *newmat = B;
386   PetscFunctionReturn(0);
387 }
388 EXTERN_C_END
389 
390 /*MC
391   MATMATLAB - MATMATLAB = "matlab" - A matrix type providing direct solvers (LU and QR) and drop tolerance
392   based ILU factorization (ILUDT) for sequential matrices via the external package Matlab.
393 
394   If Matlab is instaled (see the manual for
395   instructions on how to declare the existence of external packages),
396   a matrix type can be constructed which invokes Matlab solvers.
397   After calling MatCreate(...,A), simply call MatSetType(A,MATMATLAB).
398   This matrix type is only supported for double precision real.
399 
400   This matrix inherits from MATSEQAIJ.  As a result, MatSeqAIJSetPreallocation is
401   supported for this matrix type.  One can also call MatConvert for an inplace conversion to or from
402   the MATSEQAIJ type without data copy.
403 
404   Options Database Keys:
405 + -mat_type matlab - sets the matrix type to "matlab" during a call to MatSetFromOptions()
406 - -mat_matlab_qr   - sets the direct solver to be QR instead of LU
407 
408   Level: beginner
409 
410 .seealso: PCLU
411 M*/
412 
413 EXTERN_C_BEGIN
414 #undef __FUNCT__
415 #define __FUNCT__ "MatCreate_Matlab"
416 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_Matlab(Mat A)
417 {
418   PetscErrorCode ierr;
419 
420   PetscFunctionBegin;
421   ierr = PetscObjectChangeTypeName((PetscObject)A,MATMATLAB);CHKERRQ(ierr);
422   ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr);
423   ierr = MatConvert_SeqAIJ_Matlab(A,MATMATLAB,MAT_REUSE_MATRIX,&A);CHKERRQ(ierr);
424   PetscFunctionReturn(0);
425 }
426 EXTERN_C_END
427