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++) lusol->mnsv[i] = bb[i]; 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 /* Check the workspace allocation. */ 262 /*******************************************************************/ 263 264 nz = a->nz; 265 nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom*nz)); 266 nnz = PetscMax(nnz, 5*n); 267 268 if (nnz < lusol->luparm[12]) { 269 nnz = (int)(lusol->luroom * lusol->luparm[12]); 270 } else if ((factorizations > 0) && (lusol->luroom < 6)) { 271 lusol->luroom += 0.1; 272 } 273 274 nnz = PetscMax(nnz, (int)(lusol->luroom*(lusol->luparm[22] + lusol->luparm[23]))); 275 276 if (nnz > lusol->nnz) { 277 ierr = PetscFree3(lusol->data,lusol->indc,lusol->indr);CHKERRQ(ierr); 278 ierr = PetscMalloc3(nnz,double,&lusol->data,nnz,PetscInt,&lusol->indc,nnz,PetscInt,&lusol->indr);CHKERRQ(ierr); 279 lusol->nnz = nnz; 280 } 281 282 /*******************************************************************/ 283 /* Fill in the data for the problem. (1-based Fortran style) */ 284 /*******************************************************************/ 285 286 nz = 0; 287 for (i = 0; i < n; i++) { 288 rs = a->i[i]; 289 re = a->i[i+1]; 290 291 while (rs < re) { 292 if (a->a[rs] != 0.0) { 293 lusol->indc[nz] = i + 1; 294 lusol->indr[nz] = a->j[rs] + 1; 295 lusol->data[nz] = a->a[rs]; 296 nz++; 297 } 298 rs++; 299 } 300 } 301 302 /*******************************************************************/ 303 /* Do the factorization. */ 304 /*******************************************************************/ 305 306 LU1FAC(&m, &n, &nz, &nnz, 307 lusol->luparm, lusol->parmlu, lusol->data, 308 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 309 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, 310 lusol->iploc, lusol->iqloc, lusol->ipinv, 311 lusol->iqinv, lusol->mnsw, &status); 312 313 switch (status) { 314 case 0: /* factored */ 315 break; 316 317 case 7: /* insufficient memory */ 318 break; 319 320 case 1: 321 case -1: /* singular */ 322 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Singular matrix"); 323 324 case 3: 325 case 4: /* error conditions */ 326 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"matrix error"); 327 328 default: /* unknown condition */ 329 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"matrix unknown return code"); 330 } 331 332 factorizations++; 333 } while (status == 7); 334 F->ops->solve = MatSolve_LUSOL; 335 F->assembled = PETSC_TRUE; 336 F->preallocated = PETSC_TRUE; 337 PetscFunctionReturn(0); 338 } 339 340 #undef __FUNCT__ 341 #define __FUNCT__ "MatLUFactorSymbolic_LUSOL" 342 PetscErrorCode MatLUFactorSymbolic_LUSOL(Mat F,Mat A, IS r, IS c,const MatFactorInfo *info) 343 { 344 /************************************************************************/ 345 /* Input */ 346 /* A - matrix to factor */ 347 /* r - row permutation (ignored) */ 348 /* c - column permutation (ignored) */ 349 /* */ 350 /* Output */ 351 /* F - matrix storing the factorization; */ 352 /************************************************************************/ 353 Mat_LUSOL *lusol; 354 PetscErrorCode ierr; 355 int i, m, n, nz, nnz; 356 357 PetscFunctionBegin; 358 /************************************************************************/ 359 /* Check the arguments. */ 360 /************************************************************************/ 361 362 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 363 nz = ((Mat_SeqAIJ*)A->data)->nz; 364 365 /************************************************************************/ 366 /* Create the factorization. */ 367 /************************************************************************/ 368 369 F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 370 lusol = (Mat_LUSOL*)(F->spptr); 371 372 /************************************************************************/ 373 /* Initialize parameters */ 374 /************************************************************************/ 375 376 for (i = 0; i < 30; i++) { 377 lusol->luparm[i] = 0; 378 lusol->parmlu[i] = 0; 379 } 380 381 lusol->luparm[1] = -1; 382 lusol->luparm[2] = 5; 383 lusol->luparm[7] = 1; 384 385 lusol->parmlu[0] = 1 / Factorization_Tolerance; 386 lusol->parmlu[1] = 1 / Factorization_Tolerance; 387 lusol->parmlu[2] = Factorization_Small_Tolerance; 388 lusol->parmlu[3] = Factorization_Pivot_Tolerance; 389 lusol->parmlu[4] = Factorization_Pivot_Tolerance; 390 lusol->parmlu[5] = 3.0; 391 lusol->parmlu[6] = 0.3; 392 lusol->parmlu[7] = 0.6; 393 394 /************************************************************************/ 395 /* Allocate the workspace needed by LUSOL. */ 396 /************************************************************************/ 397 398 lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill); 399 nnz = PetscMax((int)(lusol->elbowroom*nz), 5*n); 400 401 lusol->n = n; 402 lusol->nz = nz; 403 lusol->nnz = nnz; 404 lusol->luroom = 1.75; 405 406 ierr = PetscMalloc(sizeof(int)*n,&lusol->ip); 407 ierr = PetscMalloc(sizeof(int)*n,&lusol->iq); 408 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenc); 409 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenr); 410 ierr = PetscMalloc(sizeof(int)*n,&lusol->locc); 411 ierr = PetscMalloc(sizeof(int)*n,&lusol->locr); 412 ierr = PetscMalloc(sizeof(int)*n,&lusol->iploc); 413 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqloc); 414 ierr = PetscMalloc(sizeof(int)*n,&lusol->ipinv); 415 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqinv); 416 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsw); 417 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsv); 418 419 ierr = PetscMalloc3(nnz,double,&lusol->data,nnz,PetscInt,&lusol->indc,nnz,PetscInt,&lusol->indr);CHKERRQ(ierr); 420 421 lusol->CleanUpLUSOL = PETSC_TRUE; 422 F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 423 PetscFunctionReturn(0); 424 } 425 426 EXTERN_C_BEGIN 427 #undef __FUNCT__ 428 #define __FUNCT__ "MatFactorGetSolverPackage_seqaij_lusol" 429 PetscErrorCode MatFactorGetSolverPackage_seqaij_lusol(Mat A,const MatSolverPackage *type) 430 { 431 PetscFunctionBegin; 432 *type = MATSOLVERLUSOL; 433 PetscFunctionReturn(0); 434 } 435 EXTERN_C_END 436 437 #undef __FUNCT__ 438 #define __FUNCT__ "MatGetFactor_seqaij_lusol" 439 PetscErrorCode MatGetFactor_seqaij_lusol(Mat A,MatFactorType ftype,Mat *F) 440 { 441 Mat B; 442 Mat_LUSOL *lusol; 443 PetscErrorCode ierr; 444 int m, n; 445 446 PetscFunctionBegin; 447 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 448 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 449 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 450 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 451 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 452 453 ierr = PetscNewLog(B,Mat_LUSOL,&lusol);CHKERRQ(ierr); 454 B->spptr = lusol; 455 456 B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL; 457 B->ops->destroy = MatDestroy_LUSOL; 458 459 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqaij_lusol",MatFactorGetSolverPackage_seqaij_lusol);CHKERRQ(ierr); 460 461 B->factortype = MAT_FACTOR_LU; 462 PetscFunctionReturn(0); 463 } 464 465 /*MC 466 MATSOLVERLUSOL - "lusol" - Provides direct solvers (LU) for sequential matrices 467 via the external package LUSOL. 468 469 If LUSOL is installed (see the manual for 470 instructions on how to declare the existence of external packages), 471 472 Works with MATSEQAIJ matrices 473 474 Level: beginner 475 476 .seealso: PCLU, PCFactorSetMatSolverPackage(), MatSolverPackage 477 478 M*/ 479