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