← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/detci/oldcode/h0block_ols_upd.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2006-09-10 14:01:33 UTC
  • Revision ID: james.westby@ubuntu.com-20060910140133-ib2j86trekykfsfv
Tags: upstream-3.2.3
Import upstream version 3.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/detci/oldcode/h0block_ols_upd.c
 
1
/*
 
2
** H0block_ols_upd()
 
3
**
 
4
** Parameters:
 
5
**    Ms0  =  1 if lower triangle blocks of C are used to update upper 
 
6
**            triangle parts of H0block.
 
7
*/
 
8
void H0block_ols_upd(double E, double E_est, double phase,  
 
9
      double **Cn, double **Cc, double **S, double **Hd,
 
10
      double *norm, double *ovrlap, int ac, int bc,
 
11
      int len, int Ms0)
 
12
{
 
13
   int i;
 
14
   int al,bl;
 
15
   double tval,tval2,c;
 
16
   double *ccptr, *sptr, *cnptr, *hdptr;
 
17
   double nx = 0.0, ox = 0.0;
 
18
 
 
19
   ccptr = Cc[0];
 
20
   sptr = S[0];
 
21
   cnptr = Cn[0];
 
22
   hdptr = Hd[0];
 
23
 
 
24
   for (i=0; i<len; i++) {
 
25
      tval = hdptr[i] - E;
 
26
      if (fabs(tval) < HD_MIN) tval = HD_MIN ; /* prevent /0 */
 
27
      tval = 1.0 / tval;
 
28
      tval2 = ccptr[i] + tval * (E_est * ccptr[i] - sptr[i]);
 
29
      cnptr[i] = tval2;
 
30
      nx += tval2 * tval2;
 
31
      ox += tval2 * ccptr[i];
 
32
      }         
 
33
 
 
34
 
 
35
   for (i=0; i<H0block.size; i++) {
 
36
      if (ac != H0block.alplist[i] || bc != H0block.betlist[i]) continue;
 
37
      al = H0block.alpidx[i];
 
38
      bl = H0block.betidx[i];
 
39
      tval = Hd[al][bl] - E;
 
40
      if (fabs(tval) < HD_MIN) tval = HD_MIN ; /* prevent /0 */
 
41
      tval = 1.0 / tval;
 
42
      c = Cc[al][bl];
 
43
      /* tval2 = c + tval * (E_est * c - S[al][bl]); */
 
44
      tval2 = Cn[al][bl];
 
45
      nx -= tval2 * tval2;
 
46
      ox -= tval2 * c;
 
47
      tval = c + E_est * H0block.c0bp[i];
 
48
      tval -= H0block.s0bp[i];
 
49
      Cn[al][bl] = tval;
 
50
      H0block.c0b[i] = tval; /* this should gather all of c0b. Norm later */
 
51
      nx += tval * tval;
 
52
      ox += tval * c;
 
53
      }
 
54
 
 
55
   /* get upper triangle of H0block */
 
56
   for (i=0; i<H0block.size; i++) {
 
57
      if (ac != H0block.betlist[i] || bc != H0block.alplist[i]) continue;
 
58
      al = H0block.alpidx[i];
 
59
      bl = H0block.betidx[i];
 
60
      H0block.c0b[i] = phase * Cn[bl][al]; 
 
61
      }
 
62
 
 
63
   *norm = nx;
 
64
   *ovrlap = ox;
 
65
}
 
66
 
 
67
 
 
68
/* call for Ms=0 cases.  If not Ms=0, then call H0block_ols_upd */
 
69
void H0block_diag_ols_upd(double E, double E_est, double phase, 
 
70
      double **Cn, double **Cc, double **S, double **Hd,
 
71
      double *norm, double *ovrlap, int ac, int bc, int len)
 
72
{
 
73
   int i,j;
 
74
   int al,bl;
 
75
   double tval,tval2,c;
 
76
   double nx = 0.0, ox = 0.0;
 
77
 
 
78
   for (i=0; i<len; i++) {
 
79
      tval = Hd[i][i] - E;
 
80
      if (fabs(tval) < HD_MIN) tval = HD_MIN ; /* prevent /0 */
 
81
      tval = 1.0 / tval;
 
82
      tval2 = Cc[i][i] + tval * (E_est * Cc[i][i] - S[i][i]);
 
83
      Cn[i][i] = tval2;
 
84
      nx += tval2 * tval2;
 
85
      ox += tval2 * Cc[i][i];
 
86
      
 
87
      for (j=0; j<i; j++) {
 
88
         tval = Hd[i][j] - E;
 
89
         if (fabs(tval) < HD_MIN) tval = HD_MIN ; /* prevent /0 */
 
90
         tval = 1.0 / tval;
 
91
         tval2 = Cc[i][j] + tval * (E_est * Cc[i][j] - S[i][j]);
 
92
         Cn[i][j] = tval2;
 
93
         nx += 2.0 * tval2 * tval2;
 
94
         ox += 2.0 * tval2 * Cc[i][j];
 
95
         }         
 
96
       }
 
97
 
 
98
 
 
99
   /* the code below assumes Ms=0 and that h0block members come in pairs
 
100
    * if off-diagonal
 
101
    */
 
102
   for (i=0; i<H0block.size; i++) {
 
103
      if (ac != H0block.alplist[i] || bc != H0block.betlist[i]) continue;
 
104
      al = H0block.alpidx[i];
 
105
      bl = H0block.betidx[i];
 
106
      c = Cc[al][bl];
 
107
      if (al > bl) {
 
108
         tval2 = Cn[al][bl];
 
109
         nx -= 2.0 * tval2 * tval2;
 
110
         ox -= 2.0 * tval2 * c;         
 
111
         tval = c + E_est * H0block.c0bp[i];
 
112
         tval -= H0block.s0bp[i];
 
113
         Cn[al][bl] = tval;
 
114
         nx += 2.0 * tval * tval;
 
115
         ox += 2.0 * tval * c;
 
116
         H0block.c0b[i] = tval; /* this should gather all of c0b. Norm later */
 
117
         }
 
118
      else if (al == bl) {
 
119
         tval2 = Cn[al][bl];
 
120
         nx -= tval2 * tval2;
 
121
         ox -= tval2 * c;         
 
122
         tval = c + E_est * H0block.c0bp[i];
 
123
         tval -= H0block.s0bp[i];
 
124
         Cn[al][bl] = tval;
 
125
         nx += tval * tval;
 
126
         ox += tval * c;
 
127
         H0block.c0b[i] = tval; /* this should gather all of c0b. Norm later */
 
128
         }
 
129
      }
 
130
 
 
131
   for (i=0; i<H0block.size; i++) {
 
132
      if (ac != H0block.alplist[i] || bc != H0block.betlist[i]) continue;
 
133
      al = H0block.alpidx[i];
 
134
      bl = H0block.betidx[i];
 
135
      if (al < bl) {
 
136
         H0block.c0b[i] = phase * Cn[bl][al];
 
137
         }
 
138
      }
 
139
 
 
140
   *norm = nx;
 
141
   *ovrlap = ox;
 
142
}
 
143