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