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 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 37 extern void PETSC_STDCALL LU1FAC (int *m, int *n, int *nnz, int *size, int *luparm, 38 double *parmlu, double *data, int *indc, int *indr, 39 int *rowperm, int *colperm, int *collen, int *rowlen, 40 int *colstart, int *rowstart, int *rploc, int *cploc, 41 int *rpinv, int *cpinv, double *w, int *inform); 42 43 extern void PETSC_STDCALL LU6SOL (int *mode, int *m, int *n, double *rhs, double *x, 44 int *size, int *luparm, double *parmlu, double *data, 45 int *indc, int *indr, int *rowperm, int *colperm, 46 int *collen, int *rowlen, int *colstart, int *rowstart, 47 int *inform); 48 EXTERN_C_END 49 50 extern PetscErrorCode MatDuplicate_LUSOL(Mat,MatDuplicateOption,Mat*); 51 52 typedef struct { 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 PetscBool CleanUpLUSOL; 80 81 } Mat_LUSOL; 82 83 /* LUSOL input/Output Parameters (Description uses C-style indexes 84 * 85 * Input parameters Typical value 86 * 87 * luparm(0) = nout File number for printed messages. 6 88 * luparm(1) = lprint Print level. 0 89 * < 0 suppresses output. 90 * = 0 gives error messages. 91 * = 1 gives debug output from some of the 92 * other routines in LUSOL. 93 * >= 2 gives the pivot row and column and the 94 * no. of rows and columns involved at 95 * each elimination step in lu1fac. 96 * luparm(2) = maxcol lu1fac: maximum number of columns 5 97 * searched allowed in a Markowitz-type 98 * search for the next pivot element. 99 * For some of the factorization, the 100 * number of rows searched is 101 * maxrow = maxcol - 1. 102 * 103 * 104 * Output parameters 105 * 106 * luparm(9) = inform Return code from last call to any LU routine. 107 * luparm(10) = nsing No. of singularities marked in the 108 * output array w(*). 109 * luparm(11) = jsing Column index of last singularity. 110 * luparm(12) = minlen Minimum recommended value for lena. 111 * luparm(13) = maxlen ? 112 * luparm(14) = nupdat No. of updates performed by the lu8 routines. 113 * luparm(15) = nrank No. of nonempty rows of U. 114 * luparm(16) = ndens1 No. of columns remaining when the density of 115 * the matrix being factorized reached dens1. 116 * luparm(17) = ndens2 No. of columns remaining when the density of 117 * the matrix being factorized reached dens2. 118 * luparm(18) = jumin The column index associated with dumin. 119 * luparm(19) = numl0 No. of columns in initial L. 120 * luparm(20) = lenl0 Size of initial L (no. of nonzeros). 121 * luparm(21) = lenu0 Size of initial U. 122 * luparm(22) = lenl Size of current L. 123 * luparm(23) = lenu Size of current U. 124 * luparm(24) = lrow Length of row file. 125 * luparm(25) = ncp No. of compressions of LU data structures. 126 * luparm(26) = mersum lu1fac: sum of Markowitz merit counts. 127 * luparm(27) = nutri lu1fac: triangular rows in U. 128 * luparm(28) = nltri lu1fac: triangular rows in L. 129 * luparm(29) = 130 * 131 * 132 * Input parameters Typical value 133 * 134 * parmlu(0) = elmax1 Max multiplier allowed in L 10.0 135 * during factor. 136 * parmlu(1) = elmax2 Max multiplier allowed in L 10.0 137 * during updates. 138 * parmlu(2) = small Absolute tolerance for eps**0.8 139 * treating reals as zero. IBM double: 3.0d-13 140 * parmlu(3) = utol1 Absolute tol for flagging eps**0.66667 141 * small diagonals of U. IBM double: 3.7d-11 142 * parmlu(4) = utol2 Relative tol for flagging eps**0.66667 143 * small diagonals of U. IBM double: 3.7d-11 144 * parmlu(5) = uspace Factor limiting waste space in U. 3.0 145 * In lu1fac, the row or column lists 146 * are compressed if their length 147 * exceeds uspace times the length of 148 * either file after the last compression. 149 * parmlu(6) = dens1 The density at which the Markowitz 0.3 150 * strategy should search maxcol columns 151 * and no rows. 152 * parmlu(7) = dens2 the density at which the Markowitz 0.6 153 * strategy should search only 1 column 154 * or (preferably) use a dense LU for 155 * all the remaining rows and columns. 156 * 157 * 158 * Output parameters 159 * 160 * parmlu(9) = amax Maximum element in A. 161 * parmlu(10) = elmax Maximum multiplier in current L. 162 * parmlu(11) = umax Maximum element in current U. 163 * parmlu(12) = dumax Maximum diagonal in U. 164 * parmlu(13) = dumin Minimum diagonal in U. 165 * parmlu(14) = 166 * parmlu(15) = 167 * parmlu(16) = 168 * parmlu(17) = 169 * parmlu(18) = 170 * parmlu(19) = resid lu6sol: residual after solve with U or U'. 171 * ... 172 * parmlu(29) = 173 */ 174 175 #define Factorization_Tolerance 1e-1 176 #define Factorization_Pivot_Tolerance pow(2.2204460492503131E-16, 2.0 / 3.0) 177 #define Factorization_Small_Tolerance 1e-15 /* pow(DBL_EPSILON, 0.8) */ 178 179 #undef __FUNCT__ 180 #define __FUNCT__ "MatDestroy_LUSOL" 181 PetscErrorCode MatDestroy_LUSOL(Mat A) 182 { 183 PetscErrorCode ierr; 184 Mat_LUSOL *lusol=(Mat_LUSOL *)A->spptr; 185 186 PetscFunctionBegin; 187 if (lusol && lusol->CleanUpLUSOL) { 188 ierr = PetscFree(lusol->ip);CHKERRQ(ierr); 189 ierr = PetscFree(lusol->iq);CHKERRQ(ierr); 190 ierr = PetscFree(lusol->lenc);CHKERRQ(ierr); 191 ierr = PetscFree(lusol->lenr);CHKERRQ(ierr); 192 ierr = PetscFree(lusol->locc);CHKERRQ(ierr); 193 ierr = PetscFree(lusol->locr);CHKERRQ(ierr); 194 ierr = PetscFree(lusol->iploc);CHKERRQ(ierr); 195 ierr = PetscFree(lusol->iqloc);CHKERRQ(ierr); 196 ierr = PetscFree(lusol->ipinv);CHKERRQ(ierr); 197 ierr = PetscFree(lusol->iqinv);CHKERRQ(ierr); 198 ierr = PetscFree(lusol->mnsw);CHKERRQ(ierr); 199 ierr = PetscFree(lusol->mnsv);CHKERRQ(ierr); 200 ierr = PetscFree3(lusol->data,lusol->indc,lusol->indr);CHKERRQ(ierr); 201 } 202 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 203 ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr); 204 PetscFunctionReturn(0); 205 } 206 207 #undef __FUNCT__ 208 #define __FUNCT__ "MatSolve_LUSOL" 209 PetscErrorCode MatSolve_LUSOL(Mat A,Vec b,Vec x) 210 { 211 Mat_LUSOL *lusol=(Mat_LUSOL*)A->spptr; 212 double *bb,*xx; 213 int mode=5; 214 PetscErrorCode ierr; 215 int i,m,n,nnz,status; 216 217 PetscFunctionBegin; 218 ierr = VecGetArray(x, &xx);CHKERRQ(ierr); 219 ierr = VecGetArray(b, &bb);CHKERRQ(ierr); 220 221 m = n = lusol->n; 222 nnz = lusol->nnz; 223 224 for (i = 0; i < m; i++) 225 { 226 lusol->mnsv[i] = bb[i]; 227 } 228 229 LU6SOL(&mode, &m, &n, lusol->mnsv, xx, &nnz, 230 lusol->luparm, lusol->parmlu, lusol->data, 231 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 232 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, &status); 233 234 if (status) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"solve failed, error code %d",status); 235 236 ierr = VecRestoreArray(x, &xx);CHKERRQ(ierr); 237 ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr); 238 PetscFunctionReturn(0); 239 } 240 241 #undef __FUNCT__ 242 #define __FUNCT__ "MatLUFactorNumeric_LUSOL" 243 PetscErrorCode MatLUFactorNumeric_LUSOL(Mat F,Mat A,const MatFactorInfo *info) 244 { 245 Mat_SeqAIJ *a; 246 Mat_LUSOL *lusol = (Mat_LUSOL*)F->spptr; 247 PetscErrorCode ierr; 248 int m, n, nz, nnz, status; 249 int i, rs, re; 250 int factorizations; 251 252 PetscFunctionBegin; 253 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);CHKERRQ(ierr); 254 a = (Mat_SeqAIJ *)A->data; 255 256 if (m != lusol->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"factorization struct inconsistent"); 257 258 factorizations = 0; 259 do 260 { 261 /*******************************************************************/ 262 /* Check the workspace allocation. */ 263 /*******************************************************************/ 264 265 nz = a->nz; 266 nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom*nz)); 267 nnz = PetscMax(nnz, 5*n); 268 269 if (nnz < lusol->luparm[12]){ 270 nnz = (int)(lusol->luroom * lusol->luparm[12]); 271 } else if ((factorizations > 0) && (lusol->luroom < 6)){ 272 lusol->luroom += 0.1; 273 } 274 275 nnz = PetscMax(nnz, (int)(lusol->luroom*(lusol->luparm[22] + lusol->luparm[23]))); 276 277 if (nnz > lusol->nnz){ 278 ierr = PetscFree3(lusol->data,lusol->indc,lusol->indr);CHKERRQ(ierr); 279 ierr = PetscMalloc3(nnz,double,&lusol->data,nnz,PetscInt,&lusol->indc,nnz,PetscInt,&lusol->indr);CHKERRQ(ierr); 280 lusol->nnz = nnz; 281 } 282 283 /*******************************************************************/ 284 /* Fill in the data for the problem. (1-based Fortran style) */ 285 /*******************************************************************/ 286 287 nz = 0; 288 for (i = 0; i < n; i++) 289 { 290 rs = a->i[i]; 291 re = a->i[i+1]; 292 293 while (rs < re) 294 { 295 if (a->a[rs] != 0.0) 296 { 297 lusol->indc[nz] = i + 1; 298 lusol->indr[nz] = a->j[rs] + 1; 299 lusol->data[nz] = a->a[rs]; 300 nz++; 301 } 302 rs++; 303 } 304 } 305 306 /*******************************************************************/ 307 /* Do the factorization. */ 308 /*******************************************************************/ 309 310 LU1FAC(&m, &n, &nz, &nnz, 311 lusol->luparm, lusol->parmlu, lusol->data, 312 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 313 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, 314 lusol->iploc, lusol->iqloc, lusol->ipinv, 315 lusol->iqinv, lusol->mnsw, &status); 316 317 switch(status) 318 { 319 case 0: /* factored */ 320 break; 321 322 case 7: /* insufficient memory */ 323 break; 324 325 case 1: 326 case -1: /* singular */ 327 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Singular matrix"); 328 329 case 3: 330 case 4: /* error conditions */ 331 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"matrix error"); 332 333 default: /* unknown condition */ 334 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"matrix unknown return code"); 335 } 336 337 factorizations++; 338 } while (status == 7); 339 F->ops->solve = MatSolve_LUSOL; 340 F->assembled = PETSC_TRUE; 341 F->preallocated = PETSC_TRUE; 342 PetscFunctionReturn(0); 343 } 344 345 #undef __FUNCT__ 346 #define __FUNCT__ "MatLUFactorSymbolic_LUSOL" 347 PetscErrorCode MatLUFactorSymbolic_LUSOL(Mat F,Mat A, IS r, IS c,const MatFactorInfo *info) 348 { 349 /************************************************************************/ 350 /* Input */ 351 /* A - matrix to factor */ 352 /* r - row permutation (ignored) */ 353 /* c - column permutation (ignored) */ 354 /* */ 355 /* Output */ 356 /* F - matrix storing the factorization; */ 357 /************************************************************************/ 358 Mat_LUSOL *lusol; 359 PetscErrorCode ierr; 360 int i, m, n, nz, nnz; 361 362 PetscFunctionBegin; 363 364 /************************************************************************/ 365 /* Check the arguments. */ 366 /************************************************************************/ 367 368 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 369 nz = ((Mat_SeqAIJ *)A->data)->nz; 370 371 /************************************************************************/ 372 /* Create the factorization. */ 373 /************************************************************************/ 374 375 F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 376 lusol = (Mat_LUSOL*)(F->spptr); 377 378 /************************************************************************/ 379 /* Initialize parameters */ 380 /************************************************************************/ 381 382 for (i = 0; i < 30; i++) 383 { 384 lusol->luparm[i] = 0; 385 lusol->parmlu[i] = 0; 386 } 387 388 lusol->luparm[1] = -1; 389 lusol->luparm[2] = 5; 390 lusol->luparm[7] = 1; 391 392 lusol->parmlu[0] = 1 / Factorization_Tolerance; 393 lusol->parmlu[1] = 1 / Factorization_Tolerance; 394 lusol->parmlu[2] = Factorization_Small_Tolerance; 395 lusol->parmlu[3] = Factorization_Pivot_Tolerance; 396 lusol->parmlu[4] = Factorization_Pivot_Tolerance; 397 lusol->parmlu[5] = 3.0; 398 lusol->parmlu[6] = 0.3; 399 lusol->parmlu[7] = 0.6; 400 401 /************************************************************************/ 402 /* Allocate the workspace needed by LUSOL. */ 403 /************************************************************************/ 404 405 lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill); 406 nnz = PetscMax((int)(lusol->elbowroom*nz), 5*n); 407 408 lusol->n = n; 409 lusol->nz = nz; 410 lusol->nnz = nnz; 411 lusol->luroom = 1.75; 412 413 ierr = PetscMalloc(sizeof(int)*n,&lusol->ip); 414 ierr = PetscMalloc(sizeof(int)*n,&lusol->iq); 415 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenc); 416 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenr); 417 ierr = PetscMalloc(sizeof(int)*n,&lusol->locc); 418 ierr = PetscMalloc(sizeof(int)*n,&lusol->locr); 419 ierr = PetscMalloc(sizeof(int)*n,&lusol->iploc); 420 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqloc); 421 ierr = PetscMalloc(sizeof(int)*n,&lusol->ipinv); 422 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqinv); 423 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsw); 424 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsv); 425 426 ierr = PetscMalloc3(nnz,double,&lusol->data,nnz,PetscInt,&lusol->indc,nnz,PetscInt,&lusol->indr);CHKERRQ(ierr); 427 lusol->CleanUpLUSOL = PETSC_TRUE; 428 F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 429 PetscFunctionReturn(0); 430 } 431 432 EXTERN_C_BEGIN 433 #undef __FUNCT__ 434 #define __FUNCT__ "MatFactorGetSolverPackage_seqaij_lusol" 435 PetscErrorCode MatFactorGetSolverPackage_seqaij_lusol(Mat A,const MatSolverPackage *type) 436 { 437 PetscFunctionBegin; 438 *type = MATSOLVERLUSOL; 439 PetscFunctionReturn(0); 440 } 441 EXTERN_C_END 442 443 #undef __FUNCT__ 444 #define __FUNCT__ "MatGetFactor_seqaij_lusol" 445 PetscErrorCode MatGetFactor_seqaij_lusol(Mat A,MatFactorType ftype,Mat *F) 446 { 447 Mat B; 448 Mat_LUSOL *lusol; 449 PetscErrorCode ierr; 450 int m, n; 451 452 PetscFunctionBegin; 453 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 454 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 455 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 456 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 457 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 458 459 ierr = PetscNewLog(B,Mat_LUSOL,&lusol);CHKERRQ(ierr); 460 B->spptr = lusol; 461 462 B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL; 463 B->ops->destroy = MatDestroy_LUSOL; 464 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqaij_lusol",MatFactorGetSolverPackage_seqaij_lusol);CHKERRQ(ierr); 465 B->factortype = MAT_FACTOR_LU; 466 PetscFunctionReturn(0); 467 } 468 469 /*MC 470 MATSOLVERLUSOL - "lusol" - Provides direct solvers (LU) for sequential matrices 471 via the external package LUSOL. 472 473 If LUSOL is installed (see the manual for 474 instructions on how to declare the existence of external packages), 475 476 Works with MATSEQAIJ matrices 477 478 Level: beginner 479 480 .seealso: PCLU, PCFactorSetMatSolverPackage(), MatSolverPackage 481 482 M*/ 483