1 #define PETSCMAT_DLL 2 3 /* 4 Provides an interface to the LUSOL package of .... 5 6 */ 7 #include "src/mat/impls/aij/seq/aij.h" 8 9 #if defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 10 #define LU1FAC lu1fac_ 11 #define LU6SOL lu6sol_ 12 #define M1PAGE m1page_ 13 #define M5SETX m5setx_ 14 #define M6RDEL m6rdel_ 15 #elif !defined(PETSC_HAVE_FORTRAN_CAPS) 16 #define LU1FAC lu1fac 17 #define LU6SOL lu6sol 18 #define M1PAGE m1page 19 #define M5SETX m5setx 20 #define M6RDEL m6rdel 21 #endif 22 23 EXTERN_C_BEGIN 24 /* 25 Dummy symbols that the MINOS files mi25bfac.f and mi15blas.f may require 26 */ 27 void PETSC_STDCALL M1PAGE() { 28 ; 29 } 30 void PETSC_STDCALL M5SETX() { 31 ; 32 } 33 34 void PETSC_STDCALL M6RDEL() { 35 ; 36 } 37 38 extern void PETSC_STDCALL LU1FAC (int *m, int *n, int *nnz, int *size, int *luparm, 39 double *parmlu, double *data, int *indc, int *indr, 40 int *rowperm, int *colperm, int *collen, int *rowlen, 41 int *colstart, int *rowstart, int *rploc, int *cploc, 42 int *rpinv, int *cpinv, double *w, int *inform); 43 44 extern void PETSC_STDCALL LU6SOL (int *mode, int *m, int *n, double *rhs, double *x, 45 int *size, int *luparm, double *parmlu, double *data, 46 int *indc, int *indr, int *rowperm, int *colperm, 47 int *collen, int *rowlen, int *colstart, int *rowstart, 48 int *inform); 49 EXTERN_C_END 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 PetscTruth 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->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 = PetscFree(lusol->indc);CHKERRQ(ierr); 202 } 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 != 0) SETERRQ1(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 A,MatFactorInfo *info,Mat *F) 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_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 = PetscFree(lusol->indc);CHKERRQ(ierr); 279 ierr = PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc);CHKERRQ(ierr); 280 lusol->indr = lusol->indc + nnz; 281 lusol->data = (double *)(lusol->indr + nnz); 282 lusol->nnz = nnz; 283 } 284 285 /*******************************************************************/ 286 /* Fill in the data for the problem. (1-based Fortran style) */ 287 /*******************************************************************/ 288 289 nz = 0; 290 for (i = 0; i < n; i++) 291 { 292 rs = a->i[i]; 293 re = a->i[i+1]; 294 295 while (rs < re) 296 { 297 if (a->a[rs] != 0.0) 298 { 299 lusol->indc[nz] = i + 1; 300 lusol->indr[nz] = a->j[rs] + 1; 301 lusol->data[nz] = a->a[rs]; 302 nz++; 303 } 304 rs++; 305 } 306 } 307 308 /*******************************************************************/ 309 /* Do the factorization. */ 310 /*******************************************************************/ 311 312 LU1FAC(&m, &n, &nz, &nnz, 313 lusol->luparm, lusol->parmlu, lusol->data, 314 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 315 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, 316 lusol->iploc, lusol->iqloc, lusol->ipinv, 317 lusol->iqinv, lusol->mnsw, &status); 318 319 switch(status) 320 { 321 case 0: /* factored */ 322 break; 323 324 case 7: /* insufficient memory */ 325 break; 326 327 case 1: 328 case -1: /* singular */ 329 SETERRQ(PETSC_ERR_LIB,"Singular matrix"); 330 331 case 3: 332 case 4: /* error conditions */ 333 SETERRQ(PETSC_ERR_LIB,"matrix error"); 334 335 default: /* unknown condition */ 336 SETERRQ(PETSC_ERR_LIB,"matrix unknown return code"); 337 } 338 339 factorizations++; 340 } while (status == 7); 341 (*F)->assembled = PETSC_TRUE; 342 (*F)->preallocated = PETSC_TRUE; 343 PetscFunctionReturn(0); 344 } 345 346 #undef __FUNCT__ 347 #define __FUNCT__ "MatLUFactorSymbolic_LUSOL" 348 PetscErrorCode MatLUFactorSymbolic_LUSOL(Mat A, IS r, IS c,MatFactorInfo *info, Mat *F) 349 { 350 /************************************************************************/ 351 /* Input */ 352 /* A - matrix to factor */ 353 /* r - row permutation (ignored) */ 354 /* c - column permutation (ignored) */ 355 /* */ 356 /* Output */ 357 /* F - matrix storing the factorization; */ 358 /************************************************************************/ 359 Mat B; 360 Mat_LUSOL *lusol; 361 PetscErrorCode ierr; 362 int i, m, n, nz, nnz; 363 364 PetscFunctionBegin; 365 366 /************************************************************************/ 367 /* Check the arguments. */ 368 /************************************************************************/ 369 370 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 371 nz = ((Mat_SeqAIJ *)A->data)->nz; 372 373 /************************************************************************/ 374 /* Create the factorization. */ 375 /************************************************************************/ 376 377 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 378 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 379 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 380 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 381 382 B->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 383 B->ops->solve = MatSolve_LUSOL; 384 B->factor = MAT_FACTOR_LU; 385 lusol = (Mat_LUSOL*)(B->spptr); 386 387 /************************************************************************/ 388 /* Initialize parameters */ 389 /************************************************************************/ 390 391 for (i = 0; i < 30; i++) 392 { 393 lusol->luparm[i] = 0; 394 lusol->parmlu[i] = 0; 395 } 396 397 lusol->luparm[1] = -1; 398 lusol->luparm[2] = 5; 399 lusol->luparm[7] = 1; 400 401 lusol->parmlu[0] = 1 / Factorization_Tolerance; 402 lusol->parmlu[1] = 1 / Factorization_Tolerance; 403 lusol->parmlu[2] = Factorization_Small_Tolerance; 404 lusol->parmlu[3] = Factorization_Pivot_Tolerance; 405 lusol->parmlu[4] = Factorization_Pivot_Tolerance; 406 lusol->parmlu[5] = 3.0; 407 lusol->parmlu[6] = 0.3; 408 lusol->parmlu[7] = 0.6; 409 410 /************************************************************************/ 411 /* Allocate the workspace needed by LUSOL. */ 412 /************************************************************************/ 413 414 lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill); 415 nnz = PetscMax((int)(lusol->elbowroom*nz), 5*n); 416 417 lusol->n = n; 418 lusol->nz = nz; 419 lusol->nnz = nnz; 420 lusol->luroom = 1.75; 421 422 ierr = PetscMalloc(sizeof(int)*n,&lusol->ip); 423 ierr = PetscMalloc(sizeof(int)*n,&lusol->iq); 424 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenc); 425 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenr); 426 ierr = PetscMalloc(sizeof(int)*n,&lusol->locc); 427 ierr = PetscMalloc(sizeof(int)*n,&lusol->locr); 428 ierr = PetscMalloc(sizeof(int)*n,&lusol->iploc); 429 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqloc); 430 ierr = PetscMalloc(sizeof(int)*n,&lusol->ipinv); 431 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqinv); 432 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsw); 433 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsv); 434 435 ierr = PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc); 436 lusol->indr = lusol->indc + nnz; 437 lusol->data = (double *)(lusol->indr + nnz); 438 lusol->CleanUpLUSOL = PETSC_TRUE; 439 B->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 440 B->ops->solve = MatSolve_LUSOL; 441 *F = B; 442 PetscFunctionReturn(0); 443 } 444 445 #undef __FUNCT__ 446 #define __FUNCT__ "MatGetFactor_seqaij_lusol" 447 PetscErrorCode MatGetFactor_seqaij_lusol(Mat A,MatFactorType ftype,Mat *F) 448 { 449 Mat B; 450 Mat_LUSOL *lusol; 451 PetscErrorCode ierr; 452 int i, m, n, nz, nnz; 453 454 PetscFunctionBegin; 455 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 456 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 457 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 458 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 459 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 460 461 ierr = PetscNewLog(B,Mat_LUSOL,&lusol);CHKERRQ(ierr); 462 B->spptr = lusol; 463 464 B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL; 465 B->ops->destroy = MatDestroy_LUSOL; 466 B->factor = MAT_FACTOR_LU; 467 PetscFunctionReturn(0); 468 } 469 470 /*MC 471 MATLUSOL - MATLUSOL = "lusol" - A matrix type providing direct solvers (LU) for sequential matrices 472 via the external package LUSOL. 473 474 If LUSOL is installed (see the manual for 475 instructions on how to declare the existence of external packages), 476 a matrix type can be constructed which invokes LUSOL solvers. 477 After calling MatCreate(...,A), simply call MatSetType(A,MATLUSOL). 478 This matrix type is only supported for double precision real. 479 480 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 481 supported for this matrix type. MatConvert can be called for a fast inplace conversion 482 to and from the MATSEQAIJ matrix type. 483 484 Options Database Keys: 485 . -mat_type lusol - sets the matrix type to "lusol" during a call to MatSetFromOptions() 486 487 Level: beginner 488 489 .seealso: PCLU 490 M*/ 491