1 #include <../src/mat/impls/aij/seq/aij.h> 2 #include <../src/mat/impls/aij/mpi/mpiaij.h> 3 #include <StrumpackSparseSolver.h> 4 5 /* 6 GLOBAL - The sparse matrix and right hand side are all stored initially on process 0. Should be called centralized 7 DISTRIBUTED - The sparse matrix and right hand size are initially stored across the entire MPI communicator. 8 */ 9 typedef enum {GLOBAL, DISTRIBUTED} STRUMPACK_MatInputMode; 10 const char *STRUMPACK_MatInputModes[] = {"GLOBAL","DISTRIBUTED","STRUMPACK_MatInputMode","PETSC_",0}; 11 12 typedef struct { 13 STRUMPACK_SparseSolver S; 14 STRUMPACK_MatInputMode MatInputMode; 15 } STRUMPACK_data; 16 17 extern PetscErrorCode MatFactorInfo_STRUMPACK_MPI(Mat,PetscViewer); 18 extern PetscErrorCode MatLUFactorNumeric_STRUMPACK_MPI(Mat,Mat,const MatFactorInfo*); 19 extern PetscErrorCode MatDestroy_STRUMPACK_MPI(Mat); 20 extern PetscErrorCode MatView_STRUMPACK_MPI(Mat,PetscViewer); 21 extern PetscErrorCode MatSolve_STRUMPACK_MPI(Mat,Vec,Vec); 22 extern PetscErrorCode MatLUFactorSymbolic_STRUMPACK_MPI(Mat,Mat,IS,IS,const MatFactorInfo*); 23 extern PetscErrorCode MatDestroy_MPIAIJ(Mat); 24 25 #undef __FUNCT__ 26 #define __FUNCT__ "MatGetDiagonal_STRUMPACK_MPI" 27 PetscErrorCode MatGetDiagonal_STRUMPACK_MPI(Mat A,Vec v) 28 { 29 PetscFunctionBegin; 30 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Mat type: STRUMPACK_MPI factor"); 31 PetscFunctionReturn(0); 32 } 33 34 #undef __FUNCT__ 35 #define __FUNCT__ "MatDestroy_STRUMPACK_MPI" 36 PetscErrorCode MatDestroy_STRUMPACK_MPI(Mat A) 37 { 38 STRUMPACK_data *sp = (STRUMPACK_data*)A->spptr; 39 PetscErrorCode ierr; 40 PetscBool flg; 41 42 PetscFunctionBegin; 43 /* Deallocate STRUMPACK_MPI storage */ 44 PetscStackCall("STRUMPACK_destroy",STRUMPACK_destroy(&(sp->S))); 45 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 46 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 47 if (flg) { 48 ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr); 49 } else { 50 ierr = MatDestroy_MPIAIJ(A);CHKERRQ(ierr); 51 } 52 PetscFunctionReturn(0); 53 } 54 55 #undef __FUNCT__ 56 #define __FUNCT__ "MatSolve_STRUMPACK_MPI" 57 PetscErrorCode MatSolve_STRUMPACK_MPI(Mat A,Vec b_mpi,Vec x) 58 { 59 STRUMPACK_data *sp = (STRUMPACK_data*)A->spptr; 60 STRUMPACK_RETURN_CODE sp_err; 61 PetscErrorCode ierr; 62 PetscMPIInt size; 63 PetscInt N=A->cmap->N; 64 const PetscScalar *bptr; 65 PetscScalar *xptr; 66 Vec x_seq,b_seq; 67 IS iden; 68 VecScatter scat; 69 70 PetscFunctionBegin; 71 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 72 if (size > 1 && sp->MatInputMode == GLOBAL) { 73 ierr = VecCreateSeq(PETSC_COMM_SELF,N,&x_seq);CHKERRQ(ierr); 74 ierr = VecGetArray(x_seq,&xptr);CHKERRQ(ierr); 75 /* global mat input, convert b to b_seq */ 76 ierr = VecCreateSeq(PETSC_COMM_SELF,N,&b_seq);CHKERRQ(ierr); 77 ierr = ISCreateStride(PETSC_COMM_SELF,N,0,1,&iden);CHKERRQ(ierr); 78 ierr = VecScatterCreate(b_mpi,iden,b_seq,iden,&scat);CHKERRQ(ierr); 79 ierr = ISDestroy(&iden);CHKERRQ(ierr); 80 ierr = VecScatterBegin(scat,b_mpi,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 81 ierr = VecScatterEnd(scat,b_mpi,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 82 ierr = VecGetArrayRead(b_seq,&bptr);CHKERRQ(ierr); 83 } else { /* size==1 || distributed mat input */ 84 ierr = VecGetArray(x,&xptr);CHKERRQ(ierr); 85 ierr = VecGetArrayRead(b_mpi,&bptr);CHKERRQ(ierr); 86 } 87 88 PetscStackCall("STRUMPACK_solve",sp_err = STRUMPACK_solve(sp->S,(PetscScalar*)bptr,xptr,0)); 89 90 if (sp_err != STRUMPACK_SUCCESS) { 91 if (sp_err == STRUMPACK_MATRIX_NOT_SET) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); 92 else if (sp_err == STRUMPACK_REORDERING_ERROR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); 93 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: solve failed"); 94 } 95 } 96 97 if (size > 1 && sp->MatInputMode == GLOBAL) { 98 ierr = VecRestoreArrayRead(b_seq,&bptr);CHKERRQ(ierr); 99 ierr = VecDestroy(&b_seq);CHKERRQ(ierr); 100 /* convert seq x to mpi x */ 101 ierr = VecRestoreArray(x_seq,&xptr);CHKERRQ(ierr); 102 ierr = VecScatterBegin(scat,x_seq,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 103 ierr = VecScatterEnd(scat,x_seq,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 104 ierr = VecScatterDestroy(&scat);CHKERRQ(ierr); 105 ierr = VecDestroy(&x_seq);CHKERRQ(ierr); 106 } else { 107 ierr = VecRestoreArray(x,&xptr);CHKERRQ(ierr); 108 ierr = VecRestoreArrayRead(b_mpi,&bptr);CHKERRQ(ierr); 109 } 110 111 PetscFunctionReturn(0); 112 } 113 114 #undef __FUNCT__ 115 #define __FUNCT__ "MatMatSolve_STRUMPACK_MPI" 116 PetscErrorCode MatMatSolve_STRUMPACK_MPI(Mat A,Mat B_mpi,Mat X) 117 { 118 PetscErrorCode ierr; 119 PetscBool flg; 120 121 PetscFunctionBegin; 122 ierr = PetscObjectTypeCompareAny((PetscObject)B_mpi,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 123 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix"); 124 ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 125 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix"); 126 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_STRUMPACK_MPI() is not implemented yet"); 127 PetscFunctionReturn(0); 128 } 129 130 131 #undef __FUNCT__ 132 #define __FUNCT__ "MatLUFactorNumeric_STRUMPACK_MPI" 133 PetscErrorCode MatLUFactorNumeric_STRUMPACK_MPI(Mat F,Mat A,const MatFactorInfo *info) 134 { 135 STRUMPACK_data *sp = (STRUMPACK_data*)F->spptr; 136 STRUMPACK_RETURN_CODE sp_err; 137 Mat *tseq,A_seq = NULL; 138 Mat_SeqAIJ *A_d,*A_o; 139 PetscErrorCode ierr; 140 PetscInt M=A->rmap->N,m=A->rmap->n,N=A->cmap->N; 141 PetscMPIInt size; 142 IS isrow; 143 144 PetscFunctionBegin; 145 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 146 147 if (sp->MatInputMode == GLOBAL) { /* global mat input */ 148 if (size > 1) { /* convert mpi A to seq mat A */ 149 ierr = ISCreateStride(PETSC_COMM_SELF,M,0,1,&isrow);CHKERRQ(ierr); 150 ierr = MatGetSubMatrices(A,1,&isrow,&isrow,MAT_INITIAL_MATRIX,&tseq);CHKERRQ(ierr); 151 ierr = ISDestroy(&isrow);CHKERRQ(ierr); 152 A_seq = *tseq; 153 ierr = PetscFree(tseq);CHKERRQ(ierr); 154 A_d = (Mat_SeqAIJ*)A_seq->data; 155 } else { 156 PetscBool flg; 157 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr); 158 if (flg) { 159 Mat_MPIAIJ *At = (Mat_MPIAIJ*)A->data; 160 A = At->A; 161 } 162 A_d = (Mat_SeqAIJ*)A->data; 163 } 164 PetscStackCall("STRUMPACK_set_csr_matrix",STRUMPACK_set_csr_matrix(sp->S,&N,A_d->i,A_d->j,A_d->a,0)); 165 } else { /* distributed mat input */ 166 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 167 A_d = (Mat_SeqAIJ*)(mat->A)->data; 168 A_o = (Mat_SeqAIJ*)(mat->B)->data; 169 PetscStackCall("STRUMPACK_set_MPIAIJ_matrix",STRUMPACK_set_MPIAIJ_matrix(sp->S,&m,A_d->i,A_d->j,A_d->a,A_o->i,A_o->j,A_o->a,mat->garray)); 170 } 171 172 /* Reorder and Factor the matrix. */ 173 /* TODO figure out how to avoid reorder if the matrix values changed, but the pattern remains the same. */ 174 PetscStackCall("STRUMPACK_reorder",sp_err = STRUMPACK_reorder(sp->S)); 175 PetscStackCall("STRUMPACK_factor",sp_err = STRUMPACK_factor(sp->S)); 176 if (sp_err != STRUMPACK_SUCCESS) { 177 if (sp_err == STRUMPACK_MATRIX_NOT_SET) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); 178 else if (sp_err == STRUMPACK_REORDERING_ERROR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); 179 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: factorization failed"); 180 } 181 182 if (sp->MatInputMode == GLOBAL && size > 1) { 183 ierr = MatDestroy(&A_seq);CHKERRQ(ierr); 184 } 185 186 PetscFunctionReturn(0); 187 } 188 189 #undef __FUNCT__ 190 #define __FUNCT__ "MatLUFactorSymbolic_STRUMPACK_MPI" 191 PetscErrorCode MatLUFactorSymbolic_STRUMPACK_MPI(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 192 { 193 PetscFunctionBegin; 194 F->ops->lufactornumeric = MatLUFactorNumeric_STRUMPACK_MPI; 195 F->ops->solve = MatSolve_STRUMPACK_MPI; 196 F->ops->matsolve = MatMatSolve_STRUMPACK_MPI; 197 PetscFunctionReturn(0); 198 } 199 200 #undef __FUNCT__ 201 #define __FUNCT__ "MatFactorGetSolverPackage_aij_strumpack_mpi" 202 PetscErrorCode MatFactorGetSolverPackage_aij_strumpack_mpi(Mat A,const MatSolverPackage *type) 203 { 204 PetscFunctionBegin; 205 *type = MATSOLVERSTRUMPACK_MPI; 206 PetscFunctionReturn(0); 207 } 208 209 #undef __FUNCT__ 210 #define __FUNCT__ "MatGetFactor_aij_strumpack_mpi" 211 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_strumpack_mpi(Mat A,MatFactorType ftype,Mat *F) 212 { 213 Mat B; 214 STRUMPACK_data *sp; 215 PetscErrorCode ierr; 216 PetscInt M=A->rmap->N,N=A->cmap->N; 217 PetscMPIInt size; 218 STRUMPACK_INTERFACE iface; 219 PetscBool verb; 220 PetscInt argc; 221 char **args; 222 223 PetscFunctionBegin; 224 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 225 226 /* Create the factorization matrix */ 227 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 228 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,M,N);CHKERRQ(ierr); 229 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 230 ierr = MatSeqAIJSetPreallocation(B,0,NULL); 231 ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr); 232 233 B->ops->lufactorsymbolic = MatLUFactorSymbolic_STRUMPACK_MPI; 234 B->ops->view = MatView_STRUMPACK_MPI; 235 B->ops->destroy = MatDestroy_STRUMPACK_MPI; 236 B->ops->getdiagonal = MatGetDiagonal_STRUMPACK_MPI; 237 238 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_aij_strumpack_mpi);CHKERRQ(ierr); 239 240 B->factortype = MAT_FACTOR_LU; 241 242 ierr = PetscNewLog(B,&sp);CHKERRQ(ierr); 243 B->spptr = sp; 244 245 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"STRUMPACK_MPI Options","Mat");CHKERRQ(ierr); 246 247 sp->MatInputMode = DISTRIBUTED; 248 iface = STRUMPACK_MPI_DIST; 249 250 ierr = PetscOptionsEnum("-mat_strumpack_mpi_matinput","Matrix input mode (global or distributed)","None",STRUMPACK_MatInputModes, 251 (PetscEnum)sp->MatInputMode,(PetscEnum*)&sp->MatInputMode,NULL);CHKERRQ(ierr); 252 if (sp->MatInputMode == DISTRIBUTED && size == 1) sp->MatInputMode = GLOBAL; 253 if (sp->MatInputMode == DISTRIBUTED) { 254 iface = STRUMPACK_MPI_DIST; 255 } else if (sp->MatInputMode == GLOBAL) { 256 iface = STRUMPACK_MPI; 257 } 258 if (PetscLogPrintInfo) verb = PETSC_TRUE; 259 else verb = PETSC_FALSE; 260 ierr = PetscOptionsBool("-mat_strumpack_verbose","Print STRUMPACK information","None",verb,&verb,NULL);CHKERRQ(ierr); 261 262 PetscOptionsEnd(); 263 264 ierr = PetscGetArgs(&argc,&args);CHKERRQ(ierr); 265 266 #if defined(PETSC_USE_64BIT_INDICES) 267 #if defined(PETSC_USE_COMPLEX) 268 #if defined(PETSC_USE_REAL_SINGLE) 269 PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_FLOATCOMPLEX_64,iface,argc,args,verb)); 270 #else 271 PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_DOUBLECOMPLEX_64,iface,argc,args,verb)); 272 #endif 273 #else 274 #if defined(PETSC_USE_REAL_SINGLE) 275 PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_FLOAT_64,iface,argc,args,verb)); 276 #else 277 PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_DOUBLE_64,iface,argc,args,verb)); 278 #endif 279 #endif 280 #else 281 #if defined(PETSC_USE_COMPLEX) 282 #if defined(PETSC_USE_REAL_SINGLE) 283 PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_FLOATCOMPLEX,iface,argc,args,verb)); 284 #else 285 PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_DOUBLECOMPLEX,iface,argc,args,verb)); 286 #endif 287 #else 288 #if defined(PETSC_USE_REAL_SINGLE) 289 PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_FLOAT,iface,argc,args,verb)); 290 #else 291 PetscStackCall("STRUMPACK_init",STRUMPACK_init(&(sp->S),PetscObjectComm((PetscObject)A),STRUMPACK_DOUBLE,iface,argc,args,verb)); 292 #endif 293 #endif 294 #endif 295 PetscStackCall("STRUMPACK_set_from_options",STRUMPACK_set_from_options(sp->S)); 296 297 *F = B; 298 PetscFunctionReturn(0); 299 } 300 301 #undef __FUNCT__ 302 #define __FUNCT__ "MatSolverPackageRegister_STRUMPACK_MPI" 303 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_STRUMPACK_MPI(void) 304 { 305 PetscErrorCode ierr; 306 307 PetscFunctionBegin; 308 ierr = MatSolverPackageRegister(MATSOLVERSTRUMPACK_MPI,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack_mpi);CHKERRQ(ierr); 309 ierr = MatSolverPackageRegister(MATSOLVERSTRUMPACK_MPI,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack_mpi);CHKERRQ(ierr); 310 PetscFunctionReturn(0); 311 } 312 313 #undef __FUNCT__ 314 #define __FUNCT__ "MatFactorInfo_STRUMPACK_MPI" 315 PetscErrorCode MatFactorInfo_STRUMPACK_MPI(Mat A,PetscViewer viewer) 316 { 317 PetscErrorCode ierr; 318 319 PetscFunctionBegin; 320 /* check if matrix is strumpack type */ 321 if (A->ops->solve != MatSolve_STRUMPACK_MPI) PetscFunctionReturn(0); 322 ierr = PetscViewerASCIIPrintf(viewer,"STRUMPACK_MPI sparse solver!\n");CHKERRQ(ierr); 323 PetscFunctionReturn(0); 324 } 325 326 #undef __FUNCT__ 327 #define __FUNCT__ "MatView_STRUMPACK_MPI" 328 PetscErrorCode MatView_STRUMPACK_MPI(Mat A,PetscViewer viewer) 329 { 330 PetscErrorCode ierr; 331 PetscBool iascii; 332 PetscViewerFormat format; 333 334 PetscFunctionBegin; 335 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 336 if (iascii) { 337 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 338 if (format == PETSC_VIEWER_ASCII_INFO) { 339 ierr = MatFactorInfo_STRUMPACK_MPI(A,viewer);CHKERRQ(ierr); 340 } 341 } 342 PetscFunctionReturn(0); 343 } 344 345 346 /*MC 347 MATSOLVERSSTRUMPACK_MPI - Parallel direct solver package for LU factorization 348 349 Use ./configure --download-strumpack --download-parmetis --download-metis --download-ptscotch (and ScaLAPACK?/BLACS?) to have PETSc installed with SuperLU_DIST 350 351 Use -pc_type lu -pc_factor_mat_solver_package strumpack to us this direct solver 352 353 Works with AIJ matrices 354 355 .seealso: PCLU 356 357 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 358 359 M*/ 360 361