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