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