casacore
SquareMatrix.h
Go to the documentation of this file.
1 //# SquareMatrix.h: Fast Square Matrix class with fixed (templated) size
2 //# Copyright (C) 1996,1997,1999,2001
3 //# Associated Universities, Inc. Washington DC, USA.
4 //#
5 //# This library is free software; you can redistribute it and/or modify it
6 //# under the terms of the GNU Library General Public License as published by
7 //# the Free Software Foundation; either version 2 of the License, or (at your
8 //# option) any later version.
9 //#
10 //# This library is distributed in the hope that it will be useful, but WITHOUT
11 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 //# License for more details.
14 //#
15 //# You should have received a copy of the GNU Library General Public License
16 //# along with this library; if not, write to the Free Software Foundation,
17 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 //#
19 //# Correspondence concerning AIPS++ should be addressed as follows:
20 //# Internet email: aips2-request@nrao.edu.
21 //# Postal address: AIPS++ Project Office
22 //# National Radio Astronomy Observatory
23 //# 520 Edgemont Road
24 //# Charlottesville, VA 22903-2475 USA
25 //#
26 //#
27 //# $Id$
28 
29 #ifndef SCIMATH_SQUAREMATRIX_H
30 #define SCIMATH_SQUAREMATRIX_H
31 
32 #include <casacore/casa/aips.h>
33 #include <casacore/casa/BasicSL/Complex.h>
34 #include <casacore/casa/Arrays/Matrix.h>
35 #include <casacore/casa/iosfwd.h>
36 
37 namespace casacore { //# NAMESPACE CASACORE - BEGIN
38 
39 //# forward declarations
40 template <class T, Int n> class RigidVector;
41 
42 // <summary>
43 // Fast Square Matrix class with fixed (templated) size
44 // </summary>
45 // <use visibility=export>
46 // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
47 // </reviewed>
48 // <prerequisite>
49 // <li> Complex
50 // <li> Matrix
51 // </prerequisite>
52 //
53 // <etymology>
54 // SquareMatrix is a specialized class for small (<5x5) square matrices.
55 // </etymology>
56 //
57 // <synopsis>
58 // SquareMatrix provides operations similar to the Matrix class, but it is
59 // much faster for small arrays. One important difference is that operators *=
60 // and * do matrix products for SquareMatrices instead of element by element
61 // multiplication. SquareMatrices also optimize operations internally for
62 // scalar identity matrices (diagonal matrix with all elements equal) and
63 // diagonal matrices. The different types of SquareMatrix are created by
64 // constructors and operator= taking either a scalar, a vector or a full
65 // matrix.
66 // </synopsis>
67 //
68 // <example>
69 // <srcblock>
70 // // create two SquareMatrices
71 // SquareMatrix<Float,2> sm1(3.0); // a scalar identity matrix
72 // Vector<Float> vec(2); vec(0)=2.0; vec(1)=3.0;
73 // SquareMatrix<Float,2> sm2(vec); // a diagonal matrix
74 // // multiply the matrices
75 // // Note: A*=B is equivalent to A=A*B where '*' is matrix multiplication
76 // sm1*=sm2; // sm1 now diagonal
77 // </srcblock>
78 // </example>
79 //
80 // <motivation>
81 // The basic Matrix classes are rather inefficient for small sizes,
82 // new and delete tend to dominate the execution time for computationally
83 // intensive code. The SquareMatrix classes circumvent this by having a
84 // compile-time fixed size c-array internally. The SquareMatrix class have
85 // fixed zero origin and no increments, this allows fast indexing,
86 // copying and math operations. As mentioned in the synopsis, the SquareMatrix
87 // classes also avoid unnecessary operations for simple matrices
88 // (scalar-identity and diagonal).
89 // </motivation>
90 //
91 // <templating arg=T>
92 // <li> real() is called for T=Complex/DComplex
93 // </templating>
94 //
95 // <thrown>
96 // <li> No exceptions
97 // </thrown>
98 //
99 // <todo asof="1996/11/06">
100 // <li> when the Sun native compiler improves, some explicit instantiations
101 // can be replaced by their templated equivalent and two constructor
102 // calls with for loops can be moved out of line.
103 // <li> not all operators and math functions available for Matrix are
104 // implemented yet, add on as-needed basis.
105 // </todo>
106 
107 template <class T, Int n> class SquareMatrix {
108  // Friends currently need to be explicit (non templated) type to work.
109  friend class RigidVector<T,n>;
110  //# friend class SquareMatrix<Complex,n>; // for real()
111  // friend class SquareMatrix<T,n*n>;// Sun native does not accept this
112  // friend class SquareMatrix<Complex,4>; // for directProduct of 2x2
113  // Global friend function for product of Complex matrix and Float 4-vector
115  const RigidVector<Float,4>& v);
116  // Global friend function to calculate direct product
117  friend SquareMatrix<Complex,4>&
119  const SquareMatrix<Complex,2>& left,
120  const SquareMatrix<Complex,2>& right);
121 public:
122  // Enum used internally to optimize operations.
124  // Destructor
126  // Default constructor - creates a unity matrix at present, this may not
127  // be what we want (non-intuitive?)
128  SquareMatrix() : type_p(ScalarId) {a_p[0][0]=T(1);}
129  // Create a matrix of a given type, no initialization
130  SquareMatrix(int itype) : type_p(itype) {}
131  // Copy construct a SquareMatrix, a true copy is made.
133  // Construct from c-style matrix (by copying elements).
134  SquareMatrix(const T a[n][n]) {operator=(a);}
135  // Construct from Matrix.
136  SquareMatrix(const Matrix<T>& mat) {operator=(mat);}
137  // Construct from c-style vector, creates a diagonal matrix.
138  SquareMatrix(const T vec[n]){operator=(vec);}
139  // Construct from Vector, creates a diagonal matrix.
140  SquareMatrix(const Vector<T>& vec) {operator=(vec);}
141  // Construct from scalar, creates a scalar-identity matrix
142  SquareMatrix(const T& scalar) : type_p(ScalarId) { a_p[0][0]=scalar; }
143  // Assignment, uses copy semantics.
145  // Assign a c-style matrix, creates a general matrix.
146  SquareMatrix<T,n>& operator=(const T a[n][n]) {
147  type_p=General;
148  const T* pa=&a[0][0];
149  T* pa_p=&a_p[0][0];
150  for (Int i=0; i<n*n; i++) *pa_p++=*pa++;
151  return *this;
152  }
153  // Assign a Matrix, creates a general matrix.
155  // Assign a c-style vector, creates a diagonal matrix
156  SquareMatrix<T,n>& operator=(const T vec[n]) {
158  for (Int i=0; i<n; i++) a_p[i][i]=vec[i];
159  return *this;
160  }
161  // Assign a Vector, creates a diagonal matrix
163  // Assign a scalar, creates a scalar-identity matrix
165  type_p=ScalarId; a_p[0][0]=val; return *this;
166  }
167  // Add two SquareMatrices, element by element.
169  // Matrix product of 'this' SquareMatrix with other,
170  // i.e., A*=B; is equivalent with A=A*B where '*' is matrix multiplication.
172  // Scalar multiplication
174  // Indexing, only const indexing is allowed. You cannot change the
175  // matrix via indexing. No bounds checking.
176  T operator()(Int i, Int j) const {
177  switch (type_p) {
178  case ScalarId: return (i==j) ? a_p[0][0] : T();
179  break;
180  case Diagonal: return (i==j) ? a_p[i][i] : T();
181  break;
182  }
183  return a_p[i][j];
184  }
185  // Non const indexing, throws exception if you try to change an element
186  // which would require a type change of the matrix
187  T& operator()(Int i, Int j) {
188  switch (type_p) {
189  case ScalarId: return (i==j) ? a_p[0][0] : throwInvAccess();
190  break;
191  case Diagonal: return (i==j) ? a_p[i][i] : throwInvAccess();
192  break;
193  }
194  return a_p[i][j];
195  }
196 
197  //# The following does not compile with Sun native, replaced by explicit
198  //# global function.
199  //# direct product : dp= this (directproduct) other
200  //# SquareMatrix<T,n*n>&
201  //# directProduct(SquareMatrix<T,n*n>& dp,
202  //# const SquareMatrix<T,n>& other) const;
203  // For a <src>SquareMatrix<Complex,n></src>:
204  // set the argument result to the real part of the matrix
205  // (and return result by reference to allow use in
206  // expressions without creating temporary).
207  //# SquareMatrix<Float,n>& real(SquareMatrix<Float,n>& result) const;
208  // For a <src>SquareMatrix<Complex,n></src>:
209  // return the real part of the matrix.
210  //# SquareMatrix<Float,n> real() const {
211  //# SquareMatrix<Float,n> result;
212  //# return real(result);
213  //# }
214  // Conjugate the matrix in place(!).
216  // Tranpose and conjugate the matrix in place(!).
218  // Conjugate the matrix, return it in result (and by ref)
220  // Tranpose and conjugate the matrix, return it in result (and by ref)
222  // Compute the inverse of the matrix and return it in result (also
223  // returns result by reference).
225  // Return the inverse of the matrix by value.
227  { SquareMatrix<T,n> result; return inverse(result);}
228  // Assign 'this' to the Matrix result, also return result by reference.
229  Matrix<T>& matrix(Matrix<T>& result) const;
230  // Convert the SquareMatrix to a Matrix.
232  { Matrix<T> result(n,n); return matrix(result);}
233 
234 private:
236  T a_p[n][n];
238 };
239 
240 
241 //# the following does not compile with Sun native but should...
242 //# expanded by hand for types and sizes needed
243 //#template<class T, Int n>
244 //#ostream& operator<<(ostream& os, const SquareMatrix<T,n>& m) {
245 //# return os<<m.matrix();
246 //#}
247 //#
248 //#template<class T, Int n> inline SquareMatrix<T,n> operator+(const SquareMatrix<T,n>& left,
249 //# const SquareMatrix<T,n>& right) {
250 //# SquareMatrix<T,n> result(left);
251 //# return result+=right;
252 //#}
253 //#template<class T, Int n> inline SquareMatrix<T,n> operator*(const SquareMatrix<T,n>& left,
254 //# const SquareMatrix<T,n>& right)
255 //#{
256 //# SquareMatrix<T,n> result(left);
257 //# return result*=right;
258 //#}
259 //#
260 //#template<class T, Int n> inline SquareMatrix<T,n*n> directProduct(const SquareMatrix<T,n>& left,
261 //# const SquareMatrix<T,n>& right)
262 //#{
263 //# SquareMatrix<T,n*n> result;
264 //# return left.directProduct(result,right);
265 //#}
266 //#
267 //#template<class T, Int n> inline SquareMatrix<T,n*n>&
268 //#directProduct(SquareMatrix<T,n*n>& result,
269 //# const SquareMatrix<T,n>& left,
270 //# const SquareMatrix<T,n>& right)
271 //#{
272 //# return left.directProduct(result,right);
273 //#}
274 //#template<class T, Int n> inline SquareMatrix<T,n> conj(
275 //# const SquareMatrix<T,n>& m) {
276 //# SquareMatrix<T,n> result(m);
277 //# return result.conj();
278 //#}
279 //#
280 //#template<class T, Int n> inline SquareMatrix<T,n> adjoint(
281 //# const SquareMatrix<T,n>& m) {
282 //# SquareMatrix<T,n> result(m);
283 //# return result.adjoint();
284 //#}
285 //#
286 
287 // <summary>
288 // Various global math and IO functions.
289 // </summary>
290 // <group name=SqM_global_functions>
291 
292 // Calculate direct product of two SquareMatrices.
295  const SquareMatrix<Complex,2>& right);
296 
297 // Return conjugate of SquareMatrix.
299 
300 // Return conjugate of SquareMatrix.
302 
303 // Return adjoint of SquareMatrix.
305 
306 // Return adjoint of SquareMatrix.
308 
309 // Write SquareMatrix to output, uses Matrix to do the work.
310 ostream& operator<<(ostream& os, const SquareMatrix<Complex,2>& m);
311 ostream& operator<<(ostream& os, const SquareMatrix<Complex,4>& m);
312 ostream& operator<<(ostream& os, const SquareMatrix<Float,2>& m);
313 ostream& operator<<(ostream& os, const SquareMatrix<Float,4>& m);
314 // </group>
315 
316 
317 
318 } //# NAMESPACE CASACORE - END
319 
320 #ifndef CASACORE_NO_AUTO_TEMPLATES
321 #include <casacore/scimath/Mathematics/SquareMatrix.tcc>
322 #endif //# CASACORE_NO_AUTO_TEMPLATES
323 #endif
friend SquareMatrix< Complex, 4 > & directProduct(SquareMatrix< Complex, 4 > &result, const SquareMatrix< Complex, 2 > &left, const SquareMatrix< Complex, 2 > &right)
Global friend function to calculate direct product.
SquareMatrix< T, n > & conj()
For a SquareMatrix<Complex,n>: set the argument result to the real part of the matrix (and return res...
SquareMatrix(const Vector< T > &vec)
Construct from Vector, creates a diagonal matrix.
Definition: SquareMatrix.h:140
T operator()(Int i, Int j) const
Indexing, only const indexing is allowed.
Definition: SquareMatrix.h:176
SquareMatrix< T, n > & conj(SquareMatrix< T, n > &result)
Conjugate the matrix, return it in result (and by ref)
SquareMatrix< T, n > & operator*=(Float f)
Scalar multiplication.
SquareMatrix< T, n > & operator=(const Vector< T > &v)
Assign a Vector, creates a diagonal matrix.
SquareMatrix< T, n > & adjoint()
Tranpose and conjugate the matrix in place(!).
SquareMatrix< T, n > & operator=(const Matrix< T > &m)
Assign a Matrix, creates a general matrix.
T & operator()(Int i, Int j)
Non const indexing, throws exception if you try to change an element which would require a type chang...
Definition: SquareMatrix.h:187
SquareMatrix< T, n > & operator+=(const SquareMatrix< T, n > &other)
Add two SquareMatrices, element by element.
Matrix< T > & matrix(Matrix< T > &result) const
Assign 'this' to the Matrix result, also return result by reference.
friend RigidVector< Complex, 4 > operator*(const SquareMatrix< Complex, 4 > &m, const RigidVector< Float, 4 > &v)
friend class SquareMatrix<T,n*n>;// Sun native does not accept this friend class SquareMatrix<Complex...
SquareMatrix< T, n > & inverse(SquareMatrix< T, n > &result) const
Compute the inverse of the matrix and return it in result (also returns result by reference).
SquareMatrix(const T vec[n])
Construct from c-style vector, creates a diagonal matrix.
Definition: SquareMatrix.h:138
~SquareMatrix()
Destructor.
Definition: SquareMatrix.h:125
SquareMatrix< T, n > inverse() const
Return the inverse of the matrix by value.
Definition: SquareMatrix.h:226
SquareMatrix< T, n > & operator=(const T a[n][n])
Assign a c-style matrix, creates a general matrix.
Definition: SquareMatrix.h:146
SquareMatrix(const SquareMatrix< T, n > &m)
Copy construct a SquareMatrix, a true copy is made.
Definition: SquareMatrix.h:132
SquareMatrix< T, n > & adjoint(SquareMatrix< T, n > &result)
Tranpose and conjugate the matrix, return it in result (and by ref)
SquareMatrix(const T a[n][n])
Construct from c-style matrix (by copying elements).
Definition: SquareMatrix.h:134
Matrix< T > matrix() const
Convert the SquareMatrix to a Matrix.
Definition: SquareMatrix.h:231
SquareMatrix< T, n > & operator=(const T vec[n])
Assign a c-style vector, creates a diagonal matrix.
Definition: SquareMatrix.h:156
SquareMatrix(const Matrix< T > &mat)
Construct from Matrix.
Definition: SquareMatrix.h:136
SquareMatrix(const T &scalar)
Construct from scalar, creates a scalar-identity matrix.
Definition: SquareMatrix.h:142
SquareMatrix()
Default constructor - creates a unity matrix at present, this may not be what we want (non-intuitive?...
Definition: SquareMatrix.h:128
SquareMatrix< T, n > & operator=(const SquareMatrix< T, n > &m)
Assignment, uses copy semantics.
SquareMatrix(int itype)
Create a matrix of a given type, no initialization.
Definition: SquareMatrix.h:130
SquareMatrix< T, n > & operator*=(const SquareMatrix< T, n > &other)
Matrix product of 'this' SquareMatrix with other, i.e., A*=B; is equivalent with A=A*B where '*' is m...
SquareMatrix< T, n > & operator=(T val)
Assign a scalar, creates a scalar-identity matrix.
Definition: SquareMatrix.h:164
this file contains all the compiler specific defines
Definition: mainpage.dox:28
LatticeExprNode pa(const LatticeExprNode &left, const LatticeExprNode &right)
This function finds 180/pi*atan2(left,right)/2.
float Float
Definition: aipstype.h:54
int Int
Definition: aipstype.h:50
ostream & operator<<(ostream &os, const SquareMatrix< Complex, 2 > &m)
Write SquareMatrix to output, uses Matrix to do the work.
ostream & operator<<(ostream &os, const SquareMatrix< Complex, 4 > &m)
SquareMatrix< Complex, 4 > conj(const SquareMatrix< Complex, 4 > &m)
Return conjugate of SquareMatrix.
ostream & operator<<(ostream &os, const SquareMatrix< Float, 4 > &m)
SquareMatrix< Complex, 2 > adjoint(const SquareMatrix< Complex, 2 > &m)
Return adjoint of SquareMatrix.
SquareMatrix< Complex, 2 > conj(const SquareMatrix< Complex, 2 > &m)
Return conjugate of SquareMatrix.
ostream & operator<<(ostream &os, const SquareMatrix< Float, 2 > &m)
SquareMatrix< Complex, 4 > adjoint(const SquareMatrix< Complex, 4 > &m)
Return adjoint of SquareMatrix.
SquareMatrix< Complex, 4 > directProduct(const SquareMatrix< Complex, 2 > &left, const SquareMatrix< Complex, 2 > &right)
Calculate direct product of two SquareMatrices.