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