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