ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/diagc.c
Revision: 27
Committed: Tue Jul 8 19:10:25 2008 UTC (11 years, 10 months ago) by tjod
File size: 4643 byte(s)
Log Message:
move mengine src to reflect new freemol directory structure

Line User Rev File contents
1 tjod 3 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "pcmod.h"
5    
6     void tred2(double **,int , double *, double *);
7     void tqli(double *, double *, int , double **);
8     void message_alert(char *, char *);
9     double sign(double ,double );
10    
11     void tred2(double **a,int n, double *d, double *e)
12     {
13    
14     int i,j,k,l, nn;
15     double scale, hh, h, g, f;
16    
17     nn = n-1;
18     for (i= nn; i >= 1; i--)
19     {
20     l = i - 1;
21     h = scale = 0.0;
22     if (l > 0)
23     {
24     for (k=0; k <= l; k++)
25     scale += fabs(a[i][k]);
26     if (scale == 0.0)
27     e[i] = a[i][l];
28     else
29     {
30     for (k=0; k <= l; k++)
31     {
32     a[i][k] /= scale;
33     h += a[i][k]*a[i][k];
34     }
35     f = a[i][l];
36     g = f>0 ? -sqrt(h) : sqrt(h);
37     e[i] = scale*g;
38     h -= f*g;
39     a[i][l] = f-g;
40     f = 0.0;
41     for (j=0; j <= l; j++)
42     {
43     a[j][i] = a[i][j]/h;
44     g = 0.0;
45     for (k=0; k <= j; k++)
46     g += a[j][k]*a[i][k];
47     for (k = j+1; k <= l; k++)
48     g += a[k][j]*a[i][k];
49     e[j] = g/h;
50     f += e[j]*a[i][j];
51     }
52     hh = f/(h+h);
53     for (j=0; j <= l; j++)
54     {
55     f = a[i][j];
56     e[j] = g = e[j] - hh*f;
57     for (k=0; k <= j; k++)
58     a[j][k] -= (f*e[k]+g*a[i][k]);
59     }
60     }
61     }else
62     e[i] = a[i][l];
63     d[i] = h;
64     }
65     d[0] = 0.0;
66     e[0] = 0.0;
67     for (i=0; i < n; i++)
68     {
69     l = i - 1;
70     if (i > 0)
71     {
72     for (j=0; j <= l; j++)
73     {
74     g = 0.0;
75     for (k=0; k <= l; k++)
76     g += a[i][k]*a[k][j];
77     for (k=0; k <= l; k++)
78     a[k][j] -= g*a[k][i];
79     }
80     }
81     d[i] = a[i][i];
82     a[i][i] = 1.0;
83     for (j=0; j <= l; j++)
84     {
85     a[j][i] = a[i][j] = 0.0;
86     }
87     }
88     }
89     /* =========================================== */
90     void tqli(double *d, double *e, int n, double **z)
91     {
92     int m,l,iter,i,k,j;
93     double s,r,p,g,f,dd,c,b;
94    
95     for (i=1; i < n; i++)
96     e[i-1] = e[i];
97     e[n-1] = 0.0;
98    
99     for (l=0; l < n; l++)
100     {
101     iter = 0;
102     do
103     {
104     for (m=l; m < n-1; m++)
105     {
106     dd = fabs(d[m]) + fabs(d[m+1]);
107     if (fabs(e[m])+dd == dd) break;
108     }
109     if (m!= l)
110     {
111     if (iter++ == 30)
112     {
113     fprintf(pcmoutfile,"Error in tqli\n");
114     message_alert("Error in Tqli","Error");
115     return;
116     }
117     g = (d[l+1]-d[l])/(2.0*e[l]);
118     r = sqrt((g*g)+1.0);
119     g = d[m]-d[l]+e[l]/(g+sign(r,g));
120     s = c = 1.0;
121     p = 0.0;
122     for (i=m-1; i >= l; i--)
123     {
124     f = s*e[i];
125     b = c*e[i];
126     if (fabs(f) >= fabs(g))
127     {
128     c = g/f;
129     r = sqrt((c*c)+1.0);
130     e[i+1] = f*r;
131     c *= (s=1.0/r);
132     } else
133     {
134     s = f/g;
135     r = sqrt((s*s)+1.0);
136     e[i+1] = g*r;
137     s *= (c=1.0/r);
138     }
139     g = d[i+1]-p;
140     r = (d[i]-g)*s+2.0*c*b;
141     p = s*r;
142     d[i+1] = g+p;
143     g = c*r-b;
144     for (k=0; k < n; k++)
145     {
146     f = z[k][i+1];
147     z[k][i+1] = s*z[k][i]+c*f;
148     z[k][i] = c*z[k][i]-s*f;
149     }
150     }
151     d[l] = d[l]-p;
152     e[l] = g;
153     e[m] = 0.0;
154     }
155     } while (m!= l);
156     }
157     // sort d and z
158     for (i=0; i < n; i++)
159     {
160     p = d[k=i];
161     for (j=i+1; j < n; j++)
162     if (d[j] >= p) p = d[k=j];
163     if (k != i)
164     {
165     d[k] = d[i];
166     d[i] = p;
167     for (j=0; j < n; j++)
168     {
169     p = z[j][i];
170     z[j][i] = z[j][k];
171     z[j][k] = p;
172     }
173     }
174     }
175     }