Actual source code: gqt.c

      1: #include <petscsys.h>
      2: #include <petscblaslapack.h>

      4: static PetscErrorCode estsv(PetscInt n, PetscReal *r, PetscInt ldr, PetscReal *svmin, PetscReal *z)
      5: {
      6:   PetscBLASInt   blas1=1, blasn, blasnmi, blasj, blasldr;
      7:   PetscInt       i,j;
      8:   PetscReal      e,temp,w,wm,ynorm,znorm,s,sm;

     12:   PetscBLASIntCast(n,&blasn);
     13:   PetscBLASIntCast(ldr,&blasldr);
     14:   for (i=0;i<n;i++) {
     15:     z[i]=0.0;
     16:   }
     17:   e = PetscAbs(r[0]);
     18:   if (e == 0.0) {
     19:     *svmin = 0.0;
     20:     z[0] = 1.0;
     21:   } else {
     22:     /* Solve R'*y = e */
     23:     for (i=0;i<n;i++) {
     24:       /* Scale y. The scaling factor (0.01) reduces the number of scalings */
     25:       if (z[i] >= 0.0) e =-PetscAbs(e);
     26:       else             e = PetscAbs(e);

     28:       if (PetscAbs(e - z[i]) > PetscAbs(r[i + ldr*i])) {
     29:         temp = PetscMin(0.01,PetscAbs(r[i + ldr*i]))/PetscAbs(e-z[i]);
     30:         PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &temp, z, &blas1));
     31:         e = temp*e;
     32:       }

     34:       /* Determine the two possible choices of y[i] */
     35:       if (r[i + ldr*i] == 0.0) {
     36:         w = wm = 1.0;
     37:       } else {
     38:         w  = (e - z[i]) / r[i + ldr*i];
     39:         wm = - (e + z[i]) / r[i + ldr*i];
     40:       }

     42:       /*  Chose y[i] based on the predicted value of y[j] for j>i */
     43:       s  = PetscAbs(e - z[i]);
     44:       sm = PetscAbs(e + z[i]);
     45:       for (j=i+1;j<n;j++) {
     46:         sm += PetscAbs(z[j] + wm * r[i + ldr*j]);
     47:       }
     48:       if (i < n-1) {
     49:         PetscBLASIntCast(n-i-1,&blasnmi);
     50:         PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasnmi, &w, &r[i + ldr*(i+1)], &blasldr, &z[i+1], &blas1));
     51:         PetscStackCallBLAS("BLASasum",s += BLASasum_(&blasnmi, &z[i+1], &blas1));
     52:       }
     53:       if (s < sm) {
     54:         temp = wm - w;
     55:         w = wm;
     56:         if (i < n-1) {
     57:           PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasnmi, &temp, &r[i + ldr*(i+1)], &blasldr, &z[i+1], &blas1));
     58:         }
     59:       }
     60:       z[i] = w;
     61:     }

     63:     PetscStackCallBLAS("BLASnrm2",ynorm = BLASnrm2_(&blasn, z, &blas1));

     65:     /* Solve R*z = y */
     66:     for (j=n-1; j>=0; j--) {
     67:       /* Scale z */
     68:       if (PetscAbs(z[j]) > PetscAbs(r[j + ldr*j])) {
     69:         temp = PetscMin(0.01, PetscAbs(r[j + ldr*j] / z[j]));
     70:         PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &temp, z, &blas1));
     71:         ynorm *=temp;
     72:       }
     73:       if (r[j + ldr*j] == 0) {
     74:         z[j] = 1.0;
     75:       } else {
     76:         z[j] = z[j] / r[j + ldr*j];
     77:       }
     78:       temp = -z[j];
     79:       PetscBLASIntCast(j,&blasj);
     80:       PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasj,&temp,&r[0+ldr*j],&blas1,z,&blas1));
     81:     }

     83:     /* Compute svmin and normalize z */
     84:     PetscStackCallBLAS("BLASnrm2",znorm = 1.0 / BLASnrm2_(&blasn, z, &blas1));
     85:     *svmin = ynorm*znorm;
     86:     PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &znorm, z, &blas1));
     87:   }
     88:   return(0);
     89: }

     91: /*
     92: c     ***********
     93: c
     94: c     Subroutine gqt
     95: c
     96: c     Given an n by n symmetric matrix A, an n-vector b, and a
     97: c     positive number delta, this subroutine determines a vector
     98: c     x which approximately minimizes the quadratic function
     99: c
    100: c           f(x) = (1/2)*x'*A*x + b'*x
    101: c
    102: c     subject to the Euclidean norm constraint
    103: c
    104: c           norm(x) <= delta.
    105: c
    106: c     This subroutine computes an approximation x and a Lagrange
    107: c     multiplier par such that either par is zero and
    108: c
    109: c            norm(x) <= (1+rtol)*delta,
    110: c
    111: c     or par is positive and
    112: c
    113: c            abs(norm(x) - delta) <= rtol*delta.
    114: c
    115: c     If xsol is the solution to the problem, the approximation x
    116: c     satisfies
    117: c
    118: c            f(x) <= ((1 - rtol)**2)*f(xsol)
    119: c
    120: c     The subroutine statement is
    121: c
    122: c       subroutine gqt(n,a,lda,b,delta,rtol,atol,itmax,
    123: c                        par,f,x,info,z,wa1,wa2)
    124: c
    125: c     where
    126: c
    127: c       n is an integer variable.
    128: c         On entry n is the order of A.
    129: c         On exit n is unchanged.
    130: c
    131: c       a is a double precision array of dimension (lda,n).
    132: c         On entry the full upper triangle of a must contain the
    133: c            full upper triangle of the symmetric matrix A.
    134: c         On exit the array contains the matrix A.
    135: c
    136: c       lda is an integer variable.
    137: c         On entry lda is the leading dimension of the array a.
    138: c         On exit lda is unchanged.
    139: c
    140: c       b is an double precision array of dimension n.
    141: c         On entry b specifies the linear term in the quadratic.
    142: c         On exit b is unchanged.
    143: c
    144: c       delta is a double precision variable.
    145: c         On entry delta is a bound on the Euclidean norm of x.
    146: c         On exit delta is unchanged.
    147: c
    148: c       rtol is a double precision variable.
    149: c         On entry rtol is the relative accuracy desired in the
    150: c            solution. Convergence occurs if
    151: c
    152: c              f(x) <= ((1 - rtol)**2)*f(xsol)
    153: c
    154: c         On exit rtol is unchanged.
    155: c
    156: c       atol is a double precision variable.
    157: c         On entry atol is the absolute accuracy desired in the
    158: c            solution. Convergence occurs when
    159: c
    160: c              norm(x) <= (1 + rtol)*delta
    161: c
    162: c              max(-f(x),-f(xsol)) <= atol
    163: c
    164: c         On exit atol is unchanged.
    165: c
    166: c       itmax is an integer variable.
    167: c         On entry itmax specifies the maximum number of iterations.
    168: c         On exit itmax is unchanged.
    169: c
    170: c       par is a double precision variable.
    171: c         On entry par is an initial estimate of the Lagrange
    172: c            multiplier for the constraint norm(x) <= delta.
    173: c         On exit par contains the final estimate of the multiplier.
    174: c
    175: c       f is a double precision variable.
    176: c         On entry f need not be specified.
    177: c         On exit f is set to f(x) at the output x.
    178: c
    179: c       x is a double precision array of dimension n.
    180: c         On entry x need not be specified.
    181: c         On exit x is set to the final estimate of the solution.
    182: c
    183: c       info is an integer variable.
    184: c         On entry info need not be specified.
    185: c         On exit info is set as follows:
    186: c
    187: c            info = 1  The function value f(x) has the relative
    188: c                      accuracy specified by rtol.
    189: c
    190: c            info = 2  The function value f(x) has the absolute
    191: c                      accuracy specified by atol.
    192: c
    193: c            info = 3  Rounding errors prevent further progress.
    194: c                      On exit x is the best available approximation.
    195: c
    196: c            info = 4  Failure to converge after itmax iterations.
    197: c                      On exit x is the best available approximation.
    198: c
    199: c       z is a double precision work array of dimension n.
    200: c
    201: c       wa1 is a double precision work array of dimension n.
    202: c
    203: c       wa2 is a double precision work array of dimension n.
    204: c
    205: c     Subprograms called
    206: c
    207: c       MINPACK-2  ......  destsv
    208: c
    209: c       LAPACK  .........  dpotrf
    210: c
    211: c       Level 1 BLAS  ...  daxpy, dcopy, ddot, dnrm2, dscal
    212: c
    213: c       Level 2 BLAS  ...  dtrmv, dtrsv
    214: c
    215: c     MINPACK-2 Project. October 1993.
    216: c     Argonne National Laboratory and University of Minnesota.
    217: c     Brett M. Averick, Richard Carter, and Jorge J. More'
    218: c
    219: c     ***********
    220: */
    221: PetscErrorCode gqt(PetscInt n, PetscReal *a, PetscInt lda, PetscReal *b,
    222:                    PetscReal delta, PetscReal rtol, PetscReal atol,
    223:                    PetscInt itmax, PetscReal *retpar, PetscReal *retf,
    224:                    PetscReal *x, PetscInt *retinfo, PetscInt *retits,
    225:                    PetscReal *z, PetscReal *wa1, PetscReal *wa2)
    226: {
    228:   PetscReal      f=0.0,p001=0.001,p5=0.5,minusone=-1,delta2=delta*delta;
    229:   PetscInt       iter, j, rednc,info;
    230:   PetscBLASInt   indef;
    231:   PetscBLASInt   blas1=1, blasn, iblas, blaslda, blasldap1, blasinfo;
    232:   PetscReal      alpha, anorm, bnorm, parc, parf, parl, pars, par=*retpar,paru, prod, rxnorm, rznorm=0.0, temp, xnorm;

    235:   PetscBLASIntCast(n,&blasn);
    236:   PetscBLASIntCast(lda,&blaslda);
    237:   PetscBLASIntCast(lda+1,&blasldap1);
    238:   parf   = 0.0;
    239:   xnorm  = 0.0;
    240:   rxnorm = 0.0;
    241:   rednc  = 0;
    242:   for (j=0; j<n; j++) {
    243:     x[j] = 0.0;
    244:     z[j] = 0.0;
    245:   }

    247:   /* Copy the diagonal and save A in its lower triangle */
    248:   PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn,a,&blasldap1, wa1, &blas1));
    249:   for (j=0;j<n-1;j++) {
    250:     PetscBLASIntCast(n - j - 1,&iblas);
    251:     PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[j + lda*(j+1)], &blaslda, &a[j+1 + lda*j], &blas1));
    252:   }

    254:   /* Calculate the l1-norm of A, the Gershgorin row sums, and the
    255:    l2-norm of b */
    256:   anorm = 0.0;
    257:   for (j=0;j<n;j++) {
    258:     PetscStackCallBLAS("BLASasum",wa2[j] = BLASasum_(&blasn, &a[0 + lda*j], &blas1));CHKMEMQ;
    259:     anorm = PetscMax(anorm,wa2[j]);
    260:   }
    261:   for (j=0;j<n;j++) {
    262:     wa2[j] = wa2[j] - PetscAbs(wa1[j]);
    263:   }
    264:   PetscStackCallBLAS("BLASnrm2",bnorm = BLASnrm2_(&blasn,b,&blas1));CHKMEMQ;
    265:   /* Calculate a lower bound, pars, for the domain of the problem.
    266:    Also calculate an upper bound, paru, and a lower bound, parl,
    267:    for the Lagrange multiplier. */
    268:   pars = parl = paru = -anorm;
    269:   for (j=0;j<n;j++) {
    270:     pars = PetscMax(pars, -wa1[j]);
    271:     parl = PetscMax(parl, wa1[j] + wa2[j]);
    272:     paru = PetscMax(paru, -wa1[j] + wa2[j]);
    273:   }
    274:   parl = PetscMax(bnorm/delta - parl,pars);
    275:   parl = PetscMax(0.0,parl);
    276:   paru = PetscMax(0.0, bnorm/delta + paru);

    278:   /* If the input par lies outside of the interval (parl, paru),
    279:    set par to the closer endpoint. */

    281:   par = PetscMax(par,parl);
    282:   par = PetscMin(par,paru);

    284:   /* Special case: parl == paru */
    285:   paru = PetscMax(paru, (1.0 + rtol)*parl);

    287:   /* Beginning of an iteration */

    289:   info = 0;
    290:   for (iter=1;iter<=itmax;iter++) {
    291:     /* Safeguard par */
    292:     if (par <= pars && paru > 0) {
    293:       par = PetscMax(p001, PetscSqrtScalar(parl/paru)) * paru;
    294:     }

    296:     /* Copy the lower triangle of A into its upper triangle and  compute A + par*I */

    298:     for (j=0;j<n-1;j++) {
    299:       PetscBLASIntCast(n - j - 1,&iblas);
    300:       PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[j+1 + j*lda], &blas1,&a[j + (j+1)*lda], &blaslda));
    301:     }
    302:     for (j=0;j<n;j++) {
    303:       a[j + j*lda] = wa1[j] + par;
    304:     }

    306:     /* Attempt the Cholesky factorization of A without referencing the lower triangular part. */
    307:     PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("U",&blasn,a,&blaslda,&indef));

    309:     /* Case 1: A + par*I is pos. def. */
    310:     if (indef == 0) {

    312:       /* Compute an approximate solution x and save the last value of par with A + par*I pos. def. */

    314:       parf = par;
    315:       PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn, b, &blas1, wa2, &blas1));
    316:       PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&blasn,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
    317:       if (blasinfo) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKtrtrs() returned info %d",blasinfo);
    318:       PetscStackCallBLAS("BLASnrm2",rxnorm = BLASnrm2_(&blasn, wa2, &blas1));
    319:       PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","N","N",&blasn,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
    320:       if (blasinfo) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKtrtrs() returned info %d",blasinfo);

    322:       PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn, wa2, &blas1, x, &blas1));
    323:       PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &minusone, x, &blas1));
    324:       PetscStackCallBLAS("BLASnrm2",xnorm = BLASnrm2_(&blasn, x, &blas1));CHKMEMQ;

    326:       /* Test for convergence */
    327:       if (PetscAbs(xnorm - delta) <= rtol*delta || (par == 0  && xnorm <= (1.0+rtol)*delta)) {
    328:         info = 1;
    329:       }

    331:       /* Compute a direction of negative curvature and use this information to improve pars. */
    332:       estsv(n,a,lda,&rznorm,z);CHKMEMQ;
    333:       pars = PetscMax(pars, par-rznorm*rznorm);

    335:       /* Compute a negative curvature solution of the form x + alpha*z,  where norm(x+alpha*z)==delta */

    337:       rednc = 0;
    338:       if (xnorm < delta) {
    339:         /* Compute alpha */
    340:         PetscStackCallBLAS("BLASdot",prod = BLASdot_(&blasn, z, &blas1, x, &blas1)/delta);
    341:         temp = (delta - xnorm)*((delta + xnorm)/delta);
    342:         alpha = temp/(PetscAbs(prod) + PetscSqrtScalar(prod*prod + temp/delta));
    343:         if (prod >= 0) alpha = PetscAbs(alpha);
    344:         else alpha =-PetscAbs(alpha);

    346:         /* Test to decide if the negative curvature step produces a larger reduction than with z=0 */
    347:         rznorm = PetscAbs(alpha) * rznorm;
    348:         if ((rznorm*rznorm + par*xnorm*xnorm)/(delta2) <= par) {
    349:           rednc = 1;
    350:         }
    351:         /* Test for convergence */
    352:         if (p5 * rznorm*rznorm / delta2 <= rtol*(1.0-p5*rtol)*(par + rxnorm*rxnorm/delta2)) {
    353:           info = 1;
    354:         } else if (info == 0 && (p5*(par + rxnorm*rxnorm/delta2) <= atol/delta2)) {
    355:           info = 2;
    356:         }
    357:       }

    359:       /* Compute the Newton correction parc to par. */
    360:       if (xnorm == 0) {
    361:         parc = -par;
    362:       } else {
    363:         PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn, x, &blas1, wa2, &blas1));
    364:         temp = 1.0/xnorm;
    365:         PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &temp, wa2, &blas1));
    366:         PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&blasn, &blas1, a, &blaslda, wa2, &blasn, &blasinfo));
    367:         if (blasinfo) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKtrtrs() returned info %d",blasinfo);
    368:         PetscStackCallBLAS("BLASnrm2",temp = BLASnrm2_(&blasn, wa2, &blas1));
    369:         parc = (xnorm - delta)/(delta*temp*temp);
    370:       }

    372:       /* update parl or paru */
    373:       if (xnorm > delta) {
    374:         parl = PetscMax(parl, par);
    375:       } else if (xnorm < delta) {
    376:         paru = PetscMin(paru, par);
    377:       }
    378:     } else {
    379:       /* Case 2: A + par*I is not pos. def. */

    381:       /* Use the rank information from the Cholesky decomposition to update par. */

    383:       if (indef > 1) {
    384:         /* Restore column indef to A + par*I. */
    385:         iblas = indef - 1;
    386:         PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[indef-1 + 0*lda],&blaslda,&a[0 + (indef-1)*lda],&blas1));
    387:         a[indef-1 + (indef-1)*lda] = wa1[indef-1] + par;

    389:         /* compute parc. */
    390:         PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[0 + (indef-1)*lda], &blas1, wa2, &blas1));
    391:         PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&iblas,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
    392:         if (blasinfo) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKtrtrs() returned info %d",blasinfo);
    393:         PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,wa2,&blas1,&a[0 + (indef-1)*lda],&blas1));
    394:         PetscStackCallBLAS("BLASnrm2",temp = BLASnrm2_(&iblas,&a[0 + (indef-1)*lda],&blas1));CHKMEMQ;
    395:         a[indef-1 + (indef-1)*lda] -= temp*temp;
    396:         PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","N","N",&iblas,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
    397:         if (blasinfo) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKtrtrs() returned info %d",blasinfo);
    398:       }

    400:       wa2[indef-1] = -1.0;
    401:       iblas = indef;
    402:       PetscStackCallBLAS("BLASnrm2",temp = BLASnrm2_(&iblas,wa2,&blas1));
    403:       parc = - a[indef-1 + (indef-1)*lda]/(temp*temp);
    404:       pars = PetscMax(pars,par+parc);

    406:       /* If necessary, increase paru slightly.
    407:        This is needed because in some exceptional situations
    408:        paru is the optimal value of par. */

    410:       paru = PetscMax(paru, (1.0+rtol)*pars);
    411:     }

    413:     /* Use pars to update parl */
    414:     parl = PetscMax(parl,pars);

    416:     /* Test for converged. */
    417:     if (info == 0) {
    418:       if (iter == itmax) info=4;
    419:       if (paru <= (1.0+p5*rtol)*pars) info=3;
    420:       if (paru == 0.0) info = 2;
    421:     }

    423:     /* If exiting, store the best approximation and restore
    424:      the upper triangle of A. */

    426:     if (info != 0) {
    427:       /* Compute the best current estimates for x and f. */
    428:       par = parf;
    429:       f = -p5 * (rxnorm*rxnorm + par*xnorm*xnorm);
    430:       if (rednc) {
    431:         f = -p5 * (rxnorm*rxnorm + par*delta*delta - rznorm*rznorm);
    432:         PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasn, &alpha, z, &blas1, x, &blas1));
    433:       }
    434:       /* Restore the upper triangle of A */
    435:       for (j = 0; j<n; j++) {
    436:         PetscBLASIntCast(n - j - 1,&iblas);
    437:         PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[j+1 + j*lda],&blas1, &a[j + (j+1)*lda],&blaslda));
    438:       }
    439:       PetscBLASIntCast(lda+1,&iblas);
    440:       PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn,wa1,&blas1,a,&iblas));
    441:       break;
    442:     }
    443:     par = PetscMax(parl,par+parc);
    444:   }
    445:   *retpar  = par;
    446:   *retf    = f;
    447:   *retinfo = info;
    448:   *retits  = iter;
    449:   CHKMEMQ;
    450:   return(0);
    451: }