1 /* 2 Provides an interface to the LUSOL package of .... 3 4 */ 5 #include <../src/mat/impls/aij/seq/aij.h> 6 7 #if defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 8 #define LU1FAC lu1fac_ 9 #define LU6SOL lu6sol_ 10 #define M1PAGE m1page_ 11 #define M5SETX m5setx_ 12 #define M6RDEL m6rdel_ 13 #elif !defined(PETSC_HAVE_FORTRAN_CAPS) 14 #define LU1FAC lu1fac 15 #define LU6SOL lu6sol 16 #define M1PAGE m1page 17 #define M5SETX m5setx 18 #define M6RDEL m6rdel 19 #endif 20 21 /* 22 Dummy symbols that the MINOS files mi25bfac.f and mi15blas.f may require 23 */ 24 PETSC_EXTERN void M1PAGE() 25 { 26 ; 27 } 28 PETSC_EXTERN void M5SETX() 29 { 30 ; 31 } 32 33 PETSC_EXTERN void M6RDEL() 34 { 35 ; 36 } 37 38 PETSC_EXTERN void LU1FAC(int *m, int *n, int *nnz, int *size, int *luparm, double *parmlu, double *data, int *indc, int *indr, int *rowperm, int *colperm, int *collen, int *rowlen, int *colstart, int *rowstart, int *rploc, int *cploc, int *rpinv, int *cpinv, double *w, int *inform); 39 40 PETSC_EXTERN void LU6SOL(int *mode, int *m, int *n, double *rhs, double *x, int *size, int *luparm, double *parmlu, double *data, int *indc, int *indr, int *rowperm, int *colperm, int *collen, int *rowlen, int *colstart, int *rowstart, int *inform); 41 42 extern PetscErrorCode MatDuplicate_LUSOL(Mat, MatDuplicateOption, Mat *); 43 44 typedef struct { 45 double *data; 46 int *indc; 47 int *indr; 48 49 int *ip; 50 int *iq; 51 int *lenc; 52 int *lenr; 53 int *locc; 54 int *locr; 55 int *iploc; 56 int *iqloc; 57 int *ipinv; 58 int *iqinv; 59 double *mnsw; 60 double *mnsv; 61 62 double elbowroom; 63 double luroom; /* Extra space allocated when factor fails */ 64 double parmlu[30]; /* Input/output to LUSOL */ 65 66 int n; /* Number of rows/columns in matrix */ 67 int nz; /* Number of nonzeros */ 68 int nnz; /* Number of nonzeros allocated for factors */ 69 int luparm[30]; /* Input/output to LUSOL */ 70 71 PetscBool CleanUpLUSOL; 72 73 } Mat_LUSOL; 74 75 /* 76 LUSOL input/Output Parameters (Description uses C-style indexes 77 78 Input parameters Typical value 79 luparm(0) = nout File number for printed messages. 6 80 luparm(1) = lprint Print level. 0 81 < 0 suppresses output. 82 = 0 gives error messages. 83 = 1 gives debug output from some of the 84 other routines in LUSOL. 85 >= 2 gives the pivot row and column and the 86 no. of rows and columns involved at 87 each elimination step in lu1fac. 88 luparm(2) = maxcol lu1fac: maximum number of columns 5 89 searched allowed in a Markowitz-type 90 search for the next pivot element. 91 For some of the factorization, the 92 number of rows searched is 93 maxrow = maxcol - 1. 94 95 Output parameters: 96 97 luparm(9) = inform Return code from last call to any LU routine. 98 luparm(10) = nsing No. of singularities marked in the 99 output array w(*). 100 luparm(11) = jsing Column index of last singularity. 101 luparm(12) = minlen Minimum recommended value for lena. 102 luparm(13) = maxlen ? 103 luparm(14) = nupdat No. of updates performed by the lu8 routines. 104 luparm(15) = nrank No. of nonempty rows of U. 105 luparm(16) = ndens1 No. of columns remaining when the density of 106 the matrix being factorized reached dens1. 107 luparm(17) = ndens2 No. of columns remaining when the density of 108 the matrix being factorized reached dens2. 109 luparm(18) = jumin The column index associated with dumin. 110 luparm(19) = numl0 No. of columns in initial L. 111 luparm(20) = lenl0 Size of initial L (no. of nonzeros). 112 luparm(21) = lenu0 Size of initial U. 113 luparm(22) = lenl Size of current L. 114 luparm(23) = lenu Size of current U. 115 luparm(24) = lrow Length of row file. 116 luparm(25) = ncp No. of compressions of LU data structures. 117 luparm(26) = mersum lu1fac: sum of Markowitz merit counts. 118 luparm(27) = nutri lu1fac: triangular rows in U. 119 luparm(28) = nltri lu1fac: triangular rows in L. 120 luparm(29) = 121 122 Input parameters Typical value 123 parmlu(0) = elmax1 Max multiplier allowed in L 10.0 124 during factor. 125 parmlu(1) = elmax2 Max multiplier allowed in L 10.0 126 during updates. 127 parmlu(2) = small Absolute tolerance for eps**0.8 128 treating reals as zero. IBM double: 3.0d-13 129 parmlu(3) = utol1 Absolute tol for flagging eps**0.66667 130 small diagonals of U. IBM double: 3.7d-11 131 parmlu(4) = utol2 Relative tol for flagging eps**0.66667 132 small diagonals of U. IBM double: 3.7d-11 133 parmlu(5) = uspace Factor limiting waste space in U. 3.0 134 In lu1fac, the row or column lists 135 are compressed if their length 136 exceeds uspace times the length of 137 either file after the last compression. 138 parmlu(6) = dens1 The density at which the Markowitz 0.3 139 strategy should search maxcol columns 140 and no rows. 141 parmlu(7) = dens2 the density at which the Markowitz 0.6 142 strategy should search only 1 column 143 or (preferably) use a dense LU for 144 all the remaining rows and columns. 145 146 Output parameters: 147 parmlu(9) = amax Maximum element in A. 148 parmlu(10) = elmax Maximum multiplier in current L. 149 parmlu(11) = umax Maximum element in current U. 150 parmlu(12) = dumax Maximum diagonal in U. 151 parmlu(13) = dumin Minimum diagonal in U. 152 parmlu(14) = 153 parmlu(15) = 154 parmlu(16) = 155 parmlu(17) = 156 parmlu(18) = 157 parmlu(19) = resid lu6sol: residual after solve with U or U'. 158 ... 159 parmlu(29) = 160 */ 161 162 #define Factorization_Tolerance 1e-1 163 #define Factorization_Pivot_Tolerance pow(2.2204460492503131E-16, 2.0 / 3.0) 164 #define Factorization_Small_Tolerance 1e-15 /* pow(DBL_EPSILON, 0.8) */ 165 166 static PetscErrorCode MatDestroy_LUSOL(Mat A) 167 { 168 Mat_LUSOL *lusol = (Mat_LUSOL *)A->spptr; 169 170 PetscFunctionBegin; 171 if (lusol && lusol->CleanUpLUSOL) { 172 PetscCall(PetscFree(lusol->ip)); 173 PetscCall(PetscFree(lusol->iq)); 174 PetscCall(PetscFree(lusol->lenc)); 175 PetscCall(PetscFree(lusol->lenr)); 176 PetscCall(PetscFree(lusol->locc)); 177 PetscCall(PetscFree(lusol->locr)); 178 PetscCall(PetscFree(lusol->iploc)); 179 PetscCall(PetscFree(lusol->iqloc)); 180 PetscCall(PetscFree(lusol->ipinv)); 181 PetscCall(PetscFree(lusol->iqinv)); 182 PetscCall(PetscFree(lusol->mnsw)); 183 PetscCall(PetscFree(lusol->mnsv)); 184 PetscCall(PetscFree3(lusol->data, lusol->indc, lusol->indr)); 185 } 186 PetscCall(PetscFree(A->spptr)); 187 PetscCall(MatDestroy_SeqAIJ(A)); 188 PetscFunctionReturn(PETSC_SUCCESS); 189 } 190 191 static PetscErrorCode MatSolve_LUSOL(Mat A, Vec b, Vec x) 192 { 193 Mat_LUSOL *lusol = (Mat_LUSOL *)A->spptr; 194 double *xx; 195 const double *bb; 196 int mode = 5; 197 int i, m, n, nnz, status; 198 199 PetscFunctionBegin; 200 PetscCall(VecGetArray(x, &xx)); 201 PetscCall(VecGetArrayRead(b, &bb)); 202 203 m = n = lusol->n; 204 nnz = lusol->nnz; 205 206 for (i = 0; i < m; i++) lusol->mnsv[i] = bb[i]; 207 208 LU6SOL(&mode, &m, &n, lusol->mnsv, xx, &nnz, lusol->luparm, lusol->parmlu, lusol->data, lusol->indc, lusol->indr, lusol->ip, lusol->iq, lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, &status); 209 210 PetscCheck(!status, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "solve failed, error code %d", status); 211 212 PetscCall(VecRestoreArray(x, &xx)); 213 PetscCall(VecRestoreArrayRead(b, &bb)); 214 PetscFunctionReturn(PETSC_SUCCESS); 215 } 216 217 static PetscErrorCode MatLUFactorNumeric_LUSOL(Mat F, Mat A, const MatFactorInfo *info) 218 { 219 Mat_SeqAIJ *a; 220 Mat_LUSOL *lusol = (Mat_LUSOL *)F->spptr; 221 int m, n, nz, nnz, status; 222 int i, rs, re; 223 int factorizations; 224 225 PetscFunctionBegin; 226 PetscCall(MatGetSize(A, &m, &n)); 227 a = (Mat_SeqAIJ *)A->data; 228 229 PetscCheck(m == lusol->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "factorization struct inconsistent"); 230 231 factorizations = 0; 232 do { 233 /*******************************************************************/ 234 /* Check the workspace allocation. */ 235 /*******************************************************************/ 236 237 nz = a->nz; 238 nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom * nz)); 239 nnz = PetscMax(nnz, 5 * n); 240 241 if (nnz < lusol->luparm[12]) { 242 nnz = (int)(lusol->luroom * lusol->luparm[12]); 243 } else if ((factorizations > 0) && (lusol->luroom < 6)) { 244 lusol->luroom += 0.1; 245 } 246 247 nnz = PetscMax(nnz, (int)(lusol->luroom * (lusol->luparm[22] + lusol->luparm[23]))); 248 249 if (nnz > lusol->nnz) { 250 PetscCall(PetscFree3(lusol->data, lusol->indc, lusol->indr)); 251 PetscCall(PetscMalloc3(nnz, &lusol->data, nnz, &lusol->indc, nnz, &lusol->indr)); 252 lusol->nnz = nnz; 253 } 254 255 /* Fill in the data for the problem. (1-based Fortran style) */ 256 nz = 0; 257 for (i = 0; i < n; i++) { 258 rs = a->i[i]; 259 re = a->i[i + 1]; 260 261 while (rs < re) { 262 if (a->a[rs] != 0.0) { 263 lusol->indc[nz] = i + 1; 264 lusol->indr[nz] = a->j[rs] + 1; 265 lusol->data[nz] = a->a[rs]; 266 nz++; 267 } 268 rs++; 269 } 270 } 271 272 /* Do the factorization. */ 273 LU1FAC(&m, &n, &nz, &nnz, lusol->luparm, lusol->parmlu, lusol->data, lusol->indc, lusol->indr, lusol->ip, lusol->iq, lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, lusol->iploc, lusol->iqloc, lusol->ipinv, lusol->iqinv, lusol->mnsw, &status); 274 275 switch (status) { 276 case 0: /* factored */ 277 break; 278 279 case 7: /* insufficient memory */ 280 break; 281 282 case 1: 283 case -1: /* singular */ 284 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Singular matrix"); 285 286 case 3: 287 case 4: /* error conditions */ 288 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "matrix error"); 289 290 default: /* unknown condition */ 291 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "matrix unknown return code"); 292 } 293 294 factorizations++; 295 } while (status == 7); 296 F->ops->solve = MatSolve_LUSOL; 297 F->assembled = PETSC_TRUE; 298 F->preallocated = PETSC_TRUE; 299 PetscFunctionReturn(PETSC_SUCCESS); 300 } 301 302 static PetscErrorCode MatLUFactorSymbolic_LUSOL(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info) 303 { 304 /* 305 Input 306 A - matrix to factor 307 r - row permutation (ignored) 308 c - column permutation (ignored) 309 310 Output 311 F - matrix storing the factorization; 312 */ 313 Mat_LUSOL *lusol; 314 int i, m, n, nz, nnz; 315 316 PetscFunctionBegin; 317 /* Check the arguments. */ 318 PetscCall(MatGetSize(A, &m, &n)); 319 nz = ((Mat_SeqAIJ *)A->data)->nz; 320 321 /* Create the factorization. */ 322 F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 323 lusol = (Mat_LUSOL *)F->spptr; 324 325 /* Initialize parameters */ 326 for (i = 0; i < 30; i++) { 327 lusol->luparm[i] = 0; 328 lusol->parmlu[i] = 0; 329 } 330 331 lusol->luparm[1] = -1; 332 lusol->luparm[2] = 5; 333 lusol->luparm[7] = 1; 334 335 lusol->parmlu[0] = 1 / Factorization_Tolerance; 336 lusol->parmlu[1] = 1 / Factorization_Tolerance; 337 lusol->parmlu[2] = Factorization_Small_Tolerance; 338 lusol->parmlu[3] = Factorization_Pivot_Tolerance; 339 lusol->parmlu[4] = Factorization_Pivot_Tolerance; 340 lusol->parmlu[5] = 3.0; 341 lusol->parmlu[6] = 0.3; 342 lusol->parmlu[7] = 0.6; 343 344 /* Allocate the workspace needed by LUSOL. */ 345 lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill); 346 nnz = PetscMax((int)(lusol->elbowroom * nz), 5 * n); 347 348 lusol->n = n; 349 lusol->nz = nz; 350 lusol->nnz = nnz; 351 lusol->luroom = 1.75; 352 353 PetscCall(PetscMalloc(sizeof(int) * n, &lusol->ip)); 354 PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iq)); 355 PetscCall(PetscMalloc(sizeof(int) * n, &lusol->lenc)); 356 PetscCall(PetscMalloc(sizeof(int) * n, &lusol->lenr)); 357 PetscCall(PetscMalloc(sizeof(int) * n, &lusol->locc)); 358 PetscCall(PetscMalloc(sizeof(int) * n, &lusol->locr)); 359 PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iploc)); 360 PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iqloc)); 361 PetscCall(PetscMalloc(sizeof(int) * n, &lusol->ipinv)); 362 PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iqinv)); 363 PetscCall(PetscMalloc(sizeof(double) * n, &lusol->mnsw)); 364 PetscCall(PetscMalloc(sizeof(double) * n, &lusol->mnsv)); 365 PetscCall(PetscMalloc3(nnz, &lusol->data, nnz, &lusol->indc, nnz, &lusol->indr)); 366 367 lusol->CleanUpLUSOL = PETSC_TRUE; 368 F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 369 PetscFunctionReturn(PETSC_SUCCESS); 370 } 371 372 static PetscErrorCode MatFactorGetSolverType_seqaij_lusol(Mat A, MatSolverType *type) 373 { 374 PetscFunctionBegin; 375 *type = MATSOLVERLUSOL; 376 PetscFunctionReturn(PETSC_SUCCESS); 377 } 378 379 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_lusol(Mat A, MatFactorType ftype, Mat *F) 380 { 381 Mat B; 382 Mat_LUSOL *lusol; 383 int m, n; 384 385 PetscFunctionBegin; 386 PetscCall(MatGetSize(A, &m, &n)); 387 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 388 PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n)); 389 PetscCall(MatSetType(B, ((PetscObject)A)->type_name)); 390 PetscCall(MatSeqAIJSetPreallocation(B, 0, NULL)); 391 392 PetscCall(PetscNew(&lusol)); 393 B->spptr = lusol; 394 395 B->trivialsymbolic = PETSC_TRUE; 396 B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL; 397 B->ops->destroy = MatDestroy_LUSOL; 398 399 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_lusol)); 400 401 B->factortype = MAT_FACTOR_LU; 402 PetscCall(PetscFree(B->solvertype)); 403 PetscCall(PetscStrallocpy(MATSOLVERLUSOL, &B->solvertype)); 404 PetscFunctionReturn(PETSC_SUCCESS); 405 } 406 407 PETSC_INTERN PetscErrorCode MatSolverTypeRegister_Lusol(void) 408 { 409 PetscFunctionBegin; 410 PetscCall(MatSolverTypeRegister(MATSOLVERLUSOL, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_lusol)); 411 PetscFunctionReturn(PETSC_SUCCESS); 412 } 413 414 /*MC 415 MATSOLVERLUSOL - "lusol" - Provides direct solvers, LU, for sequential matrices 416 via the external package LUSOL. 417 418 Works with `MATSEQAIJ` matrices 419 420 Level: beginner 421 422 .seealso: [](ch_matrices), `Mat`, `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType` 423 M*/ 424