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 If LUSOL is installed (see the manual for 444 instructions on how to declare the existence of external packages), 445 446 Works with `MATSEQAIJ` matrices 447 448 Level: beginner 449 450 .seealso: `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType` 451 M*/ 452