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