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