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

Comparing:
trunk/src/math/RectMatrix.hpp (file contents), Revision 74 by tim, Wed Oct 13 23:53:40 2004 UTC vs.
branches/development/src/math/RectMatrix.hpp (file contents), Revision 1787 by gezelter, Wed Aug 29 18:13:11 2012 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, 24107 (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 <
26 <
42 >
43   /**
44   * @file RectMatrix.hpp
45   * @author Teng Lin
# Line 33 | Line 49
49  
50   #ifndef MATH_RECTMATRIX_HPP
51   #define MATH_RECTMATRIX_HPP
52 <
52 > #include <math.h>
53   #include <cmath>
54   #include "Vector.hpp"
55  
56 < namespace oopse {
41 <    const double epsilon = 0.000001;
56 > namespace OpenMD {
57  
58 <    template<typename T>
59 <    inline bool equal(T e1, T e2) {
60 <        return e1 == e2;
58 >  /**
59 >   * @class RectMatrix RectMatrix.hpp "math/RectMatrix.hpp"
60 >   * @brief rectangular matrix class
61 >   */
62 >  template<typename Real, unsigned int Row, unsigned int Col>
63 >  class RectMatrix {
64 >  public:
65 >    typedef Real ElemType;
66 >    typedef Real* ElemPoinerType;
67 >            
68 >    /** default constructor */
69 >    RectMatrix() {
70 >      for (unsigned int i = 0; i < Row; i++)
71 >        for (unsigned int j = 0; j < Col; j++)
72 >          this->data_[i][j] = 0.0;
73      }
74  
75 <    template<>
76 <    inline bool equal(float e1, float e2) {
77 <        return fabs(e1 - e2) < epsilon;
75 >    /** Constructs and initializes every element of this matrix to a scalar */
76 >    RectMatrix(Real s) {
77 >      for (unsigned int i = 0; i < Row; i++)
78 >        for (unsigned int j = 0; j < Col; j++)
79 >          this->data_[i][j] = s;
80      }
81  
82 <    template<>
83 <    inline bool equal(double e1, double e2) {
84 <        return fabs(e1 - e2) < epsilon;
82 >    RectMatrix(Real* array) {
83 >      for (unsigned int i = 0; i < Row; i++)
84 >        for (unsigned int j = 0; j < Col; j++)
85 >          this->data_[i][j] = array[i * Row + j];
86      }
87  
88 +    /** copy constructor */
89 +    RectMatrix(const RectMatrix<Real, Row, Col>& m) {
90 +      *this = m;
91 +    }
92 +            
93 +    /** destructor*/
94 +    ~RectMatrix() {}
95 +
96 +    /** copy assignment operator */
97 +    RectMatrix<Real, Row, Col>& operator =(const RectMatrix<Real, Row, Col>& m) {
98 +      if (this == &m)
99 +        return *this;
100 +                
101 +      for (unsigned int i = 0; i < Row; i++)
102 +        for (unsigned int j = 0; j < Col; j++)
103 +          this->data_[i][j] = m.data_[i][j];
104 +      return *this;
105 +    }
106 +            
107      /**
108 <     * @class RectMatrix RectMatrix.hpp "math/RectMatrix.hpp"
109 <     * @brief rectangular matrix class
108 >     * Return the reference of a single element of this matrix.
109 >     * @return the reference of a single element of this matrix
110 >     * @param i row index
111 >     * @param j Column index
112       */
113 <    template<typename Real, unsigned int Row, unsigned int Col>
114 <    class RectMatrix {
115 <        public:
113 >    Real& operator()(unsigned int i, unsigned int j) {
114 >      //assert( i < Row && j < Col);
115 >      return this->data_[i][j];
116 >    }
117  
118 <        /** default constructor */
119 <        RectMatrix() {
120 <            for (unsigned int i = 0; i < Row; i++)
121 <                for (unsigned int j = 0; j < Col; j++)
122 <                    data_[i][j] = 0.0;
123 <         }
118 >    /**
119 >     * Return the value of a single element of this matrix.
120 >     * @return the value of a single element of this matrix
121 >     * @param i row index
122 >     * @param j Column index
123 >     */        
124 >    Real operator()(unsigned int i, unsigned int j) const  {
125 >                
126 >      return this->data_[i][j];  
127 >    }
128  
129 <        /** Constructs and initializes every element of this matrix to a scalar */
130 <        RectMatrix(Real s) {
131 <            for (unsigned int i = 0; i < Row; i++)
132 <                for (unsigned int j = 0; j < Col; j++)
133 <                    data_[i][j] = s;
134 <        }
129 >    /**
130 >     * Copy the internal data to an array
131 >     * @param array the pointer of destination array
132 >     */
133 >    void getArray(Real* array) {
134 >      for (unsigned int i = 0; i < Row; i++) {
135 >        for (unsigned int j = 0; j < Col; j++) {
136 >          array[i * Row + j] = this->data_[i][j];
137 >        }
138 >      }
139 >    }
140  
80        /** copy constructor */
81        RectMatrix(const RectMatrix<Real, Row, Col>& m) {
82            *this = m;
83        }
84        
85        /** destructor*/
86        ~RectMatrix() {}
141  
142 <        /** copy assignment operator */
143 <        RectMatrix<Real, Row, Col>& operator =(const RectMatrix<Real, Row, Col>& m) {
144 <            if (this == &m)
145 <                return *this;
92 <            
93 <            for (unsigned int i = 0; i < Row; i++)
94 <                for (unsigned int j = 0; j < Col; j++)
95 <                    data_[i][j] = m.data_[i][j];
96 <            return *this;
97 <        }
98 <        
99 <        /**
100 <         * Return the reference of a single element of this matrix.
101 <         * @return the reference of a single element of this matrix
102 <         * @param i row index
103 <         * @param j colum index
104 <         */
105 <        double& operator()(unsigned int i, unsigned int j) {
106 <            //assert( i < Row && j < Col);
107 <            return data_[i][j];
108 <        }
142 >    /** Returns the pointer of internal array */
143 >    Real* getArrayPointer() {
144 >      return &this->data_[0][0];
145 >    }
146  
147 <        /**
148 <         * Return the value of a single element of this matrix.
149 <         * @return the value of a single element of this matrix
150 <         * @param i row index
151 <         * @param j colum index
152 <         */        
153 <        double operator()(unsigned int i, unsigned int j) const  {
117 <            
118 <            return data_[i][j];  
119 <        }
147 >    /**
148 >     * Returns a row of  this matrix as a vector.
149 >     * @return a row of  this matrix as a vector
150 >     * @param row the row index
151 >     */                
152 >    Vector<Real, Row> getRow(unsigned int row) {
153 >      Vector<Real, Row> v;
154  
155 <        /**
156 <         * Returns a row of  this matrix as a vector.
123 <         * @return a row of  this matrix as a vector
124 <         * @param row the row index
125 <         */                
126 <        Vector<Real, Row> getRow(unsigned int row) {
127 <            Vector<Real, Row> v;
155 >      for (unsigned int i = 0; i < Col; i++)
156 >        v[i] = this->data_[row][i];
157  
158 <            for (unsigned int i = 0; i < Row; i++)
159 <                v[i] = data_[row][i];
158 >      return v;
159 >    }
160  
161 <            return v;
162 <        }
161 >    /**
162 >     * Sets a row of  this matrix
163 >     * @param row the row index
164 >     * @param v the vector to be set
165 >     */                
166 >    void setRow(unsigned int row, const Vector<Real, Row>& v) {
167  
168 <        /**
169 <         * Sets a row of  this matrix
170 <         * @param row the row index
138 <         * @param v the vector to be set
139 <         */                
140 <         void setRow(unsigned int row, const Vector<Real, Row>& v) {
168 >      for (unsigned int i = 0; i < Col; i++)
169 >        this->data_[row][i] = v[i];
170 >    }
171  
172 <            for (unsigned int i = 0; i < Row; i++)
173 <                data_[row][i] = v[i];
174 <         }
175 <
176 <        /**
177 <         * Returns a column of  this matrix as a vector.
178 <         * @return a column of  this matrix as a vector
149 <         * @param col the column index
150 <         */                
151 <        Vector<Real, Col> getColum(unsigned int col) {
152 <            Vector<Real, Col> v;
172 >    /**
173 >     * Returns a column of  this matrix as a vector.
174 >     * @return a column of  this matrix as a vector
175 >     * @param col the column index
176 >     */                
177 >    Vector<Real, Col> getColumn(unsigned int col) {
178 >      Vector<Real, Col> v;
179  
180 <            for (unsigned int j = 0; j < Col; j++)
181 <                v[j] = data_[j][col];
180 >      for (unsigned int j = 0; j < Row; j++)
181 >        v[j] = this->data_[j][col];
182  
183 <            return v;
184 <        }
183 >      return v;
184 >    }
185  
186 <        /**
187 <         * Sets a column of  this matrix
188 <         * @param col the column index
189 <         * @param v the vector to be set
190 <         */                
191 <         void setColum(unsigned int col, const Vector<Real, Col>& v){
186 >    /**
187 >     * Sets a column of  this matrix
188 >     * @param col the column index
189 >     * @param v the vector to be set
190 >     */                
191 >    void setColumn(unsigned int col, const Vector<Real, Col>& v){
192  
193 <            for (unsigned int j = 0; j < Col; j++)
194 <                data_[j][col] = v[j];
195 <         }        
193 >      for (unsigned int j = 0; j < Row; j++)
194 >        this->data_[j][col] = v[j];
195 >    }        
196  
197 <        /**
198 <         * Tests if this matrix is identical to matrix m
199 <         * @return true if this matrix is equal to the matrix m, return false otherwise
200 <         * @m matrix to be compared
201 <         *
202 <         * @todo replace operator == by template function equal
203 <         */
178 <        bool operator ==(const RectMatrix<Real, Row, Col>& m) {
179 <            for (unsigned int i = 0; i < Row; i++)
180 <                for (unsigned int j = 0; j < Col; j++)
181 <                    if (!equal(data_[i][j], m.data_[i][j]))
182 <                        return false;
197 >    /**
198 >     * swap two rows of this matrix
199 >     * @param i the first row
200 >     * @param j the second row
201 >     */
202 >    void swapRow(unsigned int i, unsigned int j){
203 >      assert(i < Row && j < Row);
204  
205 <            return true;
206 <        }
205 >      for (unsigned int k = 0; k < Col; k++)
206 >        std::swap(this->data_[i][k], this->data_[j][k]);
207 >    }
208  
209 <        /**
210 <         * Tests if this matrix is not equal to matrix m
211 <         * @return true if this matrix is not equal to the matrix m, return false otherwise
212 <         * @m matrix to be compared
213 <         */
214 <        bool operator !=(const RectMatrix<Real, Row, Col>& m) {
215 <            return !(*this == m);
216 <        }
209 >    /**
210 >     * swap two Columns of this matrix
211 >     * @param i the first Column
212 >     * @param j the second Column
213 >     */
214 >    void swapColumn(unsigned int i, unsigned int j){
215 >      assert(i < Col && j < Col);
216 >                    
217 >      for (unsigned int k = 0; k < Row; k++)
218 >        std::swap(this->data_[k][i], this->data_[k][j]);
219 >    }
220  
221 <        /** Negates the value of this matrix in place. */          
222 <        inline void negate() {
223 <            for (unsigned int i = 0; i < Row; i++)
224 <                for (unsigned int j = 0; j < Col; j++)
225 <                    data_[i][j] = -data_[i][j];
226 <        }
227 <        
228 <        /**
229 <        * Sets the value of this matrix to the negation of matrix m.
230 <        * @param m the source matrix
231 <        */
232 <        inline void negate(const RectMatrix<Real, Row, Col>& m) {
208 <            for (unsigned int i = 0; i < Row; i++)
209 <                for (unsigned int j = 0; j < Col; j++)
210 <                    data_[i][j] = -m.data_[i][j];        
211 <        }
212 <        
213 <        /**
214 <        * Sets the value of this matrix to the sum of itself and m (*this += m).
215 <        * @param m the other matrix
216 <        */
217 <        inline void add( const RectMatrix<Real, Row, Col>& m ) {
218 <            for (unsigned int i = 0; i < Row; i++)
219 <                for (unsigned int j = 0; j < Col; j++)        
220 <                data_[i][j] += m.data_[i][j];
221 <        }
222 <        
223 <        /**
224 <        * Sets the value of this matrix to the sum of m1 and m2 (*this = m1 + m2).
225 <        * @param m1 the first matrix
226 <        * @param m2 the second matrix
227 <        */
228 <        inline void add( const RectMatrix<Real, Row, Col>& m1, const RectMatrix<Real, Row, Col>& m2 ) {
229 <            for (unsigned int i = 0; i < Row; i++)
230 <                for (unsigned int j = 0; j < Col; j++)        
231 <                data_[i][j] = m1.data_[i][j] + m2.data_[i][j];
232 <        }
233 <        
234 <        /**
235 <        * Sets the value of this matrix to the difference  of itself and m (*this -= m).
236 <        * @param m the other matrix
237 <        */
238 <        inline void sub( const RectMatrix<Real, Row, Col>& m ) {
239 <            for (unsigned int i = 0; i < Row; i++)
240 <                for (unsigned int j = 0; j < Col; j++)        
241 <                    data_[i][j] -= m.data_[i][j];
242 <        }
243 <        
244 <        /**
245 <        * Sets the value of this matrix to the difference of matrix m1 and m2 (*this = m1 - m2).
246 <        * @param m1 the first matrix
247 <        * @param m2 the second matrix
248 <        */
249 <        inline void sub( const RectMatrix<Real, Row, Col>& m1, const RectMatrix<Real, Row, Col>& m2){
250 <            for (unsigned int i = 0; i < Row; i++)
251 <                for (unsigned int j = 0; j < Col; j++)        
252 <                    data_[i][j] = m1.data_[i][j] - m2.data_[i][j];
253 <        }
221 >    /**
222 >     * Tests if this matrix is identical to matrix m
223 >     * @return true if this matrix is equal to the matrix m, return false otherwise
224 >     * @m matrix to be compared
225 >     *
226 >     * @todo replace operator == by template function equal
227 >     */
228 >    bool operator ==(const RectMatrix<Real, Row, Col>& m) {
229 >      for (unsigned int i = 0; i < Row; i++)
230 >        for (unsigned int j = 0; j < Col; j++)
231 >          if (!equal(this->data_[i][j], m.data_[i][j]))
232 >            return false;
233  
234 <        /**
235 <        * Sets the value of this matrix to the scalar multiplication of itself (*this *= s).
257 <        * @param s the scalar value
258 <        */
259 <        inline void mul( double s ) {
260 <            for (unsigned int i = 0; i < Row; i++)
261 <                for (unsigned int j = 0; j < Col; j++)  
262 <                    data_[i][j] *= s;
263 <        }
234 >      return true;
235 >    }
236  
237 <        /**
238 <        * Sets the value of this matrix to the scalar multiplication of matrix m  (*this = s * m).
239 <        * @param s the scalar value
240 <        * @param m the matrix
241 <        */
242 <        inline void mul( double s, const RectMatrix<Real, Row, Col>& m ) {
243 <            for (unsigned int i = 0; i < Row; i++)
244 <                for (unsigned int j = 0; j < Col; j++)  
273 <                    data_[i][j] = s * m.data_[i][j];
274 <        }
237 >    /**
238 >     * Tests if this matrix is not equal to matrix m
239 >     * @return true if this matrix is not equal to the matrix m, return false otherwise
240 >     * @m matrix to be compared
241 >     */
242 >    bool operator !=(const RectMatrix<Real, Row, Col>& m) {
243 >      return !(*this == m);
244 >    }
245  
246 <        /**
247 <        * Sets the value of this matrix to the scalar division of itself  (*this /= s ).
248 <        * @param s the scalar value
249 <        */            
250 <        inline void div( double s) {
251 <            for (unsigned int i = 0; i < Row; i++)
252 <                for (unsigned int j = 0; j < Col; j++)  
253 <                    data_[i][j] /= s;
254 <        }
255 <
256 <        /**
257 <        * Sets the value of this matrix to the scalar division of matrix m  (*this = m /s).
258 <        * @param s the scalar value
259 <        * @param m the matrix
260 <        */
261 <        inline void div( double s, const RectMatrix<Real, Row, Col>& m ) {
262 <            for (unsigned int i = 0; i < Row; i++)
263 <                for (unsigned int j = 0; j < Col; j++)  
264 <                    data_[i][j] = m.data_[i][j] / s;
265 <        }
266 <
267 <        /**
268 <         *  Multiples a scalar into every element of this matrix.
269 <         * @param s the scalar value
270 <         */
271 <        RectMatrix<Real, Row, Col>& operator *=(const double s) {
272 <            this->mul(s);
273 <            return *this;
274 <        }
275 <
276 <        /**
277 <         *  Divides every element of this matrix by a scalar.
278 <         * @param s the scalar value
279 <         */
280 <        RectMatrix<Real, Row, Col>& operator /=(const double s) {
281 <            this->div(s);
282 <            return *this;
283 <        }
284 <
285 <        /**
286 <         * Sets the value of this matrix to the sum of the other matrix and itself (*this += m).
287 <         * @param m the other matrix
288 <         */
289 <        RectMatrix<Real, Row, Col>& operator += (const RectMatrix<Real, Row, Col>& m) {
290 <            add(m);
291 <            return *this;
292 <         }
293 <
294 <        /**
295 <         * Sets the value of this matrix to the differerence of itself and the other matrix (*this -= m)
296 <         * @param m the other matrix
297 <         */
298 <        RectMatrix<Real, Row, Col>& operator -= (const RectMatrix<Real, Row, Col>& m){
299 <            sub(m);
300 <            return *this;
301 <        }
302 <
303 <        /** Return the transpose of this matrix */
334 <        RectMatrix<Real,  Col, Row> transpose(){
335 <            RectMatrix<Real,  Col, Row> result;
336 <            
337 <            for (unsigned int i = 0; i < Row; i++)
338 <                for (unsigned int j = 0; j < Col; j++)              
339 <                    result(j, i) = data_[i][j];
246 >    /** Negates the value of this matrix in place. */          
247 >    inline void negate() {
248 >      for (unsigned int i = 0; i < Row; i++)
249 >        for (unsigned int j = 0; j < Col; j++)
250 >          this->data_[i][j] = -this->data_[i][j];
251 >    }
252 >            
253 >    /**
254 >     * Sets the value of this matrix to the negation of matrix m.
255 >     * @param m the source matrix
256 >     */
257 >    inline void negate(const RectMatrix<Real, Row, Col>& m) {
258 >      for (unsigned int i = 0; i < Row; i++)
259 >        for (unsigned int j = 0; j < Col; j++)
260 >          this->data_[i][j] = -m.data_[i][j];        
261 >    }
262 >            
263 >    /**
264 >     * Sets the value of this matrix to the sum of itself and m (*this += m).
265 >     * @param m the other matrix
266 >     */
267 >    inline void add( const RectMatrix<Real, Row, Col>& m ) {
268 >      for (unsigned int i = 0; i < Row; i++)
269 >        for (unsigned int j = 0; j < Col; j++)        
270 >          this->data_[i][j] += m.data_[i][j];
271 >    }
272 >            
273 >    /**
274 >     * Sets the value of this matrix to the sum of m1 and m2 (*this = m1 + m2).
275 >     * @param m1 the first matrix
276 >     * @param m2 the second matrix
277 >     */
278 >    inline void add( const RectMatrix<Real, Row, Col>& m1, const RectMatrix<Real, Row, Col>& m2 ) {
279 >      for (unsigned int i = 0; i < Row; i++)
280 >        for (unsigned int j = 0; j < Col; j++)        
281 >          this->data_[i][j] = m1.data_[i][j] + m2.data_[i][j];
282 >    }
283 >            
284 >    /**
285 >     * Sets the value of this matrix to the difference  of itself and m (*this -= m).
286 >     * @param m the other matrix
287 >     */
288 >    inline void sub( const RectMatrix<Real, Row, Col>& m ) {
289 >      for (unsigned int i = 0; i < Row; i++)
290 >        for (unsigned int j = 0; j < Col; j++)        
291 >          this->data_[i][j] -= m.data_[i][j];
292 >    }
293 >            
294 >    /**
295 >     * Sets the value of this matrix to the difference of matrix m1 and m2 (*this = m1 - m2).
296 >     * @param m1 the first matrix
297 >     * @param m2 the second matrix
298 >     */
299 >    inline void sub( const RectMatrix<Real, Row, Col>& m1, const RectMatrix<Real, Row, Col>& m2){
300 >      for (unsigned int i = 0; i < Row; i++)
301 >        for (unsigned int j = 0; j < Col; j++)        
302 >          this->data_[i][j] = m1.data_[i][j] - m2.data_[i][j];
303 >    }
304  
305 <            return result;
306 <        }
307 <        
308 <        protected:
309 <            Real data_[Row][Col];
310 <    };
305 >    /**
306 >     * Sets the value of this matrix to the scalar multiplication of itself (*this *= s).
307 >     * @param s the scalar value
308 >     */
309 >    inline void mul( Real s ) {
310 >      for (unsigned int i = 0; i < Row; i++)
311 >        for (unsigned int j = 0; j < Col; j++)  
312 >          this->data_[i][j] *= s;
313 >    }
314  
315 <    /** Negate the value of every element of this matrix. */
316 <    template<typename Real, unsigned int Row, unsigned int Col>
317 <    inline RectMatrix<Real, Row, Col> operator -(const RectMatrix<Real, Row, Col>& m) {
318 <        RectMatrix<Real, Row, Col> result(m);
315 >    /**
316 >     * Sets the value of this matrix to the scalar multiplication of matrix m  (*this = s * m).
317 >     * @param s the scalar value
318 >     * @param m the matrix
319 >     */
320 >    inline void mul( Real s, const RectMatrix<Real, Row, Col>& m ) {
321 >      for (unsigned int i = 0; i < Row; i++)
322 >        for (unsigned int j = 0; j < Col; j++)  
323 >          this->data_[i][j] = s * m.data_[i][j];
324 >    }
325  
326 <        result.negate();
326 >    /**
327 >     * Sets the value of this matrix to the scalar division of itself  (*this /= s ).
328 >     * @param s the scalar value
329 >     */            
330 >    inline void div( Real s) {
331 >      for (unsigned int i = 0; i < Row; i++)
332 >        for (unsigned int j = 0; j < Col; j++)  
333 >          this->data_[i][j] /= s;
334 >    }
335  
336 <        return result;
336 >    /**
337 >     * Sets the value of this matrix to the scalar division of matrix m  (*this = m /s).
338 >     * @param s the scalar value
339 >     * @param m the matrix
340 >     */
341 >    inline void div( Real s, const RectMatrix<Real, Row, Col>& m ) {
342 >      for (unsigned int i = 0; i < Row; i++)
343 >        for (unsigned int j = 0; j < Col; j++)  
344 >          this->data_[i][j] = m.data_[i][j] / s;
345      }
346 <    
346 >
347      /**
348 <    * Return the sum of two matrixes  (m1 + m2).
349 <    * @return the sum of two matrixes
350 <    * @param m1 the first matrix
351 <    * @param m2 the second matrix
352 <    */
353 <    template<typename Real, unsigned int Row, unsigned int Col>
354 <    inline RectMatrix<Real, Row, Col> operator + (const RectMatrix<Real, Row, Col>& m1,const RectMatrix<Real, Row, Col>& m2) {
366 <        RectMatrix<Real, Row, Col> result;
348 >     *  Multiples a scalar into every element of this matrix.
349 >     * @param s the scalar value
350 >     */
351 >    RectMatrix<Real, Row, Col>& operator *=(const Real s) {
352 >      this->mul(s);
353 >      return *this;
354 >    }
355  
356 <        result.add(m1, m2);
356 >    /**
357 >     *  Divides every element of this matrix by a scalar.
358 >     * @param s the scalar value
359 >     */
360 >    RectMatrix<Real, Row, Col>& operator /=(const Real s) {
361 >      this->div(s);
362 >      return *this;
363 >    }
364  
365 <        return result;
365 >    /**
366 >     * Sets the value of this matrix to the sum of the other matrix and itself (*this += m).
367 >     * @param m the other matrix
368 >     */
369 >    RectMatrix<Real, Row, Col>& operator += (const RectMatrix<Real, Row, Col>& m) {
370 >      add(m);
371 >      return *this;
372      }
373 <    
373 >
374      /**
375 <    * Return the difference of two matrixes  (m1 - m2).
376 <    * @return the sum of two matrixes
377 <    * @param m1 the first matrix
378 <    * @param m2 the second matrix
379 <    */
380 <    template<typename Real, unsigned int Row, unsigned int Col>
381 <    inline RectMatrix<Real, Row, Col> operator - (const RectMatrix<Real, Row, Col>& m1, const RectMatrix<Real, Row, Col>& m2) {
381 <        RectMatrix<Real, Row, Col> result;
375 >     * Sets the value of this matrix to the differerence of itself and the other matrix (*this -= m)
376 >     * @param m the other matrix
377 >     */
378 >    RectMatrix<Real, Row, Col>& operator -= (const RectMatrix<Real, Row, Col>& m){
379 >      sub(m);
380 >      return *this;
381 >    }
382  
383 <        result.sub(m1, m2);
383 >    /** Return the transpose of this matrix */
384 >    RectMatrix<Real,  Col, Row> transpose() const{
385 >      RectMatrix<Real,  Col, Row> result;
386 >                
387 >      for (unsigned int i = 0; i < Row; i++)
388 >        for (unsigned int j = 0; j < Col; j++)              
389 >          result(j, i) = this->data_[i][j];
390  
391 <        return result;
391 >      return result;
392      }
393  
394 <    /**
395 <    * Return the multiplication of scalra and  matrix  (m * s).
396 <    * @return the multiplication of a scalra and  a matrix
397 <    * @param m the matrix
392 <    * @param s the scalar
393 <    */
394 <    template<typename Real, unsigned int Row, unsigned int Col>
395 <    inline RectMatrix<Real, Row, Col> operator *(const RectMatrix<Real, Row, Col>& m, Real s) {
396 <        RectMatrix<Real, Row, Col> result;
394 >    template<class MatrixType>
395 >    void setSubMatrix(unsigned int beginRow, unsigned int beginCol, const MatrixType& m) {
396 >        assert(beginRow + m.getNRow() -1 <= getNRow());
397 >        assert(beginCol + m.getNCol() -1 <= getNCol());
398  
399 <        result.mul(s, m);
399 >        for (unsigned int i = 0; i < m.getNRow(); ++i)
400 >            for (unsigned int j = 0; j < m.getNCol(); ++j)
401 >                this->data_[beginRow+i][beginCol+j] = m(i, j);
402 >    }
403  
404 <        return result;
404 >    template<class MatrixType>
405 >    void getSubMatrix(unsigned int beginRow, unsigned int beginCol, MatrixType& m) {
406 >        assert(beginRow + m.getNRow() -1 <= getNRow());
407 >        assert(beginCol + m.getNCol() - 1 <= getNCol());
408 >
409 >        for (unsigned int i = 0; i < m.getNRow(); ++i)
410 >            for (unsigned int j = 0; j < m.getNCol(); ++j)
411 >                m(i, j) = this->data_[beginRow+i][beginCol+j];
412      }
413 +    
414 +    unsigned int getNRow() const {return Row;}
415 +    unsigned int getNCol() const {return Col;}        
416  
417 <    /**
418 <    * Return the multiplication of a scalra and  a matrix  (s * m).
419 <    * @return the multiplication of a scalra and  a matrix
406 <    * @param s the scalar
407 <    * @param m the matrix
408 <    */
409 <    template<typename Real, unsigned int Row, unsigned int Col>
410 <    inline RectMatrix<Real, Row, Col> operator *(Real s, const RectMatrix<Real, Row, Col>& m) {
411 <        RectMatrix<Real, Row, Col> result;
417 >  protected:
418 >    Real data_[Row][Col];
419 >  };
420  
421 <        result.mul(s, m);
421 >  /** Negate the value of every element of this matrix. */
422 >  template<typename Real, unsigned int Row, unsigned int Col>
423 >  inline RectMatrix<Real, Row, Col> operator -(const RectMatrix<Real, Row, Col>& m) {
424 >    RectMatrix<Real, Row, Col> result(m);
425  
426 <        return result;
427 <    }
426 >    result.negate();
427 >
428 >    return result;
429 >  }
430      
431 <    /**
432 <    * Return the multiplication of two matrixes  (m1 * m2).
433 <    * @return the multiplication of two matrixes
434 <    * @param m1 the first matrix
435 <    * @param m2 the second matrix
436 <    */
437 <    template<typename Real, unsigned int Row, unsigned int Col, unsigned int SameDim>
438 <    inline RectMatrix<Real, Row, Col> operator *(const RectMatrix<Real, Row, SameDim>& m1, const RectMatrix<Real, SameDim, Col>& m2) {
439 <        RectMatrix<Real, Row, Col> result;
431 >  /**
432 >   * Return the sum of two matrixes  (m1 + m2).
433 >   * @return the sum of two matrixes
434 >   * @param m1 the first matrix
435 >   * @param m2 the second matrix
436 >   */
437 >  template<typename Real, unsigned int Row, unsigned int Col>
438 >  inline RectMatrix<Real, Row, Col> operator + (const RectMatrix<Real, Row, Col>& m1,const RectMatrix<Real, Row, Col>& m2) {
439 >    RectMatrix<Real, Row, Col> result;
440  
441 <            for (unsigned int i = 0; i < Row; i++)
429 <                for (unsigned int j = 0; j < Col; j++)
430 <                    for (unsigned int k = 0; k < SameDim; k++)
431 <                        result(i, j)  = m1(i, k) * m2(k, j);                
441 >    result.add(m1, m2);
442  
443 <        return result;
444 <    }
443 >    return result;
444 >  }
445      
446 <    /**
447 <    * Return the multiplication of  a matrix and a vector  (m * v).
448 <    * @return the multiplication of a matrix and a vector
449 <    * @param m the matrix
450 <    * @param v the vector
451 <    */
452 <    template<typename Real, unsigned int Row, unsigned int Col>
453 <    inline Vector<Real, Row> operator *(const RectMatrix<Real, Row, Col>& m, const Vector<Real, Col>& v) {
454 <        Vector<Real, Row> result;
446 >  /**
447 >   * Return the difference of two matrixes  (m1 - m2).
448 >   * @return the sum of two matrixes
449 >   * @param m1 the first matrix
450 >   * @param m2 the second matrix
451 >   */
452 >  template<typename Real, unsigned int Row, unsigned int Col>
453 >  inline RectMatrix<Real, Row, Col> operator - (const RectMatrix<Real, Row, Col>& m1, const RectMatrix<Real, Row, Col>& m2) {
454 >    RectMatrix<Real, Row, Col> result;
455  
456 <        for (unsigned int i = 0; i < Row ; i++)
457 <            for (unsigned int j = 0; j < Col ; j++)            
458 <                result[i] += m(i, j) * v[j];
456 >    result.sub(m1, m2);
457 >
458 >    return result;
459 >  }
460 >
461 >  /**
462 >   * Return the multiplication of scalra and  matrix  (m * s).
463 >   * @return the multiplication of a scalra and  a matrix
464 >   * @param m the matrix
465 >   * @param s the scalar
466 >   */
467 >  template<typename Real, unsigned int Row, unsigned int Col>
468 >  inline RectMatrix<Real, Row, Col> operator *(const RectMatrix<Real, Row, Col>& m, Real s) {
469 >    RectMatrix<Real, Row, Col> result;
470 >
471 >    result.mul(s, m);
472 >
473 >    return result;
474 >  }
475 >
476 >  /**
477 >   * Return the multiplication of a scalra and  a matrix  (s * m).
478 >   * @return the multiplication of a scalra and  a matrix
479 >   * @param s the scalar
480 >   * @param m the matrix
481 >   */
482 >  template<typename Real, unsigned int Row, unsigned int Col>
483 >  inline RectMatrix<Real, Row, Col> operator *(Real s, const RectMatrix<Real, Row, Col>& m) {
484 >    RectMatrix<Real, Row, Col> result;
485 >
486 >    result.mul(s, m);
487 >
488 >    return result;
489 >  }
490 >    
491 >  /**
492 >   * Return the multiplication of two matrixes  (m1 * m2).
493 >   * @return the multiplication of two matrixes
494 >   * @param m1 the first matrix
495 >   * @param m2 the second matrix
496 >   */
497 >  template<typename Real, unsigned int Row, unsigned int Col, unsigned int SameDim>
498 >  inline RectMatrix<Real, Row, Col> operator *(const RectMatrix<Real, Row, SameDim>& m1, const RectMatrix<Real, SameDim, Col>& m2) {
499 >    RectMatrix<Real, Row, Col> result;
500 >
501 >    for (unsigned int i = 0; i < Row; i++)
502 >      for (unsigned int j = 0; j < Col; j++)
503 >        for (unsigned int k = 0; k < SameDim; k++)
504 >          result(i, j)  += m1(i, k) * m2(k, j);                
505 >
506 >    return result;
507 >  }
508 >    
509 >  /**
510 >   * Returns the multiplication of  a matrix and a vector  (m * v).
511 >   * @return the multiplication of a matrix and a vector
512 >   * @param m the matrix
513 >   * @param v the vector
514 >   */
515 >  template<typename Real, unsigned int Row, unsigned int Col>
516 >  inline Vector<Real, Row> operator *(const RectMatrix<Real, Row, Col>& m, const Vector<Real, Col>& v) {
517 >    Vector<Real, Row> result;
518 >
519 >    for (unsigned int i = 0; i < Row ; i++)
520 >      for (unsigned int j = 0; j < Col ; j++)            
521 >        result[i] += m(i, j) * v[j];
522              
523 <        return result;                                                                
524 <    }
523 >    return result;                                                                
524 >  }
525  
526 <    /**
527 <    * Return the scalar division of matrix   (m / s).
528 <    * @return the scalar division of matrix  
529 <    * @param m the matrix
530 <    * @param s the scalar
531 <    */
532 <    template<typename Real, unsigned int Row, unsigned int Col>
533 <    inline RectMatrix<Real, Row, Col> operator /(const RectMatrix<Real, Row, Col>& m, Real s) {
534 <        RectMatrix<Real, Row, Col> result;
526 >  /**
527 >   * Returns the multiplication of a vector transpose and a matrix  (v^T * m).
528 >   * @return the multiplication of a vector transpose and a matrix
529 >   * @param v the vector
530 >   * @param m the matrix
531 >   */
532 >  template<typename Real, unsigned int Row, unsigned int Col>
533 >  inline Vector<Real, Col> operator *(const Vector<Real, Row>& v, const RectMatrix<Real, Row, Col>& m) {
534 >    Vector<Real, Row> result;
535 >    
536 >    for (unsigned int i = 0; i < Col ; i++)
537 >      for (unsigned int j = 0; j < Row ; j++)            
538 >        result[i] += v[j] * m(j, i);
539 >            
540 >    return result;                                                                
541 >  }
542  
543 <        result.div(s, m);
543 >  /**
544 >   * Return the scalar division of matrix   (m / s).
545 >   * @return the scalar division of matrix  
546 >   * @param m the matrix
547 >   * @param s the scalar
548 >   */
549 >  template<typename Real, unsigned int Row, unsigned int Col>
550 >  inline RectMatrix<Real, Row, Col> operator /(const RectMatrix<Real, Row, Col>& m, Real s) {
551 >    RectMatrix<Real, Row, Col> result;
552  
553 <        return result;
554 <    }    
553 >    result.div(s, m);
554 >
555 >    return result;
556 >  }    
557 >
558 >  
559 >  /**
560 >   * Returns the vector (cross) product of two matrices.  This
561 >   * operation is defined in:
562 >   *
563 >   * W. Smith, "Point Multipoles in the Ewald Summation (Revisited),"
564 >   * CCP5 Newsletter No 46., pp. 18-30.
565 >   *
566 >   * Equation 21 defines:
567 >   * V_alpha = \sum_\beta [ A_{\alpha+1,\beta} * B_{\alpha+2,\beta}
568 >                           -A_{\alpha+2,\beta} * B_{\alpha+2,\beta} ]
569 >   * where \alpha+1 and \alpha+2 are regarded as cyclic permuations of the
570 >   * matrix indices (i.e. for a 3x3 matrix, when \alpha = 2, \alpha + 1 = 3,
571 >   * and \alpha + 2 = 1).
572 >   *
573 >   * @param t1 first matrix
574 >   * @param t2 second matrix
575 >   * @return the cross product (vector product) of t1 and t2
576 >   */
577 >  template<typename Real, unsigned int Row, unsigned int Col>
578 >  inline Vector<Real, Row> cross( const RectMatrix<Real, Row, Col>& t1, const RectMatrix<Real, Row, Col>& t2 ) {
579 >    Vector<Real, Row> result;
580 >    unsigned int i1;
581 >    unsigned int i2;
582 >    
583 >    for (unsigned int i = 0; i < Row; i++) {
584 >      i1 = (i+1)%Row;
585 >      i2 = (i+2)%Row;
586 >      
587 >      for (unsigned int j =0; j < Col; j++) {        
588 >        result[i] = t1(i1,j) * t2(i2,j) - t1(i2,j) * t2(i1,j);
589 >      }
590 >    }
591 >    
592 >    return result;
593 >  }
594 >  
595 >  
596 >  /**
597 >   * Write to an output stream
598 >   */
599 >  template<typename Real,  unsigned int Row, unsigned int Col>
600 >  std::ostream &operator<< ( std::ostream& o, const RectMatrix<Real, Row, Col>& m) {
601 >    for (unsigned int i = 0; i < Row ; i++) {
602 >      o << "(";
603 >      for (unsigned int j = 0; j < Col ; j++) {
604 >        o << m(i, j);
605 >        if (j != Col -1)
606 >          o << "\t";
607 >      }
608 >      o << ")" << std::endl;
609 >    }
610 >    return o;        
611 >  }    
612   }
613   #endif //MATH_RECTMATRIX_HPP

Comparing:
trunk/src/math/RectMatrix.hpp (property svn:keywords), Revision 74 by tim, Wed Oct 13 23:53:40 2004 UTC vs.
branches/development/src/math/RectMatrix.hpp (property svn:keywords), Revision 1787 by gezelter, Wed Aug 29 18:13:11 2012 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines