casacore
Loading...
Searching...
No Matches
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
37namespace casacore { //# NAMESPACE CASACORE - BEGIN
38
39//# forward declarations
40template <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
107template <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
119 const SquareMatrix<Complex,2>& left,
120 const SquareMatrix<Complex,2>& right);
121public:
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]) {
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
234private:
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>
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.
310ostream& operator<<(ostream& os, const SquareMatrix<Complex,2>& m);
311ostream& operator<<(ostream& os, const SquareMatrix<Complex,4>& m);
312ostream& operator<<(ostream& os, const SquareMatrix<Float,2>& m);
313ostream& 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
T & operator()(Int i, Int j)
Non const indexing, throws exception if you try to change an element which would require a type chang...
Matrix< T > matrix() const
Convert the SquareMatrix to a Matrix.
SquareMatrix< T, n > & operator=(const Vector< T > &v)
Assign a Vector, creates a diagonal matrix.
SquareMatrix< T, n > & operator+=(const SquareMatrix< T, n > &other)
Add two SquareMatrices, element by element.
SquareMatrix< T, n > inverse() const
Return the inverse of the matrix by value.
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(SquareMatrix< T, n > &result)
Conjugate the matrix, return it in result (and by ref)
SquareMatrix< T, n > & operator=(T val)
Assign a scalar, creates a scalar-identity matrix.
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 > & adjoint(SquareMatrix< T, n > &result)
Tranpose and conjugate the matrix, return it in result (and by ref)
SquareMatrix(const Vector< T > &vec)
Construct from Vector, creates a diagonal matrix.
T operator()(Int i, Int j) const
Indexing, only const indexing is allowed.
SquareMatrix< T, n > & operator=(const T vec[n])
Assign a c-style vector, creates a diagonal matrix.
SquareMatrix< T, n > & operator=(const Matrix< T > &m)
Assign a Matrix, creates a general matrix.
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 > & operator*=(Float f)
Scalar multiplication.
SquareMatrix< T, n > & adjoint()
Tranpose and conjugate the matrix in place(!).
Matrix< T > & matrix(Matrix< T > &result) const
Assign 'this' to the Matrix result, also return result by reference.
SquareMatrix(const T vec[n])
Construct from c-style vector, creates a diagonal matrix.
~SquareMatrix()
Destructor.
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 SquareMatrix< T, n > &m)
Copy construct a SquareMatrix, a true copy is made.
SquareMatrix(const T a[n][n])
Construct from c-style matrix (by copying elements).
SquareMatrix(const Matrix< T > &mat)
Construct from Matrix.
SquareMatrix(const T &scalar)
Construct from scalar, creates a scalar-identity matrix.
SquareMatrix()
Default constructor - creates a unity matrix at present, this may not be what we want (non-intuitive?...
SquareMatrix(int itype)
Create a matrix of a given type, no initialization.
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< T, n > & operator=(const T a[n][n])
Assign a c-style matrix, creates a general matrix.
SquareMatrix< T, n > & operator=(const SquareMatrix< T, n > &m)
Assignment, uses copy semantics.
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
SquareMatrix< Complex, 2 > conj(const SquareMatrix< Complex, 2 > &m)
Return conjugate of SquareMatrix.
SquareMatrix< Complex, 4 > directProduct(const SquareMatrix< Complex, 2 > &left, const SquareMatrix< Complex, 2 > &right)
Calculate direct product of two SquareMatrices.
SquareMatrix< Complex, 2 > adjoint(const SquareMatrix< Complex, 2 > &m)
Return adjoint of SquareMatrix.
ostream & operator<<(ostream &os, const SquareMatrix< Complex, 4 > &m)
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< Float, 2 > &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, 4 > adjoint(const SquareMatrix< Complex, 4 > &m)
Return adjoint of SquareMatrix.