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