1 /*$Id: mumps.c,v 1.10 2001/08/15 15:56:50 bsmith Exp $*/ 2 /* 3 Provides an interface to the MUMPS_4.2_beta sparse solver 4 */ 5 6 #include "src/mat/impls/aij/seq/aij.h" 7 #include "src/mat/impls/aij/mpi/mpiaij.h" 8 #include "src/mat/impls/sbaij/seq/sbaij.h" 9 #include "src/mat/impls/sbaij/mpi/mpisbaij.h" 10 11 EXTERN_C_BEGIN 12 #if defined(PETSC_USE_COMPLEX) 13 #include "zmumps_c.h" 14 #else 15 #include "dmumps_c.h" 16 #endif 17 EXTERN_C_END 18 #define JOB_INIT -1 19 #define JOB_END -2 20 /* macros s.t. indices match MUMPS documentation */ 21 #define ICNTL(I) icntl[(I)-1] 22 #define CNTL(I) cntl[(I)-1] 23 #define INFOG(I) infog[(I)-1] 24 #define RINFOG(I) rinfog[(I)-1] 25 26 typedef struct { 27 #if defined(PETSC_USE_COMPLEX) 28 ZMUMPS_STRUC_C id; 29 #else 30 DMUMPS_STRUC_C id; 31 #endif 32 MatStructure matstruc; 33 int myid,size,*irn,*jcn,sym; 34 PetscScalar *val; 35 MPI_Comm comm_mumps; 36 PetscTruth isAIJ,CleanUpMUMPS; 37 int (*MatDestroy)(Mat); 38 int (*MatAssemblyEnd)(Mat,MatAssemblyType); 39 int (*MatView)(Mat,PetscViewer); 40 } Mat_AIJ_MUMPS; 41 42 EXTERN int MatFactorInfo_MUMPS(Mat,PetscViewer); 43 44 /* convert Petsc mpiaij matrix to triples: row[nz], col[nz], val[nz] */ 45 /* 46 input: 47 A - matrix in mpiaij format 48 shift - 0: C style output triple; 1: Fortran style output triple. 49 valOnly - FALSE: spaces are allocated and values are set for the triple 50 TRUE: only the values in v array are updated 51 output: 52 nnz - dim of r, c, and v (number of local nonzero entries of A) 53 r, c, v - row and col index, matrix values (matrix triples) 54 */ 55 int MatConvertToTriples(Mat A,int shift,PetscTruth valOnly,int *nnz,int **r, int **c, PetscScalar **v) 56 { 57 int *ai, *aj, *bi, *bj, rstart,nz, *garray; 58 int ierr,i,j,jj,jB,irow,m=A->m,*ajj,*bjj,countA,countB,colA_start,jcol; 59 int *row,*col,*diagA; 60 PetscScalar *av, *bv,*val; 61 Mat_AIJ_MUMPS *mumps = (Mat_AIJ_MUMPS *)A->spptr; 62 PetscTruth isAIJ; 63 64 PetscFunctionBegin; 65 66 if (mumps->isAIJ){ 67 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 68 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)(mat->A)->data; 69 Mat_SeqAIJ *bb=(Mat_SeqAIJ*)(mat->B)->data; 70 nz = aa->nz + bb->nz; 71 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= mat->rstart; 72 garray = mat->garray; 73 av=aa->a; bv=bb->a; 74 75 } else { 76 Mat_MPISBAIJ *mat = (Mat_MPISBAIJ*)A->data; 77 if (mat->bs > 1) SETERRQ1(PETSC_ERR_SUP," bs=%d is not supported yet\n", mat->bs); 78 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)(mat->A)->data; 79 Mat_SeqBAIJ *bb=(Mat_SeqBAIJ*)(mat->B)->data; 80 nz = aa->s_nz + bb->nz; 81 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= mat->rstart; 82 garray = mat->garray; 83 av=aa->a; bv=bb->a; 84 } 85 86 if (!valOnly){ 87 ierr = PetscMalloc(nz*sizeof(int),&row);CHKERRQ(ierr); 88 ierr = PetscMalloc(nz*sizeof(int),&col);CHKERRQ(ierr); 89 ierr = PetscMalloc(nz*sizeof(PetscScalar),&val);CHKERRQ(ierr); 90 *r = row; *c = col; *v = val; 91 } else { 92 row = *r; col = *c; val = *v; 93 } 94 *nnz = nz; 95 96 jj = 0; jB = 0; irow = rstart; 97 for ( i=0; i<m; i++ ) { 98 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 99 countA = ai[i+1] - ai[i]; 100 countB = bi[i+1] - bi[i]; 101 bjj = bj + bi[i]; 102 103 /* get jB, the starting local col index for the 2nd B-part */ 104 colA_start = rstart + ajj[0]; /* the smallest col index for A */ 105 for (j=0; j<countB; j++){ 106 jcol = garray[bjj[j]]; 107 if (jcol > colA_start) { jB = j; break; } 108 if (j==countB-1) jB = countB; 109 } 110 111 /* B-part, smaller col index */ 112 colA_start = rstart + ajj[0]; /* the smallest col index for A */ 113 for (j=0; j<jB; j++){ 114 jcol = garray[bjj[j]]; 115 if (!valOnly){ 116 row[jj] = irow + shift; col[jj] = jcol + shift; 117 } 118 val[jj++] = *bv++; 119 } 120 /* A-part */ 121 for (j=0; j<countA; j++){ 122 if (!valOnly){ 123 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 124 } 125 val[jj++] = *av++; 126 } 127 /* B-part, larger col index */ 128 for (j=jB; j<countB; j++){ 129 if (!valOnly){ 130 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 131 } 132 val[jj++] = *bv++; 133 } 134 irow++; 135 } 136 137 PetscFunctionReturn(0); 138 } 139 140 #undef __FUNCT__ 141 #define __FUNCT__ "MatDestroy_AIJ_MUMPS" 142 int MatDestroy_AIJ_MUMPS(Mat A) 143 { 144 Mat_AIJ_MUMPS *lu = (Mat_AIJ_MUMPS*)A->spptr; 145 int ierr,size=lu->size,(*destroy)(Mat)=lu->MatDestroy; 146 147 PetscFunctionBegin; 148 149 if (lu->CleanUpMUMPS) { 150 /* Terminate instance, deallocate memories */ 151 lu->id.job=JOB_END; 152 #if defined(PETSC_USE_COMPLEX) 153 zmumps_c(&lu->id); 154 #else 155 dmumps_c(&lu->id); 156 #endif 157 if (lu->irn) { ierr = PetscFree(lu->irn);CHKERRQ(ierr);} 158 if (lu->jcn) { ierr = PetscFree(lu->jcn);CHKERRQ(ierr);} 159 if (size>1 && lu->val) { ierr = PetscFree(lu->val);CHKERRQ(ierr);} 160 161 ierr = MPI_Comm_free(&(lu->comm_mumps));CHKERRQ(ierr); 162 } 163 164 ierr = PetscFree(lu);CHKERRQ(ierr); 165 ierr = (*destroy)(A);CHKERRQ(ierr); 166 PetscFunctionReturn(0); 167 } 168 169 #undef __FUNCT__ 170 #define __FUNCT__ "MatAssemblyEnd_AIJ_MUMPS" 171 int MatAssemblyEnd_AIJ_MUMPS(Mat A,MatAssemblyType mode) { 172 int ierr; 173 Mat_AIJ_MUMPS *mumps=(Mat_AIJ_MUMPS*)A->spptr; 174 175 PetscFunctionBegin; 176 ierr = (*mumps->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 177 ierr = MatUseMUMPS_MPIAIJ(A);CHKERRQ(ierr); 178 PetscFunctionReturn(0); 179 } 180 181 #undef __FUNCT__ 182 #define __FUNCT__ "MatView_AIJ_MUMPS" 183 int MatView_AIJ_MUMPS(Mat A,PetscViewer viewer) { 184 int ierr; 185 PetscTruth isascii; 186 PetscViewerFormat format; 187 Mat_AIJ_MUMPS *mumps=(Mat_AIJ_MUMPS*)(A->spptr); 188 189 PetscFunctionBegin; 190 ierr = (*mumps->MatView)(A,viewer);CHKERRQ(ierr); 191 192 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 193 if (isascii) { 194 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 195 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 196 ierr = MatFactorInfo_MUMPS(A,viewer);CHKERRQ(ierr); 197 } 198 } 199 PetscFunctionReturn(0); 200 } 201 202 #undef __FUNCT__ 203 #define __FUNCT__ "MatSolve_AIJ_MUMPS" 204 int MatSolve_AIJ_MUMPS(Mat A,Vec b,Vec x) 205 { 206 Mat_AIJ_MUMPS *lu = (Mat_AIJ_MUMPS*)A->spptr; 207 PetscScalar *rhs,*array; 208 Vec x_seq; 209 IS iden; 210 VecScatter scat; 211 int ierr; 212 213 PetscFunctionBegin; 214 if (lu->size > 1){ 215 if (!lu->myid){ 216 ierr = VecCreateSeq(PETSC_COMM_SELF,A->N,&x_seq);CHKERRQ(ierr); 217 ierr = ISCreateStride(PETSC_COMM_SELF,A->N,0,1,&iden);CHKERRQ(ierr); 218 } else { 219 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&x_seq);CHKERRQ(ierr); 220 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&iden);CHKERRQ(ierr); 221 } 222 ierr = VecScatterCreate(b,iden,x_seq,iden,&scat);CHKERRQ(ierr); 223 ierr = ISDestroy(iden);CHKERRQ(ierr); 224 225 ierr = VecScatterBegin(b,x_seq,INSERT_VALUES,SCATTER_FORWARD,scat);CHKERRQ(ierr); 226 ierr = VecScatterEnd(b,x_seq,INSERT_VALUES,SCATTER_FORWARD,scat);CHKERRQ(ierr); 227 if (!lu->myid) {ierr = VecGetArray(x_seq,&array);CHKERRQ(ierr);} 228 } else { /* size == 1 */ 229 ierr = VecCopy(b,x);CHKERRQ(ierr); 230 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 231 } 232 if (!lu->myid) { /* define rhs on the host */ 233 #if defined(PETSC_USE_COMPLEX) 234 lu->id.rhs = (mumps_double_complex*)array; 235 #else 236 lu->id.rhs = array; 237 #endif 238 } 239 240 /* solve phase */ 241 lu->id.job=3; 242 #if defined(PETSC_USE_COMPLEX) 243 zmumps_c(&lu->id); 244 #else 245 dmumps_c(&lu->id); 246 #endif 247 if (lu->id.INFOG(1) < 0) { 248 SETERRQ1(1,"Error reported by MUMPS in solve phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 249 } 250 251 /* convert mumps solution x_seq to petsc mpi x */ 252 if (lu->size > 1) { 253 if (!lu->myid){ 254 ierr = VecRestoreArray(x_seq,&array);CHKERRQ(ierr); 255 } 256 ierr = VecScatterBegin(x_seq,x,INSERT_VALUES,SCATTER_REVERSE,scat);CHKERRQ(ierr); 257 ierr = VecScatterEnd(x_seq,x,INSERT_VALUES,SCATTER_REVERSE,scat);CHKERRQ(ierr); 258 ierr = VecScatterDestroy(scat);CHKERRQ(ierr); 259 ierr = VecDestroy(x_seq);CHKERRQ(ierr); 260 } else { 261 ierr = VecRestoreArray(x,&array);CHKERRQ(ierr); 262 } 263 264 PetscFunctionReturn(0); 265 } 266 267 #undef __FUNCT__ 268 #define __FUNCT__ "MatFactorNumeric_MPIAIJ_MUMPS" 269 int MatFactorNumeric_AIJ_MUMPS(Mat A,Mat *F) 270 { 271 Mat_AIJ_MUMPS *lu = (Mat_AIJ_MUMPS*)(*F)->spptr; 272 int rnz,nnz,ierr,nz,i,M=A->M,*ai,*aj,icntl; 273 PetscTruth valOnly,flg; 274 275 PetscFunctionBegin; 276 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 277 (*F)->ops->solve = MatSolve_AIJ_MUMPS; 278 279 /* Initialize a MUMPS instance */ 280 ierr = MPI_Comm_rank(A->comm, &lu->myid); 281 ierr = MPI_Comm_size(A->comm,&lu->size);CHKERRQ(ierr); 282 lu->id.job = JOB_INIT; 283 ierr = MPI_Comm_dup(A->comm,&(lu->comm_mumps));CHKERRQ(ierr); 284 lu->id.comm_fortran = lu->comm_mumps; 285 286 /* Set mumps options */ 287 ierr = PetscOptionsBegin(A->comm,A->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 288 lu->id.par=1; /* host participates factorizaton and solve */ 289 lu->id.sym=lu->sym; 290 if (lu->sym == 2){ 291 ierr = PetscOptionsInt("-mat_mumps_sym","SYM: (1,2)","None",lu->id.sym,&icntl,&flg);CHKERRQ(ierr); 292 if (flg && icntl == 1) lu->id.sym=icntl; /* matrix is spd */ 293 } 294 #if defined(PETSC_USE_COMPLEX) 295 zmumps_c(&lu->id); 296 #else 297 dmumps_c(&lu->id); 298 #endif 299 300 if (lu->size == 1){ 301 lu->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 302 } else { 303 lu->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 304 } 305 306 icntl=-1; 307 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",lu->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 308 if (flg && icntl > 0) { 309 lu->id.ICNTL(4)=icntl; /* and use mumps default icntl(i), i=1,2,3 */ 310 } else { /* no output */ 311 lu->id.ICNTL(1) = 0; /* error message, default= 6 */ 312 lu->id.ICNTL(2) = -1; /* output stream for diagnostic printing, statistics, and warning. default=0 */ 313 lu->id.ICNTL(3) = -1; /* output stream for global information, default=6 */ 314 lu->id.ICNTL(4) = 0; /* level of printing, 0,1,2,3,4, default=2 */ 315 } 316 ierr = PetscOptionsInt("-mat_mumps_icntl_6","ICNTL(6): matrix prescaling (0 to 7)","None",lu->id.ICNTL(6),&lu->id.ICNTL(6),PETSC_NULL);CHKERRQ(ierr); 317 icntl=-1; 318 ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): matrix ordering (0 to 7)","None",lu->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr); 319 if (flg) { 320 if (icntl== 1){ 321 SETERRQ(PETSC_ERR_SUP,"pivot order be set by the user in PERM_IN -- not supported by the PETSc/MUMPS interface\n"); 322 } else { 323 lu->id.ICNTL(7) = icntl; 324 } 325 } 326 ierr = PetscOptionsInt("-mat_mumps_icntl_9","ICNTL(9): A or A^T x=b to be solved. 1: A; otherwise: A^T","None",lu->id.ICNTL(9),&lu->id.ICNTL(9),PETSC_NULL);CHKERRQ(ierr); 327 ierr = PetscOptionsInt("-mat_mumps_icntl_10","ICNTL(10): max num of refinements","None",lu->id.ICNTL(10),&lu->id.ICNTL(10),PETSC_NULL);CHKERRQ(ierr); 328 ierr = PetscOptionsInt("-mat_mumps_icntl_11","ICNTL(11): error analysis, a positive value returns statistics (by -sles_view)","None",lu->id.ICNTL(11),&lu->id.ICNTL(11),PETSC_NULL);CHKERRQ(ierr); 329 ierr = PetscOptionsInt("-mat_mumps_icntl_12","ICNTL(12): efficiency control","None",lu->id.ICNTL(12),&lu->id.ICNTL(12),PETSC_NULL);CHKERRQ(ierr); 330 ierr = PetscOptionsInt("-mat_mumps_icntl_13","ICNTL(13): efficiency control","None",lu->id.ICNTL(13),&lu->id.ICNTL(13),PETSC_NULL);CHKERRQ(ierr); 331 ierr = PetscOptionsInt("-mat_mumps_icntl_14","ICNTL(14): efficiency control","None",lu->id.ICNTL(14),&lu->id.ICNTL(14),PETSC_NULL);CHKERRQ(ierr); 332 ierr = PetscOptionsInt("-mat_mumps_icntl_15","ICNTL(15): efficiency control","None",lu->id.ICNTL(15),&lu->id.ICNTL(15),PETSC_NULL);CHKERRQ(ierr); 333 334 /* 335 ierr = PetscOptionsInt("-mat_mumps_icntl_16","ICNTL(16): 1: rank detection; 2: rank detection and nullspace","None",lu->id.ICNTL(16),&icntl,&flg);CHKERRQ(ierr); 336 if (flg){ 337 if (icntl >-1 && icntl <3 ){ 338 if (lu->myid==0) lu->id.ICNTL(16) = icntl; 339 } else { 340 SETERRQ1(PETSC_ERR_SUP,"ICNTL(16)=%d -- not supported\n",icntl); 341 } 342 } 343 */ 344 345 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",lu->id.CNTL(1),&lu->id.CNTL(1),PETSC_NULL);CHKERRQ(ierr); 346 ierr = PetscOptionsReal("-mat_mumps_cntl_2","CNTL(2): stopping criterion of refinement","None",lu->id.CNTL(2),&lu->id.CNTL(2),PETSC_NULL);CHKERRQ(ierr); 347 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",lu->id.CNTL(3),&lu->id.CNTL(3),PETSC_NULL);CHKERRQ(ierr); 348 PetscOptionsEnd(); 349 } 350 351 /* define matrix A */ 352 switch (lu->id.ICNTL(18)){ 353 case 0: /* centralized assembled matrix input (size=1) */ 354 if (!lu->myid) { 355 if (lu->isAIJ){ 356 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data; 357 nz = aa->nz; 358 ai = aa->i; aj = aa->j; lu->val = aa->a; 359 } else { 360 Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)A->data; 361 nz = aa->s_nz; 362 ai = aa->i; aj = aa->j; lu->val = aa->a; 363 } 364 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ /* first numeric factorization, get irn and jcn */ 365 ierr = PetscMalloc(nz*sizeof(int),&lu->irn);CHKERRQ(ierr); 366 ierr = PetscMalloc(nz*sizeof(int),&lu->jcn);CHKERRQ(ierr); 367 nz = 0; 368 for (i=0; i<M; i++){ 369 rnz = ai[i+1] - ai[i]; 370 while (rnz--) { /* Fortran row/col index! */ 371 lu->irn[nz] = i+1; lu->jcn[nz] = (*aj)+1; aj++; nz++; 372 } 373 } 374 } 375 } 376 break; 377 case 3: /* distributed assembled matrix input (size>1) */ 378 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 379 valOnly = PETSC_FALSE; 380 } else { 381 valOnly = PETSC_TRUE; /* only update mat values, not row and col index */ 382 } 383 ierr = MatConvertToTriples(A,1,valOnly, &nnz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 384 break; 385 default: SETERRQ(PETSC_ERR_SUP,"Matrix input format is not supported by MUMPS."); 386 } 387 388 /* analysis phase */ 389 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 390 lu->id.n = M; 391 switch (lu->id.ICNTL(18)){ 392 case 0: /* centralized assembled matrix input */ 393 if (!lu->myid) { 394 lu->id.nz =nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 395 if (lu->id.ICNTL(6)>1){ 396 #if defined(PETSC_USE_COMPLEX) 397 lu->id.a = (mumps_double_complex*)lu->val; 398 #else 399 lu->id.a = lu->val; 400 #endif 401 } 402 } 403 break; 404 case 3: /* distributed assembled matrix input (size>1) */ 405 lu->id.nz_loc = nnz; 406 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 407 if (lu->id.ICNTL(6)>1) { 408 #if defined(PETSC_USE_COMPLEX) 409 lu->id.a_loc = (mumps_double_complex*)lu->val; 410 #else 411 lu->id.a_loc = lu->val; 412 #endif 413 } 414 break; 415 } 416 lu->id.job=1; 417 #if defined(PETSC_USE_COMPLEX) 418 zmumps_c(&lu->id); 419 #else 420 dmumps_c(&lu->id); 421 #endif 422 if (lu->id.INFOG(1) < 0) { 423 SETERRQ1(1,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 424 } 425 } 426 427 /* numerical factorization phase */ 428 if(lu->id.ICNTL(18) == 0) { 429 if (lu->myid == 0) { 430 #if defined(PETSC_USE_COMPLEX) 431 lu->id.a = (mumps_double_complex*)lu->val; 432 #else 433 lu->id.a = lu->val; 434 #endif 435 } 436 } else { 437 #if defined(PETSC_USE_COMPLEX) 438 lu->id.a_loc = (mumps_double_complex*)lu->val; 439 #else 440 lu->id.a_loc = lu->val; 441 #endif 442 } 443 lu->id.job=2; 444 #if defined(PETSC_USE_COMPLEX) 445 zmumps_c(&lu->id); 446 #else 447 dmumps_c(&lu->id); 448 #endif 449 if (lu->id.INFOG(1) < 0) { 450 SETERRQ1(1,"1, Error reported by MUMPS in numerical factorization phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 451 } 452 453 if (lu->myid==0 && lu->id.ICNTL(16) > 0){ 454 SETERRQ1(1," lu->id.ICNTL(16):=%d\n",lu->id.INFOG(16)); 455 } 456 457 (*F)->assembled = PETSC_TRUE; 458 lu->matstruc = SAME_NONZERO_PATTERN; 459 PetscFunctionReturn(0); 460 } 461 462 /* Note the Petsc r and c permutations are ignored */ 463 #undef __FUNCT__ 464 #define __FUNCT__ "MatLUFactorSymbolic_AIJ_MUMPS" 465 int MatLUFactorSymbolic_AIJ_MUMPS(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 466 { 467 Mat B; 468 Mat_AIJ_MUMPS *lu; 469 int ierr; 470 471 PetscFunctionBegin; 472 473 /* Create the factorization matrix */ 474 ierr = MatCreate(A->comm,A->m,A->n,A->M,A->N,&B);CHKERRQ(ierr); 475 ierr = MatSetType(B,MATAIJMUMPS);CHKERRQ(ierr); 476 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 477 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 478 479 B->ops->lufactornumeric = MatFactorNumeric_AIJ_MUMPS; 480 B->factor = FACTOR_LU; 481 lu = (Mat_AIJ_MUMPS*)B->spptr; 482 lu->sym = 0; 483 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 484 485 *F = B; 486 PetscFunctionReturn(0); 487 } 488 489 /* Note the Petsc r permutation is ignored */ 490 #undef __FUNCT__ 491 #define __FUNCT__ "MatCholeskyFactorSymbolic_AIJ_MUMPS" 492 int MatCholeskyFactorSymbolic_AIJ_MUMPS(Mat A,IS r,MatFactorInfo *info,Mat *F) 493 { 494 Mat B; 495 Mat_AIJ_MUMPS *lu; 496 int ierr; 497 498 PetscFunctionBegin; 499 500 /* Create the factorization matrix */ 501 ierr = MatCreate(A->comm,A->m,A->n,A->M,A->N,&B);CHKERRQ(ierr); 502 ierr = MatSetType(B,MATAIJMUMPS);CHKERRQ(ierr); 503 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 504 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 505 506 B->ops->choleskyfactornumeric = MatFactorNumeric_AIJ_MUMPS; 507 B->factor = FACTOR_CHOLESKY; 508 lu = (Mat_AIJ_MUMPS *)B->spptr; 509 lu->sym = 2; 510 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 511 512 *F = B; 513 PetscFunctionReturn(0); 514 } 515 516 #undef __FUNCT__ 517 #define __FUNCT__ "MatUseMUMPS_AIJ" 518 int MatUseMUMPS_AIJ(Mat A) 519 { 520 PetscFunctionBegin; 521 A->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJ_MUMPS; 522 A->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJ_MUMPS; 523 A->ops->lufactornumeric = MatFactorNumeric_AIJ_MUMPS; 524 525 PetscFunctionReturn(0); 526 } 527 528 int MatFactorInfo_MUMPS(Mat A,PetscViewer viewer) 529 { 530 Mat_AIJ_MUMPS *lu= (Mat_AIJ_MUMPS*)A->spptr; 531 int ierr; 532 533 PetscFunctionBegin; 534 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 535 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",lu->id.sym);CHKERRQ(ierr); 536 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",lu->id.par);CHKERRQ(ierr); 537 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",lu->id.ICNTL(4));CHKERRQ(ierr); 538 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",lu->id.ICNTL(5));CHKERRQ(ierr); 539 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",lu->id.ICNTL(6));CHKERRQ(ierr); 540 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (matrix ordering): %d \n",lu->id.ICNTL(7));CHKERRQ(ierr); 541 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(9) (A/A^T x=b is solved): %d \n",lu->id.ICNTL(9));CHKERRQ(ierr); 542 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr); 543 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",lu->id.ICNTL(11));CHKERRQ(ierr); 544 if (lu->myid == 0 && lu->id.ICNTL(11)>0) { 545 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(4) (inf norm of input mat): %g\n",lu->id.RINFOG(4));CHKERRQ(ierr); 546 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(5) (inf norm of solution): %g\n",lu->id.RINFOG(5));CHKERRQ(ierr); 547 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(6) (inf norm of residual): %g\n",lu->id.RINFOG(6));CHKERRQ(ierr); 548 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(7),RINFOG(8) (backward error est): %g\n",lu->id.RINFOG(7),lu->id.RINFOG(8));CHKERRQ(ierr); 549 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(9) (error estimate): %g \n",lu->id.RINFOG(9));CHKERRQ(ierr); 550 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr); 551 552 } 553 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",lu->id.ICNTL(12));CHKERRQ(ierr); 554 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",lu->id.ICNTL(13));CHKERRQ(ierr); 555 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (efficiency control): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr); 556 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(15) (efficiency control): %d \n",lu->id.ICNTL(15));CHKERRQ(ierr); 557 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",lu->id.ICNTL(18));CHKERRQ(ierr); 558 559 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",lu->id.CNTL(1));CHKERRQ(ierr); 560 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr); 561 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absolute pivoting threshold): %g \n",lu->id.CNTL(3));CHKERRQ(ierr); 562 PetscFunctionReturn(0); 563 } 564 565 EXTERN_C_BEGIN 566 #undef __FUNCT__ 567 #define __FUNCT__ "MatCreate_AIJ_MUMPS" 568 int MatCreate_AIJ_MUMPS(Mat A) { 569 int ierr,size; 570 MPI_Comm comm; 571 Mat_AIJ_MUMPS *mumps; 572 573 PetscFunctionBegin; 574 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 575 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 576 if (size == 1) { 577 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 578 } else { 579 ierr = MatSetType(A,MATMPIAIJ);CHKERRQ(ierr); 580 } 581 ierr = MatUseMUMPS_AIJ(A); 582 583 ierr = PetscNew(Mat_AIJ_MUMPS,&mumps);CHKERRQ(ierr); 584 mumps->MatDestroy = A->ops->destroy; 585 mumps->MatAssemblyEnd = A->ops->assemblyend; 586 mumps->MatView = A->ops->view; 587 mumps->CleanUpMUMPS = PETSC_FALSE; 588 A->spptr = (void *)mumps; 589 A->ops->destroy = MatDestroy_AIJ_MUMPS; 590 A->ops->assemblyend = MatAssemblyEnd_AIJ_MUMPS; 591 A->ops->view = MatView_AIJ_MUMPS; 592 PetscFunctionReturn(0); 593 } 594 EXTERN_C_END 595 596 EXTERN_C_BEGIN 597 #undef __FUNCT__ 598 #define __FUNCT__ "MatCreate_SBAIJ_MUMPS" 599 int MatCreate_SBAIJ_MUMPS(Mat A) { 600 int ierr,size; 601 MPI_Comm comm; 602 Mat_AIJ_MUMPS *mumps; 603 604 PetscFunctionBegin; 605 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 606 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 607 if (size == 1) { 608 ierr = MatSetType(A,MATSEQSBAIJ);CHKERRQ(ierr); 609 } else { 610 ierr = MatSetType(A,MATMPISBAIJ);CHKERRQ(ierr); 611 } 612 ierr=MatUseMUMPS_AIJ(A); 613 614 ierr = PetscNew(Mat_AIJ_MUMPS,&mumps);CHKERRQ(ierr); 615 mumps->MatDestroy = A->ops->destroy; 616 mumps->MatAssemblyEnd = A->ops->assemblyend; 617 mumps->CleanUpMUMPS = PETSC_FALSE; 618 A->spptr = (void *)mumps; 619 A->ops->destroy = MatDestroy_AIJ_MUMPS; 620 A->ops->assemblyend = MatAssemblyEnd_AIJ_MUMPS; 621 622 PetscFunctionReturn(0); 623 } 624 EXTERN_C_END 625 626 EXTERN_C_BEGIN 627 #undef __FUNCT__ 628 #define __FUNCT__ "MatLoad_AIJ_MUMPS" 629 int MatLoad_AIJ_MUMPS(PetscViewer viewer,MatType type,Mat *A) { 630 int ierr,size,(*r)(PetscViewer,MatType,Mat*); 631 MPI_Comm comm; 632 633 PetscFunctionBegin; 634 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 635 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 636 if (size == 1) { 637 ierr = PetscFListFind(comm,MatLoadList,MATSEQAIJ,(void(**)(void))&r);CHKERRQ(ierr); 638 } else { 639 ierr = PetscFListFind(comm,MatLoadList,MATMPIAIJ,(void(**)(void))&r);CHKERRQ(ierr); 640 } 641 ierr = (*r)(viewer,type,A);CHKERRQ(ierr); 642 PetscFunctionReturn(0); 643 } 644 EXTERN_C_END 645 646 EXTERN_C_BEGIN 647 #undef __FUNCT__ 648 #define __FUNCT__ "MatLoad_SBAIJ_MUMPS" 649 int MatLoad_SBAIJ_MUMPS(PetscViewer viewer,MatType type,Mat *A) { 650 int ierr,size,(*r)(PetscViewer,MatType,Mat*); 651 MPI_Comm comm; 652 653 PetscFunctionBegin; 654 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 655 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 656 if (size == 1) { 657 ierr = PetscFListFind(comm,MatLoadList,MATSEQSBAIJ,(void(**)(void))&r);CHKERRQ(ierr); 658 } else { 659 ierr = PetscFListFind(comm,MatLoadList,MATMPISBAIJ,(void(**)(void))&r);CHKERRQ(ierr); 660 } 661 ierr = (*r)(viewer,type,A);CHKERRQ(ierr); 662 PetscFunctionReturn(0); 663 } 664 EXTERN_C_END 665