1 #define PETSCMAT_DLL 2 3 /* 4 Provides an interface to the PaStiX sparse solver 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 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat,PetscReal,PetscTruth*); 11 12 EXTERN_C_BEGIN 13 #include "mpi.h" 14 #include "pastix.h" 15 static int iter_print = 0; 16 static int csr_iter_print = 0; 17 EXTERN_C_END 18 19 typedef struct Mat_Pastix_ { 20 pastix_data_t *pastix_data; /* Pastix data storage structure */ 21 MatStructure matstruc; 22 PetscInt n; /* Number of columns in the matrix */ 23 PetscInt *colptr; /* Index of first element of each column in row and val */ 24 PetscInt *row; /* Row of each element of the matrix */ 25 PetscScalar *val; /* Value of each element of the matrix */ 26 PetscInt *perm; /* Permutation tabular */ 27 PetscInt *invp; /* Reverse permutation tabular */ 28 PetscScalar *rhs; /* Rhight-hand-side member */ 29 PetscInt rhsnbr; /* Rhight-hand-side number (must be 1) */ 30 PetscInt iparm[64]; /* Integer parameters */ 31 double dparm[64]; /* Floating point parameters */ 32 MPI_Comm pastix_comm; /* PaStiX MPI communicator */ 33 PetscMPIInt commRank; /* MPI rank */ 34 PetscMPIInt commSize; /* MPI communicator size */ 35 PetscTruth CleanUpPastix; /* Boolean indicating if we call PaStiX clean step */ 36 VecScatter scat_rhs; 37 VecScatter scat_sol; 38 Vec b_seq,x_seq; 39 PetscTruth isAIJ; 40 PetscInt nSolve; /* Number of consecutive solve */ 41 PetscErrorCode (*MatDestroy)(Mat); 42 } Mat_Pastix; 43 44 EXTERN PetscErrorCode MatDuplicate_Pastix(Mat,MatDuplicateOption,Mat*); 45 46 /* 47 convert Petsc seqaij matrix to CSC: colptr[n], row[nz], val[nz] 48 49 input: 50 A - matrix in seqaij or mpisbaij (bs=1) format 51 valOnly - FALSE: spaces are allocated and values are set for the CSC 52 TRUE: Only fill values 53 output: 54 n - Size of the matrix 55 colptr - Index of first element of each column in row and val 56 row - Row of each element of the matrix 57 values - Value of each element of the matrix 58 */ 59 PetscErrorCode MatConvertToCSC(Mat A, 60 PetscTruth valOnly, 61 PetscInt *n, 62 PetscInt **colptr, 63 PetscInt **row, 64 PetscScalar **values) { 65 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data; 66 PetscInt *rowptr = aa->i; 67 PetscInt *col = aa->j; 68 PetscScalar *rvalues = aa->a; 69 PetscInt m = A->rmap->N; 70 PetscInt nnz; 71 PetscInt i,j, k; 72 PetscInt base = 1; 73 PetscInt idx; 74 PetscErrorCode ierr; 75 PetscInt colidx; 76 PetscInt *colcount; 77 PetscTruth isSym; 78 79 80 PetscFunctionBegin; 81 /* Allocate the CSC */ 82 83 84 MatIsSymmetric_SeqAIJ(A,0.0,&isSym); 85 *n = A->cmap->N; 86 87 /* PaStiX only needs triangular matrix if matrix is symmetric 88 */ 89 if (isSym) 90 { 91 nnz = (aa->nz - *n)/2 + *n; 92 } 93 else 94 { 95 nnz = aa->nz; 96 } 97 98 ierr = PetscMalloc((*n)*sizeof(PetscInt) ,&colcount);CHKERRQ(ierr); 99 if (!valOnly){ 100 ierr = PetscMalloc(((*n)+1) *sizeof(PetscInt) ,colptr );CHKERRQ(ierr); 101 ierr = PetscMalloc( nnz *sizeof(PetscInt) ,row);CHKERRQ(ierr); 102 ierr = PetscMalloc( nnz *sizeof(PetscScalar),values);CHKERRQ(ierr); 103 104 for (i = 0; i < m; i++) 105 colcount[i] = 0; 106 /* Fill-in colptr */ 107 for (i = 0; i < m; i++) 108 for (j = rowptr[i]; j < rowptr[i+1]; j++) 109 if (!isSym || col[j] <= i) 110 colcount[col[j]]++; 111 112 (*colptr)[0] = base; 113 for (j = 0; j < *n; j++) { 114 (*colptr)[j+1] = (*colptr)[j] + colcount[j]; 115 /* in next loop we fill starting from (*colptr)[colidx] - base */ 116 colcount[j] = -base; 117 } 118 119 /* Fill-in rows and values */ 120 for (i = 0; i < m; i++) { 121 for (j = rowptr[i]; j < rowptr[i+1]; j++) { 122 if (!isSym || col[j] <= i) 123 { 124 colidx = col[j]; 125 idx = (*colptr)[colidx] + colcount[colidx]; 126 (*row)[idx] = i + base; 127 (*values)[idx] = rvalues[j]; 128 colcount[colidx]++; 129 } 130 } 131 } 132 } 133 else { 134 /* Fill-in only values */ 135 for (i = 0; i < m; i++) { 136 for (j = rowptr[i]; j < rowptr[i+1]; j++) { 137 colidx = col[j]; 138 if (!isSym || col[j] <= i) 139 { 140 /* look for the value to fill */ 141 for (k = (*colptr)[colidx] - base; 142 k < (*colptr)[colidx + 1] - base; 143 k++) { 144 if ((*row)[k] == i) { 145 (*values)[k] = rvalues[j]; 146 break; 147 } 148 } 149 /* shouldn't happen, overflow */ 150 if (k == (*colptr)[colidx + 1] - base) 151 PetscFunctionReturn(1); 152 } 153 } 154 } 155 } 156 ierr = PetscFree(colcount);CHKERRQ(ierr); 157 158 PetscFunctionReturn(0); 159 } 160 161 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "MatDestroy_Pastix" 165 /* 166 Call clean step of PaStiX if lu->CleanUpPastix == true. 167 Free the CSC matrix. 168 */ 169 PetscErrorCode MatDestroy_Pastix(Mat A) 170 { 171 Mat_Pastix *lu=(Mat_Pastix*)A->spptr; 172 PetscErrorCode ierr; 173 PetscMPIInt size=lu->commSize; 174 175 PetscFunctionBegin; 176 if (lu->CleanUpPastix) { 177 /* Terminate instance, deallocate memories */ 178 if (size > 1){ 179 ierr = VecScatterDestroy(lu->scat_rhs);CHKERRQ(ierr); 180 ierr = VecDestroy(lu->b_seq);CHKERRQ(ierr); 181 if (lu->nSolve && lu->scat_sol){ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr);} 182 if (lu->nSolve && lu->x_seq){ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr);} 183 } 184 185 lu->iparm[IPARM_START_TASK]=API_TASK_CLEAN; 186 lu->iparm[IPARM_END_TASK]=API_TASK_CLEAN; 187 188 pastix((pastix_data_t **)&(lu->pastix_data), 189 lu->pastix_comm, 190 (pastix_int_t) lu->n, 191 (pastix_int_t*) lu->colptr, 192 (pastix_int_t*) lu->row, 193 (pastix_float_t*) lu->val, 194 (pastix_int_t*) lu->perm, 195 (pastix_int_t*) lu->invp, 196 (pastix_float_t*) lu->rhs, 197 (pastix_int_t) lu->rhsnbr, 198 (pastix_int_t*) lu->iparm, 199 lu->dparm); 200 201 ierr = PetscFree(lu->colptr);CHKERRQ(ierr); 202 ierr = PetscFree(lu->row); CHKERRQ(ierr); 203 ierr = PetscFree(lu->val); CHKERRQ(ierr); 204 ierr = PetscFree(lu->perm); CHKERRQ(ierr); 205 ierr = PetscFree(lu->invp); CHKERRQ(ierr); 206 /* ierr = PetscFree(lu->rhs); CHKERRQ(ierr); */ 207 ierr = MPI_Comm_free(&(lu->pastix_comm));CHKERRQ(ierr); 208 209 } 210 ierr = (lu->MatDestroy)(A);CHKERRQ(ierr); 211 PetscFunctionReturn(0); 212 } 213 214 #undef __FUNCT__ 215 #define __FUNCT__ "MatSolve_PaStiX" 216 /* 217 Gather right-hand-side. 218 Call for Solve step. 219 Scatter solution. 220 */ 221 PetscErrorCode MatSolve_PaStiX(Mat A,Vec b,Vec x) 222 { 223 Mat_Pastix *lu=(Mat_Pastix*)A->spptr; 224 PetscScalar *array; 225 Vec x_seq; 226 IS is_iden,is_petsc; 227 PetscErrorCode ierr; 228 PetscInt i,j; 229 230 PetscFunctionBegin; 231 lu->rhsnbr = 1; 232 x_seq = lu->b_seq; 233 if (lu->commSize > 1){ 234 /* PaStiX only supports centralized rhs. Scatter b into a seqential rhs vector */ 235 ierr = VecScatterBegin(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 236 ierr = VecScatterEnd(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 237 ierr = VecGetArray(x_seq,&array); 238 CHKERRQ(ierr); 239 } 240 else { /* size == 1 */ 241 ierr = VecCopy(b,x);CHKERRQ(ierr); 242 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 243 } 244 lu->rhs = array; 245 if (lu->commSize == 1){ 246 ierr = VecRestoreArray(x,&array);CHKERRQ(ierr); 247 } else { 248 ierr = VecRestoreArray(x_seq,&array);CHKERRQ(ierr); 249 } 250 251 /* solve phase */ 252 /*-------------*/ 253 lu->iparm[IPARM_START_TASK] = API_TASK_SOLVE; 254 lu->iparm[IPARM_END_TASK] = API_TASK_REFINE; 255 lu->iparm[IPARM_RHS_MAKING] = API_RHS_B; 256 257 pastix((pastix_data_t **)&(lu->pastix_data), 258 (MPI_Comm) lu->pastix_comm, 259 (pastix_int_t) lu->n, 260 (pastix_int_t*) lu->colptr, 261 (pastix_int_t*) lu->row, 262 (pastix_float_t*) lu->val, 263 (pastix_int_t*) lu->perm, 264 (pastix_int_t*) lu->invp, 265 (pastix_float_t*) lu->rhs, 266 (pastix_int_t) lu->rhsnbr, 267 (pastix_int_t*) lu->iparm, 268 (double*) lu->dparm); 269 270 if (lu->iparm[IPARM_ERROR_NUMBER] < 0) { 271 SETERRQ1(PETSC_ERR_LIB,"Error reported by PaStiX in solve phase: lu->iparm[IPARM_ERROR_NUMBER] = %d\n", 272 lu->iparm[IPARM_ERROR_NUMBER] ); 273 } 274 275 if (lu->commSize == 1){ 276 ierr = VecRestoreArray(x,&(lu->rhs));CHKERRQ(ierr); 277 } else { 278 ierr = VecRestoreArray(x_seq,&(lu->rhs));CHKERRQ(ierr); 279 } 280 281 if (lu->commSize > 1) { /* convert PaStiX centralized solution to petsc mpi x */ 282 ierr = VecScatterBegin(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 283 ierr = VecScatterEnd(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 284 } 285 lu->nSolve++; 286 PetscFunctionReturn(0); 287 } 288 289 #if !defined(PETSC_USE_COMPLEX) 290 /* 291 TODO: Fill this function 292 I didn't fill this function 293 because I didn't understood its goal. 294 */ 295 296 /* 297 input: 298 F: numeric factor 299 output: 300 nneg: total number of pivots 301 nzero: 0 302 npos: (global dimension of F) - nneg 303 */ 304 305 #undef __FUNCT__ 306 #define __FUNCT__ "MatGetInertia_SBAIJPASTIX" 307 PetscErrorCode MatGetInertia_SBAIJPASTIX(Mat F,int *nneg,int *nzero,int *npos) 308 { 309 Mat_Pastix *lu =(Mat_Pastix*)F->spptr; 310 PetscErrorCode ierr; 311 PetscMPIInt size; 312 313 PetscFunctionBegin; 314 /* ierr = MPI_Comm_size(((PetscObject)F)->comm,&size);CHKERRQ(ierr); */ 315 /* /\* PASTIX 4.3.1 calls ScaLAPACK when ICNTL(13)=0 (default), which does not offer the possibility to compute the inertia of a dense matrix. Set ICNTL(13)=1 to skip ScaLAPACK *\/ */ 316 /* if (size > 1 && lu->id.ICNTL(13) != 1){ */ 317 /* SETERRQ1(PETSC_ERR_ARG_WRONG,"ICNTL(13)=%d. -mat_pastix_icntl_13 must be set as 1 for correct global matrix inertia\n",lu->id.INFOG(13)); */ 318 /* } */ 319 /* if (nneg){ */ 320 /* if (!lu->commSize){ */ 321 /* *nneg = lu->id.INFOG(12); */ 322 /* } */ 323 /* ierr = MPI_Bcast(nneg,1,MPI_INT,0,lu->comm_pastix);CHKERRQ(ierr); */ 324 /* } */ 325 /* if (nzero) *nzero = lu->iparm[IPARM_NNZEROS]; */ 326 /* if (npos) *npos = F->rmap->N - (*nneg); */ 327 PetscFunctionReturn(0); 328 } 329 #endif /* !defined(PETSC_USE_COMPLEX) */ 330 331 /* 332 Numeric factorisation using PaStiX solver. 333 334 */ 335 #undef __FUNCT__ 336 #define __FUNCT__ "MatFactorNumeric_PASTIX" 337 PetscErrorCode MatFactorNumeric_PaStiX(Mat F,Mat A,const MatFactorInfo *info) 338 { 339 Mat_Pastix *lu =(Mat_Pastix*)(F)->spptr; 340 Mat *tseq,A_seq = PETSC_NULL; 341 PetscErrorCode ierr = 0; 342 PetscInt rnz,nnz,nz=0,i,*ai,*aj,icntl; 343 PetscInt M=A->rmap->N,N=A->cmap->N; 344 PetscTruth valOnly,flg, isSym; 345 Mat F_diag; 346 IS is_iden; 347 Vec b; 348 IS isrow; 349 PetscTruth isSeqAIJ,isSeqSBAIJ; 350 Mat_SeqAIJ *aa,*bb; 351 352 PetscFunctionBegin; 353 ierr = PetscTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 354 ierr = PetscTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 355 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 356 (F)->ops->solve = MatSolve_PaStiX; 357 358 /* Initialize a PASTIX instance */ 359 ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(lu->pastix_comm));CHKERRQ(ierr); 360 ierr = MPI_Comm_rank(lu->pastix_comm, &lu->commRank); CHKERRQ(ierr); 361 ierr = MPI_Comm_size(lu->pastix_comm, &lu->commSize); CHKERRQ(ierr); 362 363 /* Set pastix options */ 364 lu->iparm[IPARM_MODIFY_PARAMETER] = API_NO; 365 lu->iparm[IPARM_START_TASK] = API_TASK_INIT; 366 lu->iparm[IPARM_END_TASK] = API_TASK_INIT; 367 lu->rhsnbr = 1; 368 369 /* Call to set default pastix options */ 370 pastix((pastix_data_t **)&(lu->pastix_data), 371 (MPI_Comm) lu->pastix_comm, 372 (pastix_int_t) lu->n, 373 (pastix_int_t*) lu->colptr, 374 (pastix_int_t*) lu->row, 375 (pastix_float_t*) lu->val, 376 (pastix_int_t*) lu->perm, 377 (pastix_int_t*) lu->invp, 378 (pastix_float_t*) lu->rhs, 379 (pastix_int_t) lu->rhsnbr, 380 (pastix_int_t*) lu->iparm, 381 (double*) lu->dparm); 382 383 ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"PaStiX Options","Mat");CHKERRQ(ierr); 384 385 icntl=-1; 386 lu->iparm[IPARM_VERBOSE] = API_VERBOSE_NO; 387 ierr = PetscOptionsInt("-mat_pastix_verbose","iparm[IPARM_VERBOSE] : level of printing (0 to 2)","None", 388 lu->iparm[IPARM_VERBOSE],&icntl,&flg);CHKERRQ(ierr); 389 if ((flg && icntl > 0) || PetscLogPrintInfo) { 390 lu->iparm[IPARM_VERBOSE] = icntl; 391 } 392 icntl=-1; 393 ierr = PetscOptionsInt("-mat_pastix_threadnbr","iparm[IPARM_THREAD_NBR] : Number of thread by MPI node", 394 "None",lu->iparm[IPARM_THREAD_NBR],&icntl,PETSC_NULL);CHKERRQ(ierr); 395 if ((flg && icntl > 0)) { 396 lu->iparm[IPARM_THREAD_NBR] = icntl; 397 } 398 PetscOptionsEnd(); 399 valOnly = PETSC_FALSE; 400 } 401 else { 402 valOnly = PETSC_TRUE; 403 } 404 405 lu->iparm[IPARM_MATRIX_VERIFICATION] = API_YES; 406 407 /* convert mpi A to seq mat A */ 408 ierr = ISCreateStride(PETSC_COMM_SELF,M,0,1,&isrow);CHKERRQ(ierr); 409 ierr = MatGetSubMatrices(A,1,&isrow,&isrow,MAT_INITIAL_MATRIX,&tseq);CHKERRQ(ierr); 410 ierr = ISDestroy(isrow);CHKERRQ(ierr); 411 A_seq = *tseq; 412 ierr = PetscFree(tseq);CHKERRQ(ierr); 413 414 ierr = MatConvertToCSC(A_seq,valOnly, &lu->n, &lu->colptr, &lu->row, &lu->val); CHKERRQ(ierr); 415 ierr = PetscMalloc((lu->n)*sizeof(PetscInt) ,&(lu->perm));CHKERRQ(ierr); 416 ierr = PetscMalloc((lu->n)*sizeof(PetscInt) ,&(lu->invp));CHKERRQ(ierr); 417 418 MatIsSymmetric_SeqAIJ(A_seq,0.0,&isSym); 419 420 if (isSym) { 421 /* On symmetric matrix, LLT */ 422 lu->iparm[IPARM_SYM] = API_SYM_YES; 423 lu->iparm[IPARM_FACTORIZATION] = API_FACT_LLT; 424 } 425 else { 426 /* On unsymmetric matrix, LU */ 427 lu->iparm[IPARM_SYM] = API_SYM_NO; 428 lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU; 429 } 430 431 /*----------------*/ 432 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 433 if (!(isSeqAIJ || isSeqSBAIJ)) { 434 /* PaStiX only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 435 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 436 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 437 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 438 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 439 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 440 441 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 442 ierr = VecScatterCreate(lu->b_seq,is_iden,b,is_iden,&lu->scat_sol);CHKERRQ(ierr); 443 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 444 ierr = VecDestroy(b);CHKERRQ(ierr); 445 } 446 lu->iparm[IPARM_START_TASK] = API_TASK_ORDERING; 447 lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT; 448 449 pastix((pastix_data_t **)&(lu->pastix_data), 450 (MPI_Comm) lu->pastix_comm, 451 (pastix_int_t) lu->n, 452 (pastix_int_t*) lu->colptr, 453 (pastix_int_t*) lu->row, 454 (pastix_float_t*) lu->val, 455 (pastix_int_t*) lu->perm, 456 (pastix_int_t*) lu->invp, 457 (pastix_float_t*) lu->rhs, 458 (pastix_int_t) lu->rhsnbr, 459 (pastix_int_t*) lu->iparm, 460 (double*) lu->dparm); 461 if (lu->iparm[IPARM_ERROR_NUMBER] < 0) { 462 SETERRQ1(PETSC_ERR_LIB,"Error reported by PaStiX in analysis phase: ipparm(IPARM_ERROR_NUMBER)=%d\n", 463 lu->iparm[IPARM_ERROR_NUMBER]); 464 } 465 } 466 else { 467 lu->iparm[IPARM_START_TASK] = API_TASK_NUMFACT; 468 lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT; 469 pastix((pastix_data_t **)&(lu->pastix_data), 470 (MPI_Comm) lu->pastix_comm, 471 (pastix_int_t) lu->n, 472 (pastix_int_t*) lu->colptr, 473 (pastix_int_t*) lu->row, 474 (pastix_float_t*) lu->val, 475 (pastix_int_t*) lu->perm, 476 (pastix_int_t*) lu->invp, 477 (pastix_float_t*) lu->rhs, 478 (pastix_int_t) lu->rhsnbr, 479 (pastix_int_t*) lu->iparm, 480 (double*) lu->dparm); 481 482 if (lu->iparm[IPARM_ERROR_NUMBER] < 0) { 483 SETERRQ1(PETSC_ERR_LIB,"Error reported by PaStiX in analysis phase: ipparm(IPARM_ERROR_NUMBER)=%d\n", 484 lu->iparm[IPARM_ERROR_NUMBER]); 485 } 486 } 487 488 if (lu->commSize > 1){ 489 if ((F)->factor == MAT_FACTOR_LU){ 490 F_diag = ((Mat_MPIAIJ *)(F)->data)->A; 491 } else { 492 F_diag = ((Mat_MPISBAIJ *)(F)->data)->A; 493 } 494 F_diag->assembled = PETSC_TRUE; 495 if (lu->nSolve){ 496 ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr); 497 ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr); 498 } 499 } 500 (F)->assembled = PETSC_TRUE; 501 lu->matstruc = SAME_NONZERO_PATTERN; 502 lu->CleanUpPastix = PETSC_TRUE; 503 lu->nSolve = 0; 504 PetscFunctionReturn(0); 505 } 506 507 508 /* Note the Petsc r and c permutations are ignored */ 509 #undef __FUNCT__ 510 #define __FUNCT__ "MatLUFactorSymbolic_AIJPASTIX" 511 PetscErrorCode MatLUFactorSymbolic_AIJPASTIX(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 512 { 513 Mat_Pastix *lu = (Mat_Pastix*)F->spptr; 514 PetscErrorCode ierr = 0; 515 516 PetscFunctionBegin; 517 lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU; 518 lu->iparm[IPARM_SYM] = API_SYM_YES; 519 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 520 F->ops->lufactornumeric = MatFactorNumeric_PaStiX; 521 PetscFunctionReturn(0); 522 } 523 524 525 /* Note the Petsc r permutation is ignored */ 526 #undef __FUNCT__ 527 #define __FUNCT__ "MatCholeskyFactorSymbolic_SBAIJPASTIX" 528 PetscErrorCode MatCholeskyFactorSymbolic_SBAIJPASTIX(Mat F,Mat A,IS r,const MatFactorInfo *info) 529 { 530 Mat_Pastix *lu = (Mat_Pastix*)(F)->spptr; 531 PetscErrorCode ierr = 0; 532 533 PetscFunctionBegin; 534 lu->iparm[IPARM_FACTORIZATION] = API_FACT_LLT; 535 lu->iparm[IPARM_SYM] = API_SYM_NO; 536 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 537 (F)->ops->choleskyfactornumeric = MatFactorNumeric_PaStiX; 538 #if !defined(PETSC_USE_COMPLEX) 539 (F)->ops->getinertia = MatGetInertia_SBAIJPASTIX; 540 #endif 541 PetscFunctionReturn(0); 542 } 543 544 #undef __FUNCT__ 545 #define __FUNCT__ "MatFactorInfo_PaStiX" 546 PetscErrorCode MatFactorInfo_PaStiX(Mat A,PetscViewer viewer) { 547 Mat_Pastix *lu=(Mat_Pastix*)A->spptr; 548 PetscErrorCode ierr; 549 550 PetscFunctionBegin; 551 /* check if matrix is pastix type */ 552 if (A->ops->solve != MatSolve_PaStiX) PetscFunctionReturn(0); 553 554 ierr = PetscViewerASCIIPrintf(viewer,"PaStiX run parameters:\n");CHKERRQ(ierr); 555 ierr = PetscViewerASCIIPrintf(viewer," Matrix type : %s \n", 556 ((lu->iparm[IPARM_SYM] == API_SYM_YES)?"Symmetric":"Unsymmetric"));CHKERRQ(ierr); 557 ierr = PetscViewerASCIIPrintf(viewer," Level of printing (0,1,2): %d \n", 558 lu->iparm[IPARM_VERBOSE]);CHKERRQ(ierr); 559 ierr = PetscViewerASCIIPrintf(viewer," Number of refinements iterations : %d \n", 560 lu->iparm[IPARM_NBITER]);CHKERRQ(ierr); 561 ierr = PetscPrintf(PETSC_COMM_SELF," Error : %g \n", 562 lu->dparm[DPARM_RELATIVE_ERROR]);CHKERRQ(ierr); 563 PetscFunctionReturn(0); 564 } 565 566 #undef __FUNCT__ 567 #define __FUNCT__ "MatView_PaStiX" 568 PetscErrorCode MatView_PaStiX(Mat A,PetscViewer viewer) 569 { 570 PetscErrorCode ierr; 571 PetscTruth iascii; 572 PetscViewerFormat format; 573 574 PetscFunctionBegin; 575 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 576 if (iascii) { 577 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 578 if (format == PETSC_VIEWER_ASCII_INFO){ 579 ierr = MatFactorInfo_PaStiX(A,viewer);CHKERRQ(ierr); 580 } 581 } 582 PetscFunctionReturn(0); 583 } 584 585 586 /*MC 587 MATAIJPASTIX - MATAIJPASTIX = "aijpastix" - A matrix type providing direct solvers (LU) for distributed 588 and sequential matrices via the external package PaStiX. 589 590 If PaStiX is installed (see the manual for instructions 591 on how to declare the existence of external packages), 592 a matrix type can be constructed which invokes PaStiX solvers. 593 After calling MatCreate(...,A), simply call MatSetType(A,MATAIJPASTIX), then 594 optionally call MatSeqAIJSetPreallocation() or MatMPIAIJSetPreallocation() etc DO NOT 595 call MatCreateSeqAIJ/MPIAIJ() directly or the preallocation information will be LOST! 596 597 If created with a single process communicator, this matrix type inherits from MATSEQAIJ. 598 Otherwise, this matrix type inherits from MATMPIAIJ. Hence for single process communicators, 599 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 600 for communicators controlling multiple processes. It is recommended that you call both of 601 the above preallocation routines for simplicity. One can also call MatConvert() for an inplace 602 conversion to or from the MATSEQAIJ or MATMPIAIJ type (depending on the communicator size) 603 without data copy AFTER the matrix values are set. 604 605 Options Database Keys: 606 + -mat_type sbaijpastix - sets the matrix type to "sbaijpastix" during a call to MatSetFromOptions() 607 . -mat_pastix_sym <0,1> - 0 the matrix is symmetric, 1 unsymmetric 608 . -mat_pastix_verbose <0,1,2> - print level 609 . -mat_pastix_threadnbr <integer> - Set the thread number by MPI task. 610 611 Level: beginner 612 613 .seealso: MATSBAIJPASTIX 614 M*/ 615 616 617 #undef __FUNCT__ 618 #define __FUNCT__ "MatGetInfo_PaStiX" 619 PetscErrorCode MatGetInfo_PaStiX(Mat A,MatInfoType flag,MatInfo *info) 620 { 621 Mat_Pastix *lu =(Mat_Pastix*)A->spptr; 622 PetscErrorCode ierr; 623 624 PetscFunctionBegin; 625 info->block_size = 1.0; 626 info->nz_allocated = lu->iparm[IPARM_NNZEROS]; 627 info->nz_used = lu->iparm[IPARM_NNZEROS]; 628 info->nz_unneeded = 0.0; 629 info->assemblies = 0.0; 630 info->mallocs = 0.0; 631 info->memory = 0.0; 632 info->fill_ratio_given = 0; 633 info->fill_ratio_needed = 0; 634 info->factor_mallocs = 0; 635 PetscFunctionReturn(0); 636 } 637 638 /*MC 639 MATSBAIJPASTIX - MATSBAIJPASTIX = "sbaijpastix" - A symmetric matrix type providing direct solvers (Cholesky) for 640 distributed and sequential matrices via the external package PaStiX. 641 642 If PaStiX is installed (see the manual for instructions 643 on how to declare the existence of external packages), 644 a matrix type can be constructed which invokes PaStiX solvers. 645 After calling MatCreate(...,A), simply call MatSetType(A,MATSBAIJPASTIX), then 646 optionally call MatSeqSBAIJSetPreallocation() or MatMPISBAIJSetPreallocation() DO NOT 647 call MatCreateSeqSBAIJ/MPISBAIJ() directly or the preallocation information will be LOST! 648 649 If created with a single process communicator, this matrix type inherits from MATSEQSBAIJ. 650 Otherwise, this matrix type inherits from MATMPISBAIJ. Hence for single process communicators, 651 MatSeqSBAIJSetPreallocation() is supported, and similarly MatMPISBAIJSetPreallocation() is supported 652 for communicators controlling multiple processes. It is recommended that you call both of 653 the above preallocation routines for simplicity. One can also call MatConvert() for an inplace 654 conversion to or from the MATSEQSBAIJ or MATMPISBAIJ type (depending on the communicator size) 655 without data copy AFTER the matrix values have been set. 656 657 Options Database Keys: 658 + -mat_type aijpastix - sets the matrix type to "aijpastix" during a call to MatSetFromOptions() 659 . -mat_pastix_sym <0,1> - 0 the matrix is symmetric, 1 unsymmetric 660 . -mat_pastix_verbose <0,1,2> - print level 661 . -mat_pastix_threadnbr <integer> - Set the thread number by MPI task. 662 Level: beginner 663 664 .seealso: MATAIJPASTIX 665 M*/ 666 667 EXTERN_C_BEGIN 668 #undef __FUNCT__ 669 #define __FUNCT__ "MatFactorGetSolverPackage_pastix" 670 PetscErrorCode MatFactorGetSolverPackage_pastix(Mat A,const MatSolverPackage *type) 671 { 672 PetscErrorCode ierr = 0; 673 674 PetscFunctionBegin; 675 *type = MAT_SOLVER_PASTIX; 676 PetscFunctionReturn(0); 677 } 678 EXTERN_C_END 679 680 EXTERN_C_BEGIN 681 /* 682 The seq and mpi versions of this function are the same 683 */ 684 #undef __FUNCT__ 685 #define __FUNCT__ "MatGetFactor_seqaij_pastix" 686 PetscErrorCode MatGetFactor_seqaij_pastix(Mat A,MatFactorType ftype,Mat *F) 687 { 688 Mat B; 689 PetscErrorCode ierr; 690 Mat_Pastix *pastix; 691 692 PetscFunctionBegin; 693 if (ftype != MAT_FACTOR_LU) { 694 SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix"); 695 } 696 /* Create the factorization matrix */ 697 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 698 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 699 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 700 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 701 702 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX; 703 B->ops->view = MatView_PaStiX; 704 B->ops->getinfo = MatGetInfo_PaStiX; 705 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C", 706 "MatFactorGetSolverPackage_pastix", 707 MatFactorGetSolverPackage_pastix);CHKERRQ(ierr); 708 B->factor = MAT_FACTOR_LU; 709 710 ierr = PetscNewLog(B,Mat_Pastix,&pastix);CHKERRQ(ierr); 711 pastix->CleanUpPastix = PETSC_FALSE; 712 pastix->isAIJ = PETSC_TRUE; 713 pastix->scat_rhs = PETSC_NULL; 714 pastix->scat_sol = PETSC_NULL; 715 pastix->nSolve = 0; 716 pastix->MatDestroy = B->ops->destroy; 717 B->ops->destroy = MatDestroy_Pastix; 718 B->spptr = (void*)pastix; 719 720 *F = B; 721 PetscFunctionReturn(0); 722 } 723 EXTERN_C_END 724 725 EXTERN_C_BEGIN 726 #undef __FUNCT__ 727 #define __FUNCT__ "MatGetFactor_mpiaij_pastix" 728 PetscErrorCode MatGetFactor_mpiaij_pastix(Mat A,MatFactorType ftype,Mat *F) 729 { 730 Mat B; 731 PetscErrorCode ierr; 732 Mat_Pastix *pastix; 733 734 PetscFunctionBegin; 735 if (ftype != MAT_FACTOR_LU) { 736 SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix"); 737 } 738 /* Create the factorization matrix */ 739 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 740 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 741 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 742 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 743 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 744 745 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX; 746 B->ops->view = MatView_PaStiX; 747 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B, 748 "MatFactorGetSolverPackage_C", 749 "MatFactorGetSolverPackage_pastix", 750 MatFactorGetSolverPackage_pastix);CHKERRQ(ierr); 751 B->factor = MAT_FACTOR_LU; 752 753 ierr = PetscNewLog(B,Mat_Pastix,&pastix);CHKERRQ(ierr); 754 pastix->CleanUpPastix = PETSC_FALSE; 755 pastix->isAIJ = PETSC_TRUE; 756 pastix->scat_rhs = PETSC_NULL; 757 pastix->scat_sol = PETSC_NULL; 758 pastix->nSolve = 0; 759 pastix->MatDestroy = B->ops->destroy; 760 B->ops->destroy = MatDestroy_Pastix; 761 B->spptr = (void*)pastix; 762 763 *F = B; 764 PetscFunctionReturn(0); 765 } 766 EXTERN_C_END 767 768 EXTERN_C_BEGIN 769 #undef __FUNCT__ 770 #define __FUNCT__ "MatGetFactor_seqsbaij_pastix" 771 PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat A,MatFactorType ftype,Mat *F) 772 { 773 Mat B; 774 PetscErrorCode ierr; 775 Mat_Pastix *pastix; 776 777 PetscFunctionBegin; 778 if (ftype != MAT_FACTOR_CHOLESKY) { 779 SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix"); 780 } 781 /* Create the factorization matrix */ 782 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 783 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 784 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 785 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 786 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 787 788 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX; 789 B->ops->view = MatView_PaStiX; 790 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B, 791 "MatFactorGetSolverPackage_C", 792 "MatFactorGetSolverPackage_pastix", 793 MatFactorGetSolverPackage_pastix);CHKERRQ(ierr); 794 795 B->factor = MAT_FACTOR_CHOLESKY; 796 797 ierr = PetscNewLog(B,Mat_Pastix,&pastix);CHKERRQ(ierr); 798 pastix->CleanUpPastix = PETSC_FALSE; 799 pastix->isAIJ = PETSC_TRUE; 800 pastix->scat_rhs = PETSC_NULL; 801 pastix->scat_sol = PETSC_NULL; 802 pastix->nSolve = 0; 803 pastix->MatDestroy = B->ops->destroy; 804 B->ops->destroy = MatDestroy_Pastix; 805 B->spptr = (void*)pastix; 806 807 *F = B; 808 PetscFunctionReturn(0); 809 } 810 EXTERN_C_END 811 812 EXTERN_C_BEGIN 813 #undef __FUNCT__ 814 #define __FUNCT__ "MatGetFactor_mpisbaij_pastix" 815 PetscErrorCode MatGetFactor_mpisbaij_pastix(Mat A,MatFactorType ftype,Mat *F) 816 { 817 Mat B; 818 PetscErrorCode ierr; 819 Mat_Pastix *pastix; 820 821 PetscFunctionBegin; 822 if (ftype != MAT_FACTOR_CHOLESKY) { 823 SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix"); 824 } 825 /* Create the factorization matrix */ 826 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 827 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 828 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 829 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 830 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 831 832 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX; 833 B->ops->view = MatView_PaStiX; 834 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B, 835 "MatFactorGetSolverPackage_C", 836 "MatFactorGetSolverPackage_pastix", 837 MatFactorGetSolverPackage_pastix);CHKERRQ(ierr); 838 B->factor = MAT_FACTOR_CHOLESKY; 839 840 ierr = PetscNewLog(B,Mat_Pastix,&pastix);CHKERRQ(ierr); 841 pastix->CleanUpPastix = PETSC_FALSE; 842 pastix->isAIJ = PETSC_TRUE; 843 pastix->scat_rhs = PETSC_NULL; 844 pastix->scat_sol = PETSC_NULL; 845 pastix->nSolve = 0; 846 pastix->MatDestroy = B->ops->destroy; 847 B->ops->destroy = MatDestroy_Pastix; 848 B->spptr = (void*)pastix; 849 850 *F = B; 851 PetscFunctionReturn(0); 852 } 853 EXTERN_C_END 854