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