1 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 /* 23 Dummy symbols that the MINOS files mi25bfac.f and mi15blas.f may require 24 */ 25 PETSC_EXTERN void PETSC_STDCALL M1PAGE() 26 { 27 ; 28 } 29 PETSC_EXTERN void PETSC_STDCALL M5SETX() 30 { 31 ; 32 } 33 34 PETSC_EXTERN void PETSC_STDCALL M6RDEL() 35 { 36 ; 37 } 38 39 PETSC_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 PETSC_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 extern PetscErrorCode MatDuplicate_LUSOL(Mat,MatDuplicateOption,Mat*); 52 53 typedef struct { 54 double *data; 55 int *indc; 56 int *indr; 57 58 int *ip; 59 int *iq; 60 int *lenc; 61 int *lenr; 62 int *locc; 63 int *locr; 64 int *iploc; 65 int *iqloc; 66 int *ipinv; 67 int *iqinv; 68 double *mnsw; 69 double *mnsv; 70 71 double elbowroom; 72 double luroom; /* Extra space allocated when factor fails */ 73 double parmlu[30]; /* Input/output to LUSOL */ 74 75 int n; /* Number of rows/columns in matrix */ 76 int nz; /* Number of nonzeros */ 77 int nnz; /* Number of nonzeros allocated for factors */ 78 int luparm[30]; /* Input/output to LUSOL */ 79 80 PetscBool CleanUpLUSOL; 81 82 } Mat_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 #undef __FUNCT__ 181 #define __FUNCT__ "MatDestroy_LUSOL" 182 PetscErrorCode MatDestroy_LUSOL(Mat A) 183 { 184 PetscErrorCode ierr; 185 Mat_LUSOL *lusol=(Mat_LUSOL*)A->spptr; 186 187 PetscFunctionBegin; 188 if (lusol && lusol->CleanUpLUSOL) { 189 ierr = PetscFree(lusol->ip);CHKERRQ(ierr); 190 ierr = PetscFree(lusol->iq);CHKERRQ(ierr); 191 ierr = PetscFree(lusol->lenc);CHKERRQ(ierr); 192 ierr = PetscFree(lusol->lenr);CHKERRQ(ierr); 193 ierr = PetscFree(lusol->locc);CHKERRQ(ierr); 194 ierr = PetscFree(lusol->locr);CHKERRQ(ierr); 195 ierr = PetscFree(lusol->iploc);CHKERRQ(ierr); 196 ierr = PetscFree(lusol->iqloc);CHKERRQ(ierr); 197 ierr = PetscFree(lusol->ipinv);CHKERRQ(ierr); 198 ierr = PetscFree(lusol->iqinv);CHKERRQ(ierr); 199 ierr = PetscFree(lusol->mnsw);CHKERRQ(ierr); 200 ierr = PetscFree(lusol->mnsv);CHKERRQ(ierr); 201 ierr = PetscFree3(lusol->data,lusol->indc,lusol->indr);CHKERRQ(ierr); 202 } 203 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 204 ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr); 205 PetscFunctionReturn(0); 206 } 207 208 #undef __FUNCT__ 209 #define __FUNCT__ "MatSolve_LUSOL" 210 PetscErrorCode MatSolve_LUSOL(Mat A,Vec b,Vec x) 211 { 212 Mat_LUSOL *lusol=(Mat_LUSOL*)A->spptr; 213 double *xx; 214 const double *bb; 215 int mode=5; 216 PetscErrorCode ierr; 217 int i,m,n,nnz,status; 218 219 PetscFunctionBegin; 220 ierr = VecGetArray(x, &xx);CHKERRQ(ierr); 221 ierr = VecGetArrayRead(b, &bb);CHKERRQ(ierr); 222 223 m = n = lusol->n; 224 nnz = lusol->nnz; 225 226 for (i = 0; i < m; i++) lusol->mnsv[i] = bb[i]; 227 228 LU6SOL(&mode, &m, &n, lusol->mnsv, xx, &nnz, 229 lusol->luparm, lusol->parmlu, lusol->data, 230 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 231 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, &status); 232 233 if (status) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"solve failed, error code %d",status); 234 235 ierr = VecRestoreArray(x, &xx);CHKERRQ(ierr); 236 ierr = VecRestoreArrayRead(b, &bb);CHKERRQ(ierr); 237 PetscFunctionReturn(0); 238 } 239 240 #undef __FUNCT__ 241 #define __FUNCT__ "MatLUFactorNumeric_LUSOL" 242 PetscErrorCode MatLUFactorNumeric_LUSOL(Mat F,Mat A,const MatFactorInfo *info) 243 { 244 Mat_SeqAIJ *a; 245 Mat_LUSOL *lusol = (Mat_LUSOL*)F->spptr; 246 PetscErrorCode ierr; 247 int m, n, nz, nnz, status; 248 int i, rs, re; 249 int factorizations; 250 251 PetscFunctionBegin; 252 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);CHKERRQ(ierr); 253 a = (Mat_SeqAIJ*)A->data; 254 255 if (m != lusol->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"factorization struct inconsistent"); 256 257 factorizations = 0; 258 do { 259 /*******************************************************************/ 260 /* Check the workspace allocation. */ 261 /*******************************************************************/ 262 263 nz = a->nz; 264 nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom*nz)); 265 nnz = PetscMax(nnz, 5*n); 266 267 if (nnz < lusol->luparm[12]) { 268 nnz = (int)(lusol->luroom * lusol->luparm[12]); 269 } else if ((factorizations > 0) && (lusol->luroom < 6)) { 270 lusol->luroom += 0.1; 271 } 272 273 nnz = PetscMax(nnz, (int)(lusol->luroom*(lusol->luparm[22] + lusol->luparm[23]))); 274 275 if (nnz > lusol->nnz) { 276 ierr = PetscFree3(lusol->data,lusol->indc,lusol->indr);CHKERRQ(ierr); 277 ierr = PetscMalloc3(nnz,&lusol->data,nnz,&lusol->indc,nnz,&lusol->indr);CHKERRQ(ierr); 278 lusol->nnz = nnz; 279 } 280 281 /*******************************************************************/ 282 /* Fill in the data for the problem. (1-based Fortran style) */ 283 /*******************************************************************/ 284 285 nz = 0; 286 for (i = 0; i < n; i++) { 287 rs = a->i[i]; 288 re = a->i[i+1]; 289 290 while (rs < re) { 291 if (a->a[rs] != 0.0) { 292 lusol->indc[nz] = i + 1; 293 lusol->indr[nz] = a->j[rs] + 1; 294 lusol->data[nz] = a->a[rs]; 295 nz++; 296 } 297 rs++; 298 } 299 } 300 301 /*******************************************************************/ 302 /* Do the factorization. */ 303 /*******************************************************************/ 304 305 LU1FAC(&m, &n, &nz, &nnz, 306 lusol->luparm, lusol->parmlu, lusol->data, 307 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 308 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, 309 lusol->iploc, lusol->iqloc, lusol->ipinv, 310 lusol->iqinv, lusol->mnsw, &status); 311 312 switch (status) { 313 case 0: /* factored */ 314 break; 315 316 case 7: /* insufficient memory */ 317 break; 318 319 case 1: 320 case -1: /* singular */ 321 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Singular matrix"); 322 323 case 3: 324 case 4: /* error conditions */ 325 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"matrix error"); 326 327 default: /* unknown condition */ 328 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"matrix unknown return code"); 329 } 330 331 factorizations++; 332 } while (status == 7); 333 F->ops->solve = MatSolve_LUSOL; 334 F->assembled = PETSC_TRUE; 335 F->preallocated = PETSC_TRUE; 336 PetscFunctionReturn(0); 337 } 338 339 #undef __FUNCT__ 340 #define __FUNCT__ "MatLUFactorSymbolic_LUSOL" 341 PetscErrorCode MatLUFactorSymbolic_LUSOL(Mat F,Mat A, IS r, IS c,const MatFactorInfo *info) 342 { 343 /************************************************************************/ 344 /* Input */ 345 /* A - matrix to factor */ 346 /* r - row permutation (ignored) */ 347 /* c - column permutation (ignored) */ 348 /* */ 349 /* Output */ 350 /* F - matrix storing the factorization; */ 351 /************************************************************************/ 352 Mat_LUSOL *lusol; 353 PetscErrorCode ierr; 354 int i, m, n, nz, nnz; 355 356 PetscFunctionBegin; 357 /************************************************************************/ 358 /* Check the arguments. */ 359 /************************************************************************/ 360 361 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 362 nz = ((Mat_SeqAIJ*)A->data)->nz; 363 364 /************************************************************************/ 365 /* Create the factorization. */ 366 /************************************************************************/ 367 368 F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 369 lusol = (Mat_LUSOL*)(F->spptr); 370 371 /************************************************************************/ 372 /* Initialize parameters */ 373 /************************************************************************/ 374 375 for (i = 0; i < 30; i++) { 376 lusol->luparm[i] = 0; 377 lusol->parmlu[i] = 0; 378 } 379 380 lusol->luparm[1] = -1; 381 lusol->luparm[2] = 5; 382 lusol->luparm[7] = 1; 383 384 lusol->parmlu[0] = 1 / Factorization_Tolerance; 385 lusol->parmlu[1] = 1 / Factorization_Tolerance; 386 lusol->parmlu[2] = Factorization_Small_Tolerance; 387 lusol->parmlu[3] = Factorization_Pivot_Tolerance; 388 lusol->parmlu[4] = Factorization_Pivot_Tolerance; 389 lusol->parmlu[5] = 3.0; 390 lusol->parmlu[6] = 0.3; 391 lusol->parmlu[7] = 0.6; 392 393 /************************************************************************/ 394 /* Allocate the workspace needed by LUSOL. */ 395 /************************************************************************/ 396 397 lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill); 398 nnz = PetscMax((int)(lusol->elbowroom*nz), 5*n); 399 400 lusol->n = n; 401 lusol->nz = nz; 402 lusol->nnz = nnz; 403 lusol->luroom = 1.75; 404 405 ierr = PetscMalloc(sizeof(int)*n,&lusol->ip); 406 ierr = PetscMalloc(sizeof(int)*n,&lusol->iq); 407 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenc); 408 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenr); 409 ierr = PetscMalloc(sizeof(int)*n,&lusol->locc); 410 ierr = PetscMalloc(sizeof(int)*n,&lusol->locr); 411 ierr = PetscMalloc(sizeof(int)*n,&lusol->iploc); 412 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqloc); 413 ierr = PetscMalloc(sizeof(int)*n,&lusol->ipinv); 414 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqinv); 415 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsw); 416 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsv); 417 418 ierr = PetscMalloc3(nnz,&lusol->data,nnz,&lusol->indc,nnz,&lusol->indr);CHKERRQ(ierr); 419 420 lusol->CleanUpLUSOL = PETSC_TRUE; 421 F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 422 PetscFunctionReturn(0); 423 } 424 425 #undef __FUNCT__ 426 #define __FUNCT__ "MatFactorGetSolverPackage_seqaij_lusol" 427 PetscErrorCode MatFactorGetSolverPackage_seqaij_lusol(Mat A,const MatSolverPackage *type) 428 { 429 PetscFunctionBegin; 430 *type = MATSOLVERLUSOL; 431 PetscFunctionReturn(0); 432 } 433 434 #undef __FUNCT__ 435 #define __FUNCT__ "MatGetFactor_seqaij_lusol" 436 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_lusol(Mat A,MatFactorType ftype,Mat *F) 437 { 438 Mat B; 439 Mat_LUSOL *lusol; 440 PetscErrorCode ierr; 441 int m, n; 442 443 PetscFunctionBegin; 444 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 445 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 446 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 447 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 448 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 449 450 ierr = PetscNewLog(B,&lusol);CHKERRQ(ierr); 451 B->spptr = lusol; 452 453 B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL; 454 B->ops->destroy = MatDestroy_LUSOL; 455 456 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqaij_lusol);CHKERRQ(ierr); 457 458 B->factortype = MAT_FACTOR_LU; 459 PetscFunctionReturn(0); 460 } 461 462 #undef __FUNCT__ 463 #define __FUNCT__ "MatSolverPackageRegister_Lusol" 464 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_Lusol(void) 465 { 466 PetscErrorCode ierr; 467 468 PetscFunctionBegin; 469 ierr = MatSolverPackageRegister(MATSOLVERLUSOL,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 470 PetscFunctionReturn(0); 471 } 472 473 /*MC 474 MATSOLVERLUSOL - "lusol" - Provides direct solvers (LU) for sequential matrices 475 via the external package LUSOL. 476 477 If LUSOL is installed (see the manual for 478 instructions on how to declare the existence of external packages), 479 480 Works with MATSEQAIJ matrices 481 482 Level: beginner 483 484 .seealso: PCLU, PCFactorSetMatSolverPackage(), MatSolverPackage 485 486 M*/ 487