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