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

Comparing trunk/OOPSE-3.0/src/math/SquareMatrix.hpp (file contents):
Revision 1569 by tim, Thu Oct 14 23:28:09 2004 UTC vs.
Revision 1594 by tim, Mon Oct 18 23:13:23 2004 UTC

# Line 78 | Line 78 | namespace oopse {
78              return m;
79          }
80  
81 <        /** Retunrs  the inversion of this matrix. */
81 >        /**
82 >         * Retunrs  the inversion of this matrix.
83 >         * @todo
84 >         */
85           SquareMatrix<Real, Dim>  inverse() {
86               SquareMatrix<Real, Dim> result;
87  
88               return result;
89          }        
90  
91 <        /** Returns the determinant of this matrix. */
91 >        /**
92 >         * Returns the determinant of this matrix.
93 >         * @todo
94 >         */
95          double determinant() const {
96              double det;
97              return det;
# Line 142 | Line 148 | namespace oopse {
148              return true;
149          }        
150  
151 +        /** @todo need implement */
152          void diagonalize() {
153 <            jacobi(m, eigenValues, ortMat);
153 >            //jacobi(m, eigenValues, ortMat);
154          }
155  
156          /**
# Line 175 | Line 182 | namespace oopse {
182           * @param w output eigenvalues
183           * @param v output eigenvectors
184           */
185 <        void jacobi(const SquareMatrix<Real, Dim>& a,
179 <                              Vector<Real, Dim>& w,
185 >        bool jacobi(const SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
186                                SquareMatrix<Real, Dim>& v);
187      };//end SquareMatrix
188  
# Line 184 | Line 190 | template<Real, int Dim>
190   #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)
191   #define MAX_ROTATIONS 60
192  
193 < template<Real, int Dim>
194 < void SquareMatrix<Real, int Dim>::jacobi(SquareMatrix<Real, Dim>& a,
195 <                                                                       Vector<Real, Dim>& w,
190 <                                                                       SquareMatrix<Real, Dim>& v) {
193 > template<typename Real, int Dim>
194 > bool SquareMatrix<Real, Dim>::jacobi(const SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
195 >                              SquareMatrix<Real, Dim>& v) {
196      const int N = Dim;                                                                      
197      int i, j, k, iq, ip;
198      double tresh, theta, tau, t, sm, s, h, g, c;
# Line 195 | Line 200 | void SquareMatrix<Real, int Dim>::jacobi(SquareMatrix<
200      Vector<Real, Dim> b, z;
201  
202      // initialize
203 <    for (ip=0; ip<N; ip++)
204 <    {
205 <        for (iq=0; iq<N; iq++) v(ip, iq) = 0.0;
206 <        v(ip, ip) = 1.0;
203 >    for (ip=0; ip<N; ip++) {
204 >        for (iq=0; iq<N; iq++)
205 >            v(ip, iq) = 0.0;
206 >        v(ip, ip) = 1.0;
207      }
208 <    for (ip=0; ip<N; ip++)
209 <    {
210 <        b(ip) = w(ip) = a(ip, ip);
211 <        z(ip) = 0.0;
208 >    
209 >    for (ip=0; ip<N; ip++) {
210 >        b(ip) = w(ip) = a(ip, ip);
211 >        z(ip) = 0.0;
212      }
213  
214      // begin rotation sequence
215 <    for (i=0; i<MAX_ROTATIONS; i++)
216 <    {
217 <        sm = 0.0;
218 <        for (ip=0; ip<2; ip++)
219 <        {
220 <            for (iq=ip+1; iq<N; iq++) sm += fabs(a(ip, iq));
221 <        }
222 <        if (sm == 0.0) break;
215 >    for (i=0; i<MAX_ROTATIONS; i++) {
216 >        sm = 0.0;
217 >        for (ip=0; ip<2; ip++) {
218 >            for (iq=ip+1; iq<N; iq++)
219 >                sm += fabs(a(ip, iq));
220 >        }
221 >        
222 >        if (sm == 0.0)
223 >            break;
224  
225 <        if (i < 4) tresh = 0.2*sm/(9);
226 <        else tresh = 0.0;
225 >        if (i < 4)
226 >            tresh = 0.2*sm/(9);
227 >        else
228 >            tresh = 0.0;
229  
230 <        for (ip=0; ip<2; ip++)
231 <        {
232 <            for (iq=ip+1; iq<N; iq++)
233 <            {
234 <                g = 100.0*fabs(a(ip, iq));
235 <                if (i > 4 && (fabs(w(ip))+g) == fabs(w(ip))
236 <                    && (fabs(w(iq))+g) == fabs(w(iq)))
237 <                {
238 <                    a(ip, iq) = 0.0;
239 <                }
240 <                else if (fabs(a(ip, iq)) > tresh)
241 <                {
242 <                    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 <        }
230 >        for (ip=0; ip<2; ip++) {
231 >            for (iq=ip+1; iq<N; iq++) {
232 >                g = 100.0*fabs(a(ip, iq));
233 >                if (i > 4 && (fabs(w(ip))+g) == fabs(w(ip))
234 >                    && (fabs(w(iq))+g) == fabs(w(iq))) {
235 >                    a(ip, iq) = 0.0;
236 >                } else if (fabs(a(ip, iq)) > tresh) {
237 >                    h = w(iq) - w(ip);
238 >                    if ( (fabs(h)+g) == fabs(h)) {
239 >                        t = (a(ip, iq)) / h;
240 >                    } else {
241 >                        theta = 0.5*h / (a(ip, iq));
242 >                        t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
243  
244 <        for (ip=0; ip<N; ip++)
245 <        {
246 <            b(ip) += z(ip);
274 <            w(ip) = b(ip);
275 <            z(ip) = 0.0;
276 <        }
277 <    }
244 >                        if (theta < 0.0)
245 >                            t = -t;
246 >                    }
247  
248 +                    c = 1.0 / sqrt(1+t*t);
249 +                    s = t*c;
250 +                    tau = s/(1.0+c);
251 +                    h = t*a(ip, iq);
252 +                    z(ip) -= h;
253 +                    z(iq) += h;
254 +                    w(ip) -= h;
255 +                    w(iq) += h;
256 +                    a(ip, iq)=0.0;
257 +                    
258 +                    for (j=0;j<ip-1;j++)
259 +                        ROT(a,j,ip,j,iq);
260 +
261 +                    for (j=ip+1;j<iq-1;j++)
262 +                        ROT(a,ip,j,j,iq);
263 +
264 +                    for (j=iq+1; j<N; j++)
265 +                        ROT(a,ip,j,iq,j);
266 +                    
267 +                    for (j=0; j<N; j++)
268 +                        ROT(v,j,ip,j,iq);
269 +                }
270 +            }
271 +        }//for (ip=0; ip<2; ip++)
272 +
273 +        for (ip=0; ip<N; ip++) {
274 +            b(ip) += z(ip);
275 +            w(ip) = b(ip);
276 +            z(ip) = 0.0;
277 +        }
278 +        
279 +    } // end for (i=0; i<MAX_ROTATIONS; i++)
280 +
281      if ( i >= MAX_ROTATIONS )
282 <        return false;
282 >        return false;
283  
284      // sort eigenfunctions
285 <    for (j=0; j<N; j++)
286 <    {
287 <        k = j;
288 <        tmp = w(k);
289 <        for (i=j; i<N; i++)
290 <        {
291 <            if (w(i) >= tmp)
292 <            {
293 <                k = i;
294 <                tmp = w(k);
295 <            }
296 <        }
297 <        if (k != j)
298 <        {
299 <            w(k) = w(j);
300 <            w(j) = tmp;
301 <            for (i=0; i<N; i++)
302 <            {
303 <                tmp = v(i, j);
302 <                v(i, j) = v(i, k);
303 <                v(i, k) = tmp;
304 <            }
305 <        }
285 >    for (j=0; j<N; j++) {
286 >        k = j;
287 >        tmp = w(k);
288 >        for (i=j; i<N; i++) {
289 >            if (w(i) >= tmp) {
290 >            k = i;
291 >            tmp = w(k);
292 >            }
293 >        }
294 >    
295 >        if (k != j) {
296 >            w(k) = w(j);
297 >            w(j) = tmp;
298 >            for (i=0; i<N; i++)  {
299 >                tmp = v(i, j);
300 >                v(i, j) = v(i, k);
301 >                v(i, k) = tmp;
302 >            }
303 >        }
304      }
305  
306      //    insure eigenvector consistency (i.e., Jacobi can compute
# Line 311 | Line 309 | void SquareMatrix<Real, int Dim>::jacobi(SquareMatrix<
309      //    hyperstreamline/other stuff. We will select the most
310      //    positive eigenvector.
311      int numPos;
312 <    for (j=0; j<N; j++)
313 <    {
314 <        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;
312 >    for (j=0; j<N; j++) {
313 >        for (numPos=0, i=0; i<N; i++) if ( v(i, j) >= 0.0 ) numPos++;
314 >        if ( numPos < 2 ) for(i=0; i<N; i++) v(i, j) *= -1.0;
315      }
316  
317      return true;
# Line 325 | Line 322 | void SquareMatrix<Real, int Dim>::jacobi(SquareMatrix<
322  
323   }
324  
328 –
329 – }
325   #endif //MATH_SQUAREMATRIX_HPP

Diff Legend

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