1 2 /* -------------------------------------------------------------------- 3 4 This file implements a subclass of the SeqAIJ matrix class that uses 5 the SuperLU sparse solver. 6 */ 7 8 /* 9 Defines the data structure for the base matrix type (SeqAIJ) 10 */ 11 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 12 13 /* 14 SuperLU include files 15 */ 16 EXTERN_C_BEGIN 17 #if defined(PETSC_USE_COMPLEX) 18 #if defined(PETSC_USE_REAL_SINGLE) 19 #include <slu_cdefs.h> 20 #else 21 #include <slu_zdefs.h> 22 #endif 23 #else 24 #if defined(PETSC_USE_REAL_SINGLE) 25 #include <slu_sdefs.h> 26 #else 27 #include <slu_ddefs.h> 28 #endif 29 #endif 30 #include <slu_util.h> 31 EXTERN_C_END 32 33 /* 34 This is the data that defines the SuperLU factored matrix type 35 */ 36 typedef struct { 37 SuperMatrix A,L,U,B,X; 38 superlu_options_t options; 39 PetscInt *perm_c; /* column permutation vector */ 40 PetscInt *perm_r; /* row permutations from partial pivoting */ 41 PetscInt *etree; 42 PetscReal *R, *C; 43 char equed[1]; 44 PetscInt lwork; 45 void *work; 46 PetscReal rpg, rcond; 47 mem_usage_t mem_usage; 48 MatStructure flg; 49 SuperLUStat_t stat; 50 Mat A_dup; 51 PetscScalar *rhs_dup; 52 GlobalLU_t Glu; 53 PetscBool needconversion; 54 55 /* Flag to clean up (non-global) SuperLU objects during Destroy */ 56 PetscBool CleanUpSuperLU; 57 } Mat_SuperLU; 58 59 /* 60 Utility function 61 */ 62 static PetscErrorCode MatView_Info_SuperLU(Mat A,PetscViewer viewer) 63 { 64 Mat_SuperLU *lu= (Mat_SuperLU*)A->data; 65 superlu_options_t options; 66 67 PetscFunctionBegin; 68 options = lu->options; 69 70 PetscCall(PetscViewerASCIIPrintf(viewer,"SuperLU run parameters:\n")); 71 PetscCall(PetscViewerASCIIPrintf(viewer," Equil: %s\n",(options.Equil != NO) ? "YES" : "NO")); 72 PetscCall(PetscViewerASCIIPrintf(viewer," ColPerm: %" PetscInt_FMT "\n",options.ColPerm)); 73 PetscCall(PetscViewerASCIIPrintf(viewer," IterRefine: %" PetscInt_FMT "\n",options.IterRefine)); 74 PetscCall(PetscViewerASCIIPrintf(viewer," SymmetricMode: %s\n",(options.SymmetricMode != NO) ? "YES" : "NO")); 75 PetscCall(PetscViewerASCIIPrintf(viewer," DiagPivotThresh: %g\n",options.DiagPivotThresh)); 76 PetscCall(PetscViewerASCIIPrintf(viewer," PivotGrowth: %s\n",(options.PivotGrowth != NO) ? "YES" : "NO")); 77 PetscCall(PetscViewerASCIIPrintf(viewer," ConditionNumber: %s\n",(options.ConditionNumber != NO) ? "YES" : "NO")); 78 PetscCall(PetscViewerASCIIPrintf(viewer," RowPerm: %" PetscInt_FMT "\n",options.RowPerm)); 79 PetscCall(PetscViewerASCIIPrintf(viewer," ReplaceTinyPivot: %s\n",(options.ReplaceTinyPivot != NO) ? "YES" : "NO")); 80 PetscCall(PetscViewerASCIIPrintf(viewer," PrintStat: %s\n",(options.PrintStat != NO) ? "YES" : "NO")); 81 PetscCall(PetscViewerASCIIPrintf(viewer," lwork: %" PetscInt_FMT "\n",lu->lwork)); 82 if (A->factortype == MAT_FACTOR_ILU) { 83 PetscCall(PetscViewerASCIIPrintf(viewer," ILU_DropTol: %g\n",options.ILU_DropTol)); 84 PetscCall(PetscViewerASCIIPrintf(viewer," ILU_FillTol: %g\n",options.ILU_FillTol)); 85 PetscCall(PetscViewerASCIIPrintf(viewer," ILU_FillFactor: %g\n",options.ILU_FillFactor)); 86 PetscCall(PetscViewerASCIIPrintf(viewer," ILU_DropRule: %" PetscInt_FMT "\n",options.ILU_DropRule)); 87 PetscCall(PetscViewerASCIIPrintf(viewer," ILU_Norm: %" PetscInt_FMT "\n",options.ILU_Norm)); 88 PetscCall(PetscViewerASCIIPrintf(viewer," ILU_MILU: %" PetscInt_FMT "\n",options.ILU_MILU)); 89 } 90 PetscFunctionReturn(0); 91 } 92 93 PetscErrorCode MatSolve_SuperLU_Private(Mat A,Vec b,Vec x) 94 { 95 Mat_SuperLU *lu = (Mat_SuperLU*)A->data; 96 const PetscScalar *barray; 97 PetscScalar *xarray; 98 PetscInt info,i,n; 99 PetscReal ferr,berr; 100 static PetscBool cite = PETSC_FALSE; 101 102 PetscFunctionBegin; 103 if (lu->lwork == -1) PetscFunctionReturn(0); 104 PetscCall(PetscCitationsRegister("@article{superlu99,\n author = {James W. Demmel and Stanley C. Eisenstat and\n John R. Gilbert and Xiaoye S. Li and Joseph W. H. Liu},\n title = {A supernodal approach to sparse partial pivoting},\n journal = {SIAM J. Matrix Analysis and Applications},\n year = {1999},\n volume = {20},\n number = {3},\n pages = {720-755}\n}\n",&cite)); 105 106 PetscCall(VecGetLocalSize(x,&n)); 107 lu->B.ncol = 1; /* Set the number of right-hand side */ 108 if (lu->options.Equil && !lu->rhs_dup) { 109 /* superlu overwrites b when Equil is used, thus create rhs_dup to keep user's b unchanged */ 110 PetscCall(PetscMalloc1(n,&lu->rhs_dup)); 111 } 112 if (lu->options.Equil) { 113 /* Copy b into rsh_dup */ 114 PetscCall(VecGetArrayRead(b,&barray)); 115 PetscCall(PetscArraycpy(lu->rhs_dup,barray,n)); 116 PetscCall(VecRestoreArrayRead(b,&barray)); 117 barray = lu->rhs_dup; 118 } else { 119 PetscCall(VecGetArrayRead(b,&barray)); 120 } 121 PetscCall(VecGetArray(x,&xarray)); 122 123 #if defined(PETSC_USE_COMPLEX) 124 #if defined(PETSC_USE_REAL_SINGLE) 125 ((DNformat*)lu->B.Store)->nzval = (singlecomplex*)barray; 126 ((DNformat*)lu->X.Store)->nzval = (singlecomplex*)xarray; 127 #else 128 ((DNformat*)lu->B.Store)->nzval = (doublecomplex*)barray; 129 ((DNformat*)lu->X.Store)->nzval = (doublecomplex*)xarray; 130 #endif 131 #else 132 ((DNformat*)lu->B.Store)->nzval = (void*)barray; 133 ((DNformat*)lu->X.Store)->nzval = xarray; 134 #endif 135 136 lu->options.Fact = FACTORED; /* Indicate the factored form of A is supplied. */ 137 if (A->factortype == MAT_FACTOR_LU) { 138 #if defined(PETSC_USE_COMPLEX) 139 #if defined(PETSC_USE_REAL_SINGLE) 140 PetscStackCallExternalVoid("SuperLU:cgssvx",cgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, 141 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, 142 &lu->Glu, &lu->mem_usage, &lu->stat, &info)); 143 #else 144 PetscStackCallExternalVoid("SuperLU:zgssvx",zgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, 145 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, 146 &lu->Glu, &lu->mem_usage, &lu->stat, &info)); 147 #endif 148 #else 149 #if defined(PETSC_USE_REAL_SINGLE) 150 PetscStackCallExternalVoid("SuperLU:sgssvx",sgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, 151 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, 152 &lu->Glu, &lu->mem_usage, &lu->stat, &info)); 153 #else 154 PetscStackCallExternalVoid("SuperLU:dgssvx",dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, 155 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, 156 &lu->Glu,&lu->mem_usage, &lu->stat, &info)); 157 #endif 158 #endif 159 } else if (A->factortype == MAT_FACTOR_ILU) { 160 #if defined(PETSC_USE_COMPLEX) 161 #if defined(PETSC_USE_REAL_SINGLE) 162 PetscStackCallExternalVoid("SuperLU:cgsisx",cgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, 163 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, 164 &lu->Glu, &lu->mem_usage, &lu->stat, &info)); 165 #else 166 PetscStackCallExternalVoid("SuperLU:zgsisx",zgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, 167 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, 168 &lu->Glu, &lu->mem_usage, &lu->stat, &info)); 169 #endif 170 #else 171 #if defined(PETSC_USE_REAL_SINGLE) 172 PetscStackCallExternalVoid("SuperLU:sgsisx",sgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, 173 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, 174 &lu->Glu, &lu->mem_usage, &lu->stat, &info)); 175 #else 176 PetscStackCallExternalVoid("SuperLU:dgsisx",dgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, 177 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, 178 &lu->Glu, &lu->mem_usage, &lu->stat, &info)); 179 #endif 180 #endif 181 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 182 if (!lu->options.Equil) { 183 PetscCall(VecRestoreArrayRead(b,&barray)); 184 } 185 PetscCall(VecRestoreArray(x,&xarray)); 186 187 if (!info || info == lu->A.ncol+1) { 188 if (lu->options.IterRefine) { 189 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Iterative Refinement:\n")); 190 PetscCall(PetscPrintf(PETSC_COMM_SELF," %8s%8s%16s%16s\n", "rhs", "Steps", "FERR", "BERR")); 191 for (i = 0; i < 1; ++i) PetscCall(PetscPrintf(PETSC_COMM_SELF," %8d%8d%16e%16e\n", i+1, lu->stat.RefineSteps, ferr, berr)); 192 } 193 } else if (info > 0) { 194 if (lu->lwork == -1) { 195 PetscCall(PetscPrintf(PETSC_COMM_SELF," ** Estimated memory: %" PetscInt_FMT " bytes\n", info - lu->A.ncol)); 196 } else { 197 PetscCall(PetscPrintf(PETSC_COMM_SELF," Warning: gssvx() returns info %" PetscInt_FMT "\n",info)); 198 } 199 } else PetscCheck(info >= 0,PETSC_COMM_SELF,PETSC_ERR_LIB, "info = %" PetscInt_FMT ", the %" PetscInt_FMT "-th argument in gssvx() had an illegal value", info,-info); 200 201 if (lu->options.PrintStat) { 202 PetscCall(PetscPrintf(PETSC_COMM_SELF,"MatSolve__SuperLU():\n")); 203 PetscStackCallExternalVoid("SuperLU:StatPrint",StatPrint(&lu->stat)); 204 } 205 PetscFunctionReturn(0); 206 } 207 208 PetscErrorCode MatSolve_SuperLU(Mat A,Vec b,Vec x) 209 { 210 Mat_SuperLU *lu = (Mat_SuperLU*)A->data; 211 212 PetscFunctionBegin; 213 if (A->factorerrortype) { 214 PetscCall(PetscInfo(A,"MatSolve is called with singular matrix factor, skip\n")); 215 PetscCall(VecSetInf(x)); 216 PetscFunctionReturn(0); 217 } 218 219 lu->options.Trans = TRANS; 220 PetscCall(MatSolve_SuperLU_Private(A,b,x)); 221 PetscFunctionReturn(0); 222 } 223 224 PetscErrorCode MatSolveTranspose_SuperLU(Mat A,Vec b,Vec x) 225 { 226 Mat_SuperLU *lu = (Mat_SuperLU*)A->data; 227 228 PetscFunctionBegin; 229 if (A->factorerrortype) { 230 PetscCall(PetscInfo(A,"MatSolve is called with singular matrix factor, skip\n")); 231 PetscCall(VecSetInf(x)); 232 PetscFunctionReturn(0); 233 } 234 235 lu->options.Trans = NOTRANS; 236 PetscCall(MatSolve_SuperLU_Private(A,b,x)); 237 PetscFunctionReturn(0); 238 } 239 240 static PetscErrorCode MatLUFactorNumeric_SuperLU(Mat F,Mat A,const MatFactorInfo *info) 241 { 242 Mat_SuperLU *lu = (Mat_SuperLU*)F->data; 243 Mat_SeqAIJ *aa; 244 PetscInt sinfo; 245 PetscReal ferr, berr; 246 NCformat *Ustore; 247 SCformat *Lstore; 248 249 PetscFunctionBegin; 250 if (lu->flg == SAME_NONZERO_PATTERN) { /* successing numerical factorization */ 251 lu->options.Fact = SamePattern; 252 /* Ref: ~SuperLU_3.0/EXAMPLE/dlinsolx2.c */ 253 Destroy_SuperMatrix_Store(&lu->A); 254 if (lu->A_dup) PetscCall(MatCopy_SeqAIJ(A,lu->A_dup,SAME_NONZERO_PATTERN)); 255 256 if (lu->lwork >= 0) { 257 PetscStackCallExternalVoid("SuperLU:Destroy_SuperNode_Matrix",Destroy_SuperNode_Matrix(&lu->L)); 258 PetscStackCallExternalVoid("SuperLU:Destroy_CompCol_Matrix",Destroy_CompCol_Matrix(&lu->U)); 259 lu->options.Fact = SamePattern; 260 } 261 } 262 263 /* Create the SuperMatrix for lu->A=A^T: 264 Since SuperLU likes column-oriented matrices,we pass it the transpose, 265 and then solve A^T X = B in MatSolve(). */ 266 if (lu->A_dup) { 267 aa = (Mat_SeqAIJ*)(lu->A_dup)->data; 268 } else { 269 aa = (Mat_SeqAIJ*)(A)->data; 270 } 271 #if defined(PETSC_USE_COMPLEX) 272 #if defined(PETSC_USE_REAL_SINGLE) 273 PetscStackCallExternalVoid("SuperLU:cCreate_CompCol_Matrix",cCreate_CompCol_Matrix(&lu->A,A->cmap->n,A->rmap->n,aa->nz,(singlecomplex*)aa->a,aa->j,aa->i,SLU_NC,SLU_C,SLU_GE)); 274 #else 275 PetscStackCallExternalVoid("SuperLU:zCreate_CompCol_Matrix",zCreate_CompCol_Matrix(&lu->A,A->cmap->n,A->rmap->n,aa->nz,(doublecomplex*)aa->a,aa->j,aa->i,SLU_NC,SLU_Z,SLU_GE)); 276 #endif 277 #else 278 #if defined(PETSC_USE_REAL_SINGLE) 279 PetscStackCallExternalVoid("SuperLU:sCreate_CompCol_Matrix",sCreate_CompCol_Matrix(&lu->A,A->cmap->n,A->rmap->n,aa->nz,aa->a,aa->j,aa->i,SLU_NC,SLU_S,SLU_GE)); 280 #else 281 PetscStackCallExternalVoid("SuperLU:dCreate_CompCol_Matrix",dCreate_CompCol_Matrix(&lu->A,A->cmap->n,A->rmap->n,aa->nz,aa->a,aa->j,aa->i,SLU_NC,SLU_D,SLU_GE)); 282 #endif 283 #endif 284 285 /* Numerical factorization */ 286 lu->B.ncol = 0; /* Indicate not to solve the system */ 287 if (F->factortype == MAT_FACTOR_LU) { 288 #if defined(PETSC_USE_COMPLEX) 289 #if defined(PETSC_USE_REAL_SINGLE) 290 PetscStackCallExternalVoid("SuperLU:cgssvx",cgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, 291 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, 292 &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo)); 293 #else 294 PetscStackCallExternalVoid("SuperLU:zgssvx",zgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, 295 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, 296 &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo)); 297 #endif 298 #else 299 #if defined(PETSC_USE_REAL_SINGLE) 300 PetscStackCallExternalVoid("SuperLU:sgssvx",sgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, 301 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, 302 &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo)); 303 #else 304 PetscStackCallExternalVoid("SuperLU:dgssvx",dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, 305 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, 306 &lu->Glu,&lu->mem_usage, &lu->stat, &sinfo)); 307 #endif 308 #endif 309 } else if (F->factortype == MAT_FACTOR_ILU) { 310 /* Compute the incomplete factorization, condition number and pivot growth */ 311 #if defined(PETSC_USE_COMPLEX) 312 #if defined(PETSC_USE_REAL_SINGLE) 313 PetscStackCallExternalVoid("SuperLU:cgsisx",cgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r,lu->etree, lu->equed, lu->R, lu->C, 314 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, 315 &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo)); 316 #else 317 PetscStackCallExternalVoid("SuperLU:zgsisx",zgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r,lu->etree, lu->equed, lu->R, lu->C, 318 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, 319 &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo)); 320 #endif 321 #else 322 #if defined(PETSC_USE_REAL_SINGLE) 323 PetscStackCallExternalVoid("SuperLU:sgsisx",sgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, 324 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, 325 &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo)); 326 #else 327 PetscStackCallExternalVoid("SuperLU:dgsisx",dgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, 328 &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, 329 &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo)); 330 #endif 331 #endif 332 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 333 if (!sinfo || sinfo == lu->A.ncol+1) { 334 if (lu->options.PivotGrowth) { 335 PetscCall(PetscPrintf(PETSC_COMM_SELF," Recip. pivot growth = %e\n", lu->rpg)); 336 } 337 if (lu->options.ConditionNumber) { 338 PetscCall(PetscPrintf(PETSC_COMM_SELF," Recip. condition number = %e\n", lu->rcond)); 339 } 340 } else if (sinfo > 0) { 341 if (A->erroriffailure) { 342 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot in row %" PetscInt_FMT,sinfo); 343 } else { 344 if (sinfo <= lu->A.ncol) { 345 if (lu->options.ILU_FillTol == 0.0) { 346 F->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 347 } 348 PetscCall(PetscInfo(F,"Number of zero pivots %" PetscInt_FMT ", ILU_FillTol %g\n",sinfo,lu->options.ILU_FillTol)); 349 } else if (sinfo == lu->A.ncol + 1) { 350 /* 351 U is nonsingular, but RCOND is less than machine 352 precision, meaning that the matrix is singular to 353 working precision. Nevertheless, the solution and 354 error bounds are computed because there are a number 355 of situations where the computed solution can be more 356 accurate than the value of RCOND would suggest. 357 */ 358 PetscCall(PetscInfo(F,"Matrix factor U is nonsingular, but is singular to working precision. The solution is computed. info %" PetscInt_FMT,sinfo)); 359 } else { /* sinfo > lu->A.ncol + 1 */ 360 F->factorerrortype = MAT_FACTOR_OUTMEMORY; 361 PetscCall(PetscInfo(F,"Number of bytes allocated when memory allocation fails %" PetscInt_FMT "\n",sinfo)); 362 } 363 } 364 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, "info = %" PetscInt_FMT ", the %" PetscInt_FMT "-th argument in gssvx() had an illegal value", sinfo,-sinfo); 365 366 if (lu->options.PrintStat) { 367 PetscCall(PetscPrintf(PETSC_COMM_SELF,"MatLUFactorNumeric_SuperLU():\n")); 368 PetscStackCallExternalVoid("SuperLU:StatPrint",StatPrint(&lu->stat)); 369 Lstore = (SCformat*) lu->L.Store; 370 Ustore = (NCformat*) lu->U.Store; 371 PetscCall(PetscPrintf(PETSC_COMM_SELF," No of nonzeros in factor L = %" PetscInt_FMT "\n", Lstore->nnz)); 372 PetscCall(PetscPrintf(PETSC_COMM_SELF," No of nonzeros in factor U = %" PetscInt_FMT "\n", Ustore->nnz)); 373 PetscCall(PetscPrintf(PETSC_COMM_SELF," No of nonzeros in L+U = %" PetscInt_FMT "\n", Lstore->nnz + Ustore->nnz - lu->A.ncol)); 374 PetscCall(PetscPrintf(PETSC_COMM_SELF," L\\U MB %.3f\ttotal MB needed %.3f\n",lu->mem_usage.for_lu/1e6, lu->mem_usage.total_needed/1e6)); 375 } 376 377 lu->flg = SAME_NONZERO_PATTERN; 378 F->ops->solve = MatSolve_SuperLU; 379 F->ops->solvetranspose = MatSolveTranspose_SuperLU; 380 F->ops->matsolve = NULL; 381 PetscFunctionReturn(0); 382 } 383 384 static PetscErrorCode MatDestroy_SuperLU(Mat A) 385 { 386 Mat_SuperLU *lu=(Mat_SuperLU*)A->data; 387 388 PetscFunctionBegin; 389 if (lu->CleanUpSuperLU) { /* Free the SuperLU datastructures */ 390 PetscStackCallExternalVoid("SuperLU:Destroy_SuperMatrix_Store",Destroy_SuperMatrix_Store(&lu->A)); 391 PetscStackCallExternalVoid("SuperLU:Destroy_SuperMatrix_Store",Destroy_SuperMatrix_Store(&lu->B)); 392 PetscStackCallExternalVoid("SuperLU:Destroy_SuperMatrix_Store",Destroy_SuperMatrix_Store(&lu->X)); 393 PetscStackCallExternalVoid("SuperLU:StatFree",StatFree(&lu->stat)); 394 if (lu->lwork >= 0) { 395 PetscStackCallExternalVoid("SuperLU:Destroy_SuperNode_Matrix",Destroy_SuperNode_Matrix(&lu->L)); 396 PetscStackCallExternalVoid("SuperLU:Destroy_CompCol_Matrix",Destroy_CompCol_Matrix(&lu->U)); 397 } 398 } 399 PetscCall(PetscFree(lu->etree)); 400 PetscCall(PetscFree(lu->perm_r)); 401 PetscCall(PetscFree(lu->perm_c)); 402 PetscCall(PetscFree(lu->R)); 403 PetscCall(PetscFree(lu->C)); 404 PetscCall(PetscFree(lu->rhs_dup)); 405 PetscCall(MatDestroy(&lu->A_dup)); 406 PetscCall(PetscFree(A->data)); 407 408 /* clear composed functions */ 409 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL)); 410 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSuperluSetILUDropTol_C",NULL)); 411 PetscFunctionReturn(0); 412 } 413 414 static PetscErrorCode MatView_SuperLU(Mat A,PetscViewer viewer) 415 { 416 PetscBool iascii; 417 PetscViewerFormat format; 418 419 PetscFunctionBegin; 420 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 421 if (iascii) { 422 PetscCall(PetscViewerGetFormat(viewer,&format)); 423 if (format == PETSC_VIEWER_ASCII_INFO) { 424 PetscCall(MatView_Info_SuperLU(A,viewer)); 425 } 426 } 427 PetscFunctionReturn(0); 428 } 429 430 PetscErrorCode MatMatSolve_SuperLU(Mat A,Mat B,Mat X) 431 { 432 Mat_SuperLU *lu = (Mat_SuperLU*)A->data; 433 PetscBool flg; 434 435 PetscFunctionBegin; 436 PetscCall(PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL)); 437 PetscCheck(flg,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix"); 438 PetscCall(PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL)); 439 PetscCheck(flg,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix"); 440 lu->options.Trans = TRANS; 441 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_SuperLU() is not implemented yet"); 442 PetscFunctionReturn(0); 443 } 444 445 static PetscErrorCode MatLUFactorSymbolic_SuperLU(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 446 { 447 Mat_SuperLU *lu = (Mat_SuperLU*)(F->data); 448 PetscInt indx; 449 PetscBool flg,set; 450 PetscReal real_input; 451 const char *colperm[] ={"NATURAL","MMD_ATA","MMD_AT_PLUS_A","COLAMD"}; /* MY_PERMC - not supported by the petsc interface yet */ 452 const char *iterrefine[]={"NOREFINE", "SINGLE", "DOUBLE", "EXTRA"}; 453 const char *rowperm[] ={"NOROWPERM", "LargeDiag"}; /* MY_PERMC - not supported by the petsc interface yet */ 454 455 PetscFunctionBegin; 456 /* Set options to F */ 457 PetscOptionsBegin(PetscObjectComm((PetscObject)F),((PetscObject)F)->prefix,"SuperLU Options","Mat"); 458 PetscCall(PetscOptionsBool("-mat_superlu_equil","Equil","None",(PetscBool)lu->options.Equil,(PetscBool*)&lu->options.Equil,NULL)); 459 PetscCall(PetscOptionsEList("-mat_superlu_colperm","ColPerm","None",colperm,4,colperm[3],&indx,&flg)); 460 if (flg) lu->options.ColPerm = (colperm_t)indx; 461 PetscCall(PetscOptionsEList("-mat_superlu_iterrefine","IterRefine","None",iterrefine,4,iterrefine[0],&indx,&flg)); 462 if (flg) lu->options.IterRefine = (IterRefine_t)indx; 463 PetscCall(PetscOptionsBool("-mat_superlu_symmetricmode","SymmetricMode","None",(PetscBool)lu->options.SymmetricMode,&flg,&set)); 464 if (set && flg) lu->options.SymmetricMode = YES; 465 PetscCall(PetscOptionsReal("-mat_superlu_diagpivotthresh","DiagPivotThresh","None",lu->options.DiagPivotThresh,&real_input,&flg)); 466 if (flg) lu->options.DiagPivotThresh = (double) real_input; 467 PetscCall(PetscOptionsBool("-mat_superlu_pivotgrowth","PivotGrowth","None",(PetscBool)lu->options.PivotGrowth,&flg,&set)); 468 if (set && flg) lu->options.PivotGrowth = YES; 469 PetscCall(PetscOptionsBool("-mat_superlu_conditionnumber","ConditionNumber","None",(PetscBool)lu->options.ConditionNumber,&flg,&set)); 470 if (set && flg) lu->options.ConditionNumber = YES; 471 PetscCall(PetscOptionsEList("-mat_superlu_rowperm","rowperm","None",rowperm,2,rowperm[lu->options.RowPerm],&indx,&flg)); 472 if (flg) lu->options.RowPerm = (rowperm_t)indx; 473 PetscCall(PetscOptionsBool("-mat_superlu_replacetinypivot","ReplaceTinyPivot","None",(PetscBool)lu->options.ReplaceTinyPivot,&flg,&set)); 474 if (set && flg) lu->options.ReplaceTinyPivot = YES; 475 PetscCall(PetscOptionsBool("-mat_superlu_printstat","PrintStat","None",(PetscBool)lu->options.PrintStat,&flg,&set)); 476 if (set && flg) lu->options.PrintStat = YES; 477 PetscCall(PetscOptionsInt("-mat_superlu_lwork","size of work array in bytes used by factorization","None",lu->lwork,&lu->lwork,NULL)); 478 if (lu->lwork > 0) { 479 /* lwork is in bytes, hence PetscMalloc() is used here, not PetscMalloc1()*/ 480 PetscCall(PetscMalloc(lu->lwork,&lu->work)); 481 } else if (lu->lwork != 0 && lu->lwork != -1) { 482 PetscCall(PetscPrintf(PETSC_COMM_SELF," Warning: lwork %" PetscInt_FMT " is not supported by SUPERLU. The default lwork=0 is used.\n",lu->lwork)); 483 lu->lwork = 0; 484 } 485 /* ilu options */ 486 PetscCall(PetscOptionsReal("-mat_superlu_ilu_droptol","ILU_DropTol","None",lu->options.ILU_DropTol,&real_input,&flg)); 487 if (flg) lu->options.ILU_DropTol = (double) real_input; 488 PetscCall(PetscOptionsReal("-mat_superlu_ilu_filltol","ILU_FillTol","None",lu->options.ILU_FillTol,&real_input,&flg)); 489 if (flg) lu->options.ILU_FillTol = (double) real_input; 490 PetscCall(PetscOptionsReal("-mat_superlu_ilu_fillfactor","ILU_FillFactor","None",lu->options.ILU_FillFactor,&real_input,&flg)); 491 if (flg) lu->options.ILU_FillFactor = (double) real_input; 492 PetscCall(PetscOptionsInt("-mat_superlu_ilu_droprull","ILU_DropRule","None",lu->options.ILU_DropRule,&lu->options.ILU_DropRule,NULL)); 493 PetscCall(PetscOptionsInt("-mat_superlu_ilu_norm","ILU_Norm","None",lu->options.ILU_Norm,&indx,&flg)); 494 if (flg) lu->options.ILU_Norm = (norm_t)indx; 495 PetscCall(PetscOptionsInt("-mat_superlu_ilu_milu","ILU_MILU","None",lu->options.ILU_MILU,&indx,&flg)); 496 if (flg) lu->options.ILU_MILU = (milu_t)indx; 497 PetscOptionsEnd(); 498 499 lu->flg = DIFFERENT_NONZERO_PATTERN; 500 lu->CleanUpSuperLU = PETSC_TRUE; 501 F->ops->lufactornumeric = MatLUFactorNumeric_SuperLU; 502 503 /* if we are here, the nonzero pattern has changed unless the user explicitly called MatLUFactorSymbolic */ 504 PetscCall(MatDestroy(&lu->A_dup)); 505 if (lu->needconversion) { 506 PetscCall(MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&lu->A_dup)); 507 } 508 if (lu->options.Equil == YES && !lu->A_dup) { /* superlu overwrites input matrix and rhs when Equil is used, thus create A_dup to keep user's A unchanged */ 509 PetscCall(MatDuplicate_SeqAIJ(A,MAT_COPY_VALUES,&lu->A_dup)); 510 } 511 PetscFunctionReturn(0); 512 } 513 514 static PetscErrorCode MatSuperluSetILUDropTol_SuperLU(Mat F,PetscReal dtol) 515 { 516 Mat_SuperLU *lu= (Mat_SuperLU*)F->data; 517 518 PetscFunctionBegin; 519 lu->options.ILU_DropTol = dtol; 520 PetscFunctionReturn(0); 521 } 522 523 /*@ 524 MatSuperluSetILUDropTol - Set SuperLU ILU drop tolerance 525 526 Logically Collective on Mat 527 528 Input Parameters: 529 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-SuperLU interface 530 - dtol - drop tolerance 531 532 Options Database: 533 . -mat_superlu_ilu_droptol <dtol> - the drop tolerance 534 535 Level: beginner 536 537 References: 538 . * - SuperLU Users' Guide 539 540 .seealso: `MatGetFactor()` 541 @*/ 542 PetscErrorCode MatSuperluSetILUDropTol(Mat F,PetscReal dtol) 543 { 544 PetscFunctionBegin; 545 PetscValidHeaderSpecific(F,MAT_CLASSID,1); 546 PetscValidLogicalCollectiveReal(F,dtol,2); 547 PetscTryMethod(F,"MatSuperluSetILUDropTol_C",(Mat,PetscReal),(F,dtol)); 548 PetscFunctionReturn(0); 549 } 550 551 PetscErrorCode MatFactorGetSolverType_seqaij_superlu(Mat A,MatSolverType *type) 552 { 553 PetscFunctionBegin; 554 *type = MATSOLVERSUPERLU; 555 PetscFunctionReturn(0); 556 } 557 558 /*MC 559 MATSOLVERSUPERLU = "superlu" - A solver package providing solvers LU and ILU for sequential matrices 560 via the external package SuperLU. 561 562 Use ./configure --download-superlu to have PETSc installed with SuperLU 563 564 Use -pc_type lu -pc_factor_mat_solver_type superlu to use this direct solver 565 566 Options Database Keys: 567 + -mat_superlu_equil <FALSE> - Equil (None) 568 . -mat_superlu_colperm <COLAMD> - (choose one of) NATURAL MMD_ATA MMD_AT_PLUS_A COLAMD 569 . -mat_superlu_iterrefine <NOREFINE> - (choose one of) NOREFINE SINGLE DOUBLE EXTRA 570 . -mat_superlu_symmetricmode: <FALSE> - SymmetricMode (None) 571 . -mat_superlu_diagpivotthresh <1> - DiagPivotThresh (None) 572 . -mat_superlu_pivotgrowth <FALSE> - PivotGrowth (None) 573 . -mat_superlu_conditionnumber <FALSE> - ConditionNumber (None) 574 . -mat_superlu_rowperm <NOROWPERM> - (choose one of) NOROWPERM LargeDiag 575 . -mat_superlu_replacetinypivot <FALSE> - ReplaceTinyPivot (None) 576 . -mat_superlu_printstat <FALSE> - PrintStat (None) 577 . -mat_superlu_lwork <0> - size of work array in bytes used by factorization (None) 578 . -mat_superlu_ilu_droptol <0> - ILU_DropTol (None) 579 . -mat_superlu_ilu_filltol <0> - ILU_FillTol (None) 580 . -mat_superlu_ilu_fillfactor <0> - ILU_FillFactor (None) 581 . -mat_superlu_ilu_droprull <0> - ILU_DropRule (None) 582 . -mat_superlu_ilu_norm <0> - ILU_Norm (None) 583 - -mat_superlu_ilu_milu <0> - ILU_MILU (None) 584 585 Notes: 586 Do not confuse this with MATSOLVERSUPERLU_DIST which is for parallel sparse solves 587 588 Cannot currently use ordering provided by PETSc. 589 590 Level: beginner 591 592 .seealso: `PCLU`, `PCILU`, `MATSOLVERSUPERLU_DIST`, `MATSOLVERMUMPS`, `PCFactorSetMatSolverType()`, `MatSolverType` 593 M*/ 594 595 static PetscErrorCode MatGetFactor_seqaij_superlu(Mat A,MatFactorType ftype,Mat *F) 596 { 597 Mat B; 598 Mat_SuperLU *lu; 599 PetscInt m=A->rmap->n,n=A->cmap->n; 600 601 PetscFunctionBegin; 602 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 603 PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,PETSC_DETERMINE,PETSC_DETERMINE)); 604 PetscCall(PetscStrallocpy("superlu",&((PetscObject)B)->type_name)); 605 PetscCall(MatSetUp(B)); 606 B->trivialsymbolic = PETSC_TRUE; 607 if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) { 608 B->ops->lufactorsymbolic = MatLUFactorSymbolic_SuperLU; 609 B->ops->ilufactorsymbolic = MatLUFactorSymbolic_SuperLU; 610 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 611 612 PetscCall(PetscFree(B->solvertype)); 613 PetscCall(PetscStrallocpy(MATSOLVERSUPERLU,&B->solvertype)); 614 615 B->ops->getinfo = MatGetInfo_External; 616 B->ops->destroy = MatDestroy_SuperLU; 617 B->ops->view = MatView_SuperLU; 618 B->factortype = ftype; 619 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 620 B->preallocated = PETSC_TRUE; 621 622 PetscCall(PetscNewLog(B,&lu)); 623 624 if (ftype == MAT_FACTOR_LU) { 625 set_default_options(&lu->options); 626 /* Comments from SuperLU_4.0/SRC/dgssvx.c: 627 "Whether or not the system will be equilibrated depends on the 628 scaling of the matrix A, but if equilibration is used, A is 629 overwritten by diag(R)*A*diag(C) and B by diag(R)*B 630 (if options->Trans=NOTRANS) or diag(C)*B (if options->Trans = TRANS or CONJ)." 631 We set 'options.Equil = NO' as default because additional space is needed for it. 632 */ 633 lu->options.Equil = NO; 634 } else if (ftype == MAT_FACTOR_ILU) { 635 /* Set the default input options of ilu: */ 636 PetscStackCallExternalVoid("SuperLU:ilu_set_default_options",ilu_set_default_options(&lu->options)); 637 } 638 lu->options.PrintStat = NO; 639 640 /* Initialize the statistics variables. */ 641 PetscStackCallExternalVoid("SuperLU:StatInit",StatInit(&lu->stat)); 642 lu->lwork = 0; /* allocate space internally by system malloc */ 643 644 /* Allocate spaces (notice sizes are for the transpose) */ 645 PetscCall(PetscMalloc1(m,&lu->etree)); 646 PetscCall(PetscMalloc1(n,&lu->perm_r)); 647 PetscCall(PetscMalloc1(m,&lu->perm_c)); 648 PetscCall(PetscMalloc1(n,&lu->R)); 649 PetscCall(PetscMalloc1(m,&lu->C)); 650 651 /* create rhs and solution x without allocate space for .Store */ 652 #if defined(PETSC_USE_COMPLEX) 653 #if defined(PETSC_USE_REAL_SINGLE) 654 PetscStackCallExternalVoid("SuperLU:cCreate_Dense_Matrix(",cCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_C, SLU_GE)); 655 PetscStackCallExternalVoid("SuperLU:cCreate_Dense_Matrix(",cCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_C, SLU_GE)); 656 #else 657 PetscStackCallExternalVoid("SuperLU:zCreate_Dense_Matrix",zCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_Z, SLU_GE)); 658 PetscStackCallExternalVoid("SuperLU:zCreate_Dense_Matrix",zCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_Z, SLU_GE)); 659 #endif 660 #else 661 #if defined(PETSC_USE_REAL_SINGLE) 662 PetscStackCallExternalVoid("SuperLU:sCreate_Dense_Matrix",sCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_S, SLU_GE)); 663 PetscStackCallExternalVoid("SuperLU:sCreate_Dense_Matrix",sCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_S, SLU_GE)); 664 #else 665 PetscStackCallExternalVoid("SuperLU:dCreate_Dense_Matrix",dCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_D, SLU_GE)); 666 PetscStackCallExternalVoid("SuperLU:dCreate_Dense_Matrix",dCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_D, SLU_GE)); 667 #endif 668 #endif 669 670 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_superlu)); 671 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSuperluSetILUDropTol_C",MatSuperluSetILUDropTol_SuperLU)); 672 B->data = lu; 673 674 *F = B; 675 PetscFunctionReturn(0); 676 } 677 678 static PetscErrorCode MatGetFactor_seqsell_superlu(Mat A,MatFactorType ftype,Mat *F) 679 { 680 Mat_SuperLU *lu; 681 682 PetscFunctionBegin; 683 PetscCall(MatGetFactor_seqaij_superlu(A,ftype,F)); 684 lu = (Mat_SuperLU*)((*F)->data); 685 lu->needconversion = PETSC_TRUE; 686 PetscFunctionReturn(0); 687 } 688 689 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_SuperLU(void) 690 { 691 PetscFunctionBegin; 692 PetscCall(MatSolverTypeRegister(MATSOLVERSUPERLU,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_seqaij_superlu)); 693 PetscCall(MatSolverTypeRegister(MATSOLVERSUPERLU,MATSEQAIJ,MAT_FACTOR_ILU,MatGetFactor_seqaij_superlu)); 694 PetscCall(MatSolverTypeRegister(MATSOLVERSUPERLU,MATSEQSELL,MAT_FACTOR_LU,MatGetFactor_seqsell_superlu)); 695 PetscCall(MatSolverTypeRegister(MATSOLVERSUPERLU,MATSEQSELL,MAT_FACTOR_ILU,MatGetFactor_seqsell_superlu)); 696 PetscFunctionReturn(0); 697 } 698