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