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;
12: PetscBLASIntCast(n,&blasn);
13: PetscBLASIntCast(ldr,&blasldr);
14: for (i=0;i<n;i++) {
15: z[i]=0.0;
16: }
17: e = PetscAbs(r[0]);
18: if (e == 0.0) {
19: *svmin = 0.0;
20: z[0] = 1.0;
21: } else {
22: /* Solve R'*y = e */
23: for (i=0;i<n;i++) {
24: /* Scale y. The scaling factor (0.01) reduces the number of scalings */
25: if (z[i] >= 0.0) e =-PetscAbs(e);
26: else e = PetscAbs(e);
28: if (PetscAbs(e - z[i]) > PetscAbs(r[i + ldr*i])) {
29: temp = PetscMin(0.01,PetscAbs(r[i + ldr*i]))/PetscAbs(e-z[i]);
30: PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &temp, z, &blas1));
31: e = temp*e;
32: }
34: /* Determine the two possible choices of y[i] */
35: if (r[i + ldr*i] == 0.0) {
36: w = wm = 1.0;
37: } else {
38: w = (e - z[i]) / r[i + ldr*i];
39: wm = - (e + z[i]) / r[i + ldr*i];
40: }
42: /* Chose y[i] based on the predicted value of y[j] for j>i */
43: s = PetscAbs(e - z[i]);
44: sm = PetscAbs(e + z[i]);
45: for (j=i+1;j<n;j++) {
46: sm += PetscAbs(z[j] + wm * r[i + ldr*j]);
47: }
48: if (i < n-1) {
49: PetscBLASIntCast(n-i-1,&blasnmi);
50: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasnmi, &w, &r[i + ldr*(i+1)], &blasldr, &z[i+1], &blas1));
51: PetscStackCallBLAS("BLASasum",s += BLASasum_(&blasnmi, &z[i+1], &blas1));
52: }
53: if (s < sm) {
54: temp = wm - w;
55: w = wm;
56: if (i < n-1) {
57: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasnmi, &temp, &r[i + ldr*(i+1)], &blasldr, &z[i+1], &blas1));
58: }
59: }
60: z[i] = w;
61: }
63: PetscStackCallBLAS("BLASnrm2",ynorm = BLASnrm2_(&blasn, z, &blas1));
65: /* Solve R*z = y */
66: for (j=n-1; j>=0; j--) {
67: /* Scale z */
68: if (PetscAbs(z[j]) > PetscAbs(r[j + ldr*j])) {
69: temp = PetscMin(0.01, PetscAbs(r[j + ldr*j] / z[j]));
70: PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &temp, z, &blas1));
71: ynorm *=temp;
72: }
73: if (r[j + ldr*j] == 0) {
74: z[j] = 1.0;
75: } else {
76: z[j] = z[j] / r[j + ldr*j];
77: }
78: temp = -z[j];
79: PetscBLASIntCast(j,&blasj);
80: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasj,&temp,&r[0+ldr*j],&blas1,z,&blas1));
81: }
83: /* Compute svmin and normalize z */
84: PetscStackCallBLAS("BLASnrm2",znorm = 1.0 / BLASnrm2_(&blasn, z, &blas1));
85: *svmin = ynorm*znorm;
86: PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &znorm, z, &blas1));
87: }
88: return(0);
89: }
91: /*
92: c ***********
93: c
94: c Subroutine gqt
95: c
96: c Given an n by n symmetric matrix A, an n-vector b, and a
97: c positive number delta, this subroutine determines a vector
98: c x which approximately minimizes the quadratic function
99: c
100: c f(x) = (1/2)*x'*A*x + b'*x
101: c
102: c subject to the Euclidean norm constraint
103: c
104: c norm(x) <= delta.
105: c
106: c This subroutine computes an approximation x and a Lagrange
107: c multiplier par such that either par is zero and
108: c
109: c norm(x) <= (1+rtol)*delta,
110: c
111: c or par is positive and
112: c
113: c abs(norm(x) - delta) <= rtol*delta.
114: c
115: c If xsol is the solution to the problem, the approximation x
116: c satisfies
117: c
118: c f(x) <= ((1 - rtol)**2)*f(xsol)
119: c
120: c The subroutine statement is
121: c
122: c subroutine gqt(n,a,lda,b,delta,rtol,atol,itmax,
123: c par,f,x,info,z,wa1,wa2)
124: c
125: c where
126: c
127: c n is an integer variable.
128: c On entry n is the order of A.
129: c On exit n is unchanged.
130: c
131: c a is a double precision array of dimension (lda,n).
132: c On entry the full upper triangle of a must contain the
133: c full upper triangle of the symmetric matrix A.
134: c On exit the array contains the matrix A.
135: c
136: c lda is an integer variable.
137: c On entry lda is the leading dimension of the array a.
138: c On exit lda is unchanged.
139: c
140: c b is an double precision array of dimension n.
141: c On entry b specifies the linear term in the quadratic.
142: c On exit b is unchanged.
143: c
144: c delta is a double precision variable.
145: c On entry delta is a bound on the Euclidean norm of x.
146: c On exit delta is unchanged.
147: c
148: c rtol is a double precision variable.
149: c On entry rtol is the relative accuracy desired in the
150: c solution. Convergence occurs if
151: c
152: c f(x) <= ((1 - rtol)**2)*f(xsol)
153: c
154: c On exit rtol is unchanged.
155: c
156: c atol is a double precision variable.
157: c On entry atol is the absolute accuracy desired in the
158: c solution. Convergence occurs when
159: c
160: c norm(x) <= (1 + rtol)*delta
161: c
162: c max(-f(x),-f(xsol)) <= atol
163: c
164: c On exit atol is unchanged.
165: c
166: c itmax is an integer variable.
167: c On entry itmax specifies the maximum number of iterations.
168: c On exit itmax is unchanged.
169: c
170: c par is a double precision variable.
171: c On entry par is an initial estimate of the Lagrange
172: c multiplier for the constraint norm(x) <= delta.
173: c On exit par contains the final estimate of the multiplier.
174: c
175: c f is a double precision variable.
176: c On entry f need not be specified.
177: c On exit f is set to f(x) at the output x.
178: c
179: c x is a double precision array of dimension n.
180: c On entry x need not be specified.
181: c On exit x is set to the final estimate of the solution.
182: c
183: c info is an integer variable.
184: c On entry info need not be specified.
185: c On exit info is set as follows:
186: c
187: c info = 1 The function value f(x) has the relative
188: c accuracy specified by rtol.
189: c
190: c info = 2 The function value f(x) has the absolute
191: c accuracy specified by atol.
192: c
193: c info = 3 Rounding errors prevent further progress.
194: c On exit x is the best available approximation.
195: c
196: c info = 4 Failure to converge after itmax iterations.
197: c On exit x is the best available approximation.
198: c
199: c z is a double precision work array of dimension n.
200: c
201: c wa1 is a double precision work array of dimension n.
202: c
203: c wa2 is a double precision work array of dimension n.
204: c
205: c Subprograms called
206: c
207: c MINPACK-2 ...... destsv
208: c
209: c LAPACK ......... dpotrf
210: c
211: c Level 1 BLAS ... daxpy, dcopy, ddot, dnrm2, dscal
212: c
213: c Level 2 BLAS ... dtrmv, dtrsv
214: c
215: c MINPACK-2 Project. October 1993.
216: c Argonne National Laboratory and University of Minnesota.
217: c Brett M. Averick, Richard Carter, and Jorge J. More'
218: c
219: c ***********
220: */
221: PetscErrorCode gqt(PetscInt n, PetscReal *a, PetscInt lda, PetscReal *b,
222: PetscReal delta, PetscReal rtol, PetscReal atol,
223: PetscInt itmax, PetscReal *retpar, PetscReal *retf,
224: PetscReal *x, PetscInt *retinfo, PetscInt *retits,
225: PetscReal *z, PetscReal *wa1, PetscReal *wa2)
226: {
228: PetscReal f=0.0,p001=0.001,p5=0.5,minusone=-1,delta2=delta*delta;
229: PetscInt iter, j, rednc,info;
230: PetscBLASInt indef;
231: PetscBLASInt blas1=1, blasn, iblas, blaslda, blasldap1, blasinfo;
232: PetscReal alpha, anorm, bnorm, parc, parf, parl, pars, par=*retpar,paru, prod, rxnorm, rznorm=0.0, temp, xnorm;
235: PetscBLASIntCast(n,&blasn);
236: PetscBLASIntCast(lda,&blaslda);
237: PetscBLASIntCast(lda+1,&blasldap1);
238: parf = 0.0;
239: xnorm = 0.0;
240: rxnorm = 0.0;
241: rednc = 0;
242: for (j=0; j<n; j++) {
243: x[j] = 0.0;
244: z[j] = 0.0;
245: }
247: /* Copy the diagonal and save A in its lower triangle */
248: PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn,a,&blasldap1, wa1, &blas1));
249: for (j=0;j<n-1;j++) {
250: PetscBLASIntCast(n - j - 1,&iblas);
251: PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[j + lda*(j+1)], &blaslda, &a[j+1 + lda*j], &blas1));
252: }
254: /* Calculate the l1-norm of A, the Gershgorin row sums, and the
255: l2-norm of b */
256: anorm = 0.0;
257: for (j=0;j<n;j++) {
258: PetscStackCallBLAS("BLASasum",wa2[j] = BLASasum_(&blasn, &a[0 + lda*j], &blas1));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: PetscStackCallBLAS("BLASnrm2",bnorm = BLASnrm2_(&blasn,b,&blas1));CHKMEMQ;
265: /* Calculate a lower bound, pars, for the domain of the problem.
266: Also calculate an upper bound, paru, and a lower bound, parl,
267: for the Lagrange multiplier. */
268: pars = parl = paru = -anorm;
269: for (j=0;j<n;j++) {
270: pars = PetscMax(pars, -wa1[j]);
271: parl = PetscMax(parl, wa1[j] + wa2[j]);
272: paru = PetscMax(paru, -wa1[j] + wa2[j]);
273: }
274: parl = PetscMax(bnorm/delta - parl,pars);
275: parl = PetscMax(0.0,parl);
276: paru = PetscMax(0.0, bnorm/delta + paru);
278: /* If the input par lies outside of the interval (parl, paru),
279: set par to the closer endpoint. */
281: par = PetscMax(par,parl);
282: par = PetscMin(par,paru);
284: /* Special case: parl == paru */
285: paru = PetscMax(paru, (1.0 + rtol)*parl);
287: /* Beginning of an iteration */
289: info = 0;
290: for (iter=1;iter<=itmax;iter++) {
291: /* Safeguard par */
292: if (par <= pars && paru > 0) {
293: par = PetscMax(p001, PetscSqrtScalar(parl/paru)) * paru;
294: }
296: /* Copy the lower triangle of A into its upper triangle and compute A + par*I */
298: for (j=0;j<n-1;j++) {
299: PetscBLASIntCast(n - j - 1,&iblas);
300: PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[j+1 + j*lda], &blas1,&a[j + (j+1)*lda], &blaslda));
301: }
302: for (j=0;j<n;j++) {
303: a[j + j*lda] = wa1[j] + par;
304: }
306: /* Attempt the Cholesky factorization of A without referencing the lower triangular part. */
307: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("U",&blasn,a,&blaslda,&indef));
309: /* Case 1: A + par*I is pos. def. */
310: if (indef == 0) {
312: /* Compute an approximate solution x and save the last value of par with A + par*I pos. def. */
314: parf = par;
315: PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn, b, &blas1, wa2, &blas1));
316: PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&blasn,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
317: if (blasinfo) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKtrtrs() returned info %d",blasinfo);
318: PetscStackCallBLAS("BLASnrm2",rxnorm = BLASnrm2_(&blasn, wa2, &blas1));
319: PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","N","N",&blasn,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
320: if (blasinfo) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKtrtrs() returned info %d",blasinfo);
322: PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn, wa2, &blas1, x, &blas1));
323: PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &minusone, x, &blas1));
324: PetscStackCallBLAS("BLASnrm2",xnorm = BLASnrm2_(&blasn, x, &blas1));CHKMEMQ;
326: /* Test for convergence */
327: if (PetscAbs(xnorm - delta) <= rtol*delta || (par == 0 && xnorm <= (1.0+rtol)*delta)) {
328: info = 1;
329: }
331: /* Compute a direction of negative curvature and use this information to improve pars. */
332: estsv(n,a,lda,&rznorm,z);CHKMEMQ;
333: pars = PetscMax(pars, par-rznorm*rznorm);
335: /* Compute a negative curvature solution of the form x + alpha*z, where norm(x+alpha*z)==delta */
337: rednc = 0;
338: if (xnorm < delta) {
339: /* Compute alpha */
340: PetscStackCallBLAS("BLASdot",prod = BLASdot_(&blasn, z, &blas1, x, &blas1)/delta);
341: temp = (delta - xnorm)*((delta + xnorm)/delta);
342: alpha = temp/(PetscAbs(prod) + PetscSqrtScalar(prod*prod + temp/delta));
343: if (prod >= 0) alpha = PetscAbs(alpha);
344: else alpha =-PetscAbs(alpha);
346: /* Test to decide if the negative curvature step produces a larger reduction than with z=0 */
347: rznorm = PetscAbs(alpha) * rznorm;
348: if ((rznorm*rznorm + par*xnorm*xnorm)/(delta2) <= par) {
349: rednc = 1;
350: }
351: /* Test for convergence */
352: if (p5 * rznorm*rznorm / delta2 <= rtol*(1.0-p5*rtol)*(par + rxnorm*rxnorm/delta2)) {
353: info = 1;
354: } else if (info == 0 && (p5*(par + rxnorm*rxnorm/delta2) <= atol/delta2)) {
355: info = 2;
356: }
357: }
359: /* Compute the Newton correction parc to par. */
360: if (xnorm == 0) {
361: parc = -par;
362: } else {
363: PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn, x, &blas1, wa2, &blas1));
364: temp = 1.0/xnorm;
365: PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &temp, wa2, &blas1));
366: PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&blasn, &blas1, a, &blaslda, wa2, &blasn, &blasinfo));
367: if (blasinfo) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKtrtrs() returned info %d",blasinfo);
368: PetscStackCallBLAS("BLASnrm2",temp = BLASnrm2_(&blasn, wa2, &blas1));
369: parc = (xnorm - delta)/(delta*temp*temp);
370: }
372: /* update parl or paru */
373: if (xnorm > delta) {
374: parl = PetscMax(parl, par);
375: } else if (xnorm < delta) {
376: paru = PetscMin(paru, par);
377: }
378: } else {
379: /* Case 2: A + par*I is not pos. def. */
381: /* Use the rank information from the Cholesky decomposition to update par. */
383: if (indef > 1) {
384: /* Restore column indef to A + par*I. */
385: iblas = indef - 1;
386: PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[indef-1 + 0*lda],&blaslda,&a[0 + (indef-1)*lda],&blas1));
387: a[indef-1 + (indef-1)*lda] = wa1[indef-1] + par;
389: /* compute parc. */
390: PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[0 + (indef-1)*lda], &blas1, wa2, &blas1));
391: PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&iblas,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
392: if (blasinfo) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKtrtrs() returned info %d",blasinfo);
393: PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,wa2,&blas1,&a[0 + (indef-1)*lda],&blas1));
394: PetscStackCallBLAS("BLASnrm2",temp = BLASnrm2_(&iblas,&a[0 + (indef-1)*lda],&blas1));CHKMEMQ;
395: a[indef-1 + (indef-1)*lda] -= temp*temp;
396: PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","N","N",&iblas,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
397: if (blasinfo) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKtrtrs() returned info %d",blasinfo);
398: }
400: wa2[indef-1] = -1.0;
401: iblas = indef;
402: PetscStackCallBLAS("BLASnrm2",temp = BLASnrm2_(&iblas,wa2,&blas1));
403: parc = - a[indef-1 + (indef-1)*lda]/(temp*temp);
404: pars = PetscMax(pars,par+parc);
406: /* If necessary, increase paru slightly.
407: This is needed because in some exceptional situations
408: paru is the optimal value of par. */
410: paru = PetscMax(paru, (1.0+rtol)*pars);
411: }
413: /* Use pars to update parl */
414: parl = PetscMax(parl,pars);
416: /* Test for converged. */
417: if (info == 0) {
418: if (iter == itmax) info=4;
419: if (paru <= (1.0+p5*rtol)*pars) info=3;
420: if (paru == 0.0) info = 2;
421: }
423: /* If exiting, store the best approximation and restore
424: the upper triangle of A. */
426: if (info != 0) {
427: /* Compute the best current estimates for x and f. */
428: par = parf;
429: f = -p5 * (rxnorm*rxnorm + par*xnorm*xnorm);
430: if (rednc) {
431: f = -p5 * (rxnorm*rxnorm + par*delta*delta - rznorm*rznorm);
432: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasn, &alpha, z, &blas1, x, &blas1));
433: }
434: /* Restore the upper triangle of A */
435: for (j = 0; j<n; j++) {
436: PetscBLASIntCast(n - j - 1,&iblas);
437: PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[j+1 + j*lda],&blas1, &a[j + (j+1)*lda],&blaslda));
438: }
439: PetscBLASIntCast(lda+1,&iblas);
440: PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn,wa1,&blas1,a,&iblas));
441: break;
442: }
443: par = PetscMax(parl,par+parc);
444: }
445: *retpar = par;
446: *retf = f;
447: *retinfo = info;
448: *retits = iter;
449: CHKMEMQ;
450: return(0);
451: }