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