1 /*$Id: superlu_DIST.c,v 1.10 2001/08/15 15:56:50 bsmith Exp $*/ 2 /* 3 Provides an interface to the SuperLU_DIST_2.0 sparse solver 4 */ 5 6 #include "src/mat/impls/aij/seq/aij.h" 7 #include "src/mat/impls/aij/mpi/mpiaij.h" 8 #if defined(PETSC_HAVE_STDLIB_H) /* This is to get arround weird problem with SuperLU on cray */ 9 #include "stdlib.h" 10 #endif 11 12 EXTERN_C_BEGIN 13 #if defined(PETSC_USE_COMPLEX) 14 #include "superlu_zdefs.h" 15 #else 16 #include "superlu_ddefs.h" 17 #endif 18 EXTERN_C_END 19 20 typedef enum { GLOBAL,DISTRIBUTED 21 } SuperLU_MatInputMode; 22 23 typedef struct { 24 int_t nprow,npcol,*row,*col; 25 gridinfo_t grid; 26 superlu_options_t options; 27 SuperMatrix A_sup; 28 ScalePermstruct_t ScalePermstruct; 29 LUstruct_t LUstruct; 30 int StatPrint; 31 int MatInputMode; 32 SOLVEstruct_t SOLVEstruct; 33 MatStructure flg; 34 MPI_Comm comm_superlu; 35 #if defined(PETSC_USE_COMPLEX) 36 doublecomplex *val; 37 #else 38 double *val; 39 #endif 40 41 /* A few function pointers for inheritance */ 42 int (*MatView)(Mat,PetscViewer); 43 int (*MatAssemblyEnd)(Mat,MatAssemblyType); 44 int (*MatDestroy)(Mat); 45 46 /* Flag to clean up (non-global) SuperLU objects during Destroy */ 47 PetscTruth CleanUpSuperLUDist; 48 } Mat_MPIAIJ_SuperLU_DIST; 49 50 #undef __FUNCT__ 51 #define __FUNCT__ "MatDestroy_MPIAIJ_SuperLU_DIST" 52 int MatDestroy_MPIAIJ_SuperLU_DIST(Mat A) 53 { 54 Mat_MPIAIJ *a = (Mat_MPIAIJ*)A->data; 55 Mat_MPIAIJ_SuperLU_DIST *lu = (Mat_MPIAIJ_SuperLU_DIST*)A->spptr; 56 int ierr, size=a->size,(*destroy)(Mat); 57 58 PetscFunctionBegin; 59 if (lu->CleanUpSuperLUDist) { 60 /* Deallocate SuperLU_DIST storage */ 61 if (lu->MatInputMode == GLOBAL) { 62 Destroy_CompCol_Matrix_dist(&lu->A_sup); 63 } else { 64 Destroy_CompRowLoc_Matrix_dist(&lu->A_sup); 65 if ( lu->options.SolveInitialized ) { 66 #if defined(PETSC_USE_COMPLEX) 67 zSolveFinalize(&lu->options, &lu->SOLVEstruct); 68 #else 69 dSolveFinalize(&lu->options, &lu->SOLVEstruct); 70 #endif 71 } 72 } 73 Destroy_LU(A->N, &lu->grid, &lu->LUstruct); 74 ScalePermstructFree(&lu->ScalePermstruct); 75 LUstructFree(&lu->LUstruct); 76 77 /* Release the SuperLU_DIST process grid. */ 78 superlu_gridexit(&lu->grid); 79 80 ierr = MPI_Comm_free(&(lu->comm_superlu));CHKERRQ(ierr); 81 } 82 destroy = lu->MatDestroy; 83 ierr = PetscFree(lu);CHKERRQ(ierr); 84 ierr = (*destroy)(A);CHKERRQ(ierr); 85 86 PetscFunctionReturn(0); 87 } 88 89 #undef __FUNCT__ 90 #define __FUNCT__ "MatView_MPIAIJ_Spooles_DIST" 91 int MatView_MPIAIJ_SuperLU_DIST(Mat A,PetscViewer viewer) 92 { 93 int ierr; 94 PetscTruth isascii; 95 PetscViewerFormat format; 96 Mat_MPIAIJ_SuperLU_DIST *lu=(Mat_MPIAIJ_SuperLU_DIST*)(A->spptr); 97 98 PetscFunctionBegin; 99 ierr = (*lu->MatView)(A,viewer);CHKERRQ(ierr); 100 101 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 102 if (isascii) { 103 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 104 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 105 ierr = MatMPIAIJFactorInfo_SuperLu(A,viewer);CHKERRQ(ierr); 106 } 107 } 108 PetscFunctionReturn(0); 109 } 110 111 #undef __FUNCT__ 112 #define __FUNCT__ "MatAssemblyEnd_MPIAIJ_SuperLU_DIST" 113 int MatAssemblyEnd_MPIAIJ_SuperLU_DIST(Mat A,MatAssemblyType mode) { 114 int ierr; 115 Mat_MPIAIJ_SuperLU_DIST *lu=(Mat_MPIAIJ_SuperLU_DIST*)(A->spptr); 116 117 PetscFunctionBegin; 118 ierr = (*lu->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 119 ierr = MatUseSuperLU_DIST_MPIAIJ(A);CHKERRQ(ierr); 120 PetscFunctionReturn(0); 121 } 122 123 #undef __FUNCT__ 124 #define __FUNCT__ "MatSolve_MPIAIJ_SuperLU_DIST" 125 int MatSolve_MPIAIJ_SuperLU_DIST(Mat A,Vec b_mpi,Vec x) 126 { 127 Mat_MPIAIJ *aa = (Mat_MPIAIJ*)A->data; 128 Mat_MPIAIJ_SuperLU_DIST *lu = (Mat_MPIAIJ_SuperLU_DIST*)A->spptr; 129 int ierr, size=aa->size; 130 int m=A->M, N=A->N; 131 SuperLUStat_t stat; 132 double berr[1]; 133 PetscScalar *bptr; 134 int info, nrhs=1; 135 Vec x_seq; 136 IS iden; 137 VecScatter scat; 138 PetscLogDouble time0,time,time_min,time_max; 139 140 PetscFunctionBegin; 141 if (size > 1) { 142 if (lu->MatInputMode == GLOBAL) { /* global mat input, convert b to x_seq */ 143 ierr = VecCreateSeq(PETSC_COMM_SELF,N,&x_seq);CHKERRQ(ierr); 144 ierr = ISCreateStride(PETSC_COMM_SELF,N,0,1,&iden);CHKERRQ(ierr); 145 ierr = VecScatterCreate(b_mpi,iden,x_seq,iden,&scat);CHKERRQ(ierr); 146 ierr = ISDestroy(iden);CHKERRQ(ierr); 147 148 ierr = VecScatterBegin(b_mpi,x_seq,INSERT_VALUES,SCATTER_FORWARD,scat);CHKERRQ(ierr); 149 ierr = VecScatterEnd(b_mpi,x_seq,INSERT_VALUES,SCATTER_FORWARD,scat);CHKERRQ(ierr); 150 ierr = VecGetArray(x_seq,&bptr);CHKERRQ(ierr); 151 } else { /* distributed mat input */ 152 ierr = VecCopy(b_mpi,x);CHKERRQ(ierr); 153 ierr = VecGetArray(x,&bptr);CHKERRQ(ierr); 154 } 155 } else { /* size == 1 */ 156 ierr = VecCopy(b_mpi,x);CHKERRQ(ierr); 157 ierr = VecGetArray(x,&bptr);CHKERRQ(ierr); 158 } 159 160 lu->options.Fact = FACTORED; /* The factored form of A is supplied. Local option used by this func. only.*/ 161 162 PStatInit(&stat); /* Initialize the statistics variables. */ 163 if (lu->StatPrint) { 164 ierr = MPI_Barrier(A->comm);CHKERRQ(ierr); /* to be removed */ 165 ierr = PetscGetTime(&time0);CHKERRQ(ierr); /* to be removed */ 166 } 167 if (lu->MatInputMode == GLOBAL) { 168 #if defined(PETSC_USE_COMPLEX) 169 pzgssvx_ABglobal(&lu->options, &lu->A_sup, &lu->ScalePermstruct,(doublecomplex*)bptr, m, nrhs, 170 &lu->grid, &lu->LUstruct, berr, &stat, &info); 171 #else 172 pdgssvx_ABglobal(&lu->options, &lu->A_sup, &lu->ScalePermstruct,bptr, m, nrhs, 173 &lu->grid, &lu->LUstruct, berr, &stat, &info); 174 #endif 175 } else { /* distributed mat input */ 176 #if defined(PETSC_USE_COMPLEX) 177 pzgssvx(&lu->options, &lu->A_sup, &lu->ScalePermstruct, (doublecomplex*)bptr, A->M, nrhs, &lu->grid, 178 &lu->LUstruct, &lu->SOLVEstruct, berr, &stat, &info); 179 if (info) SETERRQ1(1,"pzgssvx fails, info: %d\n",info); 180 #else 181 pdgssvx(&lu->options, &lu->A_sup, &lu->ScalePermstruct, bptr, A->M, nrhs, &lu->grid, 182 &lu->LUstruct, &lu->SOLVEstruct, berr, &stat, &info); 183 if (info) SETERRQ1(1,"pdgssvx fails, info: %d\n",info); 184 #endif 185 } 186 if (lu->StatPrint) { 187 ierr = PetscGetTime(&time);CHKERRQ(ierr); /* to be removed */ 188 PStatPrint(&lu->options, &stat, &lu->grid); /* Print the statistics. */ 189 } 190 PStatFree(&stat); 191 192 if (size > 1) { 193 if (lu->MatInputMode == GLOBAL){ /* convert seq x to mpi x */ 194 ierr = VecRestoreArray(x_seq,&bptr);CHKERRQ(ierr); 195 ierr = VecScatterBegin(x_seq,x,INSERT_VALUES,SCATTER_REVERSE,scat);CHKERRQ(ierr); 196 ierr = VecScatterEnd(x_seq,x,INSERT_VALUES,SCATTER_REVERSE,scat);CHKERRQ(ierr); 197 ierr = VecScatterDestroy(scat);CHKERRQ(ierr); 198 ierr = VecDestroy(x_seq);CHKERRQ(ierr); 199 } else { 200 ierr = VecRestoreArray(x,&bptr);CHKERRQ(ierr); 201 } 202 } else { 203 ierr = VecRestoreArray(x,&bptr);CHKERRQ(ierr); 204 } 205 if (lu->StatPrint) { 206 time0 = time - time0; 207 ierr = MPI_Reduce(&time0,&time_max,1,MPI_DOUBLE,MPI_MAX,0,A->comm);CHKERRQ(ierr); 208 ierr = MPI_Reduce(&time0,&time_min,1,MPI_DOUBLE,MPI_MIN,0,A->comm);CHKERRQ(ierr); 209 ierr = MPI_Reduce(&time0,&time,1,MPI_DOUBLE,MPI_SUM,0,A->comm);CHKERRQ(ierr); 210 time = time/size; /* average time */ 211 ierr = PetscPrintf(A->comm, " Time for superlu_dist solve (max/min/avg): %g / %g / %g\n\n",time_max,time_min,time);CHKERRQ(ierr); 212 } 213 PetscFunctionReturn(0); 214 } 215 216 #undef __FUNCT__ 217 #define __FUNCT__ "MatLUFactorNumeric_MPIAIJ_SuperLU_DIST" 218 int MatLUFactorNumeric_MPIAIJ_SuperLU_DIST(Mat A,Mat *F) 219 { 220 Mat_MPIAIJ *fac = (Mat_MPIAIJ*)(*F)->data,*mat; 221 Mat *tseq,A_seq = PETSC_NULL; 222 Mat_SeqAIJ *aa,*bb; 223 Mat_MPIAIJ_SuperLU_DIST *lu = (Mat_MPIAIJ_SuperLU_DIST*)(*F)->spptr; 224 int M=A->M,N=A->N,info,ierr,size=fac->size,i,*ai,*aj,*bi,*bj,nz,rstart,*garray, 225 m=A->m, irow,colA_start,j,jcol,jB,countA,countB,*bjj,*ajj; 226 SuperLUStat_t stat; 227 double *berr=0; 228 IS isrow; 229 PetscLogDouble time0[2],time[2],time_min[2],time_max[2]; 230 #if defined(PETSC_USE_COMPLEX) 231 doublecomplex *av, *bv; 232 #else 233 double *av, *bv; 234 #endif 235 236 PetscFunctionBegin; 237 if (lu->StatPrint) { 238 ierr = MPI_Barrier(A->comm);CHKERRQ(ierr); 239 ierr = PetscGetTime(&time0[0]);CHKERRQ(ierr); 240 } 241 242 if (lu->MatInputMode == GLOBAL) { /* global mat input */ 243 if (size > 1) { /* convert mpi A to seq mat A */ 244 ierr = ISCreateStride(PETSC_COMM_SELF,M,0,1,&isrow); CHKERRQ(ierr); 245 ierr = MatGetSubMatrices(A,1,&isrow,&isrow,MAT_INITIAL_MATRIX,&tseq); CHKERRQ(ierr); 246 ierr = ISDestroy(isrow);CHKERRQ(ierr); 247 248 A_seq = *tseq; 249 ierr = PetscFree(tseq);CHKERRQ(ierr); 250 aa = (Mat_SeqAIJ*)A_seq->data; 251 } else { 252 aa = (Mat_SeqAIJ*)A->data; 253 } 254 255 /* Allocate storage, then convert Petsc NR matrix to SuperLU_DIST NC */ 256 if (lu->flg == DIFFERENT_NONZERO_PATTERN) {/* first numeric factorization */ 257 #if defined(PETSC_USE_COMPLEX) 258 zallocateA_dist(N, aa->nz, &lu->val, &lu->col, &lu->row); 259 #else 260 dallocateA_dist(N, aa->nz, &lu->val, &lu->col, &lu->row); 261 #endif 262 } else { /* successive numeric factorization, sparsity pattern is reused. */ 263 Destroy_CompCol_Matrix_dist(&lu->A_sup); 264 Destroy_LU(N, &lu->grid, &lu->LUstruct); 265 lu->options.Fact = SamePattern; 266 } 267 #if defined(PETSC_USE_COMPLEX) 268 zCompRow_to_CompCol_dist(M,N,aa->nz,(doublecomplex*)aa->a,aa->j,aa->i,&lu->val,&lu->col, &lu->row); 269 #else 270 dCompRow_to_CompCol_dist(M,N,aa->nz,aa->a,aa->j,aa->i,&lu->val, &lu->col, &lu->row); 271 #endif 272 273 /* Create compressed column matrix A_sup. */ 274 #if defined(PETSC_USE_COMPLEX) 275 zCreate_CompCol_Matrix_dist(&lu->A_sup, M, N, aa->nz, lu->val, lu->col, lu->row, SLU_NC, SLU_Z, SLU_GE); 276 #else 277 dCreate_CompCol_Matrix_dist(&lu->A_sup, M, N, aa->nz, lu->val, lu->col, lu->row, SLU_NC, SLU_D, SLU_GE); 278 #endif 279 } else { /* distributed mat input */ 280 mat = (Mat_MPIAIJ*)A->data; 281 aa=(Mat_SeqAIJ*)(mat->A)->data; 282 bb=(Mat_SeqAIJ*)(mat->B)->data; 283 ai=aa->i; aj=aa->j; 284 bi=bb->i; bj=bb->j; 285 #if defined(PETSC_USE_COMPLEX) 286 av=(doublecomplex*)aa->a; 287 bv=(doublecomplex*)bb->a; 288 #else 289 av=aa->a; 290 bv=bb->a; 291 #endif 292 rstart = mat->rstart; 293 nz = aa->nz + bb->nz; 294 garray = mat->garray; 295 rstart = mat->rstart; 296 297 if (lu->flg == DIFFERENT_NONZERO_PATTERN) {/* first numeric factorization */ 298 #if defined(PETSC_USE_COMPLEX) 299 zallocateA_dist(m, nz, &lu->val, &lu->col, &lu->row); 300 #else 301 dallocateA_dist(m, nz, &lu->val, &lu->col, &lu->row); 302 #endif 303 } else { /* successive numeric factorization, sparsity pattern and perm_c are reused. */ 304 /* Destroy_CompRowLoc_Matrix_dist(&lu->A_sup); */ /* crash! */ 305 Destroy_LU(N, &lu->grid, &lu->LUstruct); 306 lu->options.Fact = SamePattern; 307 } 308 nz = 0; jB = 0; irow = mat->rstart; 309 for ( i=0; i<m; i++ ) { 310 lu->row[i] = nz; 311 countA = ai[i+1] - ai[i]; 312 countB = bi[i+1] - bi[i]; 313 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 314 bjj = bj + bi[i]; 315 316 /* B part, smaller col index */ 317 colA_start = mat->rstart + ajj[0]; /* the smallest global col index of A */ 318 for (j=0; j<countB; j++){ 319 jcol = garray[bjj[j]]; 320 if (jcol > colA_start) { 321 jB = j; 322 break; 323 } 324 lu->col[nz] = jcol; 325 lu->val[nz++] = *bv++; 326 if (j==countB-1) jB = countB; 327 } 328 329 /* A part */ 330 for (j=0; j<countA; j++){ 331 lu->col[nz] = mat->rstart + ajj[j]; 332 lu->val[nz++] = *av++; 333 } 334 335 /* B part, larger col index */ 336 for (j=jB; j<countB; j++){ 337 lu->col[nz] = garray[bjj[j]]; 338 lu->val[nz++] = *bv++; 339 } 340 } 341 lu->row[m] = nz; 342 #if defined(PETSC_USE_COMPLEX) 343 zCreate_CompRowLoc_Matrix_dist(&lu->A_sup, M, N, nz, m, rstart, 344 lu->val, lu->col, lu->row, SLU_NR_loc, SLU_Z, SLU_GE); 345 #else 346 dCreate_CompRowLoc_Matrix_dist(&lu->A_sup, M, N, nz, m, rstart, 347 lu->val, lu->col, lu->row, SLU_NR_loc, SLU_D, SLU_GE); 348 #endif 349 } 350 if (lu->StatPrint) { 351 ierr = PetscGetTime(&time[0]);CHKERRQ(ierr); 352 time0[0] = time[0] - time0[0]; 353 } 354 355 /* Factor the matrix. */ 356 PStatInit(&stat); /* Initialize the statistics variables. */ 357 358 if (lu->StatPrint) { 359 ierr = MPI_Barrier(A->comm);CHKERRQ(ierr); 360 ierr = PetscGetTime(&time0[1]);CHKERRQ(ierr); 361 } 362 363 if (lu->MatInputMode == GLOBAL) { /* global mat input */ 364 #if defined(PETSC_USE_COMPLEX) 365 pzgssvx_ABglobal(&lu->options, &lu->A_sup, &lu->ScalePermstruct, 0, M, 0, 366 &lu->grid, &lu->LUstruct, berr, &stat, &info); 367 #else 368 pdgssvx_ABglobal(&lu->options, &lu->A_sup, &lu->ScalePermstruct, 0, M, 0, 369 &lu->grid, &lu->LUstruct, berr, &stat, &info); 370 #endif 371 } else { /* distributed mat input */ 372 #if defined(PETSC_USE_COMPLEX) 373 pzgssvx(&lu->options, &lu->A_sup, &lu->ScalePermstruct, 0, M, 0, &lu->grid, 374 &lu->LUstruct, &lu->SOLVEstruct, berr, &stat, &info); 375 if (info) SETERRQ1(1,"pzgssvx fails, info: %d\n",info); 376 #else 377 pdgssvx(&lu->options, &lu->A_sup, &lu->ScalePermstruct, 0, M, 0, &lu->grid, 378 &lu->LUstruct, &lu->SOLVEstruct, berr, &stat, &info); 379 if (info) SETERRQ1(1,"pdgssvx fails, info: %d\n",info); 380 #endif 381 } 382 if (lu->StatPrint) { 383 ierr = PetscGetTime(&time[1]);CHKERRQ(ierr); /* to be removed */ 384 time0[1] = time[1] - time0[1]; 385 if (lu->StatPrint) PStatPrint(&lu->options, &stat, &lu->grid); /* Print the statistics. */ 386 } 387 PStatFree(&stat); 388 389 if (lu->MatInputMode == GLOBAL && size > 1){ 390 ierr = MatDestroy(A_seq);CHKERRQ(ierr); 391 } 392 393 if (lu->StatPrint) { 394 ierr = MPI_Reduce(time0,time_max,2,MPI_DOUBLE,MPI_MAX,0,A->comm); 395 ierr = MPI_Reduce(time0,time_min,2,MPI_DOUBLE,MPI_MIN,0,A->comm); 396 ierr = MPI_Reduce(time0,time,2,MPI_DOUBLE,MPI_SUM,0,A->comm); 397 for (i=0; i<2; i++) time[i] = time[i]/size; /* average time */ 398 ierr = PetscPrintf(A->comm, " Time for mat conversion (max/min/avg): %g / %g / %g\n",time_max[0],time_min[0],time[0]); 399 ierr = PetscPrintf(A->comm, " Time for superlu_dist fact (max/min/avg): %g / %g / %g\n\n",time_max[1],time_min[1],time[1]); 400 } 401 (*F)->assembled = PETSC_TRUE; 402 lu->flg = SAME_NONZERO_PATTERN; 403 PetscFunctionReturn(0); 404 } 405 406 /* Note the Petsc r and c permutations are ignored */ 407 #undef __FUNCT__ 408 #define __FUNCT__ "MatLUFactorSymbolic_MPIAIJ_SuperLU_DIST" 409 int MatLUFactorSymbolic_MPIAIJ_SuperLU_DIST(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 410 { 411 Mat B; 412 Mat_MPIAIJ_SuperLU_DIST *lu; 413 int ierr,M=A->M,N=A->N,size; 414 superlu_options_t options; 415 char buff[32]; 416 PetscTruth flg; 417 char *ptype[] = {"MMD_AT_PLUS_A","NATURAL","MMD_ATA","COLAMD"}; 418 char *prtype[] = {"LargeDiag","NATURAL"}; 419 PetscFunctionBegin; 420 421 /* Create the factorization matrix */ 422 ierr = MatCreate(A->comm,A->m,A->n,M,N,&B);CHKERRQ(ierr); 423 ierr = MatSetType(B,MATSUPERLUDIST);CHKERRQ(ierr); 424 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL); 425 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 426 427 B->ops->lufactornumeric = MatLUFactorNumeric_MPIAIJ_SuperLU_DIST; 428 B->ops->solve = MatSolve_MPIAIJ_SuperLU_DIST; 429 B->ops->destroy = MatDestroy_MPIAIJ_SuperLU_DIST; 430 B->factor = FACTOR_LU; 431 432 lu = (Mat_MPIAIJ_SuperLU_DIST*)(B->spptr); 433 434 /* Set the input options */ 435 set_default_options(&options); 436 lu->MatInputMode = GLOBAL; 437 ierr = MPI_Comm_dup(A->comm,&(lu->comm_superlu));CHKERRQ(ierr); 438 439 ierr = MPI_Comm_size(A->comm,&size);CHKERRQ(ierr); 440 lu->nprow = size/2; /* Default process rows. */ 441 if (lu->nprow == 0) lu->nprow = 1; 442 lu->npcol = size/lu->nprow; /* Default process columns. */ 443 444 ierr = PetscOptionsBegin(A->comm,A->prefix,"SuperLU_Dist Options","Mat");CHKERRQ(ierr); 445 446 ierr = PetscOptionsInt("-mat_aij_superlu_dist_r","Number rows in processor partition","None",lu->nprow,&lu->nprow,PETSC_NULL);CHKERRQ(ierr); 447 ierr = PetscOptionsInt("-mat_aij_superlu_dist_c","Number columns in processor partition","None",lu->npcol,&lu->npcol,PETSC_NULL);CHKERRQ(ierr); 448 if (size != lu->nprow * lu->npcol) SETERRQ(1,"Number of processes should be equal to nprow*npcol"); 449 450 ierr = PetscOptionsInt("-mat_aij_superlu_dist_matinput","Matrix input mode (0: GLOBAL; 1: DISTRIBUTED)","None",lu->MatInputMode,&lu->MatInputMode,PETSC_NULL);CHKERRQ(ierr); 451 if(lu->MatInputMode == DISTRIBUTED && size == 1) lu->MatInputMode = GLOBAL; 452 453 ierr = PetscOptionsLogical("-mat_aij_superlu_dist_equil","Equilibrate matrix","None",PETSC_TRUE,&flg,0);CHKERRQ(ierr); 454 if (!flg) { 455 options.Equil = NO; 456 } 457 458 ierr = PetscOptionsEList("-mat_aij_superlu_dist_rowperm","Row permutation","None",prtype,2,prtype[0],buff,32,&flg);CHKERRQ(ierr); 459 while (flg) { 460 ierr = PetscStrcmp(buff,"LargeDiag",&flg);CHKERRQ(ierr); 461 if (flg) { 462 options.RowPerm = LargeDiag; 463 break; 464 } 465 ierr = PetscStrcmp(buff,"NATURAL",&flg);CHKERRQ(ierr); 466 if (flg) { 467 options.RowPerm = NOROWPERM; 468 break; 469 } 470 SETERRQ1(1,"Unknown row permutation %s",buff); 471 } 472 473 ierr = PetscOptionsEList("-mat_aij_superlu_dist_colperm","Column permutation","None",ptype,4,ptype[0],buff,32,&flg);CHKERRQ(ierr); 474 while (flg) { 475 ierr = PetscStrcmp(buff,"MMD_AT_PLUS_A",&flg);CHKERRQ(ierr); 476 if (flg) { 477 options.ColPerm = MMD_AT_PLUS_A; 478 break; 479 } 480 ierr = PetscStrcmp(buff,"NATURAL",&flg);CHKERRQ(ierr); 481 if (flg) { 482 options.ColPerm = NATURAL; 483 break; 484 } 485 ierr = PetscStrcmp(buff,"MMD_ATA",&flg);CHKERRQ(ierr); 486 if (flg) { 487 options.ColPerm = MMD_ATA; 488 break; 489 } 490 ierr = PetscStrcmp(buff,"COLAMD",&flg);CHKERRQ(ierr); 491 if (flg) { 492 options.ColPerm = COLAMD; 493 break; 494 } 495 SETERRQ1(1,"Unknown column permutation %s",buff); 496 } 497 498 ierr = PetscOptionsLogical("-mat_aij_superlu_dist_replacetinypivot","Replace tiny pivots","None",PETSC_TRUE,&flg,0);CHKERRQ(ierr); 499 if (!flg) { 500 options.ReplaceTinyPivot = NO; 501 } 502 503 options.IterRefine = NOREFINE; 504 ierr = PetscOptionsLogical("-mat_aij_superlu_dist_iterrefine","Use iterative refinement","None",PETSC_FALSE,&flg,0);CHKERRQ(ierr); 505 if (flg) { 506 options.IterRefine = DOUBLE; 507 } 508 509 if (PetscLogPrintInfo) { 510 lu->StatPrint = (int)PETSC_TRUE; 511 } else { 512 lu->StatPrint = (int)PETSC_FALSE; 513 } 514 ierr = PetscOptionsLogical("-mat_aij_superlu_dist_statprint","Print factorization information","None", 515 (PetscTruth)lu->StatPrint,(PetscTruth*)&lu->StatPrint,0);CHKERRQ(ierr); 516 PetscOptionsEnd(); 517 518 /* Initialize the SuperLU process grid. */ 519 superlu_gridinit(lu->comm_superlu, lu->nprow, lu->npcol, &lu->grid); 520 521 /* Initialize ScalePermstruct and LUstruct. */ 522 ScalePermstructInit(M, N, &lu->ScalePermstruct); 523 LUstructInit(M, N, &lu->LUstruct); 524 525 lu->options = options; 526 lu->flg = DIFFERENT_NONZERO_PATTERN; 527 lu->CleanUpSuperLUDist = PETSC_TRUE; 528 *F = B; 529 PetscFunctionReturn(0); 530 } 531 532 #undef __FUNCT__ 533 #define __FUNCT__ "MatUseSuperLU_DIST_MPIAIJ" 534 int MatUseSuperLU_DIST_MPIAIJ(Mat A) 535 { 536 PetscFunctionBegin; 537 A->ops->lufactorsymbolic = MatLUFactorSymbolic_MPIAIJ_SuperLU_DIST; 538 A->ops->lufactornumeric = MatLUFactorNumeric_MPIAIJ_SuperLU_DIST; 539 PetscFunctionReturn(0); 540 } 541 542 #undef __FUNCT__ 543 #define __FUNCT__ "MatMPIAIJFactorInfo_SuperLu" 544 int MatMPIAIJFactorInfo_SuperLu(Mat A,PetscViewer viewer) 545 { 546 Mat_MPIAIJ_SuperLU_DIST *lu= (Mat_MPIAIJ_SuperLU_DIST*)A->spptr; 547 superlu_options_t options; 548 int ierr; 549 char *colperm; 550 551 PetscFunctionBegin; 552 /* check if matrix is superlu_dist type */ 553 if (A->ops->solve != MatSolve_MPIAIJ_SuperLU_DIST) PetscFunctionReturn(0); 554 555 options = lu->options; 556 ierr = PetscViewerASCIIPrintf(viewer,"SuperLU_DIST run parameters:\n");CHKERRQ(ierr); 557 ierr = PetscViewerASCIIPrintf(viewer," Equilibrate matrix %s \n",(options.Equil != NO) ? "true": "false");CHKERRQ(ierr); 558 ierr = PetscViewerASCIIPrintf(viewer," Replace tiny pivots %s \n",(options.ReplaceTinyPivot != NO) ? "true": "false");CHKERRQ(ierr); 559 ierr = PetscViewerASCIIPrintf(viewer," Use iterative refinement %s \n",(options.IterRefine == DOUBLE) ? "true": "false");CHKERRQ(ierr); 560 ierr = PetscViewerASCIIPrintf(viewer," Processors in row %d col partition %d \n",lu->nprow,lu->npcol);CHKERRQ(ierr); 561 ierr = PetscViewerASCIIPrintf(viewer," Row permutation %s \n",(options.RowPerm == NOROWPERM) ? "NATURAL": "LargeDiag");CHKERRQ(ierr); 562 if (options.ColPerm == NATURAL) { 563 colperm = "NATURAL"; 564 } else if (options.ColPerm == MMD_AT_PLUS_A) { 565 colperm = "MMD_AT_PLUS_A"; 566 } else if (options.ColPerm == MMD_ATA) { 567 colperm = "MMD_ATA"; 568 } else if (options.ColPerm == COLAMD) { 569 colperm = "COLAMD"; 570 } else { 571 SETERRQ(1,"Unknown column permutation"); 572 } 573 ierr = PetscViewerASCIIPrintf(viewer," Column permutation %s \n",colperm);CHKERRQ(ierr); 574 PetscFunctionReturn(0); 575 } 576 577 EXTERN_C_BEGIN 578 #undef __FUNCT__ 579 #define __FUNCT__ "MatCreate_MPIAIJ_SuperLU_DIST" 580 int MatCreate_MPIAIJ_SuperLU_DIST(Mat A) { 581 int ierr,size; 582 MPI_Comm comm; 583 Mat_MPIAIJ_SuperLU_DIST *lu; 584 585 PetscFunctionBegin; 586 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 587 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 588 if (size == 1) { 589 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 590 } else { 591 ierr = MatSetType(A,MATMPIAIJ);CHKERRQ(ierr); 592 } 593 ierr = MatUseSuperLU_SeqAIJ(A);CHKERRQ(ierr); 594 595 ierr = PetscNew(Mat_MPIAIJ_SuperLU_DIST,&lu);CHKERRQ(ierr); 596 lu->MatView = A->ops->view; 597 lu->MatAssemblyEnd = A->ops->assemblyend; 598 lu->MatDestroy = A->ops->destroy; 599 lu->CleanUpSuperLUDist = PETSC_FALSE; 600 A->spptr = (void*)lu; 601 A->ops->view = MatView_MPIAIJ_SuperLU_DIST; 602 A->ops->assemblyend = MatAssemblyEnd_MPIAIJ_SuperLU_DIST; 603 A->ops->destroy = MatDestroy_MPIAIJ_SuperLU_DIST; 604 PetscFunctionReturn(0); 605 } 606 EXTERN_C_END 607 608 EXTERN_C_BEGIN 609 #undef __FUNCT__ 610 #define __FUNCT__ "MatLoad_MPIAIJ_SuperLU_DIST" 611 int MatLoad_MPIAIJ_SuperLU_DIST(PetscViewer viewer,MatType type,Mat *A) { 612 int ierr,size,(*r)(PetscViewer,MatType,Mat*); 613 MPI_Comm comm; 614 615 PetscFunctionBegin; 616 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 617 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 618 if (size == 1) { 619 ierr = PetscFListFind(comm,MatLoadList,MATSEQAIJ,(void(**)(void))&r);CHKERRQ(ierr); 620 } else { 621 ierr = PetscFListFind(comm,MatLoadList,MATMPIAIJ,(void(**)(void))&r);CHKERRQ(ierr); 622 } 623 ierr = (*r)(viewer,type,A);CHKERRQ(ierr); 624 PetscFunctionReturn(0); 625 } 626 EXTERN_C_END 627