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++) z[i] = 0.0;
 13:   e = PetscAbs(r[0]);
 14:   if (e == 0.0) {
 15:     *svmin = 0.0;
 16:     z[0]   = 1.0;
 17:   } else {
 18:     /* Solve R'*y = e */
 19:     for (i = 0; i < n; i++) {
 20:       /* Scale y. The scaling factor (0.01) reduces the number of scalings */
 21:       if (z[i] >= 0.0) e = -PetscAbs(e);
 22:       else e = PetscAbs(e);

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

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

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

 55:     PetscCallBLAS("BLASnrm2", ynorm = BLASnrm2_(&blasn, z, &blas1));

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

 75:     /* Compute svmin and normalize z */
 76:     PetscCallBLAS("BLASnrm2", znorm = 1.0 / BLASnrm2_(&blasn, z, &blas1));
 77:     *svmin = ynorm * znorm;
 78:     PetscCallBLAS("BLASscal", BLASscal_(&blasn, &znorm, z, &blas1));
 79:   }
 80:   return 0;
 81: }

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

221:   PetscBLASIntCast(n, &blasn);
222:   PetscBLASIntCast(lda, &blaslda);
223:   PetscBLASIntCast(lda + 1, &blasldap1);
224:   parf   = 0.0;
225:   xnorm  = 0.0;
226:   rxnorm = 0.0;
227:   rednc  = 0;
228:   for (j = 0; j < n; j++) {
229:     x[j] = 0.0;
230:     z[j] = 0.0;
231:   }

233:   /* Copy the diagonal and save A in its lower triangle */
234:   PetscCallBLAS("BLAScopy", BLAScopy_(&blasn, a, &blasldap1, wa1, &blas1));
235:   for (j = 0; j < n - 1; j++) {
236:     PetscBLASIntCast(n - j - 1, &iblas);
237:     PetscCallBLAS("BLAScopy", BLAScopy_(&iblas, &a[j + lda * (j + 1)], &blaslda, &a[j + 1 + lda * j], &blas1));
238:   }

240:   /* Calculate the l1-norm of A, the Gershgorin row sums, and the
241:    l2-norm of b */
242:   anorm = 0.0;
243:   for (j = 0; j < n; j++) {
244:     PetscCallBLAS("BLASasum", wa2[j] = BLASasum_(&blasn, &a[0 + lda * j], &blas1));
245:     CHKMEMQ;
246:     anorm = PetscMax(anorm, wa2[j]);
247:   }
248:   for (j = 0; j < n; j++) wa2[j] = wa2[j] - PetscAbs(wa1[j]);
249:   PetscCallBLAS("BLASnrm2", bnorm = BLASnrm2_(&blasn, b, &blas1));
250:   CHKMEMQ;
251:   /* Calculate a lower bound, pars, for the domain of the problem.
252:    Also calculate an upper bound, paru, and a lower bound, parl,
253:    for the Lagrange multiplier. */
254:   pars = parl = paru = -anorm;
255:   for (j = 0; j < n; j++) {
256:     pars = PetscMax(pars, -wa1[j]);
257:     parl = PetscMax(parl, wa1[j] + wa2[j]);
258:     paru = PetscMax(paru, -wa1[j] + wa2[j]);
259:   }
260:   parl = PetscMax(bnorm / delta - parl, pars);
261:   parl = PetscMax(0.0, parl);
262:   paru = PetscMax(0.0, bnorm / delta + paru);

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

267:   par = PetscMax(par, parl);
268:   par = PetscMin(par, paru);

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

273:   /* Beginning of an iteration */

275:   info = 0;
276:   for (iter = 1; iter <= itmax; iter++) {
277:     /* Safeguard par */
278:     if (par <= pars && paru > 0) par = PetscMax(p001, PetscSqrtScalar(parl / paru)) * paru;

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

282:     for (j = 0; j < n - 1; j++) {
283:       PetscBLASIntCast(n - j - 1, &iblas);
284:       PetscCallBLAS("BLAScopy", BLAScopy_(&iblas, &a[j + 1 + j * lda], &blas1, &a[j + (j + 1) * lda], &blaslda));
285:     }
286:     for (j = 0; j < n; j++) a[j + j * lda] = wa1[j] + par;

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

291:     /* Case 1: A + par*I is pos. def. */
292:     if (indef == 0) {
293:       /* Compute an approximate solution x and save the last value of par with A + par*I pos. def. */

295:       parf = par;
296:       PetscCallBLAS("BLAScopy", BLAScopy_(&blasn, b, &blas1, wa2, &blas1));
297:       PetscCallBLAS("LAPACKtrtrs", LAPACKtrtrs_("U", "T", "N", &blasn, &blas1, a, &blaslda, wa2, &blasn, &blasinfo));
299:       PetscCallBLAS("BLASnrm2", rxnorm = BLASnrm2_(&blasn, wa2, &blas1));
300:       PetscCallBLAS("LAPACKtrtrs", LAPACKtrtrs_("U", "N", "N", &blasn, &blas1, a, &blaslda, wa2, &blasn, &blasinfo));

303:       PetscCallBLAS("BLAScopy", BLAScopy_(&blasn, wa2, &blas1, x, &blas1));
304:       PetscCallBLAS("BLASscal", BLASscal_(&blasn, &minusone, x, &blas1));
305:       PetscCallBLAS("BLASnrm2", xnorm = BLASnrm2_(&blasn, x, &blas1));
306:       CHKMEMQ;

308:       /* Test for convergence */
309:       if (PetscAbs(xnorm - delta) <= rtol * delta || (par == 0 && xnorm <= (1.0 + rtol) * delta)) info = 1;

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

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

318:       rednc = 0;
319:       if (xnorm < delta) {
320:         /* Compute alpha */
321:         PetscCallBLAS("BLASdot", prod = BLASdot_(&blasn, z, &blas1, x, &blas1) / delta);
322:         temp  = (delta - xnorm) * ((delta + xnorm) / delta);
323:         alpha = temp / (PetscAbs(prod) + PetscSqrtScalar(prod * prod + temp / delta));
324:         if (prod >= 0) alpha = PetscAbs(alpha);
325:         else alpha = -PetscAbs(alpha);

327:         /* Test to decide if the negative curvature step produces a larger reduction than with z=0 */
328:         rznorm = PetscAbs(alpha) * rznorm;
329:         if ((rznorm * rznorm + par * xnorm * xnorm) / (delta2) <= par) rednc = 1;
330:         /* Test for convergence */
331:         if (p5 * rznorm * rznorm / delta2 <= rtol * (1.0 - p5 * rtol) * (par + rxnorm * rxnorm / delta2)) {
332:           info = 1;
333:         } else if (info == 0 && (p5 * (par + rxnorm * rxnorm / delta2) <= atol / delta2)) {
334:           info = 2;
335:         }
336:       }

338:       /* Compute the Newton correction parc to par. */
339:       if (xnorm == 0) {
340:         parc = -par;
341:       } else {
342:         PetscCallBLAS("BLAScopy", BLAScopy_(&blasn, x, &blas1, wa2, &blas1));
343:         temp = 1.0 / xnorm;
344:         PetscCallBLAS("BLASscal", BLASscal_(&blasn, &temp, wa2, &blas1));
345:         PetscCallBLAS("LAPACKtrtrs", LAPACKtrtrs_("U", "T", "N", &blasn, &blas1, a, &blaslda, wa2, &blasn, &blasinfo));
347:         PetscCallBLAS("BLASnrm2", temp = BLASnrm2_(&blasn, wa2, &blas1));
348:         parc = (xnorm - delta) / (delta * temp * temp);
349:       }

351:       /* update parl or paru */
352:       if (xnorm > delta) {
353:         parl = PetscMax(parl, par);
354:       } else if (xnorm < delta) {
355:         paru = PetscMin(paru, par);
356:       }
357:     } else {
358:       /* Case 2: A + par*I is not pos. def. */

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

362:       if (indef > 1) {
363:         /* Restore column indef to A + par*I. */
364:         iblas = indef - 1;
365:         PetscCallBLAS("BLAScopy", BLAScopy_(&iblas, &a[indef - 1 + 0 * lda], &blaslda, &a[0 + (indef - 1) * lda], &blas1));
366:         a[indef - 1 + (indef - 1) * lda] = wa1[indef - 1] + par;

368:         /* compute parc. */
369:         PetscCallBLAS("BLAScopy", BLAScopy_(&iblas, &a[0 + (indef - 1) * lda], &blas1, wa2, &blas1));
370:         PetscCallBLAS("LAPACKtrtrs", LAPACKtrtrs_("U", "T", "N", &iblas, &blas1, a, &blaslda, wa2, &blasn, &blasinfo));
372:         PetscCallBLAS("BLAScopy", BLAScopy_(&iblas, wa2, &blas1, &a[0 + (indef - 1) * lda], &blas1));
373:         PetscCallBLAS("BLASnrm2", temp = BLASnrm2_(&iblas, &a[0 + (indef - 1) * lda], &blas1));
374:         CHKMEMQ;
375:         a[indef - 1 + (indef - 1) * lda] -= temp * temp;
376:         PetscCallBLAS("LAPACKtrtrs", LAPACKtrtrs_("U", "N", "N", &iblas, &blas1, a, &blaslda, wa2, &blasn, &blasinfo));
378:       }

380:       wa2[indef - 1] = -1.0;
381:       iblas          = indef;
382:       PetscCallBLAS("BLASnrm2", temp = BLASnrm2_(&iblas, wa2, &blas1));
383:       parc = -a[indef - 1 + (indef - 1) * lda] / (temp * temp);
384:       pars = PetscMax(pars, par + parc);

386:       /* If necessary, increase paru slightly.
387:        This is needed because in some exceptional situations
388:        paru is the optimal value of par. */

390:       paru = PetscMax(paru, (1.0 + rtol) * pars);
391:     }

393:     /* Use pars to update parl */
394:     parl = PetscMax(parl, pars);

396:     /* Test for converged. */
397:     if (info == 0) {
398:       if (iter == itmax) info = 4;
399:       if (paru <= (1.0 + p5 * rtol) * pars) info = 3;
400:       if (paru == 0.0) info = 2;
401:     }

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

406:     if (info != 0) {
407:       /* Compute the best current estimates for x and f. */
408:       par = parf;
409:       f   = -p5 * (rxnorm * rxnorm + par * xnorm * xnorm);
410:       if (rednc) {
411:         f = -p5 * (rxnorm * rxnorm + par * delta * delta - rznorm * rznorm);
412:         PetscCallBLAS("BLASaxpy", BLASaxpy_(&blasn, &alpha, z, &blas1, x, &blas1));
413:       }
414:       /* Restore the upper triangle of A */
415:       for (j = 0; j < n; j++) {
416:         PetscBLASIntCast(n - j - 1, &iblas);
417:         PetscCallBLAS("BLAScopy", BLAScopy_(&iblas, &a[j + 1 + j * lda], &blas1, &a[j + (j + 1) * lda], &blaslda));
418:       }
419:       PetscBLASIntCast(lda + 1, &iblas);
420:       PetscCallBLAS("BLAScopy", BLAScopy_(&blasn, wa1, &blas1, a, &iblas));
421:       break;
422:     }
423:     par = PetscMax(parl, par + parc);
424:   }
425:   *retpar  = par;
426:   *retf    = f;
427:   *retinfo = info;
428:   *retits  = iter;
429:   CHKMEMQ;
430:   return 0;
431: }