Actual source code: lusol.c

petsc-master 2015-04-17
Report Typos and Errors
  2: /*
  3:         Provides an interface to the LUSOL package of ....

  5: */
  6: #include <../src/mat/impls/aij/seq/aij.h>

  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

 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: }

 34: PETSC_EXTERN void PETSC_STDCALL M6RDEL()
 35: {
 36:   ;
 37: }

 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);

 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);

 51: extern PetscErrorCode MatDuplicate_LUSOL(Mat,MatDuplicateOption,Mat*);

 53: typedef struct  {
 54:   double *data;
 55:   int    *indc;
 56:   int    *indr;

 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;

 71:   double elbowroom;
 72:   double luroom;                /* Extra space allocated when factor fails   */
 73:   double parmlu[30];            /* Input/output to LUSOL                     */

 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                     */

 80:   PetscBool CleanUpLUSOL;

 82: } Mat_LUSOL;

 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:  */

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) */

182: PetscErrorCode MatDestroy_LUSOL(Mat A)
183: {
185:   Mat_LUSOL      *lusol=(Mat_LUSOL*)A->spptr;

188:   if (lusol && lusol->CleanUpLUSOL) {
189:     PetscFree(lusol->ip);
190:     PetscFree(lusol->iq);
191:     PetscFree(lusol->lenc);
192:     PetscFree(lusol->lenr);
193:     PetscFree(lusol->locc);
194:     PetscFree(lusol->locr);
195:     PetscFree(lusol->iploc);
196:     PetscFree(lusol->iqloc);
197:     PetscFree(lusol->ipinv);
198:     PetscFree(lusol->iqinv);
199:     PetscFree(lusol->mnsw);
200:     PetscFree(lusol->mnsv);
201:     PetscFree3(lusol->data,lusol->indc,lusol->indr);
202:   }
203:   PetscFree(A->spptr);
204:   MatDestroy_SeqAIJ(A);
205:   return(0);
206: }

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;
217:   int            i,m,n,nnz,status;

220:   VecGetArray(x, &xx);
221:   VecGetArrayRead(b, &bb);

223:   m   = n = lusol->n;
224:   nnz = lusol->nnz;

226:   for (i = 0; i < m; i++) lusol->mnsv[i] = bb[i];

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);

233:   if (status) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"solve failed, error code %d",status);

235:   VecRestoreArray(x, &xx);
236:   VecRestoreArrayRead(b, &bb);
237:   return(0);
238: }

242: PetscErrorCode MatLUFactorNumeric_LUSOL(Mat F,Mat A,const MatFactorInfo *info)
243: {
244:   Mat_SeqAIJ     *a;
245:   Mat_LUSOL      *lusol = (Mat_LUSOL*)F->spptr;
247:   int            m, n, nz, nnz, status;
248:   int            i, rs, re;
249:   int            factorizations;

252:   MatGetSize(A,&m,&n);
253:   a    = (Mat_SeqAIJ*)A->data;

255:   if (m != lusol->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"factorization struct inconsistent");

257:   factorizations = 0;
258:   do {
259:     /*******************************************************************/
260:     /* Check the workspace allocation.                                 */
261:     /*******************************************************************/

263:     nz  = a->nz;
264:     nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom*nz));
265:     nnz = PetscMax(nnz, 5*n);

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:     }

273:     nnz = PetscMax(nnz, (int)(lusol->luroom*(lusol->luparm[22] + lusol->luparm[23])));

275:     if (nnz > lusol->nnz) {
276:       PetscFree3(lusol->data,lusol->indc,lusol->indr);
277:       PetscMalloc3(nnz,&lusol->data,nnz,&lusol->indc,nnz,&lusol->indr);
278:       lusol->nnz = nnz;
279:     }

281:     /*******************************************************************/
282:     /* Fill in the data for the problem.      (1-based Fortran style)  */
283:     /*******************************************************************/

285:     nz = 0;
286:     for (i = 0; i < n; i++) {
287:       rs = a->i[i];
288:       re = a->i[i+1];

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:     }

301:     /*******************************************************************/
302:     /* Do the factorization.                                           */
303:     /*******************************************************************/

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);

312:     switch (status) {
313:     case 0:         /* factored */
314:       break;

316:     case 7:         /* insufficient memory */
317:       break;

319:     case 1:
320:     case -1:        /* singular */
321:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Singular matrix");

323:     case 3:
324:     case 4:         /* error conditions */
325:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"matrix error");

327:     default:        /* unknown condition */
328:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"matrix unknown return code");
329:     }

331:     factorizations++;
332:   } while (status == 7);
333:   F->ops->solve   = MatSolve_LUSOL;
334:   F->assembled    = PETSC_TRUE;
335:   F->preallocated = PETSC_TRUE;
336:   return(0);
337: }

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;
354:   int            i, m, n, nz, nnz;

357:   /************************************************************************/
358:   /* Check the arguments.                                                 */
359:   /************************************************************************/

361:   MatGetSize(A, &m, &n);
362:   nz   = ((Mat_SeqAIJ*)A->data)->nz;

364:   /************************************************************************/
365:   /* Create the factorization.                                            */
366:   /************************************************************************/

368:   F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL;
369:   lusol                   = (Mat_LUSOL*)(F->spptr);

371:   /************************************************************************/
372:   /* Initialize parameters                                                */
373:   /************************************************************************/

375:   for (i = 0; i < 30; i++) {
376:     lusol->luparm[i] = 0;
377:     lusol->parmlu[i] = 0;
378:   }

380:   lusol->luparm[1] = -1;
381:   lusol->luparm[2] = 5;
382:   lusol->luparm[7] = 1;

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;

393:   /************************************************************************/
394:   /* Allocate the workspace needed by LUSOL.                              */
395:   /************************************************************************/

397:   lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill);
398:   nnz              = PetscMax((int)(lusol->elbowroom*nz), 5*n);

400:   lusol->n      = n;
401:   lusol->nz     = nz;
402:   lusol->nnz    = nnz;
403:   lusol->luroom = 1.75;

405:   PetscMalloc(sizeof(int)*n,&lusol->ip);
406:   PetscMalloc(sizeof(int)*n,&lusol->iq);
407:   PetscMalloc(sizeof(int)*n,&lusol->lenc);
408:   PetscMalloc(sizeof(int)*n,&lusol->lenr);
409:   PetscMalloc(sizeof(int)*n,&lusol->locc);
410:   PetscMalloc(sizeof(int)*n,&lusol->locr);
411:   PetscMalloc(sizeof(int)*n,&lusol->iploc);
412:   PetscMalloc(sizeof(int)*n,&lusol->iqloc);
413:   PetscMalloc(sizeof(int)*n,&lusol->ipinv);
414:   PetscMalloc(sizeof(int)*n,&lusol->iqinv);
415:   PetscMalloc(sizeof(double)*n,&lusol->mnsw);
416:   PetscMalloc(sizeof(double)*n,&lusol->mnsv);

418:   PetscMalloc3(nnz,&lusol->data,nnz,&lusol->indc,nnz,&lusol->indr);

420:   lusol->CleanUpLUSOL     = PETSC_TRUE;
421:   F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL;
422:   return(0);
423: }

427: PetscErrorCode MatFactorGetSolverPackage_seqaij_lusol(Mat A,const MatSolverPackage *type)
428: {
430:   *type = MATSOLVERLUSOL;
431:   return(0);
432: }

436: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_lusol(Mat A,MatFactorType ftype,Mat *F)
437: {
438:   Mat            B;
439:   Mat_LUSOL      *lusol;
441:   int            m, n;

444:   MatGetSize(A, &m, &n);
445:   MatCreate(PetscObjectComm((PetscObject)A),&B);
446:   MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
447:   MatSetType(B,((PetscObject)A)->type_name);
448:   MatSeqAIJSetPreallocation(B,0,NULL);

450:   PetscNewLog(B,&lusol);
451:   B->spptr = lusol;

453:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL;
454:   B->ops->destroy          = MatDestroy_LUSOL;

456:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqaij_lusol);

458:   B->factortype = MAT_FACTOR_LU;
459:   return(0);
460: }

464: PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_Lusol(void)
465: {

469:   MatSolverPackageRegister(MATSOLVERLUSOL,MATSEQAIJ,        MAT_FACTOR_LU,MatGetFactor_seqaij_lusol);
470:   return(0);
471: }

473: /*MC
474:   MATSOLVERLUSOL - "lusol" - Provides direct solvers (LU) for sequential matrices
475:                          via the external package LUSOL.

477:   If LUSOL is installed (see the manual for
478:   instructions on how to declare the existence of external packages),

480:   Works with MATSEQAIJ matrices

482:    Level: beginner

484: .seealso: PCLU, PCFactorSetMatSolverPackage(), MatSolverPackage

486: M*/