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