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