📄 cg.c
字号:
/* napack/cg.f -- translated by f2c (version 20050501).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#ifdef __cplusplus
extern "C" {
#endif
#include "v3p_netlib.h"
#include <assert.h>
#include <stdio.h>
/* ________________________________________________________ */
/* | | */
/* | MINIMIZE A FUNCTION USING THE FLETCHER-REEVES FORM | */
/* | OF THE CONJUGATE GRADIENT METHOD | */
/* | WITH (OR WITHOUT) PRECONDITIONING | */
/* | | */
/* | INPUT: | */
/* | | */
/* | X --ARRAY CONTAINING STARTING GUESS | */
/* | | */
/* | STEP --STARTING GUESS FOR MINIMIZER IN DIREC- | */
/* | TION OF NEGATIVE GRADIENT DURING FIRST | */
/* | ITERATION (E. G. STEP=1) WHEN STEP=0, | */
/* | THE PROGRAM SELECTS A STARTING GUESS | */
/* | | */
/* | T --COMPUTING TOLERANCE (ITERATIONS STOP | */
/* | WHEN MAX-NORM OF GRADIENT .LE. T) | */
/* | | */
/* | LIMIT --MAXIMUM NUMBER OF ITERATIONS | */
/* | | */
/* | N --NUMBER OF UNKNOWNS | */
/* | | */
/* | M --NUMBER OF ITERATIONS UNTIL THE SEARCH | */
/* | DIRECTIONS ARE RENORMALIZED ALONG THE | */
/* | NEGATIVE GRADIENT (TYPICALLY, M = N) | */
/* | | */
/* | VALUE --NAME OF COST EVALUATION FUNC. ROUTINE | */
/* | (EXTERNAL IN MAIN PROGRAM) | */
/* | VALUE(X) IS VALUE OF COST AT X | */
/* | | */
/* | GRAD --NAME OF GRADIENT EVALUATION SUBROUTINE | */
/* | (EXTERNAL IN MAIN PROGRAM) | */
/* | GRAD(G,X) PUTS IN G THE GRADIENT AT X | */
/* | | */
/* | BOTH --NAME SUBROUTINE TO EVALUATE BOTH COST | */
/* | AND ITS GRADIENT (EXTERNAL IN MAIN | */
/* | PROGRAM) BOTH(V,G,X) PUTS THE VALUE IN | */
/* | V AND THE GRADIENT IN G FOR THE POINT X| */
/* | | */
/* | PRE --NAME OF PRECONDITIONING SUBROUTINE | */
/* | (EXTERNAL IN MAIN PROGRAM) | */
/* | PRE(Y,Z) APPLIES THE PRECONDITIONER TO | */
/* | Z, STORING THE RESULT IN Y. | */
/* | IF PRECONDITIONING NOT USED SET Y = Z | */
/* | | */
/* | H --WORK ARRAY (LENGTH AT LEAST 3N) | */
/* | | */
/* | OUTPUT: | */
/* | | */
/* | X --MINIMIZER | */
/* | | */
/* | E --MAX-NORM OF GRADIENT | */
/* | | */
/* | IT --NUMBER OF ITERATIONS PERFORMED | */
/* | | */
/* | STEP --STEP SIZE ALONG SEARCH DIRECTION FOR | */
/* | FINAL ITERATION | */
/* | | */
/* | BUILTIN FUNCTIONS: DABS,DEXP,IDINT,DLOG,DSQRT,DMAX1,| */
/* | DMIN1,DSIGN | */
/* | PACKAGE ROUTINES: CUB,FD,FV,FVD,INS | */
/* |________________________________________________________| */
/*< SUBROUTINE CG(X,E,IT,STEP,T,LIMIT,N,M,VALUE,GRAD,BOTH,PRE,H) >*/
/* Subroutine */ int cg_(doublereal *x, doublereal *e, integer *it,
doublereal *step, doublereal *t, integer *limit, integer *n, integer *
m,
double (*value)(double*,void*),
void (*grad)(double*,double*,void*),
void (*both)(double*,double*,double*,void*),
void (*pre)(double*,double*,void*),
doublereal *h__,
void* userdata,
integer* error_code)
{
/* Initialized data */
static doublereal a1 = .1; /* constant */
static doublereal a2 = .9; /* constant */
static doublereal a3 = 5.; /* constant */
static doublereal a4 = .2; /* constant */
static doublereal a5 = 10.; /* constant */
static doublereal a6 = .9; /* constant */
static doublereal a7 = .3; /* constant */
/* System generated locals */
integer h_dim1, h_offset, i__1;
doublereal d__1, d__2;
/* Builtin functions */
double log(doublereal), exp(doublereal), d_sign(doublereal *, doublereal *
), sqrt(doublereal);
integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
e_wsle();
/* Subroutine */ int s_stop(char *, ftnlen);
/* Local variables */
doublereal a, b, c__, d__, f, g;
integer i__, j, k, l=0;
doublereal p, q, r__, s, v=0, w=0, y[50], z__[50], a8,
c0, c1=0, d0, f0, f1, l3, da, db, fa, fb, fc;
extern doublereal fd_(doublereal *, doublereal *, doublereal *, integer *,
void (*grad)(double*,double*,void*), void*);
integer na=0, nb, nc, nd, iq=0;
extern doublereal fv_(doublereal *, doublereal *, doublereal *, integer *,
double (*value)(double*,void*), void*);
extern /* Subroutine */ int cub_(doublereal *, doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *), fvd_(
doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, integer *, void (*)(double*,double*,double*,void*), void*),
ins_(doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, doublereal *, integer *, doublereal *, doublereal *)
;
/*< INTEGER I,IT,J,K,L,LIMIT,M,N,NA,NB,NC,ND >*/
/*< REAL*8 H(N,1),X(1),Y(50),Z(50),A1,A2,A3,A4,A5,A6,A7,A8,A,B,C,C0,C1 >*/
/*< REAL*8 D,D0,DA,DB,E,F,F0,F1,FA,FB,FC,G,L3,P,Q,R,S,STEP,T,V,W >*/
/*< REAL*8 FV,FD,VALUE >*/
/*< EXTERNAL BOTH,GRAD,PRE,VALUE >*/
/*< DATA A1/.1D0/,A2/.9D0/,A3/5.D0/,A4/.2D0/,A5/10.D0/,A6/.9D0/ >*/
/* Parameter adjustments */
--x;
h_dim1 = *n;
h_offset = 1 + h_dim1;
h__ -= h_offset;
/* Initialize to no error. */
if(error_code)
{
*error_code = 0;
}
/* Function Body */
/*< DATA A7/.3D0/ >*/
/*< A8 = A3 + .01D0 >*/
a8 = a3 + .01;
/*< IT = 0 >*/
*it = 0;
/*< CALL BOTH(F,H(1,3),X) >*/
(*both)(&f, &h__[h_dim1 * 3 + 1], &x[1], userdata);
/*< E = 0. >*/
*e = (float)0.;
/*< DO 10 I = 1,N >*/
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< 10 IF ( DABS(H(I,3)) .GT. E ) E = DABS(H(I,3)) >*/
/* L10: */
if ((d__1 = h__[i__ + h_dim1 * 3], abs(d__1)) > *e) {
*e = (d__2 = h__[i__ + h_dim1 * 3], abs(d__2));
}
}
/*< IF ( E .LE. T ) RETURN >*/
if (*e <= *t) {
return 0;
}
/*< L3 = 1./DLOG(A3) >*/
l3 = (float)1. / log(a3);
/*< CALL PRE(H(1,2),H(1,3)) >*/
(*pre)(&h__[(h_dim1 << 1) + 1], &h__[h_dim1 * 3 + 1], userdata);
/*< A = STEP >*/
a = *step;
/*< IF ( A .GT. 0. ) GOTO 30 >*/
if (a > (float)0.) {
goto L30;
}
/*< DO 20 I = 1,N >*/
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< 20 IF ( DABS(X(I)) .GT. A ) A = DABS(X(I)) >*/
/* L20: */
if ((d__1 = x[i__], abs(d__1)) > a) {
a = (d__2 = x[i__], abs(d__2));
}
}
/*< A = .01*A/E >*/
a = a * (float).01 / *e;
/*< IF ( A .EQ. 0. ) A = 1. >*/
if (a == (float)0.) {
a = (float)1.;
}
/*< 30 G = 0. >*/
L30:
g = (float)0.;
/*< DO 40 I = 1,N >*/
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< 40 G = G + H(I,2)*H(I,3) >*/
/* L40: */
g += h__[i__ + (h_dim1 << 1)] * h__[i__ + h_dim1 * 3];
}
/*< IF ( G .LT. 0. ) GOTO 620 >*/
if (g < (float)0.) {
goto L620;
}
/*< 50 L = 0 >*/
L50:
l = 0;
/*< DO 60 I = 1,N >*/
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< 60 H(I,1) = -H(I,2) >*/
/* L60: */
h__[i__ + h_dim1] = -h__[i__ + (h_dim1 << 1)];
}
/*< D = -G >*/
d__ = -g;
/*< 70 FA = FV(A,X,H,N,VALUE) >*/
L70:
fa = fv_(&a, &x[1], &h__[h_offset], n, value, userdata);
/*< C0 = A >*/
c0 = a;
/*< F0 = FA >*/
f0 = fa;
/*< J = 2 >*/
j = 2;
/*< Y(1) = 0. >*/
y[0] = (float)0.;
/*< Z(1) = F >*/
z__[0] = f;
/*< Y(2) = A >*/
y[1] = a;
/*< Z(2) = FA >*/
z__[1] = fa;
/*< V = A1*D >*/
v = a1 * d__;
/*< W = A2*D >*/
w = a2 * d__;
/*< IQ = 0 >*/
iq = 0;
/*< IF ( FA .LE. F ) GOTO 80 >*/
if (fa <= f) {
goto L80;
}
/*< C = A >*/
c__ = a;
/*< B = 0. >*/
b = (float)0.;
/*< A = 0. >*/
a = (float)0.;
/*< FC = FA >*/
fc = fa;
/*< FB = F >*/
fb = f;
/*< FA = F >*/
fa = f;
/*< GOTO 90 >*/
goto L90;
/*< 80 C = 0. >*/
L80:
c__ = (float)0.;
/*< B = 0. >*/
b = (float)0.;
/*< FC = F >*/
fc = f;
/*< FB = F >*/
fb = f;
/*< IQ = 1 >*/
iq = 1;
/*< 90 NA = 0 >*/
L90:
na = 0;
/*< NB = 0 >*/
nb = 0;
/*< NC = 0 >*/
nc = 0;
/*< ND = 0 >*/
nd = 0;
/*< Q = (D+(F-F0)/C0)/C0 >*/
q = (d__ + (f - f0) / c0) / c0;
/*< IF ( Q .LT. 0. ) GOTO 110 >*/
if (q < (float)0.) {
goto L110;
}
/*< Q = A >*/
q = a;
/*< 100 ND = ND + 1 >*/
L100:
++nd;
/*< IF ( ND .GT. 25 ) GOTO 610 >*/
if (nd > 25) {
goto L610;
}
/*< Q = A3*Q >*/
q = a3 * q;
/*< P = FV(Q,X,H,N,VALUE) >*/
p = fv_(&q, &x[1], &h__[h_offset], n, value, userdata);
/*< CALL INS(Q,P,A,B,C,FA,FB,FC,J,Y,Z) >*/
ins_(&q, &p, &a, &b, &c__, &fa, &fb, &fc, &j, y, z__);
/*< IF ( P-F .LT. W*Q ) GOTO 100 >*/
if (p - f < w * q) {
goto L100;
}
/*< GOTO 260 >*/
goto L260;
/*< 110 Q = .5*D/Q >*/
L110:
q = d__ * (float).5 / q;
/*< IF ( Q .LT. .01*C0 ) Q = .01*C0 >*/
if (q < c0 * (float).01) {
q = c0 * (float).01;
}
/*< P = FV(Q,X,H,N,VALUE) >*/
p = fv_(&q, &x[1], &h__[h_offset], n, value, userdata);
/*< IF ( P .LE. F0 ) GOTO 120 >*/
if (p <= f0) {
goto L120;
}
/*< F1 = F0 >*/
f1 = f0;
/*< C1 = C0 >*/
c1 = c0;
/*< F0 = P >*/
f0 = p;
/*< C0 = Q >*/
c0 = q;
/*< GOTO 130 >*/
goto L130;
/*< 120 F1 = P >*/
L120:
f1 = p;
/*< C1 = Q >*/
c1 = q;
/*< 130 CALL INS(Q,P,A,B,C,FA,FB,FC,J,Y,Z) >*/
L130:
ins_(&q, &p, &a, &b, &c__, &fa, &fb, &fc, &j, y, z__);
/*< 135 IF ( A .EQ. 0. ) GOTO 140 >*/
L135:
if (a == (float)0.) {
goto L140;
}
/*< IF ( FA-F .GE. V*A ) GOTO 160 >*/
if (fa - f >= v * a) {
goto L160;
}
/*< IF ( FA-F .LT. W*A ) GOTO 210 >*/
if (fa - f < w * a) {
goto L210;
}
/*< GOTO 280 >*/
goto L280;
/*< 140 Q = C0 >*/
L140:
q = c0;
/*< IF ( C1 .LT. Q ) Q = C1 >*/
if (c1 < q) {
q = c1;
}
/*< 150 NA = NA + 1 >*/
L150:
++na;
/*< IF ( NA .GT. 25 ) GOTO 630 >*/
if (na > 25) {
goto L630;
}
/*< Q = A4*Q >*/
q = a4 * q;
/*< P = FV(Q,X,H,N,VALUE) >*/
p = fv_(&q, &x[1], &h__[h_offset], n, value, userdata);
/*< CALL INS(Q,P,A,B,C,FA,FB,FC,J,Y,Z) >*/
ins_(&q, &p, &a, &b, &c__, &fa, &fb, &fc, &j, y, z__);
/*< IF ( P-F .GE. V*Q ) GOTO 150 >*/
if (p - f >= v * q) {
goto L150;
}
/*< GOTO 250 >*/
goto L250;
/*< 160 IF ( C0 .GT. C1 ) GOTO 200 >*/
L160:
if (c0 > c1) {
goto L200;
}
/*< IF ( F0-F .GT. V*C0 ) GOTO 180 >*/
if (f0 - f > v * c0) {
goto L180;
}
/*< IF ( F0-F .GE. W*C0 ) GOTO 320 >*/
if (f0 - f >= w * c0) {
goto L320;
}
/*< IF ( C1 .LE. A5*C0 ) GOTO 320 >*/
if (c1 <= a5 * c0) {
goto L320;
}
/*< R = DLOG(C1/C0) >*/
r__ = log(c1 / c0);
/*< S = -IDINT(R*L3+.999) >*/
s = (doublereal) (-((integer) (r__ * l3 + (float).999)));
/*< R = .999*DEXP(R/S) >*/
r__ = exp(r__ / s) * (float).999;
/*< Q = C1 >*/
q = c1;
/*< 170 Q = Q*R >*/
L170:
q *= r__;
/*< IF ( Q .LT. C0 ) GOTO 320 >*/
if (q < c0) {
goto L320;
}
/*< P = FV(Q,X,H,N,VALUE) >*/
p = fv_(&q, &x[1], &h__[h_offset], n, value, userdata);
/*< CALL INS(Q,P,A,B,C,FA,FB,FC,J,Y,Z) >*/
ins_(&q, &p, &a, &b, &c__, &fa, &fb, &fc, &j, y, z__);
/*< NA = NA + 1 >*/
++na;
/*< IF ( P-F .GT. V*Q ) GOTO 170 >*/
if (p - f > v * q) {
goto L170;
}
/*< GOTO 320 >*/
goto L320;
/*< 180 Q = C0 >*/
L180:
q = c0;
/*< 190 NA = NA + 1 >*/
L190:
++na;
/*< IF ( NA .GT. 25 ) GOTO 630 >*/
if (na > 25) {
goto L630;
}
/*< Q = A4*Q >*/
q = a4 * q;
/*< P = FV(Q,X,H,N,VALUE) >*/
p = fv_(&q, &x[1], &h__[h_offset], n, value, userdata);
/*< CALL INS(Q,P,A,B,C,FA,FB,FC,J,Y,Z) >*/
ins_(&q, &p, &a, &b, &c__, &fa, &fb, &fc, &j, y, z__);
/*< IF ( P-F .GE. V*Q ) GOTO 190 >*/
if (p - f >= v * q) {
goto L190;
}
/*< GOTO 250 >*/
goto L250;
/*< 200 Q = A >*/
L200:
q = a;
/*< GOTO 190 >*/
goto L190;
/*< 210 IF ( C0 .LT. C1 ) GOTO 290 >*/
L210:
if (c0 < c1) {
goto L290;
}
/*< IF ( F0-F .GE. V*C0 ) GOTO 230 >*/
if (f0 - f >= v * c0) {
goto L230;
}
/*< IF ( F0-F .GE. W*C0 ) GOTO 250 >*/
if (f0 - f >= w * c0) {
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -