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