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