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 #include "engine.h" /* Matlab include file */ 10 #include "mex.h" /* Matlab include file */ 11 12 typedef struct { 13 int (*MatDuplicate)(Mat,MatDuplicateOption,Mat*); 14 int (*MatView)(Mat,PetscViewer); 15 int (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*); 16 int (*MatILUDTFactor)(Mat,MatFactorInfo*,IS,IS,Mat*); 17 int (*MatDestroy)(Mat); 18 } Mat_Matlab; 19 20 EXTERN_C_BEGIN 21 #undef __FUNCT__ 22 #define __FUNCT__ "MatConvert_Matlab_SeqAIJ" 23 int MatConvert_Matlab_SeqAIJ(Mat A,const MatType type,Mat *newmat) { 24 int ierr; 25 Mat B=*newmat; 26 Mat_Matlab *lu=(Mat_Matlab*)A->spptr; 27 28 PetscFunctionBegin; 29 if (B != A) { 30 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 31 } 32 A->ops->duplicate = lu->MatDuplicate; 33 A->ops->view = lu->MatView; 34 A->ops->lufactorsymbolic = lu->MatLUFactorSymbolic; 35 A->ops->iludtfactor = lu->MatILUDTFactor; 36 A->ops->destroy = lu->MatDestroy; 37 38 ierr = PetscFree(lu);CHKERRQ(ierr); 39 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 40 *newmat = B; 41 PetscFunctionReturn(0); 42 } 43 EXTERN_C_END 44 45 #undef __FUNCT__ 46 #define __FUNCT__ "MatDestroy_Matlab" 47 int MatDestroy_Matlab(Mat A) { 48 int ierr; 49 Mat_Matlab *lu=(Mat_Matlab*)A->spptr; 50 51 PetscFunctionBegin; 52 ierr = MatConvert_Matlab_SeqAIJ(A,MATSEQAIJ,&A);CHKERRQ(ierr); 53 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 54 PetscFunctionReturn(0); 55 } 56 57 #undef __FUNCT__ 58 #define __FUNCT__ "MatSolve_Matlab" 59 int MatSolve_Matlab(Mat A,Vec b,Vec x) 60 { 61 int ierr; 62 char *_A,*_b,*_x; 63 64 PetscFunctionBegin; 65 /* make sure objects have names; use default if not */ 66 ierr = PetscObjectName((PetscObject)b);CHKERRQ(ierr); 67 ierr = PetscObjectName((PetscObject)x);CHKERRQ(ierr); 68 69 ierr = PetscObjectGetName((PetscObject)A,&_A);CHKERRQ(ierr); 70 ierr = PetscObjectGetName((PetscObject)b,&_b);CHKERRQ(ierr); 71 ierr = PetscObjectGetName((PetscObject)x,&_x);CHKERRQ(ierr); 72 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)b);CHKERRQ(ierr); 73 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = u%s\\(l%s\\(p%s*%s));",_x,_A,_A,_A,_b);CHKERRQ(ierr); 74 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_b);CHKERRQ(ierr); 75 /* ierr = PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(A->comm),stdout);CHKERRQ(ierr); */ 76 ierr = PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)x);CHKERRQ(ierr); 77 PetscFunctionReturn(0); 78 } 79 80 #undef __FUNCT__ 81 #define __FUNCT__ "MatLUFactorNumeric_Matlab" 82 int MatLUFactorNumeric_Matlab(Mat A,Mat *F) 83 { 84 Mat_SeqAIJ *f = (Mat_SeqAIJ*)(*F)->data; 85 int ierr,len; 86 char *_A,*name; 87 88 PetscFunctionBegin; 89 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);CHKERRQ(ierr); 90 _A = A->name; 91 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); 92 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_A);CHKERRQ(ierr); 93 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 94 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 95 sprintf(name,"_%s",_A); 96 ierr = PetscObjectSetName((PetscObject)*F,name);CHKERRQ(ierr); 97 ierr = PetscFree(name);CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 101 #undef __FUNCT__ 102 #define __FUNCT__ "MatLUFactorSymbolic_Matlab" 103 int MatLUFactorSymbolic_Matlab(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 104 { 105 int ierr; 106 Mat_SeqAIJ *f; 107 108 PetscFunctionBegin; 109 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 110 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,F);CHKERRQ(ierr); 111 ierr = MatSetType(*F,A->type_name);CHKERRQ(ierr); 112 ierr = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr); 113 (*F)->ops->solve = MatSolve_Matlab; 114 (*F)->ops->lufactornumeric = MatLUFactorNumeric_Matlab; 115 (*F)->factor = FACTOR_LU; 116 f = (Mat_SeqAIJ*)(*F)->data; 117 f->lu_dtcol = info->dtcol; 118 PetscFunctionReturn(0); 119 } 120 121 /* ---------------------------------------------------------------------------------*/ 122 #undef __FUNCT__ 123 #define __FUNCT__ "MatSolve_Matlab_QR" 124 int MatSolve_Matlab_QR(Mat A,Vec b,Vec x) 125 { 126 int ierr; 127 char *_A,*_b,*_x; 128 129 PetscFunctionBegin; 130 /* make sure objects have names; use default if not */ 131 ierr = PetscObjectName((PetscObject)b);CHKERRQ(ierr); 132 ierr = PetscObjectName((PetscObject)x);CHKERRQ(ierr); 133 134 ierr = PetscObjectGetName((PetscObject)A,&_A);CHKERRQ(ierr); 135 ierr = PetscObjectGetName((PetscObject)b,&_b);CHKERRQ(ierr); 136 ierr = PetscObjectGetName((PetscObject)x,&_x);CHKERRQ(ierr); 137 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)b);CHKERRQ(ierr); 138 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = r%s\\(r%s'\\(%s*%s));",_x,_A,_A,_A+1,_b);CHKERRQ(ierr); 139 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_b);CHKERRQ(ierr); 140 /* ierr = PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(A->comm),stdout);CHKERRQ(ierr); */ 141 ierr = PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)x);CHKERRQ(ierr); 142 PetscFunctionReturn(0); 143 } 144 145 #undef __FUNCT__ 146 #define __FUNCT__ "MatLUFactorNumeric_Matlab_QR" 147 int MatLUFactorNumeric_Matlab_QR(Mat A,Mat *F) 148 { 149 int ierr,len; 150 char *_A,*name; 151 152 PetscFunctionBegin; 153 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);CHKERRQ(ierr); 154 _A = A->name; 155 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"r_%s = qr(%s');",_A,_A);CHKERRQ(ierr); 156 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 157 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 158 sprintf(name,"_%s",_A); 159 ierr = PetscObjectSetName((PetscObject)*F,name);CHKERRQ(ierr); 160 ierr = PetscFree(name);CHKERRQ(ierr); 161 PetscFunctionReturn(0); 162 } 163 164 #undef __FUNCT__ 165 #define __FUNCT__ "MatLUFactorSymbolic_Matlab_QR" 166 int MatLUFactorSymbolic_Matlab_QR(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 167 { 168 int ierr; 169 170 PetscFunctionBegin; 171 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 172 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,F);CHKERRQ(ierr); 173 ierr = MatSetType(*F,A->type_name);CHKERRQ(ierr); 174 ierr = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr); 175 (*F)->ops->solve = MatSolve_Matlab_QR; 176 (*F)->ops->lufactornumeric = MatLUFactorNumeric_Matlab_QR; 177 (*F)->factor = FACTOR_LU; 178 (*F)->assembled = PETSC_TRUE; /* required by -ksp_view */ 179 180 PetscFunctionReturn(0); 181 } 182 183 /* --------------------------------------------------------------------------------*/ 184 #undef __FUNCT__ 185 #define __FUNCT__ "MatILUDTFactor_Matlab" 186 int MatILUDTFactor_Matlab(Mat A,MatFactorInfo *info,IS isrow,IS iscol,Mat *F) 187 { 188 int ierr,len; 189 char *_A,*name; 190 191 PetscFunctionBegin; 192 if (info->dt == PETSC_DEFAULT) info->dt = .005; 193 if (info->dtcol == PETSC_DEFAULT) info->dtcol = .01; 194 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 195 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,F);CHKERRQ(ierr); 196 ierr = MatSetType(*F,A->type_name);CHKERRQ(ierr); 197 ierr = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr); 198 (*F)->ops->solve = MatSolve_SeqAIJ_Matlab; 199 (*F)->factor = FACTOR_LU; 200 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);CHKERRQ(ierr); 201 _A = A->name; 202 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,info->dtcol);CHKERRQ(ierr); 203 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"[l_%s,u_%s,p_%s] = luinc(%s',info_%s);",_A,_A,_A,_A,_A);CHKERRQ(ierr); 204 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_A);CHKERRQ(ierr); 205 206 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 207 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 208 sprintf(name,"_%s",_A); 209 ierr = PetscObjectSetName((PetscObject)*F,name);CHKERRQ(ierr); 210 ierr = PetscFree(name);CHKERRQ(ierr); 211 PetscFunctionReturn(0); 212 } 213 214 #undef __FUNCT__ 215 #define __FUNCT__ "MatFactorInfo_Matlab" 216 int MatFactorInfo_Matlab(Mat A,PetscViewer viewer) 217 { 218 int ierr; 219 220 PetscFunctionBegin; 221 ierr = PetscViewerASCIIPrintf(viewer,"Matlab run parameters: -- not written yet!\n");CHKERRQ(ierr); 222 PetscFunctionReturn(0); 223 } 224 225 #undef __FUNCT__ 226 #define __FUNCT__ "MatView_Matlab" 227 int MatView_Matlab(Mat A,PetscViewer viewer) { 228 int ierr; 229 PetscTruth isascii; 230 PetscViewerFormat format; 231 Mat_Matlab *lu=(Mat_Matlab*)(A->spptr); 232 233 PetscFunctionBegin; 234 ierr = (*lu->MatView)(A,viewer);CHKERRQ(ierr); 235 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 236 if (isascii) { 237 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 238 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 239 ierr = MatFactorInfo_Matlab(A,viewer); 240 } 241 } 242 PetscFunctionReturn(0); 243 } 244 245 #undef __FUNCT__ 246 #define __FUNCT__ "MatConvert_SeqAIJ_Matlab" 247 int MatConvert_SeqAIJ_Matlab(Mat A,const MatType type,Mat *newmat) { 248 /* This routine is only called to convert to MATMATLAB */ 249 /* from MATSEQAIJ, so we will ignore 'MatType type'. */ 250 int ierr; 251 Mat B=*newmat; 252 Mat_Matlab *lu; 253 PetscTruth qr; 254 255 PetscFunctionBegin; 256 if (B != A) { 257 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 258 } 259 260 ierr = PetscNew(Mat_Matlab,&lu);CHKERRQ(ierr); 261 lu->MatDuplicate = A->ops->duplicate; 262 lu->MatView = A->ops->view; 263 lu->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 264 lu->MatILUDTFactor = A->ops->iludtfactor; 265 lu->MatDestroy = A->ops->destroy; 266 267 B->spptr = (void*)lu; 268 B->ops->duplicate = MatDuplicate_Matlab; 269 B->ops->view = MatView_Matlab; 270 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab; 271 B->ops->iludtfactor = MatILUDTFactor_Matlab; 272 B->ops->destroy = MatDestroy_Matlab; 273 274 ierr = PetscOptionsHasName(A->prefix,"-mat_matlab_qr",&qr);CHKERRQ(ierr); 275 if (qr) { 276 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab_QR; 277 PetscLogInfo(0,"Using Matlab QR with iterative refinement for LU factorization and solves"); 278 } else { 279 PetscLogInfo(0,"Using Matlab for LU factorizations and solves."); 280 } 281 PetscLogInfo(0,"Using Matlab for ILUDT factorizations and solves."); 282 283 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_matlab_C", 284 "MatConvert_SeqAIJ_Matlab",MatConvert_SeqAIJ_Matlab);CHKERRQ(ierr); 285 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_matlab_seqaij_C", 286 "MatConvert_Matlab_SeqAIJ",MatConvert_Matlab_SeqAIJ);CHKERRQ(ierr); 287 ierr = PetscObjectChangeTypeName((PetscObject)B,MATMATLAB);CHKERRQ(ierr); 288 *newmat = B; 289 PetscFunctionReturn(0); 290 } 291 EXTERN_C_END 292 293 #undef __FUNCT__ 294 #define __FUNCT__ "MatDuplicate_Matlab" 295 int MatDuplicate_Matlab(Mat A, MatDuplicateOption op, Mat *M) { 296 int ierr; 297 Mat_Matlab *lu=(Mat_Matlab*)A->spptr; 298 299 PetscFunctionBegin; 300 ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr); 301 ierr = PetscMemcpy((*M)->spptr,lu,sizeof(Mat_Matlab));CHKERRQ(ierr); 302 PetscFunctionReturn(0); 303 } 304 305 /*MC 306 MATMATLAB - MATMATLAB = "matlab" - A matrix type providing direct solvers (LU and QR) and drop tolerance 307 based ILU factorization (ILUDT) for sequential matrices via the external package Matlab. 308 309 If Matlab is instaled (see the manual for 310 instructions on how to declare the existence of external packages), 311 a matrix type can be constructed which invokes Matlab solvers. 312 After calling MatCreate(...,A), simply call MatSetType(A,MATMATLAB). 313 This matrix type is only supported for double precision real. 314 315 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 316 supported for this matrix type. One can also call MatConvert for an inplace conversion to or from 317 the MATSEQAIJ type without data copy. 318 319 Options Database Keys: 320 + -mat_type matlab - sets the matrix type to "matlab" during a call to MatSetFromOptions() 321 - -mat_matlab_qr - sets the direct solver to be QR instead of LU 322 323 Level: beginner 324 325 .seealso: PCLU 326 M*/ 327 328 EXTERN_C_BEGIN 329 #undef __FUNCT__ 330 #define __FUNCT__ "MatCreate_Matlab" 331 int MarCreate_Matlab(Mat A) { 332 int ierr; 333 334 PetscFunctionBegin; 335 ierr = PetscObjectChangeTypeName((PetscObject)A,MATMATLAB);CHKERRQ(ierr); 336 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 337 ierr = MatConvert_SeqAIJ_Matlab(A,MATMATLAB,&A);CHKERRQ(ierr); 338 PetscFunctionReturn(0); 339 } 340 EXTERN_C_END 341