ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/diagc.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 4 months ago) by tjod
File size: 4643 byte(s)
Log Message:
test

Line File contents
1 #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 }