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 PetscCheck(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), lu->pastix_comm, lu->n, lu->colptr, lu->row, (PastixScalar *)lu->val, lu->perm, lu->invp, (PastixScalar *)lu->rhs, lu->rhsnbr, lu->iparm, lu->dparm); 188 PetscCheck(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]); 189 PetscCall(PetscFree(lu->colptr)); 190 PetscCall(PetscFree(lu->row)); 191 PetscCall(PetscFree(lu->val)); 192 PetscCall(PetscFree(lu->perm)); 193 PetscCall(PetscFree(lu->invp)); 194 PetscCallMPI(MPI_Comm_free(&(lu->pastix_comm))); 195 } 196 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL)); 197 PetscCall(PetscFree(A->data)); 198 PetscFunctionReturn(0); 199 } 200 201 /* 202 Gather right-hand-side. 203 Call for Solve step. 204 Scatter solution. 205 */ 206 PetscErrorCode MatSolve_PaStiX(Mat A, Vec b, Vec x) 207 { 208 Mat_Pastix *lu = (Mat_Pastix *)A->data; 209 PetscScalar *array; 210 Vec x_seq; 211 212 PetscFunctionBegin; 213 lu->rhsnbr = 1; 214 x_seq = lu->b_seq; 215 if (lu->commSize > 1) { 216 /* PaStiX only supports centralized rhs. Scatter b into a sequential rhs vector */ 217 PetscCall(VecScatterBegin(lu->scat_rhs, b, x_seq, INSERT_VALUES, SCATTER_FORWARD)); 218 PetscCall(VecScatterEnd(lu->scat_rhs, b, x_seq, INSERT_VALUES, SCATTER_FORWARD)); 219 PetscCall(VecGetArray(x_seq, &array)); 220 } else { /* size == 1 */ 221 PetscCall(VecCopy(b, x)); 222 PetscCall(VecGetArray(x, &array)); 223 } 224 lu->rhs = array; 225 if (lu->commSize == 1) { 226 PetscCall(VecRestoreArray(x, &array)); 227 } else { 228 PetscCall(VecRestoreArray(x_seq, &array)); 229 } 230 231 /* solve phase */ 232 /*-------------*/ 233 lu->iparm[IPARM_START_TASK] = API_TASK_SOLVE; 234 lu->iparm[IPARM_END_TASK] = API_TASK_REFINE; 235 lu->iparm[IPARM_RHS_MAKING] = API_RHS_B; 236 237 PASTIX_CALL(&(lu->pastix_data), lu->pastix_comm, lu->n, lu->colptr, lu->row, (PastixScalar *)lu->val, lu->perm, lu->invp, (PastixScalar *)lu->rhs, lu->rhsnbr, lu->iparm, lu->dparm); 238 PetscCheck(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]); 239 240 if (lu->commSize == 1) { 241 PetscCall(VecRestoreArray(x, &(lu->rhs))); 242 } else { 243 PetscCall(VecRestoreArray(x_seq, &(lu->rhs))); 244 } 245 246 if (lu->commSize > 1) { /* convert PaStiX centralized solution to petsc mpi x */ 247 PetscCall(VecScatterBegin(lu->scat_sol, x_seq, x, INSERT_VALUES, SCATTER_FORWARD)); 248 PetscCall(VecScatterEnd(lu->scat_sol, x_seq, x, INSERT_VALUES, SCATTER_FORWARD)); 249 } 250 PetscFunctionReturn(0); 251 } 252 253 /* 254 Numeric factorisation using PaStiX solver. 255 256 */ 257 PetscErrorCode MatFactorNumeric_PaStiX(Mat F, Mat A, const MatFactorInfo *info) 258 { 259 Mat_Pastix *lu = (Mat_Pastix *)(F)->data; 260 Mat *tseq; 261 PetscInt icntl; 262 PetscInt M = A->rmap->N; 263 PetscBool valOnly, flg, isSym; 264 IS is_iden; 265 Vec b; 266 IS isrow; 267 PetscBool isSeqAIJ, isSeqSBAIJ, isMPIAIJ; 268 269 PetscFunctionBegin; 270 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &isSeqAIJ)); 271 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIAIJ, &isMPIAIJ)); 272 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQSBAIJ, &isSeqSBAIJ)); 273 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) { 274 (F)->ops->solve = MatSolve_PaStiX; 275 276 /* Initialize a PASTIX instance */ 277 PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)A), &(lu->pastix_comm))); 278 PetscCallMPI(MPI_Comm_rank(lu->pastix_comm, &lu->commRank)); 279 PetscCallMPI(MPI_Comm_size(lu->pastix_comm, &lu->commSize)); 280 281 /* Set pastix options */ 282 lu->iparm[IPARM_MODIFY_PARAMETER] = API_NO; 283 lu->iparm[IPARM_START_TASK] = API_TASK_INIT; 284 lu->iparm[IPARM_END_TASK] = API_TASK_INIT; 285 286 lu->rhsnbr = 1; 287 288 /* Call to set default pastix options */ 289 PASTIX_CALL(&(lu->pastix_data), lu->pastix_comm, lu->n, lu->colptr, lu->row, (PastixScalar *)lu->val, lu->perm, lu->invp, (PastixScalar *)lu->rhs, lu->rhsnbr, lu->iparm, lu->dparm); 290 PetscCheck(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]); 291 292 PetscOptionsBegin(PetscObjectComm((PetscObject)F), ((PetscObject)F)->prefix, "PaStiX Options", "Mat"); 293 icntl = -1; 294 lu->iparm[IPARM_VERBOSE] = API_VERBOSE_NOT; 295 PetscCall(PetscOptionsInt("-mat_pastix_verbose", "iparm[IPARM_VERBOSE] : level of printing (0 to 2)", "None", lu->iparm[IPARM_VERBOSE], &icntl, &flg)); 296 if ((flg && icntl >= 0) || PetscLogPrintInfo) lu->iparm[IPARM_VERBOSE] = icntl; 297 icntl = -1; 298 PetscCall(PetscOptionsInt("-mat_pastix_threadnbr", "iparm[IPARM_THREAD_NBR] : Number of thread by MPI node", "None", lu->iparm[IPARM_THREAD_NBR], &icntl, &flg)); 299 if ((flg && icntl > 0)) lu->iparm[IPARM_THREAD_NBR] = icntl; 300 PetscOptionsEnd(); 301 valOnly = PETSC_FALSE; 302 } else { 303 if (isSeqAIJ || isMPIAIJ) { 304 PetscCall(PetscFree(lu->colptr)); 305 PetscCall(PetscFree(lu->row)); 306 PetscCall(PetscFree(lu->val)); 307 valOnly = PETSC_FALSE; 308 } else valOnly = PETSC_TRUE; 309 } 310 311 lu->iparm[IPARM_MATRIX_VERIFICATION] = API_YES; 312 313 /* convert mpi A to seq mat A */ 314 PetscCall(ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &isrow)); 315 PetscCall(MatCreateSubMatrices(A, 1, &isrow, &isrow, MAT_INITIAL_MATRIX, &tseq)); 316 PetscCall(ISDestroy(&isrow)); 317 318 PetscCall(MatConvertToCSC(*tseq, valOnly, &lu->n, &lu->colptr, &lu->row, &lu->val)); 319 PetscCall(MatIsSymmetric(*tseq, 0.0, &isSym)); 320 PetscCall(MatDestroyMatrices(1, &tseq)); 321 322 if (!lu->perm) { 323 PetscCall(PetscMalloc1(lu->n, &(lu->perm))); 324 PetscCall(PetscMalloc1(lu->n, &(lu->invp))); 325 } 326 327 if (isSym) { 328 /* On symmetric matrix, LLT */ 329 lu->iparm[IPARM_SYM] = API_SYM_YES; 330 lu->iparm[IPARM_FACTORIZATION] = API_FACT_LDLT; 331 } else { 332 /* On unsymmetric matrix, LU */ 333 lu->iparm[IPARM_SYM] = API_SYM_NO; 334 lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU; 335 } 336 337 /*----------------*/ 338 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) { 339 if (!(isSeqAIJ || isSeqSBAIJ) && !lu->b_seq) { 340 /* PaStiX only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 341 PetscCall(VecCreateSeq(PETSC_COMM_SELF, A->cmap->N, &lu->b_seq)); 342 PetscCall(ISCreateStride(PETSC_COMM_SELF, A->cmap->N, 0, 1, &is_iden)); 343 PetscCall(MatCreateVecs(A, NULL, &b)); 344 PetscCall(VecScatterCreate(b, is_iden, lu->b_seq, is_iden, &lu->scat_rhs)); 345 PetscCall(VecScatterCreate(lu->b_seq, is_iden, b, is_iden, &lu->scat_sol)); 346 PetscCall(ISDestroy(&is_iden)); 347 PetscCall(VecDestroy(&b)); 348 } 349 lu->iparm[IPARM_START_TASK] = API_TASK_ORDERING; 350 lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT; 351 352 PASTIX_CALL(&(lu->pastix_data), lu->pastix_comm, lu->n, lu->colptr, lu->row, (PastixScalar *)lu->val, lu->perm, lu->invp, (PastixScalar *)lu->rhs, lu->rhsnbr, lu->iparm, lu->dparm); 353 PetscCheck(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]); 354 } else { 355 lu->iparm[IPARM_START_TASK] = API_TASK_NUMFACT; 356 lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT; 357 PASTIX_CALL(&(lu->pastix_data), lu->pastix_comm, lu->n, lu->colptr, lu->row, (PastixScalar *)lu->val, lu->perm, lu->invp, (PastixScalar *)lu->rhs, lu->rhsnbr, lu->iparm, lu->dparm); 358 PetscCheck(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]); 359 } 360 361 (F)->assembled = PETSC_TRUE; 362 lu->matstruc = SAME_NONZERO_PATTERN; 363 lu->CleanUpPastix = PETSC_TRUE; 364 PetscFunctionReturn(0); 365 } 366 367 /* Note the Petsc r and c permutations are ignored */ 368 PetscErrorCode MatLUFactorSymbolic_AIJPASTIX(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info) 369 { 370 Mat_Pastix *lu = (Mat_Pastix *)F->data; 371 372 PetscFunctionBegin; 373 lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU; 374 lu->iparm[IPARM_SYM] = API_SYM_YES; 375 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 376 F->ops->lufactornumeric = MatFactorNumeric_PaStiX; 377 PetscFunctionReturn(0); 378 } 379 380 PetscErrorCode MatCholeskyFactorSymbolic_SBAIJPASTIX(Mat F, Mat A, IS r, const MatFactorInfo *info) 381 { 382 Mat_Pastix *lu = (Mat_Pastix *)(F)->data; 383 384 PetscFunctionBegin; 385 lu->iparm[IPARM_FACTORIZATION] = API_FACT_LLT; 386 lu->iparm[IPARM_SYM] = API_SYM_NO; 387 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 388 (F)->ops->choleskyfactornumeric = MatFactorNumeric_PaStiX; 389 PetscFunctionReturn(0); 390 } 391 392 PetscErrorCode MatView_PaStiX(Mat A, PetscViewer viewer) 393 { 394 PetscBool iascii; 395 PetscViewerFormat format; 396 397 PetscFunctionBegin; 398 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 399 if (iascii) { 400 PetscCall(PetscViewerGetFormat(viewer, &format)); 401 if (format == PETSC_VIEWER_ASCII_INFO) { 402 Mat_Pastix *lu = (Mat_Pastix *)A->data; 403 404 PetscCall(PetscViewerASCIIPrintf(viewer, "PaStiX run parameters:\n")); 405 PetscCall(PetscViewerASCIIPrintf(viewer, " Matrix type : %s \n", ((lu->iparm[IPARM_SYM] == API_SYM_YES) ? "Symmetric" : "Unsymmetric"))); 406 PetscCall(PetscViewerASCIIPrintf(viewer, " Level of printing (0,1,2): %" PetscInt_FMT " \n", lu->iparm[IPARM_VERBOSE])); 407 PetscCall(PetscViewerASCIIPrintf(viewer, " Number of refinements iterations : %" PetscInt_FMT " \n", lu->iparm[IPARM_NBITER])); 408 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error : %g \n", lu->dparm[DPARM_RELATIVE_ERROR])); 409 } 410 } 411 PetscFunctionReturn(0); 412 } 413 414 /*MC 415 MATSOLVERPASTIX - A solver package providing direct solvers (LU) for distributed 416 and sequential matrices via the external package PaStiX. 417 418 Use ./configure --download-pastix --download-ptscotch to have PETSc installed with PasTiX 419 420 Use -pc_type lu -pc_factor_mat_solver_type pastix to use this direct solver 421 422 Options Database Keys: 423 + -mat_pastix_verbose <0,1,2> - print level 424 - -mat_pastix_threadnbr <integer> - Set the thread number by MPI task. 425 426 Notes: 427 This only works for matrices with symmetric nonzero structure, if you pass it a matrix with 428 nonsymmetric structure PasTiX and hence PETSc return with an error. 429 430 Level: beginner 431 432 .seealso: `PCFactorSetMatSolverType()`, `MatSolverType`, `MatGetFactor()` 433 M*/ 434 435 PetscErrorCode MatGetInfo_PaStiX(Mat A, MatInfoType flag, MatInfo *info) 436 { 437 Mat_Pastix *lu = (Mat_Pastix *)A->data; 438 439 PetscFunctionBegin; 440 info->block_size = 1.0; 441 info->nz_allocated = lu->iparm[IPARM_NNZEROS]; 442 info->nz_used = lu->iparm[IPARM_NNZEROS]; 443 info->nz_unneeded = 0.0; 444 info->assemblies = 0.0; 445 info->mallocs = 0.0; 446 info->memory = 0.0; 447 info->fill_ratio_given = 0; 448 info->fill_ratio_needed = 0; 449 info->factor_mallocs = 0; 450 PetscFunctionReturn(0); 451 } 452 453 static PetscErrorCode MatFactorGetSolverType_pastix(Mat A, MatSolverType *type) 454 { 455 PetscFunctionBegin; 456 *type = MATSOLVERPASTIX; 457 PetscFunctionReturn(0); 458 } 459 460 /* 461 The seq and mpi versions of this function are the same 462 */ 463 static PetscErrorCode MatGetFactor_seqaij_pastix(Mat A, MatFactorType ftype, Mat *F) 464 { 465 Mat B; 466 Mat_Pastix *pastix; 467 468 PetscFunctionBegin; 469 PetscCheck(ftype == MAT_FACTOR_LU, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix"); 470 /* Create the factorization matrix */ 471 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 472 PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 473 PetscCall(PetscStrallocpy("pastix", &((PetscObject)B)->type_name)); 474 PetscCall(MatSetUp(B)); 475 476 B->trivialsymbolic = PETSC_TRUE; 477 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX; 478 B->ops->view = MatView_PaStiX; 479 B->ops->getinfo = MatGetInfo_PaStiX; 480 481 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_pastix)); 482 483 B->factortype = MAT_FACTOR_LU; 484 485 /* set solvertype */ 486 PetscCall(PetscFree(B->solvertype)); 487 PetscCall(PetscStrallocpy(MATSOLVERPASTIX, &B->solvertype)); 488 489 PetscCall(PetscNew(&pastix)); 490 491 pastix->CleanUpPastix = PETSC_FALSE; 492 pastix->scat_rhs = NULL; 493 pastix->scat_sol = NULL; 494 B->ops->getinfo = MatGetInfo_External; 495 B->ops->destroy = MatDestroy_Pastix; 496 B->data = (void *)pastix; 497 498 *F = B; 499 PetscFunctionReturn(0); 500 } 501 502 static PetscErrorCode MatGetFactor_mpiaij_pastix(Mat A, MatFactorType ftype, Mat *F) 503 { 504 Mat B; 505 Mat_Pastix *pastix; 506 507 PetscFunctionBegin; 508 PetscCheck(ftype == MAT_FACTOR_LU, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix"); 509 /* Create the factorization matrix */ 510 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 511 PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 512 PetscCall(PetscStrallocpy("pastix", &((PetscObject)B)->type_name)); 513 PetscCall(MatSetUp(B)); 514 515 B->trivialsymbolic = PETSC_TRUE; 516 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX; 517 B->ops->view = MatView_PaStiX; 518 B->ops->getinfo = MatGetInfo_PaStiX; 519 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_pastix)); 520 521 B->factortype = MAT_FACTOR_LU; 522 523 /* set solvertype */ 524 PetscCall(PetscFree(B->solvertype)); 525 PetscCall(PetscStrallocpy(MATSOLVERPASTIX, &B->solvertype)); 526 527 PetscCall(PetscNew(&pastix)); 528 529 pastix->CleanUpPastix = PETSC_FALSE; 530 pastix->scat_rhs = NULL; 531 pastix->scat_sol = NULL; 532 B->ops->getinfo = MatGetInfo_External; 533 B->ops->destroy = MatDestroy_Pastix; 534 B->data = (void *)pastix; 535 536 *F = B; 537 PetscFunctionReturn(0); 538 } 539 540 static PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat A, MatFactorType ftype, Mat *F) 541 { 542 Mat B; 543 Mat_Pastix *pastix; 544 545 PetscFunctionBegin; 546 PetscCheck(ftype == MAT_FACTOR_CHOLESKY, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix"); 547 /* Create the factorization matrix */ 548 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 549 PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 550 PetscCall(PetscStrallocpy("pastix", &((PetscObject)B)->type_name)); 551 PetscCall(MatSetUp(B)); 552 553 B->trivialsymbolic = PETSC_TRUE; 554 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX; 555 B->ops->view = MatView_PaStiX; 556 B->ops->getinfo = MatGetInfo_PaStiX; 557 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_pastix)); 558 559 B->factortype = MAT_FACTOR_CHOLESKY; 560 561 /* set solvertype */ 562 PetscCall(PetscFree(B->solvertype)); 563 PetscCall(PetscStrallocpy(MATSOLVERPASTIX, &B->solvertype)); 564 565 PetscCall(PetscNew(&pastix)); 566 567 pastix->CleanUpPastix = PETSC_FALSE; 568 pastix->scat_rhs = NULL; 569 pastix->scat_sol = NULL; 570 B->ops->getinfo = MatGetInfo_External; 571 B->ops->destroy = MatDestroy_Pastix; 572 B->data = (void *)pastix; 573 *F = B; 574 PetscFunctionReturn(0); 575 } 576 577 static PetscErrorCode MatGetFactor_mpisbaij_pastix(Mat A, MatFactorType ftype, Mat *F) 578 { 579 Mat B; 580 Mat_Pastix *pastix; 581 582 PetscFunctionBegin; 583 PetscCheck(ftype == MAT_FACTOR_CHOLESKY, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix"); 584 585 /* Create the factorization matrix */ 586 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 587 PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 588 PetscCall(PetscStrallocpy("pastix", &((PetscObject)B)->type_name)); 589 PetscCall(MatSetUp(B)); 590 591 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX; 592 B->ops->view = MatView_PaStiX; 593 B->ops->getinfo = MatGetInfo_PaStiX; 594 B->ops->destroy = MatDestroy_Pastix; 595 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_pastix)); 596 597 B->factortype = MAT_FACTOR_CHOLESKY; 598 599 /* set solvertype */ 600 PetscCall(PetscFree(B->solvertype)); 601 PetscCall(PetscStrallocpy(MATSOLVERPASTIX, &B->solvertype)); 602 603 PetscCall(PetscNew(&pastix)); 604 605 pastix->CleanUpPastix = PETSC_FALSE; 606 pastix->scat_rhs = NULL; 607 pastix->scat_sol = NULL; 608 B->data = (void *)pastix; 609 610 *F = B; 611 PetscFunctionReturn(0); 612 } 613 614 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Pastix(void) 615 { 616 PetscFunctionBegin; 617 PetscCall(MatSolverTypeRegister(MATSOLVERPASTIX, MATMPIAIJ, MAT_FACTOR_LU, MatGetFactor_mpiaij_pastix)); 618 PetscCall(MatSolverTypeRegister(MATSOLVERPASTIX, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_pastix)); 619 PetscCall(MatSolverTypeRegister(MATSOLVERPASTIX, MATMPISBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_mpisbaij_pastix)); 620 PetscCall(MatSolverTypeRegister(MATSOLVERPASTIX, MATSEQSBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_seqsbaij_pastix)); 621 PetscFunctionReturn(0); 622 } 623