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;

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

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

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

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

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

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

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

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

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

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

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

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

277:   par = PetscMax(par,parl);
278:   par = PetscMin(par,paru);

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

283:   /* Beginning of an iteration */

285:   info = 0;
286:   for (iter=1;iter<=itmax;iter++) {
287:     /* Safeguard par */
288:     if (par <= pars && paru > 0) {
289:       par = PetscMax(p001, PetscSqrtScalar(parl/paru)) * paru;
290:     }

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

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

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

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

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

310:       parf = par;
311:       PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn, b, &blas1, wa2, &blas1));
312:       PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&blasn,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
314:       PetscStackCallBLAS("BLASnrm2",rxnorm = BLASnrm2_(&blasn, wa2, &blas1));
315:       PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","N","N",&blasn,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));

318:       PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn, wa2, &blas1, x, &blas1));
319:       PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &minusone, x, &blas1));
320:       PetscStackCallBLAS("BLASnrm2",xnorm = BLASnrm2_(&blasn, x, &blas1));CHKMEMQ;

322:       /* Test for convergence */
323:       if (PetscAbs(xnorm - delta) <= rtol*delta || (par == 0  && xnorm <= (1.0+rtol)*delta)) {
324:         info = 1;
325:       }

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

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

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

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

355:       /* Compute the Newton correction parc to par. */
356:       if (xnorm == 0) {
357:         parc = -par;
358:       } else {
359:         PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn, x, &blas1, wa2, &blas1));
360:         temp = 1.0/xnorm;
361:         PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &temp, wa2, &blas1));
362:         PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&blasn, &blas1, a, &blaslda, wa2, &blasn, &blasinfo));
364:         PetscStackCallBLAS("BLASnrm2",temp = BLASnrm2_(&blasn, wa2, &blas1));
365:         parc = (xnorm - delta)/(delta*temp*temp);
366:       }

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

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

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

385:         /* compute parc. */
386:         PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[0 + (indef-1)*lda], &blas1, wa2, &blas1));
387:         PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&iblas,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
389:         PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,wa2,&blas1,&a[0 + (indef-1)*lda],&blas1));
390:         PetscStackCallBLAS("BLASnrm2",temp = BLASnrm2_(&iblas,&a[0 + (indef-1)*lda],&blas1));CHKMEMQ;
391:         a[indef-1 + (indef-1)*lda] -= temp*temp;
392:         PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","N","N",&iblas,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
394:       }

396:       wa2[indef-1] = -1.0;
397:       iblas = indef;
398:       PetscStackCallBLAS("BLASnrm2",temp = BLASnrm2_(&iblas,wa2,&blas1));
399:       parc = - a[indef-1 + (indef-1)*lda]/(temp*temp);
400:       pars = PetscMax(pars,par+parc);

402:       /* If necessary, increase paru slightly.
403:        This is needed because in some exceptional situations
404:        paru is the optimal value of par. */

406:       paru = PetscMax(paru, (1.0+rtol)*pars);
407:     }

409:     /* Use pars to update parl */
410:     parl = PetscMax(parl,pars);

412:     /* Test for converged. */
413:     if (info == 0) {
414:       if (iter == itmax) info=4;
415:       if (paru <= (1.0+p5*rtol)*pars) info=3;
416:       if (paru == 0.0) info = 2;
417:     }

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

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