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