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

Comparing trunk/OOPSE-2.0/src/math/SquareMatrix.hpp (file contents):
Revision 1616 by tim, Wed Oct 20 18:07:08 2004 UTC vs.
Revision 1644 by tim, Mon Oct 25 22:46:19 2004 UTC

# Line 45 | Line 45 | namespace oopse {
45      template<typename Real, int Dim>
46      class SquareMatrix : public RectMatrix<Real, Dim, Dim> {
47          public:
48 +            typedef Real ElemType;
49 +            typedef Real* ElemPoinerType;
50  
51 <        /** default constructor */
52 <        SquareMatrix() {
53 <            for (unsigned int i = 0; i < Dim; i++)
54 <                for (unsigned int j = 0; j < Dim; j++)
55 <                    data_[i][j] = 0.0;
56 <         }
51 >            /** default constructor */
52 >            SquareMatrix() {
53 >                for (unsigned int i = 0; i < Dim; i++)
54 >                    for (unsigned int j = 0; j < Dim; j++)
55 >                        data_[i][j] = 0.0;
56 >             }
57  
58 <        /** copy constructor */
59 <        SquareMatrix(const RectMatrix<Real, Dim, Dim>& m)  : RectMatrix<Real, Dim, Dim>(m) {
60 <        }
59 <        
60 <        /** copy assignment operator */
61 <        SquareMatrix<Real, Dim>& operator =(const RectMatrix<Real, Dim, Dim>& m) {
62 <            RectMatrix<Real, Dim, Dim>::operator=(m);
63 <            return *this;
64 <        }
65 <                              
66 <        /** Retunrs  an identity matrix*/
58 >            /** Constructs and initializes every element of this matrix to a scalar */
59 >            SquareMatrix(Real s) : RectMatrix<Real, Dim, Dim>(s){
60 >            }
61  
62 <       static SquareMatrix<Real, Dim> identity() {
63 <            SquareMatrix<Real, Dim> m;
62 >            /** Constructs and initializes from an array */
63 >            SquareMatrix(Real* array) : RectMatrix<Real, Dim, Dim>(array){
64 >            }
65 >
66 >
67 >            /** copy constructor */
68 >            SquareMatrix(const RectMatrix<Real, Dim, Dim>& m) : RectMatrix<Real, Dim, Dim>(m) {
69 >            }
70              
71 <            for (unsigned int i = 0; i < Dim; i++)
72 <                for (unsigned int j = 0; j < Dim; j++)
73 <                    if (i == j)
74 <                        m(i, j) = 1.0;
75 <                    else
76 <                        m(i, j) = 0.0;
71 >            /** copy assignment operator */
72 >            SquareMatrix<Real, Dim>& operator =(const RectMatrix<Real, Dim, Dim>& m) {
73 >                RectMatrix<Real, Dim, Dim>::operator=(m);
74 >                return *this;
75 >            }
76 >                                  
77 >            /** Retunrs  an identity matrix*/
78  
79 <            return m;
80 <        }
79 >           static SquareMatrix<Real, Dim> identity() {
80 >                SquareMatrix<Real, Dim> m;
81 >                
82 >                for (unsigned int i = 0; i < Dim; i++)
83 >                    for (unsigned int j = 0; j < Dim; j++)
84 >                        if (i == j)
85 >                            m(i, j) = 1.0;
86 >                        else
87 >                            m(i, j) = 0.0;
88  
89 <        /**
90 <         * Retunrs  the inversion of this matrix.
83 <         * @todo need implementation
84 <         */
85 <         SquareMatrix<Real, Dim>  inverse() {
86 <             SquareMatrix<Real, Dim> result;
89 >                return m;
90 >            }
91  
92 <             return result;
93 <        }        
92 >            /**
93 >             * Retunrs  the inversion of this matrix.
94 >             * @todo need implementation
95 >             */
96 >             SquareMatrix<Real, Dim>  inverse() {
97 >                 SquareMatrix<Real, Dim> result;
98  
99 <        /**
100 <         * Returns the determinant of this matrix.
93 <         * @todo need implementation
94 <         */
95 <        Real determinant() const {
96 <            Real det;
97 <            return det;
98 <        }
99 >                 return result;
100 >            }        
101  
102 <        /** Returns the trace of this matrix. */
103 <        Real trace() const {
104 <           Real tmp = 0;
105 <          
106 <            for (unsigned int i = 0; i < Dim ; i++)
107 <                tmp += data_[i][i];
102 >            /**
103 >             * Returns the determinant of this matrix.
104 >             * @todo need implementation
105 >             */
106 >            Real determinant() const {
107 >                Real det;
108 >                return det;
109 >            }
110  
111 <            return tmp;
112 <        }
111 >            /** Returns the trace of this matrix. */
112 >            Real trace() const {
113 >               Real tmp = 0;
114 >              
115 >                for (unsigned int i = 0; i < Dim ; i++)
116 >                    tmp += data_[i][i];
117  
118 <        /** Tests if this matrix is symmetrix. */            
119 <        bool isSymmetric() const {
112 <            for (unsigned int i = 0; i < Dim - 1; i++)
113 <                for (unsigned int j = i; j < Dim; j++)
114 <                    if (fabs(data_[i][j] - data_[j][i]) > oopse::epsilon)
115 <                        return false;
116 <                    
117 <            return true;
118 <        }
118 >                return tmp;
119 >            }
120  
121 <        /** Tests if this matrix is orthogonal. */            
122 <        bool isOrthogonal() {
123 <            SquareMatrix<Real, Dim> tmp;
121 >            /** Tests if this matrix is symmetrix. */            
122 >            bool isSymmetric() const {
123 >                for (unsigned int i = 0; i < Dim - 1; i++)
124 >                    for (unsigned int j = i; j < Dim; j++)
125 >                        if (fabs(data_[i][j] - data_[j][i]) > oopse::epsilon)
126 >                            return false;
127 >                        
128 >                return true;
129 >            }
130  
131 <            tmp = *this * transpose();
131 >            /** Tests if this matrix is orthogonal. */            
132 >            bool isOrthogonal() {
133 >                SquareMatrix<Real, Dim> tmp;
134  
135 <            return tmp.isDiagonal();
127 <        }
135 >                tmp = *this * transpose();
136  
137 <        /** Tests if this matrix is diagonal. */
138 <        bool isDiagonal() const {
131 <            for (unsigned int i = 0; i < Dim ; i++)
132 <                for (unsigned int j = 0; j < Dim; j++)
133 <                    if (i !=j && fabs(data_[i][j]) > oopse::epsilon)
134 <                        return false;
135 <                    
136 <            return true;
137 <        }
137 >                return tmp.isDiagonal();
138 >            }
139  
140 <        /** Tests if this matrix is the unit matrix. */
141 <        bool isUnitMatrix() const {
142 <            if (!isDiagonal())
143 <                return false;
144 <            
145 <            for (unsigned int i = 0; i < Dim ; i++)
146 <                if (fabs(data_[i][i] - 1) > oopse::epsilon)
140 >            /** Tests if this matrix is diagonal. */
141 >            bool isDiagonal() const {
142 >                for (unsigned int i = 0; i < Dim ; i++)
143 >                    for (unsigned int j = 0; j < Dim; j++)
144 >                        if (i !=j && fabs(data_[i][j]) > oopse::epsilon)
145 >                            return false;
146 >                        
147 >                return true;
148 >            }
149 >
150 >            /** Tests if this matrix is the unit matrix. */
151 >            bool isUnitMatrix() const {
152 >                if (!isDiagonal())
153                      return false;
154                  
155 <            return true;
156 <        }        
155 >                for (unsigned int i = 0; i < Dim ; i++)
156 >                    if (fabs(data_[i][i] - 1) > oopse::epsilon)
157 >                        return false;
158 >                    
159 >                return true;
160 >            }        
161  
162 <        /** @todo need implementation */
163 <        void diagonalize() {
164 <            //jacobi(m, eigenValues, ortMat);
165 <        }
162 >            /** @todo need implementation */
163 >            void diagonalize() {
164 >                //jacobi(m, eigenValues, ortMat);
165 >            }
166  
167 <        /**
168 <         * Jacobi iteration routines for computing eigenvalues/eigenvectors of
169 <         * real symmetric matrix
170 <         *
171 <         * @return true if success, otherwise return false
172 <         * @param a symmetric matrix whose eigenvectors are to be computed. On return, the matrix is
173 <         *     overwritten
174 <         * @param w will contain the eigenvalues of the matrix On return of this function
175 <         * @param v the columns of this matrix will contain the eigenvectors. The eigenvectors are
176 <         *    normalized and mutually orthogonal.
177 <         */
178 <      
179 <        static int jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& d,
180 <                              SquareMatrix<Real, Dim>& v);
167 >            /**
168 >             * Jacobi iteration routines for computing eigenvalues/eigenvectors of
169 >             * real symmetric matrix
170 >             *
171 >             * @return true if success, otherwise return false
172 >             * @param a symmetric matrix whose eigenvectors are to be computed. On return, the matrix is
173 >             *     overwritten
174 >             * @param w will contain the eigenvalues of the matrix On return of this function
175 >             * @param v the columns of this matrix will contain the eigenvectors. The eigenvectors are
176 >             *    normalized and mutually orthogonal.
177 >             */
178 >          
179 >            static int jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& d,
180 >                                  SquareMatrix<Real, Dim>& v);
181      };//end SquareMatrix
182  
183  
# Line 197 | Line 208 | namespace oopse {
208      // normalized.
209      template<typename Real, int Dim>
210      int SquareMatrix<Real, Dim>::jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
211 <                                  SquareMatrix<Real, Dim>& v) {
212 <      const int n = Dim;  
213 <      int i, j, k, iq, ip, numPos;
214 <      Real tresh, theta, tau, t, sm, s, h, g, c, tmp;
215 <      Real bspace[4], zspace[4];
216 <      Real *b = bspace;
217 <      Real *z = zspace;
211 >                                        SquareMatrix<Real, Dim>& v) {
212 >        const int n = Dim;  
213 >        int i, j, k, iq, ip, numPos;
214 >        Real tresh, theta, tau, t, sm, s, h, g, c, tmp;
215 >        Real bspace[4], zspace[4];
216 >        Real *b = bspace;
217 >        Real *z = zspace;
218  
219 <      // only allocate memory if the matrix is large
220 <      if (n > 4)
221 <        {
222 <        b = new Real[n];
212 <        z = new Real[n];
219 >        // only allocate memory if the matrix is large
220 >        if (n > 4) {
221 >            b = new Real[n];
222 >            z = new Real[n];
223          }
224  
225 <      // initialize
226 <      for (ip=0; ip<n; ip++)
227 <        {
228 <        for (iq=0; iq<n; iq++)
229 <          {
230 <          v(ip, iq) = 0.0;
221 <          }
222 <        v(ip, ip) = 1.0;
225 >        // initialize
226 >        for (ip=0; ip<n; ip++) {
227 >            for (iq=0; iq<n; iq++) {
228 >                v(ip, iq) = 0.0;
229 >            }
230 >            v(ip, ip) = 1.0;
231          }
232 <      for (ip=0; ip<n; ip++)
233 <        {
234 <        b[ip] = w[ip] = a(ip, ip);
227 <        z[ip] = 0.0;
232 >        for (ip=0; ip<n; ip++) {
233 >            b[ip] = w[ip] = a(ip, ip);
234 >            z[ip] = 0.0;
235          }
236  
237 <      // begin rotation sequence
238 <      for (i=0; i<VTK_MAX_ROTATIONS; i++)
239 <        {
240 <        sm = 0.0;
241 <        for (ip=0; ip<n-1; ip++)
242 <          {
243 <          for (iq=ip+1; iq<n; iq++)
237 <            {
238 <            sm += fabs(a(ip, iq));
237 >        // begin rotation sequence
238 >        for (i=0; i<VTK_MAX_ROTATIONS; i++) {
239 >            sm = 0.0;
240 >            for (ip=0; ip<n-1; ip++) {
241 >                for (iq=ip+1; iq<n; iq++) {
242 >                    sm += fabs(a(ip, iq));
243 >                }
244              }
245 <          }
246 <        if (sm == 0.0)
247 <          {
243 <          break;
244 <          }
245 >            if (sm == 0.0) {
246 >                break;
247 >            }
248  
249 <        if (i < 3)                                // first 3 sweeps
250 <          {
251 <          tresh = 0.2*sm/(n*n);
252 <          }
253 <        else
251 <          {
252 <          tresh = 0.0;
253 <          }
249 >            if (i < 3) {                                // first 3 sweeps
250 >                tresh = 0.2*sm/(n*n);
251 >            } else {
252 >                tresh = 0.0;
253 >            }
254  
255 <        for (ip=0; ip<n-1; ip++)
256 <          {
257 <          for (iq=ip+1; iq<n; iq++)
258 <            {
259 <            g = 100.0*fabs(a(ip, iq));
255 >            for (ip=0; ip<n-1; ip++) {
256 >                for (iq=ip+1; iq<n; iq++) {
257 >                    g = 100.0*fabs(a(ip, iq));
258  
259 <            // after 4 sweeps
260 <            if (i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
261 <            && (fabs(w[iq])+g) == fabs(w[iq]))
262 <              {
263 <              a(ip, iq) = 0.0;
264 <              }
265 <            else if (fabs(a(ip, iq)) > tresh)
266 <              {
267 <              h = w[iq] - w[ip];
268 <              if ( (fabs(h)+g) == fabs(h))
269 <                {
270 <                t = (a(ip, iq)) / h;
271 <                }
272 <              else
273 <                {
274 <                theta = 0.5*h / (a(ip, iq));
275 <                t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
276 <                if (theta < 0.0)
277 <                  {
278 <                  t = -t;
279 <                  }
280 <                }
281 <              c = 1.0 / sqrt(1+t*t);
282 <              s = t*c;
285 <              tau = s/(1.0+c);
286 <              h = t*a(ip, iq);
287 <              z[ip] -= h;
288 <              z[iq] += h;
289 <              w[ip] -= h;
290 <              w[iq] += h;
291 <              a(ip, iq)=0.0;
259 >                    // after 4 sweeps
260 >                    if (i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
261 >                        && (fabs(w[iq])+g) == fabs(w[iq])) {
262 >                        a(ip, iq) = 0.0;
263 >                    } else if (fabs(a(ip, iq)) > tresh) {
264 >                        h = w[iq] - w[ip];
265 >                        if ( (fabs(h)+g) == fabs(h)) {
266 >                            t = (a(ip, iq)) / h;
267 >                        } else {
268 >                            theta = 0.5*h / (a(ip, iq));
269 >                            t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
270 >                            if (theta < 0.0) {
271 >                                t = -t;
272 >                            }
273 >                        }
274 >                        c = 1.0 / sqrt(1+t*t);
275 >                        s = t*c;
276 >                        tau = s/(1.0+c);
277 >                        h = t*a(ip, iq);
278 >                        z[ip] -= h;
279 >                        z[iq] += h;
280 >                        w[ip] -= h;
281 >                        w[iq] += h;
282 >                        a(ip, iq)=0.0;
283  
284 <              // ip already shifted left by 1 unit
285 <              for (j = 0;j <= ip-1;j++)
286 <                {
287 <                VTK_ROTATE(a,j,ip,j,iq);
284 >                        // ip already shifted left by 1 unit
285 >                        for (j = 0;j <= ip-1;j++) {
286 >                            VTK_ROTATE(a,j,ip,j,iq);
287 >                        }
288 >                        // ip and iq already shifted left by 1 unit
289 >                        for (j = ip+1;j <= iq-1;j++) {
290 >                            VTK_ROTATE(a,ip,j,j,iq);
291 >                        }
292 >                        // iq already shifted left by 1 unit
293 >                        for (j=iq+1; j<n; j++) {
294 >                            VTK_ROTATE(a,ip,j,iq,j);
295 >                        }
296 >                        for (j=0; j<n; j++) {
297 >                            VTK_ROTATE(v,j,ip,j,iq);
298 >                        }
299 >                    }
300                  }
298 –              // ip and iq already shifted left by 1 unit
299 –              for (j = ip+1;j <= iq-1;j++)
300 –                {
301 –                VTK_ROTATE(a,ip,j,j,iq);
302 –                }
303 –              // iq already shifted left by 1 unit
304 –              for (j=iq+1; j<n; j++)
305 –                {
306 –                VTK_ROTATE(a,ip,j,iq,j);
307 –                }
308 –              for (j=0; j<n; j++)
309 –                {
310 –                VTK_ROTATE(v,j,ip,j,iq);
311 –                }
312 –              }
301              }
314 –          }
302  
303 <        for (ip=0; ip<n; ip++)
304 <          {
305 <          b[ip] += z[ip];
306 <          w[ip] = b[ip];
307 <          z[ip] = 0.0;
321 <          }
303 >            for (ip=0; ip<n; ip++) {
304 >                b[ip] += z[ip];
305 >                w[ip] = b[ip];
306 >                z[ip] = 0.0;
307 >            }
308          }
309  
310 <      //// this is NEVER called
311 <      if ( i >= VTK_MAX_ROTATIONS )
312 <        {
313 <           std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl;
328 <           return 0;
310 >        //// this is NEVER called
311 >        if ( i >= VTK_MAX_ROTATIONS ) {
312 >            std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl;
313 >            return 0;
314          }
315  
316 <      // sort eigenfunctions                 these changes do not affect accuracy
317 <      for (j=0; j<n-1; j++)                  // boundary incorrect
318 <        {
334 <        k = j;
335 <        tmp = w[k];
336 <        for (i=j+1; i<n; i++)                // boundary incorrect, shifted already
337 <          {
338 <          if (w[i] >= tmp)                   // why exchage if same?
339 <            {
340 <            k = i;
316 >        // sort eigenfunctions                 these changes do not affect accuracy
317 >        for (j=0; j<n-1; j++) {                  // boundary incorrect
318 >            k = j;
319              tmp = w[k];
320 +            for (i=j+1; i<n; i++) {                // boundary incorrect, shifted already
321 +                if (w[i] >= tmp) {                   // why exchage if same?
322 +                    k = i;
323 +                    tmp = w[k];
324 +                }
325              }
326 <          }
327 <        if (k != j)
328 <          {
329 <          w[k] = w[j];
330 <          w[j] = tmp;
331 <          for (i=0; i<n; i++)
332 <            {
333 <            tmp = v(i, j);
351 <            v(i, j) = v(i, k);
352 <            v(i, k) = tmp;
326 >            if (k != j) {
327 >                w[k] = w[j];
328 >                w[j] = tmp;
329 >                for (i=0; i<n; i++) {
330 >                    tmp = v(i, j);
331 >                    v(i, j) = v(i, k);
332 >                    v(i, k) = tmp;
333 >                }
334              }
354 –          }
335          }
336 <      // insure eigenvector consistency (i.e., Jacobi can compute vectors that
337 <      // are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
338 <      // reek havoc in hyperstreamline/other stuff. We will select the most
339 <      // positive eigenvector.
340 <      int ceil_half_n = (n >> 1) + (n & 1);
341 <      for (j=0; j<n; j++)
342 <        {
343 <        for (numPos=0, i=0; i<n; i++)
344 <          {
345 <          if ( v(i, j) >= 0.0 )
366 <            {
367 <            numPos++;
336 >        // insure eigenvector consistency (i.e., Jacobi can compute vectors that
337 >        // are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
338 >        // reek havoc in hyperstreamline/other stuff. We will select the most
339 >        // positive eigenvector.
340 >        int ceil_half_n = (n >> 1) + (n & 1);
341 >        for (j=0; j<n; j++) {
342 >            for (numPos=0, i=0; i<n; i++) {
343 >                if ( v(i, j) >= 0.0 ) {
344 >                    numPos++;
345 >                }
346              }
347 <          }
348 <    //    if ( numPos < ceil(double(n)/double(2.0)) )
349 <        if ( numPos < ceil_half_n)
350 <          {
351 <          for(i=0; i<n; i++)
374 <            {
375 <            v(i, j) *= -1.0;
347 >            //    if ( numPos < ceil(double(n)/double(2.0)) )
348 >            if ( numPos < ceil_half_n) {
349 >                for (i=0; i<n; i++) {
350 >                    v(i, j) *= -1.0;
351 >                }
352              }
377 –          }
353          }
354  
355 <      if (n > 4)
356 <        {
357 <        delete [] b;
383 <        delete [] z;
355 >        if (n > 4) {
356 >            delete [] b;
357 >            delete [] z;
358          }
359 <      return 1;
359 >        return 1;
360      }
361  
362  

Diff Legend

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