Actual source code: gqt.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: #include <petsc.h>
  2: #include <petscblaslapack.h>

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

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

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

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

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

 62:     ynorm = BLASnrm2_(&blasn, z, &blas1);

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

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

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

236:   iter = 0;
237:   parf = 0.0;
238:   xnorm = 0.0;
239:   rxnorm = 0.0;
240:   rednc = 0;
241:   for (j=0; j<n; j++) {
242:     x[j] = 0.0;
243:     z[j] = 0.0;
244:   }

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

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

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

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

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

288:   /* Beginning of an iteration */

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

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

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

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

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

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

318:       parf = par;
319:       PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn, b, &blas1, wa2, &blas1));
320:       PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&blasn,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
321:       rxnorm = BLASnrm2_(&blasn, wa2, &blas1);
322:       PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","N","N",&blasn,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
323:       PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn, wa2, &blas1, x, &blas1));
324:       PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &minusone, x, &blas1));
325:       xnorm = BLASnrm2_(&blasn, x, &blas1);
326:       CHKMEMQ;

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

334:       /* Compute a direction of negative curvature and use this
335:        information to improve pars. */

337:       iblas=blasn*blasn;

339:       estsv(n,a,lda,&rznorm,z);
340:       CHKMEMQ;
341:       pars = PetscMax(pars, par-rznorm*rznorm);

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

346:       rednc = 0;
347:       if (xnorm < delta) {
348:         /* Compute alpha */
349:         prod = BLASdot_(&blasn, z, &blas1, x, &blas1) / delta;
350:         temp = (delta - xnorm)*((delta + xnorm)/delta);
351:         alpha = temp/(PetscAbs(prod) + PetscSqrtScalar(prod*prod + temp/delta));
352:         if (prod >= 0) alpha = PetscAbs(alpha);
353:         else alpha =-PetscAbs(alpha);

355:                 /* Test to decide if the negative curvature step
356:                    produces a larger reduction than with z=0 */
357:         rznorm = PetscAbs(alpha) * rznorm;
358:         if ((rznorm*rznorm + par*xnorm*xnorm)/(delta2) <= par) {
359:           rednc = 1;
360:         }
361:         /* Test for convergence */
362:         if (p5 * rznorm*rznorm / delta2 <= rtol*(1.0-p5*rtol)*(par + rxnorm*rxnorm/delta2)) {
363:           info = 1;
364:         } else if (info == 0 && (p5*(par + rxnorm*rxnorm/delta2) <= atol/delta2)) {
365:           info = 2;
366:         }
367:       }

369:       /* Compute the Newton correction parc to par. */
370:       if (xnorm == 0) {
371:         parc = -par;
372:       } else {
373:         PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn, x, &blas1, wa2, &blas1));
374:         temp = 1.0/xnorm;
375:         PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &temp, wa2, &blas1));
376:         PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&blasn, &blas1, a, &blaslda, wa2, &blasn, &blasinfo));
377:         temp = BLASnrm2_(&blasn, wa2, &blas1);
378:         parc = (xnorm - delta)/(delta*temp*temp);
379:       }

381:       /* update parl or paru */
382:       if (xnorm > delta) {
383:         parl = PetscMax(parl, par);
384:       } else if (xnorm < delta) {
385:         paru = PetscMin(paru, par);
386:       }
387:     } else {
388:       /* Case 2: A + par*I is not pos. def. */

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

393:       if (indef > 1) {
394:         /* Restore column indef to A + par*I. */
395:         iblas = indef - 1;
396:         PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[indef-1 + 0*lda],&blaslda,&a[0 + (indef-1)*lda],&blas1));
397:         a[indef-1 + (indef-1)*lda] = wa1[indef-1] + par;

399:                 /* compute parc. */
400:         PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[0 + (indef-1)*lda], &blas1, wa2, &blas1));
401:         PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&iblas,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
402:         PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,wa2,&blas1,&a[0 + (indef-1)*lda],&blas1));
403:         temp = BLASnrm2_(&iblas,&a[0 + (indef-1)*lda],&blas1);
404:         CHKMEMQ;
405:         a[indef-1 + (indef-1)*lda] -= temp*temp;
406:         PetscStackCallBLAS("LAPACKtrtr",LAPACKtrtrs_("U","N","N",&iblas,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
407:       }

409:       wa2[indef-1] = -1.0;
410:       iblas = indef;
411:       temp = BLASnrm2_(&iblas,wa2,&blas1);
412:       parc = - a[indef-1 + (indef-1)*lda]/(temp*temp);
413:       pars = PetscMax(pars,par+parc);

415:       /* If necessary, increase paru slightly.
416:        This is needed because in some exceptional situations
417:        paru is the optimal value of par. */

419:       paru = PetscMax(paru, (1.0+rtol)*pars);
420:     }

422:     /* Use pars to update parl */
423:     parl = PetscMax(parl,pars);

425:     /* Test for converged. */
426:     if (info == 0) {
427:       if (iter == itmax) info=4;
428:       if (paru <= (1.0+p5*rtol)*pars) info=3;
429:       if (paru == 0.0) info = 2;
430:     }

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

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