chemshell_qchem / blob
1 int i;
2 double step, lens, fac, dg1, u1, u2, u3, u4;
3 double *tmp1, *tmp2;
4 double mxlen = 1.0;
6 /* create log entry */
7 log=1;
8 desc="cg step";
9 comment=" ";
11 printf("** commence cj step\n");
13 if(!rflag){
14 printf("CG step must follow a CG restart step\n");
15 return TCL_ERROR;
16 }
18 tmp1 = (double *) ckalloc(sizeof(double) * opt->nvar);
19 tmp2 = (double *) ckalloc(sizeof(double) * opt->nvar);
21 /* calculate the restart hessian times the current gradient. */
22 /*u1 = -dot(rd,g) / yy;*/
23 u1 = -MATRIX_dot(opt->st, opt->g) / yy;
25 /* u2 = dot(rd,g)* (double) 2.0 / dy - dot(ry,g) / yy; */
26 u2 = MATRIX_dot(opt->st, opt->g) * (double) 2.0 / dy - MATRIX_dot(opt->yt, opt->g) / yy;
28 u3 = dy / yy;
30 /*for (i = 0; i < n; ++i) xb[i] = -u3 * g[i] - u1 * ry[i] - u2 * rd[i];*/
32 /* save part of new step vector */
33 for (i = 0; i < opt->nvar; ++i) {
34 tmp1[i] = -u3 * opt->g->d[i] - u1 * opt->yt->d[i] - u2 * opt->st->d[i];
35 }
37 if (debug) printf("u1 %f u2 %f u3 %f\n",u1,u2,u3);
38 if (debug) {
39 printf(" tmp1=");for(i=0;i<opt->nvar;i++)printf("%f ",tmp1[i]); printf("\n");
40 }
42 /* calculate the restart hessian times the current y. */
44 u1 = (double) 0.;
45 u2 = (double) 0.;
46 u3 = (double) 0.;
47 u4 = (double) 0.;
49 for (i = 0; i < opt->nvar; ++i) {
51 /*u1 -= (g[i] - gb[i]) * rd[i] / yy;*/
52 u1 -= (opt->g->d[i] - opt->pg->d[i]) * opt->st->d[i] / yy;
54 /* u2 = u2 - (g[i] - gb[i]) * ry[i] / yy + rd[i] * (double) 2. * (g[i] - gb[i]) / dy;*/
55 u2 = u2 - (opt->g->d[i] - opt->pg->d[i]) * opt->yt->d[i] / yy +
56 opt->st->d[i] * (double) 2. * (opt->g->d[i] - opt->pg->d[i]) / dy;
58 /* u3 += v[i] * (g[i] - gb[i]); */
59 u3 += opt->s->d[i] * (opt->g->d[i] - opt->pg->d[i]);
61 }
63 for (i = 0; i < opt->nvar; ++i) {
64 step = dy / yy * (opt->g->d[i] - opt->pg->d[i]) + u1 * opt->yt->d[i] + u2 * opt->st->d[i];
65 u4 += step * (opt->g->d[i] - opt->pg->d[i]);
66 tmp2[i] = step;
67 }
69 if (debug) printf("u1 %f u2 %f u3 %f u4 %f\n",u1,u2,u3,u4);
71 /* calculate the doubly updated hessian times the current */
72 /* gradient to obtain the search vector. */
74 u1 = (double)0.;
75 u2 = (double)0.;
76 for (i = 0; i < opt->nvar; ++i) {
78 /*u1 -= v[i] * g[i] / u3;*/
79 u1 -= opt->s->d[i] * opt->g->d[i] / u3;
81 /*u2 = u2 + (u4 / u3 + (double) 1.) * v[i] * g[i] / u3 - tmp2[i] * g[i] / u3;*/
82 u2 = u2 + (u4 / u3 + (double) 1.) * opt->s->d[i] * opt->g->d[i] / u3 - tmp2[i] * opt->g->d[i] / u3;
83 }
85 if (debug) printf("u1 %f u2 %f\n",u1,u2);
86 if (debug) {
87 printf(" tmp2=");for(i=0;i<opt->nvar;i++)printf("%f ",tmp2[i]); printf("\n");
88 printf(" tmp1=");for(i=0;i<opt->nvar;i++)printf("%f ",tmp1[i]); printf("\n");
89 printf(" step=");for(i=0;i<opt->nvar;i++)printf("%f ",opt->s->d[i]); printf("\n");
90 }
92 for (i = 0; i < opt->nvar; ++i) opt->s->d[i] = tmp1[i] - u1 * tmp2[i] - u2 * opt->s->d[i];
94 if (debug) {
95 printf(" step=");for(i=0;i<opt->nvar;i++)printf("%f ",opt->s->d[i]); printf("\n");
96 }
98 /* still need to apply non-restart alpha here */
100 /* calculate the derivative along the new search vector. */
101 /*copy(xb,v);*/
102 /* update nrst to assure at least one restart every n iterations. */
103 /*if (nrst == n) nrst = 0;*/
104 /*++nrst;*/
106 /* if the new direction is not a descent direction,stop. */
107 dg1 = MATRIX_dot(opt->s,opt->g);
108 if (dg1 > 0.) {
109 printf("CG step is uphill\n");
110 nice_err = 1;
111 /* return TCL_ERROR;*/
112 }
114 /* if |s| becomes large through numerical inaccuracies or otherwise, scale down the step */
116 lens = MATRIX_dot(opt->s,opt->s);
117 if (lens > mxlen) {
118 if (debug) printf("newopt: cg step too large, scale down to unit length");
119 fac = 1.0/lens;
120 MATRIX_scale(opt->s,fac);
121 }
122 ckfree((char *) tmp1);
123 ckfree((char *) tmp2);
