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

Comparing trunk/OOPSE-4/src/math/SquareMatrix.hpp (file contents):
Revision 1576 by tim, Fri Oct 15 18:18:12 2004 UTC vs.
Revision 1606 by tim, Tue Oct 19 23:01:03 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 need implementation
84 >         */
85           SquareMatrix<Real, Dim>  inverse() {
86               SquareMatrix<Real, Dim> result;
87  
88               return result;
89          }        
90  
91 <        /** Returns the determinant of this matrix. */
92 <        double determinant() const {
93 <            double det;
91 >        /**
92 >         * Returns the determinant of this matrix.
93 >         * @todo need implementation
94 >         */
95 >        Real determinant() const {
96 >            Real det;
97              return det;
98          }
99  
100          /** Returns the trace of this matrix. */
101 <        double trace() const {
102 <           double tmp = 0;
101 >        Real trace() const {
102 >           Real tmp = 0;
103            
104              for (unsigned int i = 0; i < Dim ; i++)
105                  tmp += data_[i][i];
# Line 142 | Line 148 | namespace oopse {
148              return true;
149          }        
150  
151 +        /** @todo need implementation */
152          void diagonalize() {
153 <            jacobi(m, eigenValues, ortMat);
153 >            //jacobi(m, eigenValues, ortMat);
154          }
155  
156          /**
# Line 158 | Line 165 | namespace oopse {
165              SquareMatrix<Real, Dim> ortMat;
166              
167              if ( !isSymmetric()){
168 <                throw();
168 >                //throw();
169              }
170              
171              SquareMatrix<Real, Dim> m(*this);
# Line 175 | Line 182 | namespace oopse {
182           * @param w output eigenvalues
183           * @param v output eigenvectors
184           */
185 <        bool jacobi(const SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
185 >        bool jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
186                                SquareMatrix<Real, Dim>& v);
187      };//end SquareMatrix
188  
# Line 184 | Line 191 | template<typename Real, int Dim>
191   #define MAX_ROTATIONS 60
192  
193   template<typename Real, int Dim>
194 < bool SquareMatrix<Real, Dim>::jacobi(const SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
194 > bool SquareMatrix<Real, Dim>::jacobi(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;
199 <    double tmp;
198 >    Real tresh, theta, tau, t, sm, s, h, g, c;
199 >    Real tmp;
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);
233 <                    if ( (fabs(h)+g) == fabs(h)) t = (a(ip, iq)) / h;
234 <                    else
235 <                    {
236 <                        theta = 0.5*h / (a(ip, iq));
237 <                        t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
238 <                        if (theta < 0.0) t = -t;
239 <                    }
240 <                    c = 1.0 / sqrt(1+t*t);
241 <                    s = t*c;
242 <                    tau = s/(1.0+c);
243 <                    h = t*a(ip, iq);
244 <                    z(ip) -= h;
245 <                    z(iq) += h;
246 <                    w(ip) -= h;
247 <                    w(iq) += h;
248 <                    a(ip, iq)=0.0;
249 <                    for (j=0;j<ip-1;j++)
250 <                    {
251 <                        ROT(a,j,ip,j,iq);
252 <                    }
253 <                    for (j=ip+1;j<iq-1;j++)
254 <                    {
255 <                        ROT(a,ip,j,j,iq);
256 <                    }
257 <                    for (j=iq+1; j<N; j++)
258 <                    {
259 <                        ROT(a,ip,j,iq,j);
260 <                    }
261 <                    for (j=0; j<N; j++)
262 <                    {
263 <                        ROT(v,j,ip,j,iq);
264 <                    }
265 <                }
266 <            }
267 <        }
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);
272 <            w(ip) = b(ip);
273 <            z(ip) = 0.0;
274 <        }
275 <    }
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);
300 <                v(i, j) = v(i, k);
301 <                v(i, k) = tmp;
302 <            }
303 <        }
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 309 | Line 309 | bool SquareMatrix<Real, Dim>::jacobi(const SquareMatri
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++;
315 <        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;

Diff Legend

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