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