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