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