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