1 /*$Id: lusol.c,v 1.11 2001/08/06 21:15:14 bsmith Exp $*/ 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 EXTERN_C_BEGIN 23 /* 24 Dummy symbols that the MINOS files mi25bfac.f and mi15blas.f may require 25 */ 26 void PETSC_STDCALL M1PAGE() { 27 ; 28 } 29 void PETSC_STDCALL M5SETX() { 30 ; 31 } 32 33 void PETSC_STDCALL M6RDEL() { 34 ; 35 } 36 37 extern void PETSC_STDCALL LU1FAC (int *m, int *n, int *nnz, int *size, int *luparm, 38 double *parmlu, double *data, int *indc, int *indr, 39 int *rowperm, int *colperm, int *collen, int *rowlen, 40 int *colstart, int *rowstart, int *rploc, int *cploc, 41 int *rpinv, int *cpinv, double *w, int *inform); 42 43 extern void PETSC_STDCALL LU6SOL (int *mode, int *m, int *n, double *rhs, double *x, 44 int *size, int *luparm, double *parmlu, double *data, 45 int *indc, int *indr, int *rowperm, int *colperm, 46 int *collen, int *rowlen, int *colstart, int *rowstart, 47 int *inform); 48 EXTERN_C_END 49 50 typedef struct 51 { 52 double *data; 53 int *indc; 54 int *indr; 55 56 int *ip; 57 int *iq; 58 int *lenc; 59 int *lenr; 60 int *locc; 61 int *locr; 62 int *iploc; 63 int *iqloc; 64 int *ipinv; 65 int *iqinv; 66 double *mnsw; 67 double *mnsv; 68 69 double elbowroom; 70 double luroom; /* Extra space allocated when factor fails */ 71 double parmlu[30]; /* Input/output to LUSOL */ 72 73 int n; /* Number of rows/columns in matrix */ 74 int nz; /* Number of nonzeros */ 75 int nnz; /* Number of nonzeros allocated for factors */ 76 int luparm[30]; /* Input/output to LUSOL */ 77 78 int (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*); 79 int (*MatDestroy)(Mat); 80 PetscTruth CleanUpLUSOL; 81 82 } Mat_SeqAIJ_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 EXTERN_C_BEGIN 181 #undef __FUNCT__ 182 #define __FUNCT__ "MatConvert_LUSOL_SeqAIJ" 183 int MatConvert_LUSOL_SeqAIJ(Mat A,MatType type,Mat *newmat) { 184 /* This routine is only called to convert an unfactored PETSc-LUSOL matrix */ 185 /* to its base PETSc type, so we will ignore 'MatType type'. */ 186 int ierr; 187 Mat B=*newmat; 188 Mat_SeqAIJ_LUSOL *lusol=(Mat_SeqAIJ_LUSOL *)A->spptr; 189 190 PetscFunctionBegin; 191 if (B != A) { 192 /* This routine was inherited from SeqAIJ. */ 193 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 194 } else { 195 B->ops->lufactorsymbolic = lusol->MatLUFactorSymbolic; 196 B->ops->destroy = lusol->MatDestroy; 197 198 ierr = PetscFree(lusol);CHKERRQ(ierr); 199 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 200 } 201 *newmat = B; 202 PetscFunctionReturn(0); 203 } 204 EXTERN_C_END 205 206 #undef __FUNCT__ 207 #define __FUNCT__ "MatDestroy_SeqAIJ_LUSOL" 208 int MatDestroy_SeqAIJ_LUSOL(Mat A) { 209 int ierr; 210 Mat_SeqAIJ_LUSOL *lusol=(Mat_SeqAIJ_LUSOL *)A->spptr; 211 212 PetscFunctionBegin; 213 if (lusol->CleanUpLUSOL) { 214 ierr = PetscFree(lusol->ip);CHKERRQ(ierr); 215 ierr = PetscFree(lusol->iq);CHKERRQ(ierr); 216 ierr = PetscFree(lusol->lenc);CHKERRQ(ierr); 217 ierr = PetscFree(lusol->lenr);CHKERRQ(ierr); 218 ierr = PetscFree(lusol->locc);CHKERRQ(ierr); 219 ierr = PetscFree(lusol->locr);CHKERRQ(ierr); 220 ierr = PetscFree(lusol->iploc);CHKERRQ(ierr); 221 ierr = PetscFree(lusol->iqloc);CHKERRQ(ierr); 222 ierr = PetscFree(lusol->ipinv);CHKERRQ(ierr); 223 ierr = PetscFree(lusol->iqinv);CHKERRQ(ierr); 224 ierr = PetscFree(lusol->mnsw);CHKERRQ(ierr); 225 ierr = PetscFree(lusol->mnsv);CHKERRQ(ierr); 226 227 ierr = PetscFree(lusol->indc);CHKERRQ(ierr); 228 } 229 230 ierr = MatConvert_LUSOL_SeqAIJ(A,MATSEQAIJ,&A); 231 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 232 PetscFunctionReturn(0); 233 } 234 235 #undef __FUNCT__ 236 #define __FUNCT__ "MatSolve_SeqAIJ_LUSOL" 237 int MatSolve_SeqAIJ_LUSOL(Mat A,Vec b,Vec x) 238 { 239 Mat_SeqAIJ_LUSOL *lusol = (Mat_SeqAIJ_LUSOL *)A->spptr; 240 double *bb, *xx; 241 int mode = 5; 242 int i, m, n, nnz, status, ierr; 243 244 PetscFunctionBegin; 245 ierr = VecGetArray(x, &xx);CHKERRQ(ierr); 246 ierr = VecGetArray(b, &bb);CHKERRQ(ierr); 247 248 m = n = lusol->n; 249 nnz = lusol->nnz; 250 251 for (i = 0; i < m; i++) 252 { 253 lusol->mnsv[i] = bb[i]; 254 } 255 256 LU6SOL(&mode, &m, &n, lusol->mnsv, xx, &nnz, 257 lusol->luparm, lusol->parmlu, lusol->data, 258 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 259 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, &status); 260 261 if (status != 0) 262 { 263 SETERRQ(PETSC_ERR_ARG_SIZ,"solve failed"); 264 } 265 266 ierr = VecRestoreArray(x, &xx);CHKERRQ(ierr); 267 ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr); 268 PetscFunctionReturn(0); 269 } 270 271 #undef __FUNCT__ 272 #define __FUNCT__ "MatLUFactorNumeric_SeqAIJ_LUSOL" 273 int MatLUFactorNumeric_SeqAIJ_LUSOL(Mat A, Mat *F) 274 { 275 Mat_SeqAIJ *a; 276 Mat_SeqAIJ_LUSOL *lusol = (Mat_SeqAIJ_LUSOL *)(*F)->spptr; 277 int m, n, nz, nnz, status; 278 int i, rs, re,ierr; 279 int factorizations; 280 281 PetscFunctionBegin; 282 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);CHKERRQ(ierr); 283 a = (Mat_SeqAIJ *)A->data; 284 285 if (m != lusol->n) { 286 SETERRQ(PETSC_ERR_ARG_SIZ,"factorization struct inconsistent"); 287 } 288 289 factorizations = 0; 290 do 291 { 292 /*******************************************************************/ 293 /* Check the workspace allocation. */ 294 /*******************************************************************/ 295 296 nz = a->nz; 297 nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom*nz)); 298 nnz = PetscMax(nnz, 5*n); 299 300 if (nnz < lusol->luparm[12]){ 301 nnz = (int)(lusol->luroom * lusol->luparm[12]); 302 } else if ((factorizations > 0) && (lusol->luroom < 6)){ 303 lusol->luroom += 0.1; 304 } 305 306 nnz = PetscMax(nnz, (int)(lusol->luroom*(lusol->luparm[22] + lusol->luparm[23]))); 307 308 if (nnz > lusol->nnz){ 309 ierr = PetscFree(lusol->indc);CHKERRQ(ierr); 310 ierr = PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc);CHKERRQ(ierr); 311 lusol->indr = lusol->indc + nnz; 312 lusol->data = (double *)(lusol->indr + nnz); 313 lusol->nnz = nnz; 314 } 315 316 /*******************************************************************/ 317 /* Fill in the data for the problem. (1-based Fortran style) */ 318 /*******************************************************************/ 319 320 nz = 0; 321 for (i = 0; i < n; i++) 322 { 323 rs = a->i[i]; 324 re = a->i[i+1]; 325 326 while (rs < re) 327 { 328 if (a->a[rs] != 0.0) 329 { 330 lusol->indc[nz] = i + 1; 331 lusol->indr[nz] = a->j[rs] + 1; 332 lusol->data[nz] = a->a[rs]; 333 nz++; 334 } 335 rs++; 336 } 337 } 338 339 /*******************************************************************/ 340 /* Do the factorization. */ 341 /*******************************************************************/ 342 343 LU1FAC(&m, &n, &nz, &nnz, 344 lusol->luparm, lusol->parmlu, lusol->data, 345 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 346 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, 347 lusol->iploc, lusol->iqloc, lusol->ipinv, 348 lusol->iqinv, lusol->mnsw, &status); 349 350 switch(status) 351 { 352 case 0: /* factored */ 353 break; 354 355 case 7: /* insufficient memory */ 356 break; 357 358 case 1: 359 case -1: /* singular */ 360 SETERRQ(1,"Singular matrix"); 361 362 case 3: 363 case 4: /* error conditions */ 364 SETERRQ(1,"matrix error"); 365 366 default: /* unknown condition */ 367 SETERRQ(1,"matrix unknown return code"); 368 } 369 370 factorizations++; 371 } while (status == 7); 372 (*F)->assembled = PETSC_TRUE; 373 PetscFunctionReturn(0); 374 } 375 376 #undef __FUNCT__ 377 #define __FUNCT__ "MatLUFactorSymbolic_SeqAIJ_LUSOL" 378 int MatLUFactorSymbolic_SeqAIJ_LUSOL(Mat A, IS r, IS c,MatFactorInfo *info, Mat *F) 379 { 380 /************************************************************************/ 381 /* Input */ 382 /* A - matrix to factor */ 383 /* r - row permutation (ignored) */ 384 /* c - column permutation (ignored) */ 385 /* */ 386 /* Output */ 387 /* F - matrix storing the factorization; */ 388 /************************************************************************/ 389 Mat B; 390 Mat_SeqAIJ_LUSOL *lusol; 391 int ierr,i, m, n, nz, nnz; 392 393 PetscFunctionBegin; 394 395 /************************************************************************/ 396 /* Check the arguments. */ 397 /************************************************************************/ 398 399 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 400 nz = ((Mat_SeqAIJ *)A->data)->nz; 401 402 /************************************************************************/ 403 /* Create the factorization. */ 404 /************************************************************************/ 405 406 ierr = MatCreate(A->comm,PETSC_DECIDE,PETSC_DECIDE,m,n,&B);CHKERRQ(ierr); 407 ierr = MatSetType(B,MATLUSOL);CHKERRQ(ierr); 408 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 409 410 B->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_LUSOL; 411 B->ops->solve = MatSolve_SeqAIJ_LUSOL; 412 B->factor = FACTOR_LU; 413 lusol = (Mat_SeqAIJ_LUSOL*)(B->spptr); 414 415 /************************************************************************/ 416 /* Initialize parameters */ 417 /************************************************************************/ 418 419 for (i = 0; i < 30; i++) 420 { 421 lusol->luparm[i] = 0; 422 lusol->parmlu[i] = 0; 423 } 424 425 lusol->luparm[1] = -1; 426 lusol->luparm[2] = 5; 427 lusol->luparm[7] = 1; 428 429 lusol->parmlu[0] = 1 / Factorization_Tolerance; 430 lusol->parmlu[1] = 1 / Factorization_Tolerance; 431 lusol->parmlu[2] = Factorization_Small_Tolerance; 432 lusol->parmlu[3] = Factorization_Pivot_Tolerance; 433 lusol->parmlu[4] = Factorization_Pivot_Tolerance; 434 lusol->parmlu[5] = 3.0; 435 lusol->parmlu[6] = 0.3; 436 lusol->parmlu[7] = 0.6; 437 438 /************************************************************************/ 439 /* Allocate the workspace needed by LUSOL. */ 440 /************************************************************************/ 441 442 lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill); 443 nnz = PetscMax((int)(lusol->elbowroom*nz), 5*n); 444 445 lusol->n = n; 446 lusol->nz = nz; 447 lusol->nnz = nnz; 448 lusol->luroom = 1.75; 449 450 ierr = PetscMalloc(sizeof(int)*n,&lusol->ip); 451 ierr = PetscMalloc(sizeof(int)*n,&lusol->iq); 452 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenc); 453 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenr); 454 ierr = PetscMalloc(sizeof(int)*n,&lusol->locc); 455 ierr = PetscMalloc(sizeof(int)*n,&lusol->locr); 456 ierr = PetscMalloc(sizeof(int)*n,&lusol->iploc); 457 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqloc); 458 ierr = PetscMalloc(sizeof(int)*n,&lusol->ipinv); 459 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqinv); 460 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsw); 461 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsv); 462 463 ierr = PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc); 464 lusol->indr = lusol->indc + nnz; 465 lusol->data = (double *)(lusol->indr + nnz); 466 lusol->CleanUpLUSOL = PETSC_TRUE; 467 *F = B; 468 PetscFunctionReturn(0); 469 } 470 471 EXTERN_C_BEGIN 472 #undef __FUNCT__ 473 #define __FUNCT__ "MatConvert_SeqAIJ_LUSOL" 474 int MatConvert_SeqAIJ_LUSOL(Mat A,MatType type,Mat *newmat) { 475 int ierr, m, n; 476 Mat_SeqAIJ_LUSOL *lusol; 477 Mat B=*newmat; 478 479 PetscFunctionBegin; 480 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 481 if (m != n) { 482 SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 483 } 484 if (B != A) { 485 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 486 } 487 488 ierr = PetscNew(Mat_SeqAIJ_LUSOL,&lusol);CHKERRQ(ierr); 489 lusol->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 490 lusol->MatDestroy = A->ops->destroy; 491 lusol->CleanUpLUSOL = PETSC_FALSE; 492 493 B->spptr = (void *)lusol; 494 B->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ_LUSOL; 495 B->ops->destroy = MatDestroy_SeqAIJ_LUSOL; 496 497 PetscLogInfo(0,"Using LUSOL for SeqAIJ LU factorization and solves."); 498 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_lusol_C", 499 "MatConvert_SeqAIJ_LUSOL",MatConvert_SeqAIJ_LUSOL);CHKERRQ(ierr); 500 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_lusol_seqaij_C", 501 "MatConvert_LUSOL_SeqAIJ",MatConvert_LUSOL_SeqAIJ);CHKERRQ(ierr); 502 ierr = PetscObjectChangeTypeName((PetscObject)B,type);CHKERRQ(ierr); 503 *newmat = B; 504 PetscFunctionReturn(0); 505 } 506 EXTERN_C_END 507 508 /*MC 509 MATLUSOL - a matrix type providing direct solvers (LU) for sequential matrices 510 via the external package LUSOL. 511 512 If LUSOL is installed (see the manual for 513 instructions on how to declare the existence of external packages), 514 a matrix type can be constructed which invokes LUSOL solvers. 515 After calling MatCreate(...,A), simply call MatSetType(A,MATLUSOL). 516 This matrix type is only supported for double precision real. 517 518 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 519 supported for this matrix type. 520 521 Options Database Keys: 522 . -mat_type lusol - sets the matrix type to lusol during a call to MatSetFromOptions() 523 524 Level: beginner 525 526 .seealso: PCLU 527 M*/ 528 529 EXTERN_C_BEGIN 530 #undef __FUNCT__ 531 #define __FUNCT__ "MatCreate_SeqAIJ_LUSOL" 532 int MatCreate_SeqAIJ_LUSOL(Mat A) 533 { 534 int ierr; 535 536 PetscFunctionBegin; 537 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 538 ierr = MatConvert_SeqAIJ_LUSOL(A,MATLUSOL,&A);CHKERRQ(ierr); 539 PetscFunctionReturn(0); 540 } 541 EXTERN_C_END 542