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=n, blasnmi, blasj, blasldr = ldr;
7: PetscInt i,j;
8: PetscReal e,temp,w,wm,ynorm,znorm,s,sm;
11: for (i=0;i<n;i++) {
12: z[i]=0.0;
13: }
14: e = PetscAbs(r[0]);
15: if (e == 0.0) {
16: *svmin = 0.0;
17: z[0] = 1.0;
18: } else {
19: /* Solve R'*y = e */
20: for (i=0;i<n;i++) {
21: /* Scale y. The scaling factor (0.01) reduces the number of scalings */
22: if (z[i] >= 0.0) e =-PetscAbs(e);
23: else e = PetscAbs(e);
25: if (PetscAbs(e - z[i]) > PetscAbs(r[i + ldr*i])) {
26: temp = PetscMin(0.01,PetscAbs(r[i + ldr*i]))/PetscAbs(e-z[i]);
27: PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &temp, z, &blas1));
28: e = temp*e;
29: }
31: /* Determine the two possible choices of y[i] */
32: if (r[i + ldr*i] == 0.0) {
33: w = wm = 1.0;
34: } else {
35: w = (e - z[i]) / r[i + ldr*i];
36: wm = - (e + z[i]) / r[i + ldr*i];
37: }
39: /* Chose y[i] based on the predicted value of y[j] for j>i */
40: s = PetscAbs(e - z[i]);
41: sm = PetscAbs(e + z[i]);
42: for (j=i+1;j<n;j++) {
43: sm += PetscAbs(z[j] + wm * r[i + ldr*j]);
44: }
45: if (i < n-1) {
46: blasnmi = n-i-1;
47: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasnmi, &w, &r[i + ldr*(i+1)], &blasldr, &z[i+1], &blas1));
48: s += BLASasum_(&blasnmi, &z[i+1], &blas1);
49: }
50: if (s < sm) {
51: temp = wm - w;
52: w = wm;
53: if (i < n-1) {
54: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasnmi, &temp, &r[i + ldr*(i+1)], &blasldr, &z[i+1], &blas1));
55: }
56: }
57: z[i] = w;
58: }
60: ynorm = BLASnrm2_(&blasn, z, &blas1);
62: /* Solve R*z = y */
63: for (j=n-1; j>=0; j--) {
64: /* Scale z */
65: if (PetscAbs(z[j]) > PetscAbs(r[j + ldr*j])) {
66: temp = PetscMin(0.01, PetscAbs(r[j + ldr*j] / z[j]));
67: PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &temp, z, &blas1));
68: ynorm *=temp;
69: }
70: if (r[j + ldr*j] == 0) {
71: z[j] = 1.0;
72: } else {
73: z[j] = z[j] / r[j + ldr*j];
74: }
75: temp = -z[j];
76: blasj=j;
77: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasj,&temp,&r[0+ldr*j],&blas1,z,&blas1));
78: }
80: /* Compute svmin and normalize z */
81: znorm = 1.0 / BLASnrm2_(&blasn, z, &blas1);
82: *svmin = ynorm*znorm;
83: PetscStackCallBLAS("BLASscal",BLASscal_(&blasn, &znorm, z, &blas1));
84: }
85: return(0);
86: }
88: /*
89: c ***********
90: c
91: c Subroutine dgqt
92: c
93: c Given an n by n symmetric matrix A, an n-vector b, and a
94: c positive number delta, this subroutine determines a vector
95: c x which approximately minimizes the quadratic function
96: c 97: c f(x) = (1/2)*x'*A*x + b'*x
98: c 99: c subject to the Euclidean norm constraint100: c101: c norm(x) <= delta.
102: c103: c This subroutine computes an approximation x and a Lagrange104: c multiplier par such that either par is zero and105: c106: c norm(x) <= (1+rtol)*delta,
107: c108: c or par is positive and109: c110: c abs(norm(x) - delta) <= rtol*delta.
111: c
112: c If xsol is the solution to the problem, the approximation x
113: c satisfies114: c115: c f(x) <= ((1 - rtol)**2)*f(xsol)
116: c117: c The subroutine statement is118: c119: c subroutine dgqt(n,a,lda,b,delta,rtol,atol,itmax,120: c par,f,x,info,z,wa1,wa2)121: c
122: c where
123: c
124: c n is an integer variable.
125: c On entry n is the order of A.
126: c On exit n is unchanged.
127: c128: c a is a double precision array of dimension (lda,n).
129: c On entry the full upper triangle of a must contain the
130: c full upper triangle of the symmetric matrix A.
131: c On exit the array contains the matrix A.
132: c
133: c lda is an integer variable.
134: c On entry lda is the leading dimension of the array a.
135: c On exit lda is unchanged.
136: c
137: c b is an double precision array of dimension n.
138: c On entry b specifies the linear term in the quadratic.
139: c On exit b is unchanged.
140: c
141: c delta is a double precision variable.
142: c On entry delta is a bound on the Euclidean norm of x.
143: c On exit delta is unchanged.
144: c
145: c rtol is a double precision variable.
146: c On entry rtol is the relative accuracy desired in the
147: c solution. Convergence occurs if
148: c149: c f(x) <= ((1 - rtol)**2)*f(xsol)
150: c
151: c On exit rtol is unchanged.
152: c
153: c atol is a double precision variable.
154: c On entry atol is the absolute accuracy desired in the
155: c solution. Convergence occurs when
156: c157: c norm(x) <= (1 + rtol)*delta
158: c
159: c max(-f(x),-f(xsol)) <= atol
160: c
161: c On exit atol is unchanged.
162: c
163: c itmax is an integer variable.
164: c On entry itmax specifies the maximum number of iterations.
165: c On exit itmax is unchanged.
166: c
167: c par is a double precision variable.
168: c On entry par is an initial estimate of the Lagrange169: c multiplier for the constraint norm(x) <= delta.
170: c On exit par contains the final estimate of the multiplier.
171: c
172: c f is a double precision variable.
173: c On entry f need not be specified.
174: c On exit f is set to f(x) at the output x.
175: c
176: c x is a double precision array of dimension n.
177: c On entry x need not be specified.
178: c On exit x is set to the final estimate of the solution.
179: c
180: c info is an integer variable.
181: c On entry info need not be specified.
182: c On exit info is set as follows:
183: c
184: c info = 1 The function value f(x) has the relative
185: c accuracy specified by rtol.
186: c
187: c info = 2 The function value f(x) has the absolute
188: c accuracy specified by atol.
189: c
190: c info = 3 Rounding errors prevent further progress.
191: c On exit x is the best available approximation.
192: c
193: c info = 4 Failure to converge after itmax iterations.
194: c On exit x is the best available approximation.
195: c
196: c z is a double precision work array of dimension n.
197: c
198: c wa1 is a double precision work array of dimension n.
199: c
200: c wa2 is a double precision work array of dimension n.
201: c
202: c Subprograms called
203: c
204: c MINPACK-2 ...... destsv
205: c
206: c LAPACK ......... dpotrf
207: c
208: c Level 1 BLAS ... daxpy, dcopy, ddot, dnrm2, dscal
209: c
210: c Level 2 BLAS ... dtrmv, dtrsv
211: c
212: c MINPACK-2 Project. October 1993.
213: c Argonne National Laboratory and University of Minnesota.
214: c Brett M. Averick, Richard Carter, and Jorge J. More'
215: c
216: c ***********
217: */
218: PetscErrorCode gqt(PetscInt n, PetscReal *a, PetscInt lda, PetscReal *b,219: PetscReal delta, PetscReal rtol, PetscReal atol,220: PetscInt itmax, PetscReal *retpar, PetscReal *retf,221: PetscReal *x, PetscInt *retinfo, PetscInt *retits,222: PetscReal *z, PetscReal *wa1, PetscReal *wa2)223: {
225: PetscReal f=0.0,p001=0.001,p5=0.5,minusone=-1,delta2=delta*delta;
226: PetscInt iter, j, rednc,info;
227: PetscBLASInt indef;
228: PetscBLASInt blas1=1, blasn=n, iblas, blaslda = lda,blasldap1=lda+1,blasinfo;
229: PetscReal alpha, anorm, bnorm, parc, parf, parl, pars, par=*retpar,paru, prod, rxnorm, rznorm=0.0, temp, xnorm;
232: parf = 0.0;
233: xnorm = 0.0;
234: rxnorm = 0.0;
235: rednc = 0;
236: for (j=0; j<n; j++) {
237: x[j] = 0.0;
238: z[j] = 0.0;
239: }
241: /* Copy the diagonal and save A in its lower triangle */
242: PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn,a,&blasldap1, wa1, &blas1));
243: for (j=0;j<n-1;j++) {
244: iblas = n - j - 1;
245: PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[j + lda*(j+1)], &blaslda, &a[j+1 + lda*j], &blas1));
246: }
248: /* Calculate the l1-norm of A, the Gershgorin row sums, and the
249: l2-norm of b */
250: anorm = 0.0;
251: for (j=0;j<n;j++) {
252: wa2[j] = BLASasum_(&blasn, &a[0 + lda*j], &blas1);
253: CHKMEMQ;
254: anorm = PetscMax(anorm,wa2[j]);
255: }
256: for (j=0;j<n;j++) {
257: wa2[j] = wa2[j] - PetscAbs(wa1[j]);
258: }
259: bnorm = BLASnrm2_(&blasn,b,&blas1);
260: CHKMEMQ;
261: /* Calculate a lower bound, pars, for the domain of the problem.
262: Also calculate an upper bound, paru, and a lower bound, parl,
263: for the Lagrange multiplier. */
264: pars = parl = paru = -anorm;
265: for (j=0;j<n;j++) {
266: pars = PetscMax(pars, -wa1[j]);
267: parl = PetscMax(parl, wa1[j] + wa2[j]);
268: paru = PetscMax(paru, -wa1[j] + wa2[j]);
269: }
270: parl = PetscMax(bnorm/delta - parl,pars);
271: parl = PetscMax(0.0,parl);
272: paru = PetscMax(0.0, bnorm/delta + paru);
274: /* If the input par lies outside of the interval (parl, paru),
275: set par to the closer endpoint. */
277: par = PetscMax(par,parl);
278: par = PetscMin(par,paru);
280: /* Special case: parl == paru */
281: paru = PetscMax(paru, (1.0 + rtol)*parl);
283: /* Beginning of an iteration */
285: info = 0;
286: for (iter=1;iter<=itmax;iter++) {
287: /* Safeguard par */
288: if (par <= pars && paru > 0) {
289: par = PetscMax(p001, PetscSqrtScalar(parl/paru)) * paru;
290: }
292: /* Copy the lower triangle of A into its upper triangle and
293: compute A + par*I */
295: for (j=0;j<n-1;j++) {
296: iblas = n - j - 1;
297: PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[j+1 + j*lda], &blas1,&a[j + (j+1)*lda], &blaslda));
298: }
299: for (j=0;j<n;j++) {
300: a[j + j*lda] = wa1[j] + par;
301: }
303: /* Attempt the Cholesky factorization of A without referencing
304: the lower triangular part. */
305: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("U",&blasn,a,&blaslda,&indef));
307: /* Case 1: A + par*I is pos. def. */
308: if (indef == 0) {
310: /* Compute an approximate solution x and save the
311: last value of par with A + par*I pos. def. */
313: parf = par;
314: PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn, b, &blas1, wa2, &blas1));
315: #if defined(PETSC_MISSING_LAPACK_TRTRS)
316: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"TRTRS - Lapack routine is unavailable.");
317: #else
318: PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&blasn,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
319: rxnorm = BLASnrm2_(&blasn, wa2, &blas1);
320: PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","N","N",&blasn,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
321: #endif
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: #if defined(PETSC_MISSING_LAPACK_TRTRS)
377: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"TRTRS - Lapack routine is unavailable.");
378: #else
379: PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&blasn, &blas1, a, &blaslda, wa2, &blasn, &blasinfo));
380: #endif
381: temp = BLASnrm2_(&blasn, wa2, &blas1);
382: parc = (xnorm - delta)/(delta*temp*temp);
383: }
385: /* update parl or paru */
386: if (xnorm > delta) {
387: parl = PetscMax(parl, par);
388: } else if (xnorm < delta) {
389: paru = PetscMin(paru, par);
390: }
391: } else {
392: /* Case 2: A + par*I is not pos. def. */
394: /* Use the rank information from the Cholesky
395: decomposition to update par. */
397: if (indef > 1) {
398: /* Restore column indef to A + par*I. */
399: iblas = indef - 1;
400: PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[indef-1 + 0*lda],&blaslda,&a[0 + (indef-1)*lda],&blas1));
401: a[indef-1 + (indef-1)*lda] = wa1[indef-1] + par;
403: /* compute parc. */
404: PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[0 + (indef-1)*lda], &blas1, wa2, &blas1));
405: #if defined(PETSC_MISSING_LAPACK_TRTRS)
406: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"TRTRS - Lapack routine is unavailable.");
407: #else
408: PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&iblas,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
409: PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,wa2,&blas1,&a[0 + (indef-1)*lda],&blas1));
410: temp = BLASnrm2_(&iblas,&a[0 + (indef-1)*lda],&blas1);
411: CHKMEMQ;
412: a[indef-1 + (indef-1)*lda] -= temp*temp;
413: PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","N","N",&iblas,&blas1,a,&blaslda,wa2,&blasn,&blasinfo));
414: #endif
415: }
417: wa2[indef-1] = -1.0;
418: iblas = indef;
419: temp = BLASnrm2_(&iblas,wa2,&blas1);
420: parc = - a[indef-1 + (indef-1)*lda]/(temp*temp);
421: pars = PetscMax(pars,par+parc);
423: /* If necessary, increase paru slightly.
424: This is needed because in some exceptional situations
425: paru is the optimal value of par. */
427: paru = PetscMax(paru, (1.0+rtol)*pars);
428: }
430: /* Use pars to update parl */
431: parl = PetscMax(parl,pars);
433: /* Test for converged. */
434: if (info == 0) {
435: if (iter == itmax) info=4;
436: if (paru <= (1.0+p5*rtol)*pars) info=3;
437: if (paru == 0.0) info = 2;
438: }
440: /* If exiting, store the best approximation and restore
441: the upper triangle of A. */
443: if (info != 0) {
444: /* Compute the best current estimates for x and f. */
445: par = parf;
446: f = -p5 * (rxnorm*rxnorm + par*xnorm*xnorm);
447: if (rednc) {
448: f = -p5 * (rxnorm*rxnorm + par*delta*delta - rznorm*rznorm);
449: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasn, &alpha, z, &blas1, x, &blas1));
450: }
451: /* Restore the upper triangle of A */
452: for (j = 0; j<n; j++) {
453: iblas = n - j - 1;
454: PetscStackCallBLAS("BLAScopy",BLAScopy_(&iblas,&a[j+1 + j*lda],&blas1, &a[j + (j+1)*lda],&blaslda));
455: }
456: iblas = lda+1;
457: PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn,wa1,&blas1,a,&iblas));
458: break;
459: }
460: par = PetscMax(parl,par+parc);
461: }
462: *retpar = par;
463: *retf = f;
464: *retinfo = info;
465: *retits = iter;
466: CHKMEMQ;
467: return(0);
468: }