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: }