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 PetscErrorCode ierr; 183 Mat_LUSOL *lusol=(Mat_LUSOL*)A->spptr; 184 185 PetscFunctionBegin; 186 if (lusol && lusol->CleanUpLUSOL) { 187 ierr = PetscFree(lusol->ip);CHKERRQ(ierr); 188 ierr = PetscFree(lusol->iq);CHKERRQ(ierr); 189 ierr = PetscFree(lusol->lenc);CHKERRQ(ierr); 190 ierr = PetscFree(lusol->lenr);CHKERRQ(ierr); 191 ierr = PetscFree(lusol->locc);CHKERRQ(ierr); 192 ierr = PetscFree(lusol->locr);CHKERRQ(ierr); 193 ierr = PetscFree(lusol->iploc);CHKERRQ(ierr); 194 ierr = PetscFree(lusol->iqloc);CHKERRQ(ierr); 195 ierr = PetscFree(lusol->ipinv);CHKERRQ(ierr); 196 ierr = PetscFree(lusol->iqinv);CHKERRQ(ierr); 197 ierr = PetscFree(lusol->mnsw);CHKERRQ(ierr); 198 ierr = PetscFree(lusol->mnsv);CHKERRQ(ierr); 199 ierr = PetscFree3(lusol->data,lusol->indc,lusol->indr);CHKERRQ(ierr); 200 } 201 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 202 ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr); 203 PetscFunctionReturn(0); 204 } 205 206 PetscErrorCode MatSolve_LUSOL(Mat A,Vec b,Vec x) 207 { 208 Mat_LUSOL *lusol=(Mat_LUSOL*)A->spptr; 209 double *xx; 210 const double *bb; 211 int mode=5; 212 PetscErrorCode ierr; 213 int i,m,n,nnz,status; 214 215 PetscFunctionBegin; 216 ierr = VecGetArray(x, &xx);CHKERRQ(ierr); 217 ierr = VecGetArrayRead(b, &bb);CHKERRQ(ierr); 218 219 m = n = lusol->n; 220 nnz = lusol->nnz; 221 222 for (i = 0; i < m; i++) lusol->mnsv[i] = bb[i]; 223 224 LU6SOL(&mode, &m, &n, lusol->mnsv, xx, &nnz, 225 lusol->luparm, lusol->parmlu, lusol->data, 226 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 227 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, &status); 228 229 if (status) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"solve failed, error code %d",status); 230 231 ierr = VecRestoreArray(x, &xx);CHKERRQ(ierr); 232 ierr = VecRestoreArrayRead(b, &bb);CHKERRQ(ierr); 233 PetscFunctionReturn(0); 234 } 235 236 PetscErrorCode MatLUFactorNumeric_LUSOL(Mat F,Mat A,const MatFactorInfo *info) 237 { 238 Mat_SeqAIJ *a; 239 Mat_LUSOL *lusol = (Mat_LUSOL*)F->spptr; 240 PetscErrorCode ierr; 241 int m, n, nz, nnz, status; 242 int i, rs, re; 243 int factorizations; 244 245 PetscFunctionBegin; 246 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 247 a = (Mat_SeqAIJ*)A->data; 248 249 if (m != lusol->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"factorization struct inconsistent"); 250 251 factorizations = 0; 252 do { 253 /*******************************************************************/ 254 /* Check the workspace allocation. */ 255 /*******************************************************************/ 256 257 nz = a->nz; 258 nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom*nz)); 259 nnz = PetscMax(nnz, 5*n); 260 261 if (nnz < lusol->luparm[12]) { 262 nnz = (int)(lusol->luroom * lusol->luparm[12]); 263 } else if ((factorizations > 0) && (lusol->luroom < 6)) { 264 lusol->luroom += 0.1; 265 } 266 267 nnz = PetscMax(nnz, (int)(lusol->luroom*(lusol->luparm[22] + lusol->luparm[23]))); 268 269 if (nnz > lusol->nnz) { 270 ierr = PetscFree3(lusol->data,lusol->indc,lusol->indr);CHKERRQ(ierr); 271 ierr = PetscMalloc3(nnz,&lusol->data,nnz,&lusol->indc,nnz,&lusol->indr);CHKERRQ(ierr); 272 lusol->nnz = nnz; 273 } 274 275 /*******************************************************************/ 276 /* Fill in the data for the problem. (1-based Fortran style) */ 277 /*******************************************************************/ 278 279 nz = 0; 280 for (i = 0; i < n; i++) { 281 rs = a->i[i]; 282 re = a->i[i+1]; 283 284 while (rs < re) { 285 if (a->a[rs] != 0.0) { 286 lusol->indc[nz] = i + 1; 287 lusol->indr[nz] = a->j[rs] + 1; 288 lusol->data[nz] = a->a[rs]; 289 nz++; 290 } 291 rs++; 292 } 293 } 294 295 /*******************************************************************/ 296 /* Do the factorization. */ 297 /*******************************************************************/ 298 299 LU1FAC(&m, &n, &nz, &nnz, 300 lusol->luparm, lusol->parmlu, lusol->data, 301 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 302 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, 303 lusol->iploc, lusol->iqloc, lusol->ipinv, 304 lusol->iqinv, lusol->mnsw, &status); 305 306 switch (status) { 307 case 0: /* factored */ 308 break; 309 310 case 7: /* insufficient memory */ 311 break; 312 313 case 1: 314 case -1: /* singular */ 315 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Singular matrix"); 316 317 case 3: 318 case 4: /* error conditions */ 319 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"matrix error"); 320 321 default: /* unknown condition */ 322 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"matrix unknown return code"); 323 } 324 325 factorizations++; 326 } while (status == 7); 327 F->ops->solve = MatSolve_LUSOL; 328 F->assembled = PETSC_TRUE; 329 F->preallocated = PETSC_TRUE; 330 PetscFunctionReturn(0); 331 } 332 333 PetscErrorCode MatLUFactorSymbolic_LUSOL(Mat F,Mat A, IS r, IS c,const MatFactorInfo *info) 334 { 335 /************************************************************************/ 336 /* Input */ 337 /* A - matrix to factor */ 338 /* r - row permutation (ignored) */ 339 /* c - column permutation (ignored) */ 340 /* */ 341 /* Output */ 342 /* F - matrix storing the factorization; */ 343 /************************************************************************/ 344 Mat_LUSOL *lusol; 345 PetscErrorCode ierr; 346 int i, m, n, nz, nnz; 347 348 PetscFunctionBegin; 349 /************************************************************************/ 350 /* Check the arguments. */ 351 /************************************************************************/ 352 353 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 354 nz = ((Mat_SeqAIJ*)A->data)->nz; 355 356 /************************************************************************/ 357 /* Create the factorization. */ 358 /************************************************************************/ 359 360 F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 361 lusol = (Mat_LUSOL*)(F->spptr); 362 363 /************************************************************************/ 364 /* Initialize parameters */ 365 /************************************************************************/ 366 367 for (i = 0; i < 30; i++) { 368 lusol->luparm[i] = 0; 369 lusol->parmlu[i] = 0; 370 } 371 372 lusol->luparm[1] = -1; 373 lusol->luparm[2] = 5; 374 lusol->luparm[7] = 1; 375 376 lusol->parmlu[0] = 1 / Factorization_Tolerance; 377 lusol->parmlu[1] = 1 / Factorization_Tolerance; 378 lusol->parmlu[2] = Factorization_Small_Tolerance; 379 lusol->parmlu[3] = Factorization_Pivot_Tolerance; 380 lusol->parmlu[4] = Factorization_Pivot_Tolerance; 381 lusol->parmlu[5] = 3.0; 382 lusol->parmlu[6] = 0.3; 383 lusol->parmlu[7] = 0.6; 384 385 /************************************************************************/ 386 /* Allocate the workspace needed by LUSOL. */ 387 /************************************************************************/ 388 389 lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill); 390 nnz = PetscMax((int)(lusol->elbowroom*nz), 5*n); 391 392 lusol->n = n; 393 lusol->nz = nz; 394 lusol->nnz = nnz; 395 lusol->luroom = 1.75; 396 397 ierr = PetscMalloc(sizeof(int)*n,&lusol->ip); 398 ierr = PetscMalloc(sizeof(int)*n,&lusol->iq); 399 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenc); 400 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenr); 401 ierr = PetscMalloc(sizeof(int)*n,&lusol->locc); 402 ierr = PetscMalloc(sizeof(int)*n,&lusol->locr); 403 ierr = PetscMalloc(sizeof(int)*n,&lusol->iploc); 404 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqloc); 405 ierr = PetscMalloc(sizeof(int)*n,&lusol->ipinv); 406 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqinv); 407 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsw); 408 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsv); 409 410 ierr = PetscMalloc3(nnz,&lusol->data,nnz,&lusol->indc,nnz,&lusol->indr);CHKERRQ(ierr); 411 412 lusol->CleanUpLUSOL = PETSC_TRUE; 413 F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 414 PetscFunctionReturn(0); 415 } 416 417 PetscErrorCode MatFactorGetSolverType_seqaij_lusol(Mat A,MatSolverType *type) 418 { 419 PetscFunctionBegin; 420 *type = MATSOLVERLUSOL; 421 PetscFunctionReturn(0); 422 } 423 424 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_lusol(Mat A,MatFactorType ftype,Mat *F) 425 { 426 Mat B; 427 Mat_LUSOL *lusol; 428 PetscErrorCode ierr; 429 int m, n; 430 431 PetscFunctionBegin; 432 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 433 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 434 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 435 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 436 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 437 438 ierr = PetscNewLog(B,&lusol);CHKERRQ(ierr); 439 B->spptr = lusol; 440 441 B->trivialsymbolic = PETSC_TRUE; 442 B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL; 443 B->ops->destroy = MatDestroy_LUSOL; 444 445 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_lusol);CHKERRQ(ierr); 446 447 B->factortype = MAT_FACTOR_LU; 448 ierr = PetscFree(B->solvertype);CHKERRQ(ierr); 449 ierr = PetscStrallocpy(MATSOLVERLUSOL,&B->solvertype);CHKERRQ(ierr); 450 451 PetscFunctionReturn(0); 452 } 453 454 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Lusol(void) 455 { 456 PetscErrorCode ierr; 457 458 PetscFunctionBegin; 459 ierr = MatSolverTypeRegister(MATSOLVERLUSOL,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 460 PetscFunctionReturn(0); 461 } 462 463 /*MC 464 MATSOLVERLUSOL - "lusol" - Provides direct solvers (LU) for sequential matrices 465 via the external package LUSOL. 466 467 If LUSOL is installed (see the manual for 468 instructions on how to declare the existence of external packages), 469 470 Works with MATSEQAIJ matrices 471 472 Level: beginner 473 474 .seealso: PCLU, PCFactorSetMatSolverType(), MatSolverType 475 476 M*/ 477