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