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 101 by tim, Mon Oct 18 23:13:23 2004 UTC vs.
branches/development/src/math/RectMatrix.hpp (file contents), Revision 1808 by gezelter, Mon Oct 22 20:42:10 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 {
56 > namespace OpenMD {
57  
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:
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 <        /** default constructor */
76 <        RectMatrix() {
77 <            for (unsigned int i = 0; i < Row; i++)
78 <                for (unsigned int j = 0; j < Col; j++)
79 <                    data_[i][j] = 0.0;
80 <         }
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 <        /** Constructs and initializes every element of this matrix to a scalar */
83 <        RectMatrix(Real s) {
84 <            for (unsigned int i = 0; i < Row; i++)
85 <                for (unsigned int j = 0; j < Col; j++)
86 <                    data_[i][j] = s;
62 <        }
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 <        }
68 <        
69 <        /** destructor*/
70 <        ~RectMatrix() {}
71 <
72 <        /** copy assignment operator */
73 <        RectMatrix<Real, Row, Col>& operator =(const RectMatrix<Real, Row, Col>& m) {
74 <            if (this == &m)
75 <                return *this;
88 >    /** copy constructor */
89 >    RectMatrix(const RectMatrix<Real, Row, Col>& m) {
90 >      *this = m;
91 >    }
92              
93 <            for (unsigned int i = 0; i < Row; i++)
94 <                for (unsigned int j = 0; j < Col; j++)
79 <                    data_[i][j] = m.data_[i][j];
80 <            return *this;
81 <        }
82 <        
83 <        /**
84 <         * Return the reference of a single element of this matrix.
85 <         * @return the reference of a single element of this matrix
86 <         * @param i row index
87 <         * @param j colum index
88 <         */
89 <        double& operator()(unsigned int i, unsigned int j) {
90 <            //assert( i < Row && j < Col);
91 <            return data_[i][j];
92 <        }
93 >    /** destructor*/
94 >    ~RectMatrix() {}
95  
96 <        /**
97 <         * Return the value of a single element of this matrix.
98 <         * @return the value of a single element of this matrix
99 <         * @param i row index
100 <         * @param j colum index
101 <         */        
102 <        double operator()(unsigned int i, unsigned int j) const  {
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 <            return data_[i][j];  
108 <        }
107 >    /**
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 >    Real& operator()(unsigned int i, unsigned int j) {
114 >      //assert( i < Row && j < Col);
115 >      return this->data_[i][j];
116 >    }
117  
118 <        /**
119 <         * Returns a row of  this matrix as a vector.
120 <         * @return a row of  this matrix as a vector
121 <         * @param row the row index
122 <         */                
123 <        Vector<Real, Row> getRow(unsigned int row) {
124 <            Vector<Real, Row> v;
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 <            for (unsigned int i = 0; i < Row; i++)
130 <                v[i] = data_[row][i];
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  
116            return v;
117        }
141  
142 <        /**
143 <         * Sets a row of  this matrix
144 <         * @param row the row index
145 <         * @param v the vector to be set
123 <         */                
124 <         void setRow(unsigned int row, const Vector<Real, Row>& v) {
142 >    /** Returns the pointer of internal array */
143 >    Real* getArrayPointer() {
144 >      return &this->data_[0][0];
145 >    }
146  
147 <            for (unsigned int i = 0; i < Row; i++)
148 <                data_[row][i] = v[i];
149 <         }
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 column of  this matrix as a vector.
132 <         * @return a column of  this matrix as a vector
133 <         * @param col the column index
134 <         */                
135 <        Vector<Real, Col> getColum(unsigned int col) {
136 <            Vector<Real, Col> v;
155 >      for (unsigned int i = 0; i < Col; i++)
156 >        v[i] = this->data_[row][i];
157  
158 <            for (unsigned int j = 0; j < Col; j++)
159 <                v[j] = data_[j][col];
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 column of  this matrix
170 <         * @param col the column index
147 <         * @param v the vector to be set
148 <         */                
149 <         void setColum(unsigned int col, const Vector<Real, Col>& v){
168 >      for (unsigned int i = 0; i < Col; i++)
169 >        this->data_[row][i] = v[i];
170 >    }
171  
172 <            for (unsigned int j = 0; j < Col; j++)
173 <                data_[j][col] = v[j];
174 <         }        
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 <        /**
181 <         * swap two rows of this matrix
157 <         * @param i the first row
158 <         * @param j the second row
159 <         */
160 <        void swapRow(unsigned int i, unsigned int j){
161 <                assert(i < Row && j < Row);
180 >      for (unsigned int j = 0; j < Row; j++)
181 >        v[j] = this->data_[j][col];
182  
183 <                for (unsigned int k = 0; k < Col; k++)
184 <                    std::swap(data_[i][k], data_[j][k]);
165 <        }
183 >      return v;
184 >    }
185  
186 <       /**
187 <         * swap two colums of this matrix
188 <         * @param i the first colum
189 <         * @param j the second colum
190 <         */
191 <        void swapColum(unsigned int i, unsigned int j){
173 <                assert(i < Col && j < Col);
174 <                
175 <                for (unsigned int k = 0; k < Row; k++)
176 <                    std::swap(data_[k][i], data_[k][j]);
177 <        }
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 <        /**
194 <         * Tests if this matrix is identical to matrix m
195 <         * @return true if this matrix is equal to the matrix m, return false otherwise
182 <         * @m matrix to be compared
183 <         *
184 <         * @todo replace operator == by template function equal
185 <         */
186 <        bool operator ==(const RectMatrix<Real, Row, Col>& m) {
187 <            for (unsigned int i = 0; i < Row; i++)
188 <                for (unsigned int j = 0; j < Col; j++)
189 <                    if (!equal(data_[i][j], m.data_[i][j]))
190 <                        return false;
193 >      for (unsigned int j = 0; j < Row; j++)
194 >        this->data_[j][col] = v[j];
195 >    }        
196  
197 <            return true;
198 <        }
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 <        /**
206 <         * Tests if this matrix is not equal to matrix m
207 <         * @return true if this matrix is not equal to the matrix m, return false otherwise
198 <         * @m matrix to be compared
199 <         */
200 <        bool operator !=(const RectMatrix<Real, Row, Col>& m) {
201 <            return !(*this == m);
202 <        }
205 >      for (unsigned int k = 0; k < Col; k++)
206 >        std::swap(this->data_[i][k], this->data_[j][k]);
207 >    }
208  
209 <        /** Negates the value of this matrix in place. */          
210 <        inline void negate() {
211 <            for (unsigned int i = 0; i < Row; i++)
212 <                for (unsigned int j = 0; j < Col; j++)
213 <                    data_[i][j] = -data_[i][j];
214 <        }
215 <        
216 <        /**
217 <        * Sets the value of this matrix to the negation of matrix m.
218 <        * @param m the source matrix
219 <        */
215 <        inline void negate(const RectMatrix<Real, Row, Col>& m) {
216 <            for (unsigned int i = 0; i < Row; i++)
217 <                for (unsigned int j = 0; j < Col; j++)
218 <                    data_[i][j] = -m.data_[i][j];        
219 <        }
220 <        
221 <        /**
222 <        * Sets the value of this matrix to the sum of itself and m (*this += m).
223 <        * @param m the other matrix
224 <        */
225 <        inline void add( const RectMatrix<Real, Row, Col>& m ) {
226 <            for (unsigned int i = 0; i < Row; i++)
227 <                for (unsigned int j = 0; j < Col; j++)        
228 <                data_[i][j] += m.data_[i][j];
229 <        }
230 <        
231 <        /**
232 <        * Sets the value of this matrix to the sum of m1 and m2 (*this = m1 + m2).
233 <        * @param m1 the first matrix
234 <        * @param m2 the second matrix
235 <        */
236 <        inline void add( const RectMatrix<Real, Row, Col>& m1, const RectMatrix<Real, Row, Col>& m2 ) {
237 <            for (unsigned int i = 0; i < Row; i++)
238 <                for (unsigned int j = 0; j < Col; j++)        
239 <                data_[i][j] = m1.data_[i][j] + m2.data_[i][j];
240 <        }
241 <        
242 <        /**
243 <        * Sets the value of this matrix to the difference  of itself and m (*this -= m).
244 <        * @param m the other matrix
245 <        */
246 <        inline void sub( const RectMatrix<Real, Row, Col>& m ) {
247 <            for (unsigned int i = 0; i < Row; i++)
248 <                for (unsigned int j = 0; j < Col; j++)        
249 <                    data_[i][j] -= m.data_[i][j];
250 <        }
251 <        
252 <        /**
253 <        * Sets the value of this matrix to the difference of matrix m1 and m2 (*this = m1 - m2).
254 <        * @param m1 the first matrix
255 <        * @param m2 the second matrix
256 <        */
257 <        inline void sub( const RectMatrix<Real, Row, Col>& m1, const RectMatrix<Real, Row, Col>& m2){
258 <            for (unsigned int i = 0; i < Row; i++)
259 <                for (unsigned int j = 0; j < Col; j++)        
260 <                    data_[i][j] = m1.data_[i][j] - m2.data_[i][j];
261 <        }
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 <        /**
222 <        * Sets the value of this matrix to the scalar multiplication of itself (*this *= s).
223 <        * @param s the scalar value
224 <        */
225 <        inline void mul( double s ) {
226 <            for (unsigned int i = 0; i < Row; i++)
227 <                for (unsigned int j = 0; j < Col; j++)  
228 <                    data_[i][j] *= s;
229 <        }
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 >     * @param 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 matrix m  (*this = s * m).
275 <        * @param s the scalar value
276 <        * @param m the matrix
277 <        */
278 <        inline void mul( double s, const RectMatrix<Real, Row, Col>& m ) {
279 <            for (unsigned int i = 0; i < Row; i++)
280 <                for (unsigned int j = 0; j < Col; j++)  
281 <                    data_[i][j] = s * m.data_[i][j];
282 <        }
234 >      return true;
235 >    }
236  
237 <        /**
238 <        * Sets the value of this matrix to the scalar division of itself  (*this /= s ).
239 <        * @param s the scalar value
240 <        */            
241 <        inline void div( double s) {
242 <            for (unsigned int i = 0; i < Row; i++)
243 <                for (unsigned int j = 0; j < Col; j++)  
244 <                    data_[i][j] /= s;
292 <        }
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 >     * @param 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 matrix m  (*this = m /s).
248 <        * @param s the scalar value
249 <        * @param m the matrix
250 <        */
299 <        inline void div( double s, const RectMatrix<Real, Row, Col>& m ) {
300 <            for (unsigned int i = 0; i < Row; i++)
301 <                for (unsigned int j = 0; j < Col; j++)  
302 <                    data_[i][j] = m.data_[i][j] / s;
303 <        }
304 <
305 <        /**
306 <         *  Multiples a scalar into every element of this matrix.
307 <         * @param s the scalar value
308 <         */
309 <        RectMatrix<Real, Row, Col>& operator *=(const double s) {
310 <            this->mul(s);
311 <            return *this;
312 <        }
313 <
314 <        /**
315 <         *  Divides every element of this matrix by a scalar.
316 <         * @param s the scalar value
317 <         */
318 <        RectMatrix<Real, Row, Col>& operator /=(const double s) {
319 <            this->div(s);
320 <            return *this;
321 <        }
322 <
323 <        /**
324 <         * Sets the value of this matrix to the sum of the other matrix and itself (*this += m).
325 <         * @param m the other matrix
326 <         */
327 <        RectMatrix<Real, Row, Col>& operator += (const RectMatrix<Real, Row, Col>& m) {
328 <            add(m);
329 <            return *this;
330 <         }
331 <
332 <        /**
333 <         * Sets the value of this matrix to the differerence of itself and the other matrix (*this -= m)
334 <         * @param m the other matrix
335 <         */
336 <        RectMatrix<Real, Row, Col>& operator -= (const RectMatrix<Real, Row, Col>& m){
337 <            sub(m);
338 <            return *this;
339 <        }
340 <
341 <        /** Return the transpose of this matrix */
342 <        RectMatrix<Real,  Col, Row> transpose(){
343 <            RectMatrix<Real,  Col, Row> result;
344 <            
345 <            for (unsigned int i = 0; i < Row; i++)
346 <                for (unsigned int j = 0; j < Col; j++)              
347 <                    result(j, i) = data_[i][j];
348 <
349 <            return result;
350 <        }
351 <        
352 <        protected:
353 <            Real data_[Row][Col];
354 <    };
355 <
356 <    /** Negate the value of every element of this matrix. */
357 <    template<typename Real, unsigned int Row, unsigned int Col>
358 <    inline RectMatrix<Real, Row, Col> operator -(const RectMatrix<Real, Row, Col>& m) {
359 <        RectMatrix<Real, Row, Col> result(m);
360 <
361 <        result.negate();
362 <
363 <        return result;
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 <    
252 >            
253      /**
254 <    * Return the sum of two matrixes  (m1 + m2).
255 <    * @return the sum of two matrixes
256 <    * @param m1 the first matrix
257 <    * @param m2 the second matrix
258 <    */
259 <    template<typename Real, unsigned int Row, unsigned int Col>
260 <    inline RectMatrix<Real, Row, Col> operator + (const RectMatrix<Real, Row, Col>& m1,const RectMatrix<Real, Row, Col>& m2) {
261 <        RectMatrix<Real, Row, Col> result;
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 <        result.add(m1, m2);
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 <        return result;
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 <    
325 >
326      /**
327 <    * Return the difference of two matrixes  (m1 - m2).
328 <    * @return the sum of two matrixes
329 <    * @param m1 the first matrix
330 <    * @param m2 the second matrix
331 <    */
332 <    template<typename Real, unsigned int Row, unsigned int Col>
333 <    inline RectMatrix<Real, Row, Col> operator - (const RectMatrix<Real, Row, Col>& m1, const RectMatrix<Real, Row, Col>& m2) {
334 <        RectMatrix<Real, Row, Col> result;
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 <        result.sub(m1, m2);
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  
347 <        return result;
347 >    /**
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      /**
357 <    * Return the multiplication of scalra and  matrix  (m * s).
358 <    * @return the multiplication of a scalra and  a matrix
359 <    * @param m the matrix
360 <    * @param s the scalar
361 <    */
362 <    template<typename Real, unsigned int Row, unsigned int Col>
363 <    inline RectMatrix<Real, Row, Col> operator *(const RectMatrix<Real, Row, Col>& m, Real s) {
404 <        RectMatrix<Real, Row, Col> result;
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 <        result.mul(s, m);
366 <
367 <        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  
374      /**
375 <    * Return the multiplication of a scalra and  a matrix  (s * m).
376 <    * @return the multiplication of a scalra and  a matrix
377 <    * @param s the scalar
378 <    * @param m the matrix
379 <    */
380 <    template<typename Real, unsigned int Row, unsigned int Col>
381 <    inline RectMatrix<Real, Row, Col> operator *(Real s, const RectMatrix<Real, Row, Col>& m) {
419 <        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.mul(s, m);
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      }
425    
426    /**
427    * Return the multiplication of two matrixes  (m1 * m2).
428    * @return the multiplication of two matrixes
429    * @param m1 the first matrix
430    * @param m2 the second matrix
431    */
432    template<typename Real, unsigned int Row, unsigned int Col, unsigned int SameDim>
433    inline RectMatrix<Real, Row, Col> operator *(const RectMatrix<Real, Row, SameDim>& m1, const RectMatrix<Real, SameDim, Col>& m2) {
434        RectMatrix<Real, Row, Col> result;
393  
394 <            for (unsigned int i = 0; i < Row; i++)
395 <                for (unsigned int j = 0; j < Col; j++)
396 <                    for (unsigned int k = 0; k < SameDim; k++)
397 <                        result(i, j)  += m1(i, k) * m2(k, j);                
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 <        return result;
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      }
443    
444    /**
445    * Return the multiplication of  a matrix and a vector  (m * v).
446    * @return the multiplication of a matrix and a vector
447    * @param m the matrix
448    * @param v the vector
449    */
450    template<typename Real, unsigned int Row, unsigned int Col>
451    inline Vector<Real, Row> operator *(const RectMatrix<Real, Row, Col>& m, const Vector<Real, Col>& v) {
452        Vector<Real, Row> result;
403  
404 <        for (unsigned int i = 0; i < Row ; i++)
405 <            for (unsigned int j = 0; j < Col ; j++)            
406 <                result[i] += m(i, j) * v[j];
407 <            
408 <        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 scalar division of matrix   (m / s).
419 <    * @return the scalar division of 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;
417 >  protected:
418 >    Real data_[Row][Col];
419 >  };
420  
421 <        result.div(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;
474 <    }    
426 >    result.negate();
427  
428 <    /**
429 <     * Write to an output stream
430 <     */
431 <    template<typename Real,  unsigned int Row, unsigned int Col>
432 <    std::ostream &operator<< ( std::ostream& o, const RectMatrix<Real, Row, Col>& m) {
433 <        for (unsigned int i = 0; i < Row ; i++) {
434 <            o << "("
435 <            for (unsigned int j = 0; j < Col ; j++) {
436 <                o << m(i, j) << "\t"
437 <            }
438 <            o << ")" << std::endl;
439 <        }
440 <        return o;        
441 <    }    
428 >    return result;
429 >  }
430 >    
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 >    result.add(m1, m2);
442 >
443 >    return result;
444 >  }
445 >    
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 >    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 >  }
525 >
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 >  /**
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 >    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 >   * \f[
568 >   * V_alpha = \sum_\beta \left[ A_{\alpha+1,\beta} * B_{\alpha+2,\beta}
569 >                           -A_{\alpha+2,\beta} * B_{\alpha+2,\beta} \right]
570 >   * \f]
571 >   * where \f[\alpha+1\f] and \f[\alpha+2\f] are regarded as cyclic permuations of the
572 >   * matrix indices (i.e. for a 3x3 matrix, when \f[\alpha = 2\f], \f[\alpha + 1 = 3 \f],
573 >   * and \f[\alpha + 2 = 1 \f] ).
574 >   *
575 >   * @param t1 first matrix
576 >   * @param t2 second matrix
577 >   * @return the cross product (vector product) of t1 and t2
578 >   */
579 >  template<typename Real, unsigned int Row, unsigned int Col>
580 >  inline Vector<Real, Row> cross( const RectMatrix<Real, Row, Col>& t1, const RectMatrix<Real, Row, Col>& t2 ) {
581 >    Vector<Real, Row> result;
582 >    unsigned int i1;
583 >    unsigned int i2;
584 >    
585 >    for (unsigned int i = 0; i < Row; i++) {
586 >      i1 = (i+1)%Row;
587 >      i2 = (i+2)%Row;
588 >      
589 >      for (unsigned int j =0; j < Col; j++) {        
590 >        result[i] = t1(i1,j) * t2(i2,j) - t1(i2,j) * t2(i1,j);
591 >      }
592 >    }
593 >    
594 >    return result;
595 >  }
596 >  
597 >  
598 >  /**
599 >   * Write to an output stream
600 >   */
601 >  template<typename Real,  unsigned int Row, unsigned int Col>
602 >  std::ostream &operator<< ( std::ostream& o, const RectMatrix<Real, Row, Col>& m) {
603 >    for (unsigned int i = 0; i < Row ; i++) {
604 >      o << "(";
605 >      for (unsigned int j = 0; j < Col ; j++) {
606 >        o << m(i, j);
607 >        if (j != Col -1)
608 >          o << "\t";
609 >      }
610 >      o << ")" << std::endl;
611 >    }
612 >    return o;        
613 >  }    
614   }
615   #endif //MATH_RECTMATRIX_HPP

Comparing:
trunk/src/math/RectMatrix.hpp (property svn:keywords), Revision 101 by tim, Mon Oct 18 23:13:23 2004 UTC vs.
branches/development/src/math/RectMatrix.hpp (property svn:keywords), Revision 1808 by gezelter, Mon Oct 22 20:42:10 2012 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines