## Revision 1264 trunk/src/linad99/eigen.cpp

eigen.cpp (revision 1264)
25 25
```    ad_exit(1);
```
26 26
```  }
```
27 27
```  dmatrix m1=symmetrize(m);
```
28
```  int n=m1.rowsize();
```
29
```  int cmin=m1.colmin();
```
30
```  int rmin=m1.rowmin();
```
31 28
```  m1.colshift(1);     // set minimum column and row indices to 1
```
32 29
```  m1.rowshift(1);
```
30
```  int n=m1.rowsize();
```
33 31
```  dvector diag(1,n);
```
34 32
```  dvector off_diag(1,n);
```
35 33

......
68 66
```    ad_exit(1);
```
69 67
```  }
```
70 68
```  int n=m.rowsize();
```
71
```  int l,k,j,i;
```
69
```  int l,j,i;
```
72 70
```  double scale,hh,h,g,f;
```
73 71

74 72
```  for (i=n;i>=2;i--)
```
......
77 75
```    h=scale=0.0;
```
78 76
```    if (l > 1)
```
79 77
```    {
```
80
```      for (k=1;k<=l;k++)
```
78
```      for (int k=1;k<=l;k++)
```
81 79
```        scale += fabs(m[i][k]);
```
82 80
```      if (scale == 0.0)
```
83 81
```        e[i]=m[i][l];
```
84 82
```      else
```
85 83
```      {
```
86
```        for (k=1;k<=l;k++)
```
84
```        for (int k=1;k<=l;k++)
```
87 85
```        {
```
88 86
```          m[i][k] /= scale;
```
89 87
```          h += m[i][k]*m[i][k];
```
......
101 99
```          m[j][i]=m[i][j]/h;
```
102 100
```        #endif
```
103 101
```          g=0.0;
```
104
```          for (k=1;k<=j;k++)
```
102
```          for (int k=1;k<=j;k++)
```
105 103
```            g += m[j][k]*m[i][k];
```
106
```          for (k=j+1;k<=l;k++)
```
104
```          for (int k=j+1;k<=l;k++)
```
107 105
```            g += m[k][j]*m[i][k];
```
108 106
```          e[j]=g/h;
```
109 107
```          f += e[j]*m[i][j];
```
......
113 111
```        {
```
114 112
```          f=m[i][j];
```
115 113
```          e[j]=g=e[j]-hh*f;
```
116
```          for (k=1;k<=j;k++)
```
114
```          for (int k=1;k<=j;k++)
```
117 115
```            m[j][k] -= (f*e[k]+g*m[i][k]);
```
118 116
```        }
```
119 117
```      }
```
......
138 136
```      for (j=1;j<=l;j++)
```
139 137
```      {
```
140 138
```        g=0.0;
```
141
```        for (k=1;k<=l;k++)
```
139
```        for (int k=1;k<=l;k++)
```
142 140
```          g += m[i][k]*m[k][j];
```
143
```        for (k=1;k<=l;k++)
```
141
```        for (int k=1;k<=l;k++)
```
144 142
```          m[k][j] -= g*m[k][i];
```
145 143
```      }
```
146 144
```    }
```
......
182 180
```{
```
183 181
```  dvector& d = (dvector&) _d;
```
184 182
```  dvector& e = (dvector&) _e;
```
183
```#ifdef EIGEN_VECTORS
```
185 184
```  dmatrix& z = (dmatrix&) _z;
```
185
```#endif
```
186 186
```  int max_iterations=30;
```
187 187
```  int n=d.size();
```
188 188
```  max_iterations+=10*(n/100);
```
189
```  int m,l,iter,i,k;
```
189
```  int m,l,iter,i;
```
190 190
```  double s,r,p,g,f,dd,c,b;
```
191 191

192 192
```  for (i=2;i<=n;i++) e[i-1]=e[i];
```
......
232 232
```          g=c*r-b;
```
233 233
```          /* Next loop can be omitted if eigenvectors not wanted */
```
234 234
```          #ifdef EIGEN_VECTORS
```
235
```            for (k=1;k<=n;k++)
```
235
```            for (int k=1;k<=n;k++)
```
236 236
```            {
```
237 237
```              f=z[k][i+1];
```
238 238
```              z[k][i+1]=s*z[k][i]+c*f;
```
......
266 266
```  int max_iterations=30;
```
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff