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