1 /*$Id: aijmatlab.c,v 1.12 2001/08/06 21:15:14 bsmith Exp $*/ 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 #if defined(PETSC_HAVE_MATLAB) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 10 #include "engine.h" /* Matlab include file */ 11 #include "mex.h" /* Matlab include file */ 12 13 #undef __FUNCT__ 14 #define __FUNCT__ "MatSolve_SeqAIJ_Matlab" 15 int MatSolve_SeqAIJ_Matlab(Mat A,Vec b,Vec x) 16 { 17 int ierr; 18 char *_A,*_b,*_x; 19 20 PetscFunctionBegin; 21 /* make sure objects have names; use default if not */ 22 ierr = PetscObjectName((PetscObject)b);CHKERRQ(ierr); 23 ierr = PetscObjectName((PetscObject)x);CHKERRQ(ierr); 24 25 ierr = PetscObjectGetName((PetscObject)A,&_A);CHKERRQ(ierr); 26 ierr = PetscObjectGetName((PetscObject)b,&_b);CHKERRQ(ierr); 27 ierr = PetscObjectGetName((PetscObject)x,&_x);CHKERRQ(ierr); 28 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)b);CHKERRQ(ierr); 29 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = u%s\\(l%s\\(p%s*%s));",_x,_A,_A,_A,_b);CHKERRQ(ierr); 30 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_b);CHKERRQ(ierr); 31 /* ierr = PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(A->comm),stdout);CHKERRQ(ierr); */ 32 ierr = PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)x);CHKERRQ(ierr); 33 PetscFunctionReturn(0); 34 } 35 36 #undef __FUNCT__ 37 #define __FUNCT__ "MatLUFactorNumeric_SeqAIJ_Matlab" 38 int MatLUFactorNumeric_SeqAIJ_Matlab(Mat A,Mat *F) 39 { 40 Mat_SeqAIJ *f = (Mat_SeqAIJ*)(*F)->data; 41 int ierr,len; 42 char *_A,*name; 43 44 PetscFunctionBegin; 45 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);CHKERRQ(ierr); 46 _A = A->name; 47 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"[l_%s,u_%s,p_%s] = lu(%s',%g);",_A,_A,_A,_A,f->lu_dtcol);CHKERRQ(ierr); 48 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_A);CHKERRQ(ierr); 49 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 50 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 51 sprintf(name,"_%s",_A); 52 ierr = PetscObjectSetName((PetscObject)*F,name);CHKERRQ(ierr); 53 ierr = PetscFree(name);CHKERRQ(ierr); 54 PetscFunctionReturn(0); 55 } 56 57 #undef __FUNCT__ 58 #define __FUNCT__ "MatLUFactorSymbolic_SeqAIJ_Matlab" 59 int MatLUFactorSymbolic_SeqAIJ_Matlab(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 60 { 61 int ierr; 62 Mat_SeqAIJ *f; 63 64 PetscFunctionBegin; 65 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 66 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,F);CHKERRQ(ierr); 67 ierr = MatSetType(*F,A->type_name);CHKERRQ(ierr); 68 ierr = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr); 69 (*F)->ops->solve = MatSolve_SeqAIJ_Matlab; 70 (*F)->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Matlab; 71 (*F)->factor = FACTOR_LU; 72 f = (Mat_SeqAIJ*)(*F)->data; 73 f->lu_dtcol = info->dtcol; 74 PetscFunctionReturn(0); 75 } 76 77 /* ---------------------------------------------------------------------------------*/ 78 #undef __FUNCT__ 79 #define __FUNCT__ "MatSolve_SeqAIJ_Matlab_QR" 80 int MatSolve_SeqAIJ_Matlab_QR(Mat A,Vec b,Vec x) 81 { 82 int ierr; 83 char *_A,*_b,*_x; 84 85 PetscFunctionBegin; 86 /* make sure objects have names; use default if not */ 87 ierr = PetscObjectName((PetscObject)b);CHKERRQ(ierr); 88 ierr = PetscObjectName((PetscObject)x);CHKERRQ(ierr); 89 90 ierr = PetscObjectGetName((PetscObject)A,&_A);CHKERRQ(ierr); 91 ierr = PetscObjectGetName((PetscObject)b,&_b);CHKERRQ(ierr); 92 ierr = PetscObjectGetName((PetscObject)x,&_x);CHKERRQ(ierr); 93 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)b);CHKERRQ(ierr); 94 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = r%s\\(r%s'\\(%s*%s));",_x,_A,_A,_A+1,_b);CHKERRQ(ierr); 95 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_b);CHKERRQ(ierr); 96 /* ierr = PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(A->comm),stdout);CHKERRQ(ierr); */ 97 ierr = PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)x);CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 101 #undef __FUNCT__ 102 #define __FUNCT__ "MatLUFactorNumeric_SeqAIJ_Matlab_QR" 103 int MatLUFactorNumeric_SeqAIJ_Matlab_QR(Mat A,Mat *F) 104 { 105 int ierr,len; 106 char *_A,*name; 107 108 PetscFunctionBegin; 109 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);CHKERRQ(ierr); 110 _A = A->name; 111 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"r_%s = qr(%s');",_A,_A);CHKERRQ(ierr); 112 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 113 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 114 sprintf(name,"_%s",_A); 115 ierr = PetscObjectSetName((PetscObject)*F,name);CHKERRQ(ierr); 116 ierr = PetscFree(name);CHKERRQ(ierr); 117 PetscFunctionReturn(0); 118 } 119 120 #undef __FUNCT__ 121 #define __FUNCT__ "MatLUFactorSymbolic_SeqAIJ_Matlab_QR" 122 int MatLUFactorSymbolic_SeqAIJ_Matlab_QR(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 123 { 124 int ierr; 125 126 PetscFunctionBegin; 127 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 128 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,F);CHKERRQ(ierr); 129 ierr = MatSetType(*F,A->type_name);CHKERRQ(ierr); 130 ierr = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr); 131 (*F)->ops->solve = MatSolve_SeqAIJ_Matlab_QR; 132 (*F)->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Matlab_QR; 133 (*F)->factor = FACTOR_LU; 134 (*F)->assembled = PETSC_TRUE; /* required by -ksp_view */ 135 136 PetscFunctionReturn(0); 137 } 138 139 /* --------------------------------------------------------------------------------*/ 140 #undef __FUNCT__ 141 #define __FUNCT__ "MatILUDTFactor_SeqAIJ_Matlab" 142 int MatILUDTFactor_SeqAIJ_Matlab(Mat A,MatFactorInfo *info,IS isrow,IS iscol,Mat *F) 143 { 144 int ierr,len; 145 char *_A,*name; 146 147 PetscFunctionBegin; 148 if (info->dt == PETSC_DEFAULT) info->dt = .005; 149 if (info->dtcol == PETSC_DEFAULT) info->dtcol = .01; 150 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 151 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,F);CHKERRQ(ierr); 152 ierr = MatSetType(*F,A->type_name);CHKERRQ(ierr); 153 ierr = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr); 154 (*F)->ops->solve = MatSolve_SeqAIJ_Matlab; 155 (*F)->factor = FACTOR_LU; 156 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);CHKERRQ(ierr); 157 _A = A->name; 158 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,info->dtcol);CHKERRQ(ierr); 159 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"[l_%s,u_%s,p_%s] = luinc(%s',info_%s);",_A,_A,_A,_A,_A);CHKERRQ(ierr); 160 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_A);CHKERRQ(ierr); 161 162 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 163 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 164 sprintf(name,"_%s",_A); 165 ierr = PetscObjectSetName((PetscObject)*F,name);CHKERRQ(ierr); 166 ierr = PetscFree(name);CHKERRQ(ierr); 167 PetscFunctionReturn(0); 168 } 169 170 int MatSeqAIJFactorInfo_Matlab(Mat A,PetscViewer viewer) 171 { 172 int ierr; 173 174 PetscFunctionBegin; 175 /* check if matrix is matlab type */ 176 if (A->ops->solve != MatSolve_SeqAIJ_Matlab) PetscFunctionReturn(0); 177 178 ierr = PetscViewerASCIIPrintf(viewer,"Matlab run parameters: -- not written yet!\n");CHKERRQ(ierr); 179 PetscFunctionReturn(0); 180 } 181 182 #undef __FUNCT__ 183 #define __FUNCT__ "MatUseMatlab_SeqAIJ" 184 int MatUseMatlab_SeqAIJ(Mat A) 185 { 186 PetscTruth qr; 187 int ierr; 188 189 PetscFunctionBegin; 190 ierr = PetscOptionsHasName(A->prefix,"-mat_aij_matlab_qr",&qr);CHKERRQ(ierr); 191 if (qr) { 192 A->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ_Matlab_QR; 193 PetscLogInfo(0,"Using Matlab QR with iterative refinement for SeqAIJ LU factorization and solves"); 194 } else { 195 A->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ_Matlab; 196 PetscLogInfo(0,"Using Matlab for SeqAIJ LU factorization and solves"); 197 } 198 A->ops->iludtfactor = MatILUDTFactor_SeqAIJ_Matlab; 199 PetscLogInfo(0,"Using Matlab for SeqAIJ ILUDT factorization and solves"); 200 PetscFunctionReturn(0); 201 } 202 203 #else 204 205 #undef __FUNCT__ 206 #define __FUNCT__ "MatUseMatlab_SeqAIJ" 207 int MatUseMatlab_SeqAIJ(Mat A) 208 { 209 PetscFunctionBegin; 210 PetscFunctionReturn(0); 211 } 212 213 #endif 214 215 216