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 74 by tim, Wed Oct 13 23:53:40 2004 UTC vs.
Revision 1924 by gezelter, Mon Aug 5 21:46:11 2013 UTC

# Line 1 | Line 1
1   /*
2 < * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
3 < *
4 < * Contact: oopse@oopse.org
5 < *
6 < * This program is free software; you can redistribute it and/or
7 < * modify it under the terms of the GNU Lesser General Public License
8 < * as published by the Free Software Foundation; either version 2.1
9 < * of the License, or (at your option) any later version.
10 < * All we ask is that proper credit is given for our work, which includes
11 < * - but is not limited to - adding the above copyright notice to the beginning
12 < * of your source code files, and to any copyright notice that you may distribute
13 < * with programs based on this work.
14 < *
15 < * This program is distributed in the hope that it will be useful,
16 < * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 < * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18 < * GNU Lesser General Public License for more details.
19 < *
20 < * You should have received a copy of the GNU Lesser General Public License
21 < * along with this program; if not, write to the Free Software
22 < * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
2 > * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4 + * The University of Notre Dame grants you ("Licensee") a
5 + * non-exclusive, royalty free, license to use, modify and
6 + * redistribute this software in source and binary code form, provided
7 + * that the following conditions are met:
8 + *
9 + * 1. Redistributions of source code must retain the above copyright
10 + *    notice, this list of conditions and the following disclaimer.
11 + *
12 + * 2. Redistributions in binary form must reproduce the above copyright
13 + *    notice, this list of conditions and the following disclaimer in the
14 + *    documentation and/or other materials provided with the
15 + *    distribution.
16 + *
17 + * This software is provided "AS IS," without a warranty of any
18 + * kind. All express or implied conditions, representations and
19 + * warranties, including any implied warranty of merchantability,
20 + * fitness for a particular purpose or non-infringement, are hereby
21 + * excluded.  The University of Notre Dame and its licensors shall not
22 + * be liable for any damages suffered by licensee as a result of
23 + * using, modifying or distributing the software or its
24 + * derivatives. In no event will the University of Notre Dame or its
25 + * licensors be liable for any lost revenue, profit or data, or for
26 + * direct, indirect, special, consequential, incidental or punitive
27 + * damages, however caused and regardless of the theory of liability,
28 + * arising out of the use of or inability to use software, even if the
29 + * University of Notre Dame has been advised of the possibility of
30 + * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39 + * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 + * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42 <
42 >
43   /**
44   * @file SquareMatrix.hpp
45   * @author Teng Lin
# Line 33 | Line 50
50   #define MATH_SQUAREMATRIX_HPP
51  
52   #include "math/RectMatrix.hpp"
53 + #include "utils/NumericConstant.hpp"
54  
55 < namespace oopse {
55 > namespace OpenMD {
56  
57 <    /**
58 <     * @class SquareMatrix SquareMatrix.hpp "math/SquareMatrix.hpp"
59 <     * @brief A square matrix class
60 <     * @template Real the element type
61 <     * @template Dim the dimension of the square matrix
62 <     */
63 <    template<typename Real, int Dim>
64 <    class SquareMatrix : public RectMatrix<Real, Dim, Dim> {
65 <        public:
57 >  /**
58 >   * @class SquareMatrix SquareMatrix.hpp "math/SquareMatrix.hpp"
59 >   * @brief A square matrix class
60 >   * \tparam Real the element type
61 >   * \tparam Dim the dimension of the square matrix
62 >   */
63 >  template<typename Real, int Dim>
64 >  class SquareMatrix : public RectMatrix<Real, Dim, Dim> {
65 >  public:
66 >    typedef Real ElemType;
67 >    typedef Real* ElemPoinerType;
68  
69 <        /** default constructor */
70 <        SquareMatrix() {
71 <            for (unsigned int i = 0; i < Dim; i++)
72 <                for (unsigned int j = 0; j < Dim; j++)
73 <                    data_[i][j] = 0.0;
74 <         }
69 >    /** default constructor */
70 >    SquareMatrix() {
71 >      for (unsigned int i = 0; i < Dim; i++)
72 >        for (unsigned int j = 0; j < Dim; j++)
73 >          this->data_[i][j] = 0.0;
74 >    }
75  
76 <        /** copy constructor */
77 <        SquareMatrix(const RectMatrix<Real, Dim, Dim>& m)  : RectMatrix<Real, Dim, Dim>(m) {
78 <        }
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*/
76 >    /** Constructs and initializes every element of this matrix to a scalar */
77 >    SquareMatrix(Real s) : RectMatrix<Real, Dim, Dim>(s){
78 >    }
79  
80 <       static SquareMatrix<Real, Dim> identity() {
81 <            SquareMatrix<Real, Dim> m;
82 <            
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;
80 >    /** Constructs and initializes from an array */
81 >    SquareMatrix(Real* array) : RectMatrix<Real, Dim, Dim>(array){
82 >    }
83  
78            return m;
79        }
84  
85 <        /** Retunrs  the inversion of this matrix. */
86 <         SquareMatrix<Real, Dim>  inverse() {
87 <             SquareMatrix<Real, Dim> result;
85 >    /** copy constructor */
86 >    SquareMatrix(const RectMatrix<Real, Dim, Dim>& m) : RectMatrix<Real, Dim, Dim>(m) {
87 >    }
88 >            
89 >    /** copy assignment operator */
90 >    SquareMatrix<Real, Dim>& operator =(const RectMatrix<Real, Dim, Dim>& m) {
91 >      RectMatrix<Real, Dim, Dim>::operator=(m);
92 >      return *this;
93 >    }
94 >                                  
95 >    /** Retunrs  an identity matrix*/
96  
97 <             return result;
98 <        }
97 >    static SquareMatrix<Real, Dim> identity() {
98 >      SquareMatrix<Real, Dim> m;
99 >                
100 >      for (unsigned int i = 0; i < Dim; i++)
101 >        for (unsigned int j = 0; j < Dim; j++)
102 >          if (i == j)
103 >            m(i, j) = 1.0;
104 >          else
105 >            m(i, j) = 0.0;
106  
107 <        
107 >      return m;
108 >    }
109  
110 <        /** Returns the determinant of this matrix. */
111 <        double determinant() const {
112 <            double det;
113 <            return det;
114 <        }
110 >    /**
111 >     * Retunrs  the inversion of this matrix.
112 >     * @todo need implementation
113 >     */
114 >    SquareMatrix<Real, Dim>  inverse() {
115 >      SquareMatrix<Real, Dim> result;
116  
117 <        /** Returns the trace of this matrix. */
118 <        double trace() const {
98 <           double tmp = 0;
99 <          
100 <            for (unsigned int i = 0; i < Dim ; i++)
101 <                tmp += data_[i][i];
117 >      return result;
118 >    }        
119  
120 <            return tmp;
121 <        }
120 >    /**
121 >     * Returns the determinant of this matrix.
122 >     * @todo need implementation
123 >     */
124 >    Real determinant() const {
125 >      Real det;
126 >      return det;
127 >    }
128 >    
129 >    /** Returns the trace of this matrix. */
130 >    Real trace() const {
131 >      Real tmp = 0;
132 >              
133 >      for (unsigned int i = 0; i < Dim ; i++)
134 >        tmp += this->data_[i][i];
135  
136 <        /** Tests if this matrix is symmetrix. */            
137 <        bool isSymmetric() const {
108 <            for (unsigned int i = 0; i < Dim - 1; i++)
109 <                for (unsigned int j = i; j < Dim; j++)
110 <                    if (fabs(data_[i][j] - data_[j][i]) > oopse::epsilon)
111 <                        return false;
112 <                    
113 <            return true;
114 <        }
136 >      return tmp;
137 >    }
138  
139 <        /** Tests if this matrix is orthogona. */            
140 <        bool isOrthogonal() {
141 <            SquareMatrix<Real, Dim> tmp;
139 >    /** Tests if this matrix is symmetrix. */            
140 >    bool isSymmetric() const {
141 >      for (unsigned int i = 0; i < Dim - 1; i++)
142 >        for (unsigned int j = i; j < Dim; j++)
143 >          if (fabs(this->data_[i][j] - this->data_[j][i]) > epsilon)
144 >            return false;
145 >                        
146 >      return true;
147 >    }
148  
149 <            tmp = *this * transpose();
149 >    /** Tests if this matrix is orthogonal. */            
150 >    bool isOrthogonal() {
151 >      SquareMatrix<Real, Dim> tmp;
152  
153 <            return tmp.isUnitMatrix();
123 <        }
153 >      tmp = *this * transpose();
154  
155 <        /** Tests if this matrix is diagonal. */
156 <        bool isDiagonal() const {
127 <            for (unsigned int i = 0; i < Dim ; i++)
128 <                for (unsigned int j = 0; j < Dim; j++)
129 <                    if (i !=j && fabs(data_[i][j]) > oopse::epsilon)
130 <                        return false;
131 <                    
132 <            return true;
133 <        }
155 >      return tmp.isDiagonal();
156 >    }
157  
158 <        /** Tests if this matrix is the unit matrix. */
159 <        bool isUnitMatrix() const {
160 <            if (!isDiagonal())
161 <                return false;
162 <            
163 <            for (unsigned int i = 0; i < Dim ; i++)
164 <                if (fabs(data_[i][i] - 1) > oopse::epsilon)
165 <                    return false;
166 <                
144 <            return true;
145 <        }        
146 <
147 <    };//end SquareMatrix
158 >    /** Tests if this matrix is diagonal. */
159 >    bool isDiagonal() const {
160 >      for (unsigned int i = 0; i < Dim ; i++)
161 >        for (unsigned int j = 0; j < Dim; j++)
162 >          if (i !=j && fabs(this->data_[i][j]) > epsilon)
163 >            return false;
164 >                        
165 >      return true;
166 >    }
167  
168 +    /**
169 +     * Returns a column vector that contains the elements from the
170 +     * diagonal of m in the order R(0) = m(0,0), R(1) = m(1,1), and so
171 +     * on.
172 +     */
173 +    Vector<Real, Dim> diagonals() const {
174 +      Vector<Real, Dim> result;
175 +      for (unsigned int i = 0; i < Dim; i++) {
176 +        result(i) = this->data_[i][i];
177 +      }
178 +      return result;
179 +    }
180 +
181 +    /** Tests if this matrix is the unit matrix. */
182 +    bool isUnitMatrix() const {
183 +      if (!isDiagonal())
184 +        return false;
185 +                
186 +      for (unsigned int i = 0; i < Dim ; i++)
187 +        if (fabs(this->data_[i][i] - 1) > epsilon)
188 +          return false;
189 +                    
190 +      return true;
191 +    }        
192 +
193 +    /** Return the transpose of this matrix */
194 +    SquareMatrix<Real,  Dim> transpose() const{
195 +      SquareMatrix<Real,  Dim> result;
196 +                
197 +      for (unsigned int i = 0; i < Dim; i++)
198 +        for (unsigned int j = 0; j < Dim; j++)              
199 +          result(j, i) = this->data_[i][j];
200 +
201 +      return result;
202 +    }
203 +            
204 +    /** @todo need implementation */
205 +    void diagonalize() {
206 +      //jacobi(m, eigenValues, ortMat);
207 +    }
208 +
209 +    /**
210 +     * Jacobi iteration routines for computing eigenvalues/eigenvectors of
211 +     * real symmetric matrix
212 +     *
213 +     * @return true if success, otherwise return false
214 +     * @param a symmetric matrix whose eigenvectors are to be computed. On return, the matrix is
215 +     *     overwritten
216 +     * @param d will contain the eigenvalues of the matrix On return of this function
217 +     * @param v the columns of this matrix will contain the eigenvectors. The eigenvectors are
218 +     *    normalized and mutually orthogonal.
219 +     */
220 +          
221 +    static int jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& d,
222 +                      SquareMatrix<Real, Dim>& v);
223 +  };//end SquareMatrix
224 +
225 +
226 +  /*=========================================================================
227 +
228 +  Program:   Visualization Toolkit
229 +  Module:    $RCSfile: SquareMatrix.hpp,v $
230 +
231 +  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
232 +  All rights reserved.
233 +  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
234 +
235 +  This software is distributed WITHOUT ANY WARRANTY; without even
236 +  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
237 +  PURPOSE.  See the above copyright notice for more information.
238 +
239 +  =========================================================================*/
240 +
241 + #define VTK_ROTATE(a,i,j,k,l) g=a(i, j);h=a(k, l);a(i, j)=g-s*(h+g*tau); \
242 +    a(k, l)=h+s*(g-h*tau)
243 +
244 + #define VTK_MAX_ROTATIONS 20
245 +
246 +  // Jacobi iteration for the solution of eigenvectors/eigenvalues of a nxn
247 +  // real symmetric matrix. Square nxn matrix a; size of matrix in n;
248 +  // output eigenvalues in w; and output eigenvectors in v. Resulting
249 +  // eigenvalues/vectors are sorted in decreasing order; eigenvectors are
250 +  // normalized.
251 +  template<typename Real, int Dim>
252 +  int SquareMatrix<Real, Dim>::jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
253 +                                      SquareMatrix<Real, Dim>& v) {
254 +    const int n = Dim;  
255 +    int i, j, k, iq, ip, numPos;
256 +    Real tresh, theta, tau, t, sm, s, h, g, c, tmp;
257 +    Real bspace[4], zspace[4];
258 +    Real *b = bspace;
259 +    Real *z = zspace;
260 +
261 +    // only allocate memory if the matrix is large
262 +    if (n > 4) {
263 +      b = new Real[n];
264 +      z = new Real[n];
265 +    }
266 +
267 +    // initialize
268 +    for (ip=0; ip<n; ip++) {
269 +      for (iq=0; iq<n; iq++) {
270 +        v(ip, iq) = 0.0;
271 +      }
272 +      v(ip, ip) = 1.0;
273 +    }
274 +    for (ip=0; ip<n; ip++) {
275 +      b[ip] = w[ip] = a(ip, ip);
276 +      z[ip] = 0.0;
277 +    }
278 +
279 +    // begin rotation sequence
280 +    for (i=0; i<VTK_MAX_ROTATIONS; i++) {
281 +      sm = 0.0;
282 +      for (ip=0; ip<n-1; ip++) {
283 +        for (iq=ip+1; iq<n; iq++) {
284 +          sm += fabs(a(ip, iq));
285 +        }
286 +      }
287 +      if (sm == 0.0) {
288 +        break;
289 +      }
290 +
291 +      if (i < 3) {                                // first 3 sweeps
292 +        tresh = 0.2*sm/(n*n);
293 +      } else {
294 +        tresh = 0.0;
295 +      }
296 +
297 +      for (ip=0; ip<n-1; ip++) {
298 +        for (iq=ip+1; iq<n; iq++) {
299 +          g = 100.0*fabs(a(ip, iq));
300 +
301 +          // after 4 sweeps
302 +          if (i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
303 +              && (fabs(w[iq])+g) == fabs(w[iq])) {
304 +            a(ip, iq) = 0.0;
305 +          } else if (fabs(a(ip, iq)) > tresh) {
306 +            h = w[iq] - w[ip];
307 +            if ( (fabs(h)+g) == fabs(h)) {
308 +              t = (a(ip, iq)) / h;
309 +            } else {
310 +              theta = 0.5*h / (a(ip, iq));
311 +              t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
312 +              if (theta < 0.0) {
313 +                t = -t;
314 +              }
315 +            }
316 +            c = 1.0 / sqrt(1+t*t);
317 +            s = t*c;
318 +            tau = s/(1.0+c);
319 +            h = t*a(ip, iq);
320 +            z[ip] -= h;
321 +            z[iq] += h;
322 +            w[ip] -= h;
323 +            w[iq] += h;
324 +            a(ip, iq)=0.0;
325 +
326 +            // ip already shifted left by 1 unit
327 +            for (j = 0;j <= ip-1;j++) {
328 +              VTK_ROTATE(a,j,ip,j,iq);
329 +            }
330 +            // ip and iq already shifted left by 1 unit
331 +            for (j = ip+1;j <= iq-1;j++) {
332 +              VTK_ROTATE(a,ip,j,j,iq);
333 +            }
334 +            // iq already shifted left by 1 unit
335 +            for (j=iq+1; j<n; j++) {
336 +              VTK_ROTATE(a,ip,j,iq,j);
337 +            }
338 +            for (j=0; j<n; j++) {
339 +              VTK_ROTATE(v,j,ip,j,iq);
340 +            }
341 +          }
342 +        }
343 +      }
344 +
345 +      for (ip=0; ip<n; ip++) {
346 +        b[ip] += z[ip];
347 +        w[ip] = b[ip];
348 +        z[ip] = 0.0;
349 +      }
350 +    }
351 +
352 +    //// this is NEVER called
353 +    if ( i >= VTK_MAX_ROTATIONS ) {
354 +      std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl;
355 +      if (n > 4) {
356 +        delete[] b;
357 +        delete[] z;
358 +      }      
359 +      return 0;
360 +    }
361 +
362 +    // sort eigenfunctions                 these changes do not affect accuracy
363 +    for (j=0; j<n-1; j++) {                  // boundary incorrect
364 +      k = j;
365 +      tmp = w[k];
366 +      for (i=j+1; i<n; i++) {                // boundary incorrect, shifted already
367 +        if (w[i] >= tmp) {                   // why exchage if same?
368 +          k = i;
369 +          tmp = w[k];
370 +        }
371 +      }
372 +      if (k != j) {
373 +        w[k] = w[j];
374 +        w[j] = tmp;
375 +        for (i=0; i<n; i++) {
376 +          tmp = v(i, j);
377 +          v(i, j) = v(i, k);
378 +          v(i, k) = tmp;
379 +        }
380 +      }
381 +    }
382 +    // insure eigenvector consistency (i.e., Jacobi can compute vectors that
383 +    // are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
384 +    // reek havoc in hyperstreamline/other stuff. We will select the most
385 +    // positive eigenvector.
386 +    int ceil_half_n = (n >> 1) + (n & 1);
387 +    for (j=0; j<n; j++) {
388 +      for (numPos=0, i=0; i<n; i++) {
389 +        if ( v(i, j) >= 0.0 ) {
390 +          numPos++;
391 +        }
392 +      }
393 +      //    if ( numPos < ceil(RealType(n)/RealType(2.0)) )
394 +      if ( numPos < ceil_half_n) {
395 +        for (i=0; i<n; i++) {
396 +          v(i, j) *= -1.0;
397 +        }
398 +      }
399 +    }
400 +
401 +    if (n > 4) {
402 +      delete [] b;
403 +      delete [] z;
404 +    }
405 +    return 1;
406 +  }
407 +
408 +
409 +  typedef SquareMatrix<RealType, 6> Mat6x6d;
410   }
411   #endif //MATH_SQUAREMATRIX_HPP
412 +

Comparing trunk/src/math/SquareMatrix.hpp (property svn:keywords):
Revision 74 by tim, Wed Oct 13 23:53:40 2004 UTC vs.
Revision 1924 by gezelter, Mon Aug 5 21:46:11 2013 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines