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