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