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