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