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 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 201 *newmat = B; 202 PetscFunctionReturn(0); 203 } 204 EXTERN_C_END 205 206 #undef __FUNCT__ 207 #define __FUNCT__ "MatDestroy_LUSOL" 208 PetscErrorCode MatDestroy_LUSOL(Mat A) 209 { 210 PetscErrorCode 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 PetscErrorCode MatSolve_LUSOL(Mat A,Vec b,Vec x) 239 { 240 Mat_LUSOL *lusol=(Mat_LUSOL*)A->spptr; 241 double *bb,*xx; 242 int mode=5; 243 PetscErrorCode ierr; 244 int i,m,n,nnz,status; 245 246 PetscFunctionBegin; 247 ierr = VecGetArray(x, &xx);CHKERRQ(ierr); 248 ierr = VecGetArray(b, &bb);CHKERRQ(ierr); 249 250 m = n = lusol->n; 251 nnz = lusol->nnz; 252 253 for (i = 0; i < m; i++) 254 { 255 lusol->mnsv[i] = bb[i]; 256 } 257 258 LU6SOL(&mode, &m, &n, lusol->mnsv, xx, &nnz, 259 lusol->luparm, lusol->parmlu, lusol->data, 260 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 261 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, &status); 262 263 if (status != 0) 264 { 265 SETERRQ(PETSC_ERR_ARG_SIZ,"solve failed"); 266 } 267 268 ierr = VecRestoreArray(x, &xx);CHKERRQ(ierr); 269 ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr); 270 PetscFunctionReturn(0); 271 } 272 273 #undef __FUNCT__ 274 #define __FUNCT__ "MatLUFactorNumeric_LUSOL" 275 PetscErrorCode MatLUFactorNumeric_LUSOL(Mat A, Mat *F) 276 { 277 Mat_SeqAIJ *a; 278 Mat_LUSOL *lusol = (Mat_LUSOL*)(*F)->spptr; 279 PetscErrorCode ierr; 280 int m, n, nz, nnz, status; 281 int i, rs, re; 282 int factorizations; 283 284 PetscFunctionBegin; 285 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);CHKERRQ(ierr); 286 a = (Mat_SeqAIJ *)A->data; 287 288 if (m != lusol->n) { 289 SETERRQ(PETSC_ERR_ARG_SIZ,"factorization struct inconsistent"); 290 } 291 292 factorizations = 0; 293 do 294 { 295 /*******************************************************************/ 296 /* Check the workspace allocation. */ 297 /*******************************************************************/ 298 299 nz = a->nz; 300 nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom*nz)); 301 nnz = PetscMax(nnz, 5*n); 302 303 if (nnz < lusol->luparm[12]){ 304 nnz = (int)(lusol->luroom * lusol->luparm[12]); 305 } else if ((factorizations > 0) && (lusol->luroom < 6)){ 306 lusol->luroom += 0.1; 307 } 308 309 nnz = PetscMax(nnz, (int)(lusol->luroom*(lusol->luparm[22] + lusol->luparm[23]))); 310 311 if (nnz > lusol->nnz){ 312 ierr = PetscFree(lusol->indc);CHKERRQ(ierr); 313 ierr = PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc);CHKERRQ(ierr); 314 lusol->indr = lusol->indc + nnz; 315 lusol->data = (double *)(lusol->indr + nnz); 316 lusol->nnz = nnz; 317 } 318 319 /*******************************************************************/ 320 /* Fill in the data for the problem. (1-based Fortran style) */ 321 /*******************************************************************/ 322 323 nz = 0; 324 for (i = 0; i < n; i++) 325 { 326 rs = a->i[i]; 327 re = a->i[i+1]; 328 329 while (rs < re) 330 { 331 if (a->a[rs] != 0.0) 332 { 333 lusol->indc[nz] = i + 1; 334 lusol->indr[nz] = a->j[rs] + 1; 335 lusol->data[nz] = a->a[rs]; 336 nz++; 337 } 338 rs++; 339 } 340 } 341 342 /*******************************************************************/ 343 /* Do the factorization. */ 344 /*******************************************************************/ 345 346 LU1FAC(&m, &n, &nz, &nnz, 347 lusol->luparm, lusol->parmlu, lusol->data, 348 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 349 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, 350 lusol->iploc, lusol->iqloc, lusol->ipinv, 351 lusol->iqinv, lusol->mnsw, &status); 352 353 switch(status) 354 { 355 case 0: /* factored */ 356 break; 357 358 case 7: /* insufficient memory */ 359 break; 360 361 case 1: 362 case -1: /* singular */ 363 SETERRQ(1,"Singular matrix"); 364 365 case 3: 366 case 4: /* error conditions */ 367 SETERRQ(1,"matrix error"); 368 369 default: /* unknown condition */ 370 SETERRQ(1,"matrix unknown return code"); 371 } 372 373 factorizations++; 374 } while (status == 7); 375 (*F)->assembled = PETSC_TRUE; 376 PetscFunctionReturn(0); 377 } 378 379 #undef __FUNCT__ 380 #define __FUNCT__ "MatLUFactorSymbolic_LUSOL" 381 PetscErrorCode MatLUFactorSymbolic_LUSOL(Mat A, IS r, IS c,MatFactorInfo *info, Mat *F) { 382 /************************************************************************/ 383 /* Input */ 384 /* A - matrix to factor */ 385 /* r - row permutation (ignored) */ 386 /* c - column permutation (ignored) */ 387 /* */ 388 /* Output */ 389 /* F - matrix storing the factorization; */ 390 /************************************************************************/ 391 Mat B; 392 Mat_LUSOL *lusol; 393 PetscErrorCode ierr; 394 int i, m, n, nz, nnz; 395 396 PetscFunctionBegin; 397 398 /************************************************************************/ 399 /* Check the arguments. */ 400 /************************************************************************/ 401 402 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 403 nz = ((Mat_SeqAIJ *)A->data)->nz; 404 405 /************************************************************************/ 406 /* Create the factorization. */ 407 /************************************************************************/ 408 409 ierr = MatCreate(A->comm,PETSC_DECIDE,PETSC_DECIDE,m,n,&B);CHKERRQ(ierr); 410 ierr = MatSetType(B,A->type_name);CHKERRQ(ierr); 411 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 412 413 B->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 414 B->ops->solve = MatSolve_LUSOL; 415 B->factor = FACTOR_LU; 416 lusol = (Mat_LUSOL*)(B->spptr); 417 418 /************************************************************************/ 419 /* Initialize parameters */ 420 /************************************************************************/ 421 422 for (i = 0; i < 30; i++) 423 { 424 lusol->luparm[i] = 0; 425 lusol->parmlu[i] = 0; 426 } 427 428 lusol->luparm[1] = -1; 429 lusol->luparm[2] = 5; 430 lusol->luparm[7] = 1; 431 432 lusol->parmlu[0] = 1 / Factorization_Tolerance; 433 lusol->parmlu[1] = 1 / Factorization_Tolerance; 434 lusol->parmlu[2] = Factorization_Small_Tolerance; 435 lusol->parmlu[3] = Factorization_Pivot_Tolerance; 436 lusol->parmlu[4] = Factorization_Pivot_Tolerance; 437 lusol->parmlu[5] = 3.0; 438 lusol->parmlu[6] = 0.3; 439 lusol->parmlu[7] = 0.6; 440 441 /************************************************************************/ 442 /* Allocate the workspace needed by LUSOL. */ 443 /************************************************************************/ 444 445 lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill); 446 nnz = PetscMax((int)(lusol->elbowroom*nz), 5*n); 447 448 lusol->n = n; 449 lusol->nz = nz; 450 lusol->nnz = nnz; 451 lusol->luroom = 1.75; 452 453 ierr = PetscMalloc(sizeof(int)*n,&lusol->ip); 454 ierr = PetscMalloc(sizeof(int)*n,&lusol->iq); 455 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenc); 456 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenr); 457 ierr = PetscMalloc(sizeof(int)*n,&lusol->locc); 458 ierr = PetscMalloc(sizeof(int)*n,&lusol->locr); 459 ierr = PetscMalloc(sizeof(int)*n,&lusol->iploc); 460 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqloc); 461 ierr = PetscMalloc(sizeof(int)*n,&lusol->ipinv); 462 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqinv); 463 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsw); 464 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsv); 465 466 ierr = PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc); 467 lusol->indr = lusol->indc + nnz; 468 lusol->data = (double *)(lusol->indr + nnz); 469 lusol->CleanUpLUSOL = PETSC_TRUE; 470 *F = B; 471 PetscFunctionReturn(0); 472 } 473 474 EXTERN_C_BEGIN 475 #undef __FUNCT__ 476 #define __FUNCT__ "MatConvert_SeqAIJ_LUSOL" 477 PetscErrorCode MatConvert_SeqAIJ_LUSOL(Mat A,const MatType type,Mat *newmat) { 478 PetscErrorCode ierr; 479 int m, n; 480 Mat_LUSOL *lusol; 481 Mat B=*newmat; 482 483 PetscFunctionBegin; 484 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 485 if (m != n) { 486 SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 487 } 488 if (B != A) { 489 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 490 } 491 492 ierr = PetscNew(Mat_LUSOL,&lusol);CHKERRQ(ierr); 493 lusol->MatDuplicate = A->ops->duplicate; 494 lusol->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 495 lusol->MatDestroy = A->ops->destroy; 496 lusol->CleanUpLUSOL = PETSC_FALSE; 497 498 B->spptr = (void*)lusol; 499 B->ops->duplicate = MatDuplicate_LUSOL; 500 B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL; 501 B->ops->destroy = MatDestroy_LUSOL; 502 503 PetscLogInfo(0,"Using LUSOL for LU factorization and solves."); 504 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_lusol_C", 505 "MatConvert_SeqAIJ_LUSOL",MatConvert_SeqAIJ_LUSOL);CHKERRQ(ierr); 506 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_lusol_seqaij_C", 507 "MatConvert_LUSOL_SeqAIJ",MatConvert_LUSOL_SeqAIJ);CHKERRQ(ierr); 508 ierr = PetscObjectChangeTypeName((PetscObject)B,type);CHKERRQ(ierr); 509 *newmat = B; 510 PetscFunctionReturn(0); 511 } 512 EXTERN_C_END 513 514 #undef __FUNCT__ 515 #define __FUNCT__ "MatDuplicate_LUSOL" 516 PetscErrorCode MatDuplicate_LUSOL(Mat A, MatDuplicateOption op, Mat *M) { 517 PetscErrorCode ierr; 518 Mat_LUSOL *lu=(Mat_LUSOL *)A->spptr; 519 PetscFunctionBegin; 520 ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr); 521 ierr = PetscMemcpy((*M)->spptr,lu,sizeof(Mat_LUSOL));CHKERRQ(ierr); 522 PetscFunctionReturn(0); 523 } 524 525 /*MC 526 MATLUSOL - MATLUSOL = "lusol" - A matrix type providing direct solvers (LU) for sequential matrices 527 via the external package LUSOL. 528 529 If LUSOL is installed (see the manual for 530 instructions on how to declare the existence of external packages), 531 a matrix type can be constructed which invokes LUSOL solvers. 532 After calling MatCreate(...,A), simply call MatSetType(A,MATLUSOL). 533 This matrix type is only supported for double precision real. 534 535 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 536 supported for this matrix type. MatConvert can be called for a fast inplace conversion 537 to and from the MATSEQAIJ matrix type. 538 539 Options Database Keys: 540 . -mat_type lusol - sets the matrix type to "lusol" during a call to MatSetFromOptions() 541 542 Level: beginner 543 544 .seealso: PCLU 545 M*/ 546 547 EXTERN_C_BEGIN 548 #undef __FUNCT__ 549 #define __FUNCT__ "MatCreate_LUSOL" 550 PetscErrorCode MatCreate_LUSOL(Mat A) 551 { 552 PetscErrorCode ierr; 553 554 PetscFunctionBegin; 555 /* Change type name before calling MatSetType to force proper construction of SeqAIJ and LUSOL types */ 556 ierr = PetscObjectChangeTypeName((PetscObject)A,MATLUSOL);CHKERRQ(ierr); 557 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 558 ierr = MatConvert_SeqAIJ_LUSOL(A,MATLUSOL,&A);CHKERRQ(ierr); 559 PetscFunctionReturn(0); 560 } 561 EXTERN_C_END 562