← Back to branch summary

~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/bin/detcas/test_bfgs.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/detcas/test_bfgs.c
1
 
/*
2
 
** test_bfgs.c
3
 
**
4
 
** Try the Numerical Recipies code 
5
 
**
6
 
** C. David Sherrill
7
 
** Georgia Institute of Technology
8
 
** March 2004
9
 
*/
10
 
 
11
 
#include <stdio.h>
12
 
#include <stdlib.h>
13
 
#include <math.h>
14
 
#include <libciomr/libciomr.h>
15
 
 
16
 
void test_bfgs(void)
17
 
{
18
 
  int ndim, i, j;
19
 
  int iter, maxiter = 20;
20
 
  double E, E_last;
21
 
  double *x_cur, *x_last, *g_cur, *g_last;
22
 
  double *dx, *dg, *hdg;
23
 
  double **hessin;
24
 
  double fac, fad, fae;
25
 
  double dfunc(double *x, double *g);
26
 
 
27
 
  ndim = 2;
28
 
  E_last = 0.0;
29
 
 
30
 
  x_cur = init_array(ndim);
31
 
  x_last = init_array(ndim);
32
 
  g_cur = init_array(ndim);
33
 
  g_last = init_array(ndim);
34
 
  dx = init_array(ndim);
35
 
  dg = init_array(ndim);
36
 
  hdg = init_array(ndim);
37
 
  hessin = block_matrix(ndim,ndim);
38
 
 
39
 
  /* guess */
40
 
  x_cur[0] = 1.0;
41
 
  x_cur[1] = 5.0;
42
 
 
43
 
  for (i=0; i<ndim; i++) {
44
 
    hessin[i][i] = 1.0;
45
 
  }
46
 
 
47
 
  /* get value and gradient */
48
 
  E = dfunc(x_cur,g_cur);
49
 
 
50
 
  for (iter=0; iter<maxiter; iter++) {
51
 
 
52
 
    printf("Iteration %d.  Energy = %12.6lf\n", iter, E);
53
 
    printf("x and gradient\n");
54
 
    for (i=0; i<ndim; i++)
55
 
      printf("%12.6lf %12.6lf\n", x_cur[i], g_cur[i]);
56
 
    printf("\n");
57
 
 
58
 
    if (fabs(E - E_last) < 0.00001) {
59
 
      printf("Converged\n");
60
 
      exit(1);
61
 
    }
62
 
 
63
 
    for (i=0; i<ndim; i++) {
64
 
      dg[i] = g_cur[i] - g_last[i];
65
 
      dx[i] = x_cur[i] - x_last[i];
66
 
    }
67
 
 
68
 
    for (i=0; i<ndim; i++) {
69
 
      hdg[i]=0.0;
70
 
      for (j=0; j<ndim; j++) hdg[i] += hessin[i][j]*dg[j];
71
 
    }
72
 
 
73
 
    fac = fae = 0.0;
74
 
    for (i=0; i<ndim; i++) {
75
 
      fac += dg[i]*dx[i];
76
 
      fae += dg[i]*hdg[i];
77
 
    }
78
 
    fac = 1.0/fac;
79
 
    fad = 1.0/fae;
80
 
    for (i=0; i<ndim; i++) dg[i] = fac*dx[i] - fad*hdg[i];
81
 
    for (i=0; i<ndim; i++) {
82
 
      for (j=0; j<ndim; j++) {
83
 
        hessin[i][j] += fac*dx[i]*dx[j] - fad*hdg[i]*hdg[j] + fae*dg[i]*dg[j];
84
 
      }
85
 
    }
86
 
 
87
 
    printf("Hessian inverse:\n");
88
 
    for (i=0; i<ndim; i++) {
89
 
      for (j=0; j<ndim; j++) {
90
 
        printf("%12.6lf ", hessin[i][j]);
91
 
      }
92
 
      printf("\n");
93
 
    }
94
 
    printf("\n");
95
 
 
96
 
    for (i=0; i<ndim; i++) {
97
 
      x_last[i] = x_cur[i];
98
 
      g_last[i] = g_cur[i];
99
 
    }
100
 
 
101
 
    for (i=0; i<ndim; i++) {
102
 
      for (j=0; j<ndim; j++) {
103
 
        x_cur[i] -= hessin[i][j] * g_cur[j];
104
 
      }
105
 
    }
106
 
 
107
 
    E_last = E;
108
 
    E = dfunc(x_cur,g_cur);
109
 
 
110
 
 
111
 
  } /* end iteration */
112
 
 
113
 
  free(x_cur);  free(x_last);  free(g_cur);  free(g_last);
114
 
  free(dx);  free(dg);  free(hdg);
115
 
  free_block(hessin);
116
 
  exit(0); 
117
 
}
118
 
 
119
 
double dfunc(double *x, double *g)
120
 
{
121
 
  double E;
122
 
  /* right now try (x-1)^2 + 10y^2, min is 0 at x=1,y=0 */
123
 
  E = (x[0] - 1.0) * (x[0] - 1.0) + 10.0 * x[1] * x[1];
124
 
  g[0] = 2.0 * (x[0] - 1.0);
125
 
  g[1] = 20.0 * x[1];
126
 
  return(E);
127
 
}
128