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