Actual source code: gqt.c

  1: #include "petsc.h"
  2: #include "petscblaslapack.h"
  3: #include "taolapack.h"



  9: static PetscErrorCode estsv(PetscInt n, PetscReal *r, PetscInt ldr, PetscReal *svmin, PetscReal *z) {

 11:     PetscBLASInt blas1=1, blasn=n, blasnmi, blasj, blasldr = ldr;
 12:     PetscInt i,j;
 13:     PetscReal e,temp,w,wm,ynorm,znorm,s,sm;
 15:     for (i=0;i<n;i++) {
 16:         z[i]=0.0;
 17:     }

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

 33:             if (PetscAbs(e - z[i]) > PetscAbs(r[i + ldr*i])) {
 34:                 temp = PetscMin(0.01,PetscAbs(r[i + ldr*i]))/
 35:                     PetscAbs(e-z[i]);
 36:                 BLASscal_(&blasn, &temp, z, &blas1);
 37:                 e = temp*e;
 38:             }
 39:         
 40:             /* Determine the two possible choices of y[i] */
 41:             if (r[i + ldr*i] == 0.0)
 42:                 w = wm = 1.0;
 43:             else {
 44:                 w = (e - z[i]) / r[i + ldr*i];
 45:                 wm = - (e + z[i]) / r[i + ldr*i];
 46:             }

 48:             /*  Chose y[i] based on the predicted value of y[j] for j>i */
 49:             s = PetscAbs(e - z[i]);
 50:             sm = PetscAbs(e + z[i]);
 51:             for (j=i+1;j<n;j++) {
 52:                 sm += PetscAbs(z[j] + wm * r[i + ldr*j]);
 53:             }
 54:             if (i < n-1) {
 55:                 blasnmi = n-i-1;
 56:                 BLASaxpy_(&blasnmi, &w, &r[i + ldr*(i+1)], &blasldr, &z[i+1], &blas1);
 57:                 s += BLASasum_(&blasnmi, &z[i+1], &blas1);
 58:             }
 59:             if (s < sm) {
 60:                 temp = wm - w;
 61:                 w = wm;
 62:                 if (i < n-1) {
 63:                     BLASaxpy_(&blasnmi, &temp, &r[i + ldr*(i+1)], &blasldr, &z[i+1], &blas1);
 64:                 }
 65:             }
 66:             z[i] = w;
 67:         }

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

 71:         /* Solve R*z = y */
 72:         for (j=n-1; j>=0; j--) {
 73:             
 74:           /* Scale z */
 75:             
 76:             if (PetscAbs(z[j]) > PetscAbs(r[j + ldr*j])) {
 77:                 temp = PetscMin(0.01, PetscAbs(r[j + ldr*j] / z[j]));
 78:                 BLASscal_(&blasn, &temp, z, &blas1);
 79:                 ynorm *=temp;
 80:             }
 81:             if (r[j + ldr*j] == 0) {
 82:                 z[j] = 1.0;
 83:             } else {
 84:                 z[j] = z[j] / r[j + ldr*j];
 85:             }
 86:             temp = -z[j];
 87:             blasj=j;
 88:             BLASaxpy_(&blasj,&temp,&r[0+ldr*j],&blas1,z,&blas1);
 89:         }
 90:                                 
 91:         /* Compute svmin and normalize z */
 92:         znorm = 1.0 / BLASnrm2_(&blasn, z, &blas1);
 93:         *svmin = ynorm*znorm;
 94:         BLASscal_(&blasn, &znorm, z, &blas1);

 96:     }
 97:                     
 98:     return(0);
 99: }
100:                      


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

236: PetscErrorCode gqt(PetscInt n, PetscReal *a, PetscInt lda, PetscReal *b, 
237:                    PetscReal delta, PetscReal rtol, PetscReal atol,
238:                    PetscInt itmax, PetscReal *retpar, PetscReal *retf,
239:                    PetscReal *x, PetscInt *retinfo, PetscInt *retits, 
240:                    PetscReal *z, PetscReal *wa1, PetscReal *wa2)
241: {
243:     PetscReal f=0.0,p001=0.001,p5=0.5,minusone=-1,delta2=delta*delta;
244:     PetscInt iter, j, rednc,info;
245:     PetscBLASInt indef;
246:     PetscBLASInt blas1=1, blasn=n, iblas, blaslda = lda,blasldap1=lda+1,blasinfo;
247:     PetscReal alpha, anorm, bnorm, parc, parf, parl, pars, par=*retpar,
248:         paru, prod, rxnorm, rznorm=0.0, temp, xnorm;
249:     
251:     iter = 0;
252:     parf = 0.0;
253:     xnorm = 0.0;
254:     rxnorm = 0.0;
255:     rednc = 0;
256:     for (j=0; j<n; j++) {
257:         x[j] = 0.0;
258:         z[j] = 0.0;
259:     }
260:     
261:     /* Copy the diagonal and save A in its lower triangle */
262:     BLAScopy_(&blasn,a,&blasldap1, wa1, &blas1);
263:     CHKMEMQ;
264:     for (j=0;j<n-1;j++) {
265:         iblas = n - j - 1;
266:         BLAScopy_(&iblas,&a[j + lda*(j+1)], &blaslda, &a[j+1 + lda*j], &blas1);
267:         CHKMEMQ;
268:     }
269:     
270:     /* Calculate the l1-norm of A, the Gershgorin row sums, and the
271:        l2-norm of b */
272:     anorm = 0.0;
273:     for (j=0;j<n;j++) {
274:         wa2[j] = BLASasum_(&blasn, &a[0 + lda*j], &blas1);
275:         CHKMEMQ;
276:         anorm = PetscMax(anorm,wa2[j]);
277:     }
278:     for (j=0;j<n;j++) {
279:         wa2[j] = wa2[j] - PetscAbs(wa1[j]);
280:     }
281:     bnorm = BLASnrm2_(&blasn,b,&blas1);
282:     CHKMEMQ;
283:     /* Calculate a lower bound, pars, for the domain of the problem.
284:        Also calculate an upper bound, paru, and a lower bound, parl,
285:        for the Lagrange multiplier. */
286:     pars = parl = paru = -anorm;
287:     for (j=0;j<n;j++) {
288:         pars = PetscMax(pars, -wa1[j]);
289:         parl = PetscMax(parl, wa1[j] + wa2[j]);
290:         paru = PetscMax(paru, -wa1[j] + wa2[j]);
291:     }
292:     parl = PetscMax(bnorm/delta - parl,pars);
293:     parl = PetscMax(0.0,parl);
294:     paru = PetscMax(0.0, bnorm/delta + paru);
295:     
296:     /* If the input par lies outside of the interval (parl, paru),
297:        set par to the closer endpoint. */

299:     par = PetscMax(par,parl);
300:     par = PetscMin(par,paru);

302:     /* Special case: parl == paru */
303:     paru = PetscMax(paru, (1.0 + rtol)*parl);
304:     
305:     /* Beginning of an iteration */

307:     info = 0;
308:     for (iter=1;iter<=itmax;iter++) {
309:         
310:       /* Safeguard par */
311:         if (par <= pars && paru > 0) {
312:             par = PetscMax(p001, PetscSqrtScalar(parl/paru)) * paru;
313:         }
314:         
315:         /* Copy the lower triangle of A into its upper triangle and 
316:            compute A + par*I */

318:         for (j=0;j<n-1;j++) {
319:             iblas = n - j - 1;
320:             BLAScopy_(&iblas,&a[j+1 + j*lda], &blas1,
321:                       &a[j + (j+1)*lda], &blaslda);
322:             CHKMEMQ;
323:         }
324:         for (j=0;j<n;j++) {
325:             a[j + j*lda] = wa1[j] + par;
326:         }

328:         /* Attempt the Cholesky factorization of A without referencing
329:            the lower triangular part. */
330:         LAPACKpotrf_("U",&blasn,a,&blaslda,&indef);
331:         CHKMEMQ;
332:         
333:         /* Case 1: A + par*I is pos. def. */
334:         if (indef == 0) {

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

339:             parf = par;
340:             BLAScopy_(&blasn, b, &blas1, wa2, &blas1);
341:             CHKMEMQ;
342:             LAPACKtrtrs_("U","T","N",&blasn,&blas1,a,&blaslda,wa2,&blasn,&blasinfo);
343:             CHKMEMQ;
344:             rxnorm = BLASnrm2_(&blasn, wa2, &blas1);
345:             LAPACKtrtrs_("U","N","N",&blasn,&blas1,a,&blaslda,wa2,&blasn,&blasinfo);
346:             CHKMEMQ;
347:             BLAScopy_(&blasn, wa2, &blas1, x, &blas1);
348:             CHKMEMQ;
349:             BLASscal_(&blasn, &minusone, x, &blas1);
350:             CHKMEMQ;
351:             xnorm = BLASnrm2_(&blasn, x, &blas1);
352:             CHKMEMQ;

354:             /* Test for convergence */
355:             if (PetscAbs(xnorm - delta) <= rtol*delta ||
356:                 (par == 0  && xnorm <= (1.0+rtol)*delta)) {
357:                 info = 1;
358:             }
359:             
360:             /* Compute a direction of negative curvature and use this
361:                information to improve pars. */

363:             iblas=blasn*blasn;
364:             
365:             estsv(n,a,lda,&rznorm,z); 
366:             CHKMEMQ;
367:             pars = PetscMax(pars, par-rznorm*rznorm);
368:             
369:             /* Compute a negative curvature solution of the form
370:                x + alpha*z,  where norm(x+alpha*z)==delta */
371:             
372:             rednc = 0;
373:             if (xnorm < delta) {
374:                 
375:                 /* Compute alpha */
376:                 prod = BLASdot_(&blasn, z, &blas1, x, &blas1) / delta;
377:                 temp = (delta - xnorm)*((delta + xnorm)/delta);
378:                 alpha = temp/(PetscAbs(prod) + PetscSqrtScalar(prod*prod + temp/delta));
379:                 if (prod >= 0) 
380:                     alpha = PetscAbs(alpha);
381:                 else 
382:                     alpha =-PetscAbs(alpha);


385:                 /* Test to decide if the negative curvature step
386:                    produces a larger reduction than with z=0 */
387:                 
388:                 rznorm = PetscAbs(alpha) * rznorm;
389:                 if ((rznorm*rznorm + par*xnorm*xnorm)/(delta2) <= par) {
390:                     rednc = 1;
391:                 }
392:                 
393:                 /* Test for convergence */
394:                 if (p5 * rznorm*rznorm / delta2 <= 
395:                     rtol*(1.0-p5*rtol)*(par + rxnorm*rxnorm/delta2)) {
396:                     info = 1;
397:                 } else if (info == 0 && 
398:                            (p5*(par + rxnorm*rxnorm/delta2) <= atol/delta2)) {
399:                     info = 2;
400:                 }
401:             }
402:             
403:             /* Compute the Newton correction parc to par. */
404:             if (xnorm == 0) {
405:                 parc = -par;
406:             } else {
407:                 BLAScopy_(&blasn, x, &blas1, wa2, &blas1);
408:                 CHKMEMQ;
409:                 temp = 1.0/xnorm;
410:                 BLASscal_(&blasn, &temp, wa2, &blas1);
411:                 CHKMEMQ;
412:                 LAPACKtrtrs_("U","T","N",&blasn, &blas1, a, &blaslda, wa2, &blasn, &blasinfo);
413:                 CHKMEMQ;
414:                 temp = BLASnrm2_(&blasn, wa2, &blas1);
415:                 parc = (xnorm - delta)/(delta*temp*temp);
416:             }
417:             
418:             /* update parl or paru */
419:             if (xnorm > delta) {
420:                 parl = PetscMax(parl, par);
421:             } else if (xnorm < delta) {
422:                 paru = PetscMin(paru, par);
423:             }
424:         } else {
425:             /* Case 2: A + par*I is not pos. def. */
426:             
427:             /* Use the rank information from the Cholesky 
428:                decomposition to update par. */

430:             if (indef > 1) {
431:                 
432:                 /* Restore column indef to A + par*I. */
433:                 iblas = indef - 1;
434:                 BLAScopy_(&iblas,&a[indef-1 + 0*lda],&blaslda,
435:                           &a[0 + (indef-1)*lda],&blas1);
436:                 CHKMEMQ;
437:                 a[indef-1 + (indef-1)*lda] = wa1[indef-1] + par;
438:                 
439:                 /* compute parc. */
440:                 
441:                 BLAScopy_(&iblas,&a[0 + (indef-1)*lda], &blas1, wa2, &blas1);
442:                 CHKMEMQ;
443:                 LAPACKtrtrs_("U","T","N",&iblas,&blas1,a,&blaslda,wa2,&blasn,&blasinfo);
444:                 CHKMEMQ;
445:                 BLAScopy_(&iblas,wa2,&blas1,&a[0 + (indef-1)*lda],&blas1); 
446:                 CHKMEMQ;
447:                 temp = BLASnrm2_(&iblas,&a[0 + (indef-1)*lda],&blas1);
448:                 CHKMEMQ;
449:                 a[indef-1 + (indef-1)*lda] -= temp*temp;
450:                 LAPACKtrtrs_("U","N","N",&iblas,&blas1,a,&blaslda,wa2,&blasn,&blasinfo);
451:                 CHKMEMQ;
452:             }
453:             
454:             wa2[indef-1] = -1.0;
455:             iblas = indef;
456:             temp = BLASnrm2_(&iblas,wa2,&blas1);
457:             parc = - a[indef-1 + (indef-1)*lda]/(temp*temp);
458:             pars = PetscMax(pars,par+parc);
459:             
460:             /* If necessary, increase paru slightly.
461:                This is needed because in some exceptional situations
462:                paru is the optimal value of par. */
463:             
464:             paru = PetscMax(paru, (1.0+rtol)*pars);
465:         }

467:         /* Use pars to update parl */
468:         parl = PetscMax(parl,pars);
469:         
470:         /* Test for termination. */
471:         if (info == 0) {
472:             if (iter == itmax) info=4;
473:             if (paru <= (1.0+p5*rtol)*pars) info=3;
474:             if (paru == 0.0) info = 2;
475:         }
476:         
477:         /* If exiting, store the best approximation and restore
478:            the upper triangle of A. */
479:         
480:         if (info != 0) {
481:             
482:           /* Compute the best current estimates for x and f. */
483:             
484:             par = parf;
485:             f = -p5 * (rxnorm*rxnorm + par*xnorm*xnorm);
486:             if (rednc) {
487:                 f = -p5 * (rxnorm*rxnorm + par*delta*delta - rznorm*rznorm);
488:                 BLASaxpy_(&blasn, &alpha, z, &blas1, x, &blas1);
489:                 CHKMEMQ;
490:             }

492:             /* Restore the upper triangle of A */
493:             
494:             for (j = 0; j<n; j++) {
495:                 iblas = n - j - 1;
496:                 BLAScopy_(&iblas,&a[j+1 + j*lda],&blas1, &a[j + (j+1)*lda],&blaslda);
497:                 CHKMEMQ;
498:             }
499:             iblas = lda+1;
500:             BLAScopy_(&blasn,wa1,&blas1,a,&iblas);
501:             CHKMEMQ;
502:             break;
503:         }
504:         par = PetscMax(parl,par+parc);
505:     }
506:     *retpar = par;
507:     *retf = f;
508:     *retinfo = info;
509:     *retits = iter;
510:     CHKMEMQ;
511:     return(0);
512: }