ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/math/SquareMatrix.hpp
(Generate patch)

Comparing trunk/src/math/SquareMatrix.hpp (file contents):
Revision 76 by tim, Thu Oct 14 23:28:09 2004 UTC vs.
Revision 93 by tim, Sun Oct 17 01:19:11 2004 UTC

# Line 175 | Line 175 | namespace oopse {
175           * @param w output eigenvalues
176           * @param v output eigenvectors
177           */
178 <        void jacobi(const SquareMatrix<Real, Dim>& a,
179 <                              Vector<Real, Dim>& w,
178 >        bool jacobi(const SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
179                                SquareMatrix<Real, Dim>& v);
180      };//end SquareMatrix
181  
# Line 184 | Line 183 | namespace oopse {
183   #define ROT(a,i,j,k,l) g=a(i, j);h=a(k, l);a(i, j)=g-s*(h+g*tau);a(k, l)=h+s*(g-h*tau)
184   #define MAX_ROTATIONS 60
185  
186 < template<Real, int Dim>
187 < void SquareMatrix<Real, int Dim>::jacobi(SquareMatrix<Real, Dim>& a,
188 <                                                                       Vector<Real, Dim>& w,
190 <                                                                       SquareMatrix<Real, Dim>& v) {
186 > template<typename Real, int Dim>
187 > bool SquareMatrix<Real, Dim>::jacobi(const SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
188 >                              SquareMatrix<Real, Dim>& v) {
189      const int N = Dim;                                                                      
190      int i, j, k, iq, ip;
191      double tresh, theta, tau, t, sm, s, h, g, c;
# Line 195 | Line 193 | void SquareMatrix<Real, int Dim>::jacobi(SquareMatrix<
193      Vector<Real, Dim> b, z;
194  
195      // initialize
196 <    for (ip=0; ip<N; ip++)
197 <    {
198 <        for (iq=0; iq<N; iq++) v(ip, iq) = 0.0;
199 <        v(ip, ip) = 1.0;
196 >    for (ip=0; ip<N; ip++) {
197 >        for (iq=0; iq<N; iq++)
198 >            v(ip, iq) = 0.0;
199 >        v(ip, ip) = 1.0;
200      }
201 <    for (ip=0; ip<N; ip++)
202 <    {
203 <        b(ip) = w(ip) = a(ip, ip);
204 <        z(ip) = 0.0;
201 >    
202 >    for (ip=0; ip<N; ip++) {
203 >        b(ip) = w(ip) = a(ip, ip);
204 >        z(ip) = 0.0;
205      }
206  
207      // begin rotation sequence
208 <    for (i=0; i<MAX_ROTATIONS; i++)
209 <    {
210 <        sm = 0.0;
211 <        for (ip=0; ip<2; ip++)
212 <        {
213 <            for (iq=ip+1; iq<N; iq++) sm += fabs(a(ip, iq));
214 <        }
215 <        if (sm == 0.0) break;
208 >    for (i=0; i<MAX_ROTATIONS; i++) {
209 >        sm = 0.0;
210 >        for (ip=0; ip<2; ip++) {
211 >            for (iq=ip+1; iq<N; iq++)
212 >                sm += fabs(a(ip, iq));
213 >        }
214 >        
215 >        if (sm == 0.0)
216 >            break;
217  
218 <        if (i < 4) tresh = 0.2*sm/(9);
219 <        else tresh = 0.0;
218 >        if (i < 4)
219 >            tresh = 0.2*sm/(9);
220 >        else
221 >            tresh = 0.0;
222  
223 <        for (ip=0; ip<2; ip++)
224 <        {
225 <            for (iq=ip+1; iq<N; iq++)
226 <            {
227 <                g = 100.0*fabs(a(ip, iq));
228 <                if (i > 4 && (fabs(w(ip))+g) == fabs(w(ip))
229 <                    && (fabs(w(iq))+g) == fabs(w(iq)))
230 <                {
231 <                    a(ip, iq) = 0.0;
232 <                }
233 <                else if (fabs(a(ip, iq)) > tresh)
234 <                {
235 <                    h = w(iq) - w(ip);
235 <                    if ( (fabs(h)+g) == fabs(h)) t = (a(ip, iq)) / h;
236 <                    else
237 <                    {
238 <                        theta = 0.5*h / (a(ip, iq));
239 <                        t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
240 <                        if (theta < 0.0) t = -t;
241 <                    }
242 <                    c = 1.0 / sqrt(1+t*t);
243 <                    s = t*c;
244 <                    tau = s/(1.0+c);
245 <                    h = t*a(ip, iq);
246 <                    z(ip) -= h;
247 <                    z(iq) += h;
248 <                    w(ip) -= h;
249 <                    w(iq) += h;
250 <                    a(ip, iq)=0.0;
251 <                    for (j=0;j<ip-1;j++)
252 <                    {
253 <                        ROT(a,j,ip,j,iq);
254 <                    }
255 <                    for (j=ip+1;j<iq-1;j++)
256 <                    {
257 <                        ROT(a,ip,j,j,iq);
258 <                    }
259 <                    for (j=iq+1; j<N; j++)
260 <                    {
261 <                        ROT(a,ip,j,iq,j);
262 <                    }
263 <                    for (j=0; j<N; j++)
264 <                    {
265 <                        ROT(v,j,ip,j,iq);
266 <                    }
267 <                }
268 <            }
269 <        }
223 >        for (ip=0; ip<2; ip++) {
224 >            for (iq=ip+1; iq<N; iq++) {
225 >                g = 100.0*fabs(a(ip, iq));
226 >                if (i > 4 && (fabs(w(ip))+g) == fabs(w(ip))
227 >                    && (fabs(w(iq))+g) == fabs(w(iq))) {
228 >                    a(ip, iq) = 0.0;
229 >                } else if (fabs(a(ip, iq)) > tresh) {
230 >                    h = w(iq) - w(ip);
231 >                    if ( (fabs(h)+g) == fabs(h)) {
232 >                        t = (a(ip, iq)) / h;
233 >                    } else {
234 >                        theta = 0.5*h / (a(ip, iq));
235 >                        t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
236  
237 <        for (ip=0; ip<N; ip++)
238 <        {
239 <            b(ip) += z(ip);
274 <            w(ip) = b(ip);
275 <            z(ip) = 0.0;
276 <        }
277 <    }
237 >                        if (theta < 0.0)
238 >                            t = -t;
239 >                    }
240  
241 +                    c = 1.0 / sqrt(1+t*t);
242 +                    s = t*c;
243 +                    tau = s/(1.0+c);
244 +                    h = t*a(ip, iq);
245 +                    z(ip) -= h;
246 +                    z(iq) += h;
247 +                    w(ip) -= h;
248 +                    w(iq) += h;
249 +                    a(ip, iq)=0.0;
250 +                    
251 +                    for (j=0;j<ip-1;j++)
252 +                        ROT(a,j,ip,j,iq);
253 +
254 +                    for (j=ip+1;j<iq-1;j++)
255 +                        ROT(a,ip,j,j,iq);
256 +
257 +                    for (j=iq+1; j<N; j++)
258 +                        ROT(a,ip,j,iq,j);
259 +                    for (j=0; j<N; j++)
260 +                        ROT(v,j,ip,j,iq);
261 +                }
262 +            }
263 +        }//for (ip=0; ip<2; ip++)
264 +
265 +        for (ip=0; ip<N; ip++) {
266 +            b(ip) += z(ip);
267 +            w(ip) = b(ip);
268 +            z(ip) = 0.0;
269 +        }
270 +        
271 +    } // end for (i=0; i<MAX_ROTATIONS; i++)
272 +
273      if ( i >= MAX_ROTATIONS )
274 <        return false;
274 >        return false;
275  
276      // sort eigenfunctions
277 <    for (j=0; j<N; j++)
278 <    {
279 <        k = j;
280 <        tmp = w(k);
281 <        for (i=j; i<N; i++)
282 <        {
283 <            if (w(i) >= tmp)
284 <            {
285 <                k = i;
286 <                tmp = w(k);
287 <            }
288 <        }
289 <        if (k != j)
290 <        {
291 <            w(k) = w(j);
292 <            w(j) = tmp;
293 <            for (i=0; i<N; i++)
294 <            {
295 <                tmp = v(i, j);
302 <                v(i, j) = v(i, k);
303 <                v(i, k) = tmp;
304 <            }
305 <        }
277 >    for (j=0; j<N; j++) {
278 >        k = j;
279 >        tmp = w(k);
280 >        for (i=j; i<N; i++) {
281 >            if (w(i) >= tmp) {
282 >            k = i;
283 >            tmp = w(k);
284 >            }
285 >        }
286 >    
287 >        if (k != j) {
288 >            w(k) = w(j);
289 >            w(j) = tmp;
290 >            for (i=0; i<N; i++)  {
291 >                tmp = v(i, j);
292 >                v(i, j) = v(i, k);
293 >                v(i, k) = tmp;
294 >            }
295 >        }
296      }
297  
298      //    insure eigenvector consistency (i.e., Jacobi can compute
# Line 311 | Line 301 | void SquareMatrix<Real, int Dim>::jacobi(SquareMatrix<
301      //    hyperstreamline/other stuff. We will select the most
302      //    positive eigenvector.
303      int numPos;
304 <    for (j=0; j<N; j++)
305 <    {
306 <        for (numPos=0, i=0; i<N; i++) if ( v(i, j) >= 0.0 ) numPos++;
317 <        if ( numPos < 2 ) for(i=0; i<N; i++) v(i, j) *= -1.0;
304 >    for (j=0; j<N; j++) {
305 >        for (numPos=0, i=0; i<N; i++) if ( v(i, j) >= 0.0 ) numPos++;
306 >        if ( numPos < 2 ) for(i=0; i<N; i++) v(i, j) *= -1.0;
307      }
308  
309      return true;
# Line 325 | Line 314 | void SquareMatrix<Real, int Dim>::jacobi(SquareMatrix<
314  
315   }
316  
328 –
329 – }
317   #endif //MATH_SQUAREMATRIX_HPP

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines