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