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 EXTERN int MatDuplicate_LUSOL(Mat,MatDuplicateOption,Mat*); 51 52 typedef struct { 53 double *data; 54 int *indc; 55 int *indr; 56 57 int *ip; 58 int *iq; 59 int *lenc; 60 int *lenr; 61 int *locc; 62 int *locr; 63 int *iploc; 64 int *iqloc; 65 int *ipinv; 66 int *iqinv; 67 double *mnsw; 68 double *mnsv; 69 70 double elbowroom; 71 double luroom; /* Extra space allocated when factor fails */ 72 double parmlu[30]; /* Input/output to LUSOL */ 73 74 int n; /* Number of rows/columns in matrix */ 75 int nz; /* Number of nonzeros */ 76 int nnz; /* Number of nonzeros allocated for factors */ 77 int luparm[30]; /* Input/output to LUSOL */ 78 79 int (*MatDuplicate)(Mat,MatDuplicateOption,Mat*); 80 int (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*); 81 int (*MatDestroy)(Mat); 82 PetscTruth CleanUpLUSOL; 83 84 } Mat_LUSOL; 85 86 /* LUSOL input/Output Parameters (Description uses C-style indexes 87 * 88 * Input parameters Typical value 89 * 90 * luparm(0) = nout File number for printed messages. 6 91 * luparm(1) = lprint Print level. 0 92 * < 0 suppresses output. 93 * = 0 gives error messages. 94 * = 1 gives debug output from some of the 95 * other routines in LUSOL. 96 * >= 2 gives the pivot row and column and the 97 * no. of rows and columns involved at 98 * each elimination step in lu1fac. 99 * luparm(2) = maxcol lu1fac: maximum number of columns 5 100 * searched allowed in a Markowitz-type 101 * search for the next pivot element. 102 * For some of the factorization, the 103 * number of rows searched is 104 * maxrow = maxcol - 1. 105 * 106 * 107 * Output parameters 108 * 109 * luparm(9) = inform Return code from last call to any LU routine. 110 * luparm(10) = nsing No. of singularities marked in the 111 * output array w(*). 112 * luparm(11) = jsing Column index of last singularity. 113 * luparm(12) = minlen Minimum recommended value for lena. 114 * luparm(13) = maxlen ? 115 * luparm(14) = nupdat No. of updates performed by the lu8 routines. 116 * luparm(15) = nrank No. of nonempty rows of U. 117 * luparm(16) = ndens1 No. of columns remaining when the density of 118 * the matrix being factorized reached dens1. 119 * luparm(17) = ndens2 No. of columns remaining when the density of 120 * the matrix being factorized reached dens2. 121 * luparm(18) = jumin The column index associated with dumin. 122 * luparm(19) = numl0 No. of columns in initial L. 123 * luparm(20) = lenl0 Size of initial L (no. of nonzeros). 124 * luparm(21) = lenu0 Size of initial U. 125 * luparm(22) = lenl Size of current L. 126 * luparm(23) = lenu Size of current U. 127 * luparm(24) = lrow Length of row file. 128 * luparm(25) = ncp No. of compressions of LU data structures. 129 * luparm(26) = mersum lu1fac: sum of Markowitz merit counts. 130 * luparm(27) = nutri lu1fac: triangular rows in U. 131 * luparm(28) = nltri lu1fac: triangular rows in L. 132 * luparm(29) = 133 * 134 * 135 * Input parameters Typical value 136 * 137 * parmlu(0) = elmax1 Max multiplier allowed in L 10.0 138 * during factor. 139 * parmlu(1) = elmax2 Max multiplier allowed in L 10.0 140 * during updates. 141 * parmlu(2) = small Absolute tolerance for eps**0.8 142 * treating reals as zero. IBM double: 3.0d-13 143 * parmlu(3) = utol1 Absolute tol for flagging eps**0.66667 144 * small diagonals of U. IBM double: 3.7d-11 145 * parmlu(4) = utol2 Relative tol for flagging eps**0.66667 146 * small diagonals of U. IBM double: 3.7d-11 147 * parmlu(5) = uspace Factor limiting waste space in U. 3.0 148 * In lu1fac, the row or column lists 149 * are compressed if their length 150 * exceeds uspace times the length of 151 * either file after the last compression. 152 * parmlu(6) = dens1 The density at which the Markowitz 0.3 153 * strategy should search maxcol columns 154 * and no rows. 155 * parmlu(7) = dens2 the density at which the Markowitz 0.6 156 * strategy should search only 1 column 157 * or (preferably) use a dense LU for 158 * all the remaining rows and columns. 159 * 160 * 161 * Output parameters 162 * 163 * parmlu(9) = amax Maximum element in A. 164 * parmlu(10) = elmax Maximum multiplier in current L. 165 * parmlu(11) = umax Maximum element in current U. 166 * parmlu(12) = dumax Maximum diagonal in U. 167 * parmlu(13) = dumin Minimum diagonal in U. 168 * parmlu(14) = 169 * parmlu(15) = 170 * parmlu(16) = 171 * parmlu(17) = 172 * parmlu(18) = 173 * parmlu(19) = resid lu6sol: residual after solve with U or U'. 174 * ... 175 * parmlu(29) = 176 */ 177 178 #define Factorization_Tolerance 1e-1 179 #define Factorization_Pivot_Tolerance pow(2.2204460492503131E-16, 2.0 / 3.0) 180 #define Factorization_Small_Tolerance 1e-15 /* pow(DBL_EPSILON, 0.8) */ 181 182 EXTERN_C_BEGIN 183 #undef __FUNCT__ 184 #define __FUNCT__ "MatConvert_LUSOL_SeqAIJ" 185 int MatConvert_LUSOL_SeqAIJ(Mat A,const MatType type,Mat *newmat) { 186 /* This routine is only called to convert an unfactored PETSc-LUSOL matrix */ 187 /* to its base PETSc type, so we will ignore 'MatType type'. */ 188 int ierr; 189 Mat B=*newmat; 190 Mat_LUSOL *lusol=(Mat_LUSOL *)A->spptr; 191 192 PetscFunctionBegin; 193 if (B != A) { 194 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 195 } 196 B->ops->duplicate = lusol->MatDuplicate; 197 B->ops->lufactorsymbolic = lusol->MatLUFactorSymbolic; 198 B->ops->destroy = lusol->MatDestroy; 199 200 ierr = PetscFree(lusol);CHKERRQ(ierr); 201 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 202 *newmat = B; 203 PetscFunctionReturn(0); 204 } 205 EXTERN_C_END 206 207 #undef __FUNCT__ 208 #define __FUNCT__ "MatDestroy_LUSOL" 209 int MatDestroy_LUSOL(Mat A) { 210 int ierr; 211 Mat_LUSOL *lusol=(Mat_LUSOL *)A->spptr; 212 213 PetscFunctionBegin; 214 if (lusol->CleanUpLUSOL) { 215 ierr = PetscFree(lusol->ip);CHKERRQ(ierr); 216 ierr = PetscFree(lusol->iq);CHKERRQ(ierr); 217 ierr = PetscFree(lusol->lenc);CHKERRQ(ierr); 218 ierr = PetscFree(lusol->lenr);CHKERRQ(ierr); 219 ierr = PetscFree(lusol->locc);CHKERRQ(ierr); 220 ierr = PetscFree(lusol->locr);CHKERRQ(ierr); 221 ierr = PetscFree(lusol->iploc);CHKERRQ(ierr); 222 ierr = PetscFree(lusol->iqloc);CHKERRQ(ierr); 223 ierr = PetscFree(lusol->ipinv);CHKERRQ(ierr); 224 ierr = PetscFree(lusol->iqinv);CHKERRQ(ierr); 225 ierr = PetscFree(lusol->mnsw);CHKERRQ(ierr); 226 ierr = PetscFree(lusol->mnsv);CHKERRQ(ierr); 227 228 ierr = PetscFree(lusol->indc);CHKERRQ(ierr); 229 } 230 231 ierr = MatConvert_LUSOL_SeqAIJ(A,MATSEQAIJ,&A); 232 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 233 PetscFunctionReturn(0); 234 } 235 236 #undef __FUNCT__ 237 #define __FUNCT__ "MatSolve_LUSOL" 238 int MatSolve_LUSOL(Mat A,Vec b,Vec x) { 239 Mat_LUSOL *lusol=(Mat_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_LUSOL" 273 int MatLUFactorNumeric_LUSOL(Mat A, Mat *F) { 274 Mat_SeqAIJ *a; 275 Mat_LUSOL *lusol = (Mat_LUSOL*)(*F)->spptr; 276 int m, n, nz, nnz, status; 277 int i, rs, re,ierr; 278 int factorizations; 279 280 PetscFunctionBegin; 281 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);CHKERRQ(ierr); 282 a = (Mat_SeqAIJ *)A->data; 283 284 if (m != lusol->n) { 285 SETERRQ(PETSC_ERR_ARG_SIZ,"factorization struct inconsistent"); 286 } 287 288 factorizations = 0; 289 do 290 { 291 /*******************************************************************/ 292 /* Check the workspace allocation. */ 293 /*******************************************************************/ 294 295 nz = a->nz; 296 nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom*nz)); 297 nnz = PetscMax(nnz, 5*n); 298 299 if (nnz < lusol->luparm[12]){ 300 nnz = (int)(lusol->luroom * lusol->luparm[12]); 301 } else if ((factorizations > 0) && (lusol->luroom < 6)){ 302 lusol->luroom += 0.1; 303 } 304 305 nnz = PetscMax(nnz, (int)(lusol->luroom*(lusol->luparm[22] + lusol->luparm[23]))); 306 307 if (nnz > lusol->nnz){ 308 ierr = PetscFree(lusol->indc);CHKERRQ(ierr); 309 ierr = PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc);CHKERRQ(ierr); 310 lusol->indr = lusol->indc + nnz; 311 lusol->data = (double *)(lusol->indr + nnz); 312 lusol->nnz = nnz; 313 } 314 315 /*******************************************************************/ 316 /* Fill in the data for the problem. (1-based Fortran style) */ 317 /*******************************************************************/ 318 319 nz = 0; 320 for (i = 0; i < n; i++) 321 { 322 rs = a->i[i]; 323 re = a->i[i+1]; 324 325 while (rs < re) 326 { 327 if (a->a[rs] != 0.0) 328 { 329 lusol->indc[nz] = i + 1; 330 lusol->indr[nz] = a->j[rs] + 1; 331 lusol->data[nz] = a->a[rs]; 332 nz++; 333 } 334 rs++; 335 } 336 } 337 338 /*******************************************************************/ 339 /* Do the factorization. */ 340 /*******************************************************************/ 341 342 LU1FAC(&m, &n, &nz, &nnz, 343 lusol->luparm, lusol->parmlu, lusol->data, 344 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 345 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, 346 lusol->iploc, lusol->iqloc, lusol->ipinv, 347 lusol->iqinv, lusol->mnsw, &status); 348 349 switch(status) 350 { 351 case 0: /* factored */ 352 break; 353 354 case 7: /* insufficient memory */ 355 break; 356 357 case 1: 358 case -1: /* singular */ 359 SETERRQ(1,"Singular matrix"); 360 361 case 3: 362 case 4: /* error conditions */ 363 SETERRQ(1,"matrix error"); 364 365 default: /* unknown condition */ 366 SETERRQ(1,"matrix unknown return code"); 367 } 368 369 factorizations++; 370 } while (status == 7); 371 (*F)->assembled = PETSC_TRUE; 372 PetscFunctionReturn(0); 373 } 374 375 #undef __FUNCT__ 376 #define __FUNCT__ "MatLUFactorSymbolic_LUSOL" 377 int MatLUFactorSymbolic_LUSOL(Mat A, IS r, IS c,MatFactorInfo *info, Mat *F) { 378 /************************************************************************/ 379 /* Input */ 380 /* A - matrix to factor */ 381 /* r - row permutation (ignored) */ 382 /* c - column permutation (ignored) */ 383 /* */ 384 /* Output */ 385 /* F - matrix storing the factorization; */ 386 /************************************************************************/ 387 Mat B; 388 Mat_LUSOL *lusol; 389 int ierr,i, m, n, nz, nnz; 390 391 PetscFunctionBegin; 392 393 /************************************************************************/ 394 /* Check the arguments. */ 395 /************************************************************************/ 396 397 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 398 nz = ((Mat_SeqAIJ *)A->data)->nz; 399 400 /************************************************************************/ 401 /* Create the factorization. */ 402 /************************************************************************/ 403 404 ierr = MatCreate(A->comm,PETSC_DECIDE,PETSC_DECIDE,m,n,&B);CHKERRQ(ierr); 405 ierr = MatSetType(B,A->type_name);CHKERRQ(ierr); 406 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 407 408 B->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 409 B->ops->solve = MatSolve_LUSOL; 410 B->factor = FACTOR_LU; 411 lusol = (Mat_LUSOL*)(B->spptr); 412 413 /************************************************************************/ 414 /* Initialize parameters */ 415 /************************************************************************/ 416 417 for (i = 0; i < 30; i++) 418 { 419 lusol->luparm[i] = 0; 420 lusol->parmlu[i] = 0; 421 } 422 423 lusol->luparm[1] = -1; 424 lusol->luparm[2] = 5; 425 lusol->luparm[7] = 1; 426 427 lusol->parmlu[0] = 1 / Factorization_Tolerance; 428 lusol->parmlu[1] = 1 / Factorization_Tolerance; 429 lusol->parmlu[2] = Factorization_Small_Tolerance; 430 lusol->parmlu[3] = Factorization_Pivot_Tolerance; 431 lusol->parmlu[4] = Factorization_Pivot_Tolerance; 432 lusol->parmlu[5] = 3.0; 433 lusol->parmlu[6] = 0.3; 434 lusol->parmlu[7] = 0.6; 435 436 /************************************************************************/ 437 /* Allocate the workspace needed by LUSOL. */ 438 /************************************************************************/ 439 440 lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill); 441 nnz = PetscMax((int)(lusol->elbowroom*nz), 5*n); 442 443 lusol->n = n; 444 lusol->nz = nz; 445 lusol->nnz = nnz; 446 lusol->luroom = 1.75; 447 448 ierr = PetscMalloc(sizeof(int)*n,&lusol->ip); 449 ierr = PetscMalloc(sizeof(int)*n,&lusol->iq); 450 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenc); 451 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenr); 452 ierr = PetscMalloc(sizeof(int)*n,&lusol->locc); 453 ierr = PetscMalloc(sizeof(int)*n,&lusol->locr); 454 ierr = PetscMalloc(sizeof(int)*n,&lusol->iploc); 455 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqloc); 456 ierr = PetscMalloc(sizeof(int)*n,&lusol->ipinv); 457 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqinv); 458 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsw); 459 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsv); 460 461 ierr = PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc); 462 lusol->indr = lusol->indc + nnz; 463 lusol->data = (double *)(lusol->indr + nnz); 464 lusol->CleanUpLUSOL = PETSC_TRUE; 465 *F = B; 466 PetscFunctionReturn(0); 467 } 468 469 EXTERN_C_BEGIN 470 #undef __FUNCT__ 471 #define __FUNCT__ "MatConvert_SeqAIJ_LUSOL" 472 int MatConvert_SeqAIJ_LUSOL(Mat A,const MatType type,Mat *newmat) { 473 int ierr, m, n; 474 Mat_LUSOL *lusol; 475 Mat B=*newmat; 476 477 PetscFunctionBegin; 478 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 479 if (m != n) { 480 SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 481 } 482 if (B != A) { 483 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 484 } 485 486 ierr = PetscNew(Mat_LUSOL,&lusol);CHKERRQ(ierr); 487 lusol->MatDuplicate = A->ops->duplicate; 488 lusol->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 489 lusol->MatDestroy = A->ops->destroy; 490 lusol->CleanUpLUSOL = PETSC_FALSE; 491 492 B->spptr = (void *)lusol; 493 B->ops->duplicate = MatDuplicate_LUSOL; 494 B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL; 495 B->ops->destroy = MatDestroy_LUSOL; 496 497 PetscLogInfo(0,"Using LUSOL for 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 #undef __FUNCT__ 509 #define __FUNCT__ "MatDuplicate_LUSOL" 510 int MatDuplicate_LUSOL(Mat A, MatDuplicateOption op, Mat *M) { 511 int ierr; 512 Mat_LUSOL *lu=(Mat_LUSOL *)A->spptr; 513 PetscFunctionBegin; 514 ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr); 515 ierr = MatConvert_SeqAIJ_LUSOL(*M,MATLUSOL,M);CHKERRQ(ierr); 516 ierr = PetscMemcpy((*M)->spptr,lu,sizeof(Mat_LUSOL));CHKERRQ(ierr); 517 PetscFunctionReturn(0); 518 } 519 520 /*MC 521 MATLUSOL - MATLUSOL = "lusol" - A matrix type providing direct solvers (LU) for sequential matrices 522 via the external package LUSOL. 523 524 If LUSOL is installed (see the manual for 525 instructions on how to declare the existence of external packages), 526 a matrix type can be constructed which invokes LUSOL solvers. 527 After calling MatCreate(...,A), simply call MatSetType(A,MATLUSOL). 528 This matrix type is only supported for double precision real. 529 530 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 531 supported for this matrix type. MatConvert can be called for a fast inplace conversion 532 to and from the MATSEQAIJ matrix type. 533 534 Options Database Keys: 535 . -mat_type lusol - sets the matrix type to "lusol" during a call to MatSetFromOptions() 536 537 Level: beginner 538 539 .seealso: PCLU 540 M*/ 541 542 EXTERN_C_BEGIN 543 #undef __FUNCT__ 544 #define __FUNCT__ "MatCreate_LUSOL" 545 int MatCreate_LUSOL(Mat A) { 546 int ierr; 547 548 PetscFunctionBegin; 549 /* Change type name before calling MatSetType to force proper construction of SeqAIJ and LUSOL types */ 550 ierr = PetscObjectChangeTypeName((PetscObject)A,MATLUSOL);CHKERRQ(ierr); 551 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 552 ierr = MatConvert_SeqAIJ_LUSOL(A,MATLUSOL,&A);CHKERRQ(ierr); 553 PetscFunctionReturn(0); 554 } 555 EXTERN_C_END 556