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