casa  $Rev:20696$
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SquareMatrix.h
Go to the documentation of this file.
00001 //# SquareMatrix.h: Fast Square Matrix class with fixed (templated) size
00002 //# Copyright (C) 1996,1997,1999,2001
00003 //# Associated Universities, Inc. Washington DC, USA.
00004 //#
00005 //# This library is free software; you can redistribute it and/or modify it
00006 //# under the terms of the GNU Library General Public License as published by
00007 //# the Free Software Foundation; either version 2 of the License, or (at your
00008 //# option) any later version.
00009 //#
00010 //# This library is distributed in the hope that it will be useful, but WITHOUT
00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00012 //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00013 //# License for more details.
00014 //#
00015 //# You should have received a copy of the GNU Library General Public License
00016 //# along with this library; if not, write to the Free Software Foundation,
00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00018 //#
00019 //# Correspondence concerning AIPS++ should be adressed as follows:
00020 //#        Internet email: aips2-request@nrao.edu.
00021 //#        Postal address: AIPS++ Project Office
00022 //#                        National Radio Astronomy Observatory
00023 //#                        520 Edgemont Road
00024 //#                        Charlottesville, VA 22903-2475 USA
00025 //#
00026 //#
00027 //# $Id: SquareMatrix.h 21024 2011-03-01 11:46:18Z gervandiepen $
00028  
00029 #ifndef SCIMATH_SQUAREMATRIX_H
00030 #define SCIMATH_SQUAREMATRIX_H
00031 
00032 #include <casa/aips.h>
00033 #include <casa/BasicSL/Complex.h>
00034 #include <casa/Arrays/Matrix.h>
00035 #include <casa/iosfwd.h>
00036 
00037 namespace casa { //# NAMESPACE CASA - BEGIN
00038 
00039 //# forward declarations
00040 template <class T, Int n> class RigidVector;
00041 
00042 // <summary>
00043 // Fast Square Matrix class with fixed (templated) size
00044 // </summary>
00045 // <use visibility=export>
00046 // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
00047 // </reviewed>
00048 // <prerequisite>
00049 //   <li> Complex
00050 //   <li> Matrix
00051 // </prerequisite>
00052 //
00053 // <etymology>
00054 // SquareMatrix is a specialized class for small (<5x5) square matrices.
00055 // </etymology>
00056 //
00057 // <synopsis>
00058 // SquareMatrix provides operations similar to the Matrix class, but it is
00059 // much faster for small arrays. One important difference is that operators *=
00060 // and * do matrix products for SquareMatrices instead of element by element
00061 // multiplication. SquareMatrices also optimize operations internally for 
00062 // scalar identity matrices (diagonal matrix with all elements equal) and
00063 // diagonal matrices. The different types of SquareMatrix are created by
00064 // constructors and operator= taking either a scalar, a vector or a full
00065 // matrix.
00066 // </synopsis>
00067 //
00068 // <example>
00069 // <srcblock>
00070 // // create two SquareMatrices
00071 // SquareMatrix<Float,2> sm1(3.0); // a scalar identity matrix
00072 // Vector<Float> vec(2); vec(0)=2.0; vec(1)=3.0;
00073 // SquareMatrix<Float,2> sm2(vec); // a diagonal matrix
00074 // // multiply the matrices
00075 // // Note: A*=B is equivalent to A=A*B where '*' is matrix multiplication
00076 // sm1*=sm2; // sm1 now diagonal
00077 // </srcblock>
00078 // </example>
00079 //
00080 // <motivation>
00081 // The basic Matrix classes are rather inefficient for small sizes,
00082 // new and delete tend to dominate the execution time for computationally
00083 // intensive code. The SquareMatrix classes circumvent this by having a
00084 // compile-time fixed size c-array internally. The SquareMatrix class have
00085 // fixed zero origin and no increments, this allows fast indexing,
00086 // copying and math operations. As mentioned in the synopsis, the SquareMatrix
00087 // classes also avoid unnecessary operations for simple matrices 
00088 // (scalar-identity and diagonal).
00089 // </motivation>
00090 //
00091 // <templating arg=T>
00092 //    <li> real() is called for T=Complex/DComplex
00093 // </templating>
00094 //
00095 // <thrown>
00096 //    <li> No exceptions
00097 // </thrown>
00098 //
00099 // <todo asof="1996/11/06">
00100 //   <li> when the Sun native compiler improves, some explicit instantiations
00101 //        can be replaced by their templated equivalent and two constructor
00102 //        calls with for loops can be moved out of line.
00103 //   <li> not all operators and math functions available for Matrix are
00104 //        implemented yet, add on as-needed basis.
00105 // </todo>
00106  
00107 template <class T, Int n> class SquareMatrix {
00108     // Friends currently need to be explicit (non templated) type to work.
00109     friend class RigidVector<T,n>;
00110     //# friend class SquareMatrix<Complex,n>; // for real()
00111     //    friend class SquareMatrix<T,n*n>;// Sun native does not accept this
00112     //    friend class SquareMatrix<Complex,4>; // for directProduct of 2x2
00113     // Global friend function for product of Complex matrix and Float 4-vector
00114     friend RigidVector<Complex,4> operator*(const SquareMatrix<Complex,4>& m,
00115     const RigidVector<Float,4>& v);
00116     // Global friend function to calculate direct product
00117     friend SquareMatrix<Complex,4>& 
00118         directProduct(SquareMatrix<Complex,4>& result,
00119                       const SquareMatrix<Complex,2>& left, 
00120                       const SquareMatrix<Complex,2>& right);
00121 public:
00122     // Enum used internally to optimize operations.
00123     enum {General, Diagonal, ScalarId};
00124     // Destructor
00125     ~SquareMatrix() {}
00126     // Default constructor - creates a unity matrix at present, this may not
00127     // be what we want (non-intuitive?)
00128     SquareMatrix() : type_p(ScalarId) {a_p[0][0]=T(1);}
00129     // Create a matrix of a given type, no initialization
00130     SquareMatrix(int itype) : type_p(itype) {}
00131     // Copy construct a SquareMatrix, a true copy is made.
00132     SquareMatrix(const SquareMatrix<T,n>& m) {operator=(m);}
00133     // Construct from c-style matrix (by copying elements).
00134     SquareMatrix(const T a[n][n]) {operator=(a);}
00135     // Construct from Matrix.
00136     SquareMatrix(const Matrix<T>& mat) {operator=(mat);}
00137     // Construct from c-style vector, creates a diagonal matrix.
00138     SquareMatrix(const T vec[n]){operator=(vec);}
00139     // Construct from Vector, creates a diagonal matrix.
00140     SquareMatrix(const Vector<T>& vec) {operator=(vec);}
00141     // Construct from scalar, creates a scalar-identity matrix
00142     SquareMatrix(const T& scalar) : type_p(ScalarId) { a_p[0][0]=scalar; }
00143     // Assignment, uses copy semantics.
00144     SquareMatrix<T,n>& operator=(const SquareMatrix<T,n>& m);
00145     // Assign a c-style matrix, creates a general matrix.
00146     SquareMatrix<T,n>& operator=(const T a[n][n]) {
00147         type_p=General;
00148         const T* pa=&a[0][0];
00149         T* pa_p=&a_p[0][0];
00150         for (Int i=0; i<n*n; i++) *pa_p++=*pa++;
00151         return *this;
00152     }   
00153     // Assign a Matrix, creates a general matrix.
00154     SquareMatrix<T,n>& operator=(const Matrix<T>& m);
00155     // Assign a c-style vector, creates a diagonal matrix
00156     SquareMatrix<T,n>& operator=(const T vec[n]) {
00157         type_p=Diagonal;
00158         for (Int i=0; i<n; i++) a_p[i][i]=vec[i];
00159         return *this;
00160     } 
00161     // Assign a Vector, creates a diagonal matrix
00162     SquareMatrix<T,n>& operator=(const Vector<T>& v);
00163     // Assign a scalar, creates a scalar-identity matrix
00164     SquareMatrix<T,n>& operator=(T val) {
00165         type_p=ScalarId; a_p[0][0]=val; return *this;
00166     }
00167     // Add two SquareMatrices, element by element.
00168     SquareMatrix<T,n>& operator+=(const SquareMatrix<T,n>& other);
00169     // Matrix product of 'this' SquareMatrix with other,
00170     // i.e., A*=B; is equivalent with A=A*B where '*' is matrix multiplication.
00171     SquareMatrix<T,n>& operator*=(const SquareMatrix<T,n>& other);
00172     // Scalar multiplication
00173     SquareMatrix<T,n>& operator*=(Float f);
00174     // Indexing, only const indexing is allowed. You cannot change the
00175     // matrix via indexing. No bounds checking.
00176     T operator()(Int i, Int j) const { 
00177       switch (type_p) {
00178       case ScalarId: return (i==j) ? a_p[0][0] : T();
00179         break;
00180       case Diagonal: return (i==j) ? a_p[i][i] : T();
00181         break;
00182       }
00183       return a_p[i][j];
00184     }
00185     // Non const indexing, throws exception if you try to change an element
00186     // which would require a type change of the matrix
00187     T& operator()(Int i, Int j) {
00188       switch (type_p) {
00189       case ScalarId: return (i==j) ? a_p[0][0] : throwInvAccess();
00190         break;
00191       case Diagonal: return (i==j) ? a_p[i][i] : throwInvAccess();
00192         break;
00193       }
00194       return a_p[i][j];
00195     }
00196 
00197     //# The following does not compile with Sun native, replaced by explicit
00198     //# global function.
00199     //# direct product : dp= this (directproduct) other
00200     //#    SquareMatrix<T,n*n>& 
00201     //# directProduct(SquareMatrix<T,n*n>& dp, 
00202     //#               const SquareMatrix<T,n>& other) const;
00203     // For a <src>SquareMatrix<Complex,n></src>: 
00204     // set the argument result to the real part of the matrix 
00205     // (and return result by reference to allow use in
00206     // expressions without creating temporary).
00207     //# SquareMatrix<Float,n>& real(SquareMatrix<Float,n>& result) const;
00208     // For a <src>SquareMatrix<Complex,n></src>: 
00209     // return the real part of the matrix.
00210     //# SquareMatrix<Float,n> real() const {
00211     //# SquareMatrix<Float,n> result;
00212     //# return real(result);
00213     //# }
00214     // Conjugate the matrix in place(!).
00215     SquareMatrix<T,n>& conj();
00216     // Tranpose and conjugate the matrix in place(!).
00217     SquareMatrix<T,n>& adjoint(); 
00218     // Conjugate the matrix, return it in result (and by ref)
00219     SquareMatrix<T,n>& conj(SquareMatrix<T,n>& result);
00220     // Tranpose and conjugate the matrix, return it in result (and by ref)
00221     SquareMatrix<T,n>& adjoint(SquareMatrix<T,n>& result);
00222     // Compute the inverse of the matrix and return it in result (also
00223     // returns result by reference).
00224     SquareMatrix<T,n>& inverse(SquareMatrix<T,n>& result) const;
00225     // Return the inverse of the matrix by value.
00226     SquareMatrix<T,n> inverse() const
00227         { SquareMatrix<T,n> result; return inverse(result);}
00228     // Assign 'this' to the Matrix result, also return result by reference.
00229     Matrix<T>& matrix(Matrix<T>& result) const;
00230     // Convert the SquareMatrix to a Matrix.
00231     Matrix<T> matrix() const
00232         { Matrix<T> result(n,n); return matrix(result);}
00233 
00234 private:
00235     T& throwInvAccess();
00236     T a_p[n][n];
00237     Int type_p;
00238 };
00239 
00240 
00241 //# the following does not compile with Sun native but should...
00242 //# expanded by hand for types and sizes needed
00243 //#template<class T, Int n> 
00244 //#ostream& operator<<(ostream& os, const SquareMatrix<T,n>& m) {
00245 //#   return os<<m.matrix();
00246 //#}
00247 //#
00248 //#template<class T, Int n> inline SquareMatrix<T,n> operator+(const SquareMatrix<T,n>& left, 
00249 //#                                                const SquareMatrix<T,n>& right) {
00250 //#    SquareMatrix<T,n> result(left);
00251 //#    return result+=right;
00252 //#}
00253 //#template<class T, Int n> inline SquareMatrix<T,n> operator*(const SquareMatrix<T,n>& left, 
00254 //#                                                const SquareMatrix<T,n>& right) 
00255 //#{
00256 //#    SquareMatrix<T,n> result(left);
00257 //#    return result*=right;
00258 //#}
00259 //#
00260 //#template<class T, Int n> inline SquareMatrix<T,n*n> directProduct(const SquareMatrix<T,n>& left, 
00261 //#                                                      const SquareMatrix<T,n>& right)
00262 //#{
00263 //#    SquareMatrix<T,n*n> result;
00264 //#    return left.directProduct(result,right);
00265 //#}
00266 //#
00267 //#template<class T, Int n> inline SquareMatrix<T,n*n>& 
00268 //#directProduct(SquareMatrix<T,n*n>& result,
00269 //#           const SquareMatrix<T,n>& left, 
00270 //#           const SquareMatrix<T,n>& right)
00271 //#{
00272 //#    return left.directProduct(result,right);
00273 //#}
00274 //#template<class T, Int n> inline SquareMatrix<T,n> conj(
00275 //#    const SquareMatrix<T,n>& m) {
00276 //#    SquareMatrix<T,n> result(m);
00277 //#    return result.conj();
00278 //#}
00279 //#
00280 //#template<class T, Int n> inline SquareMatrix<T,n> adjoint(
00281 //#    const SquareMatrix<T,n>& m) {
00282 //#    SquareMatrix<T,n> result(m);
00283 //#    return result.adjoint();
00284 //#}
00285 //#
00286 
00287 // <summary>
00288 // Various global math and IO functions.
00289 // </summary>
00290 // <group name=SqM_global_functions>
00291 
00292 // Calculate direct product of two SquareMatrices.
00293 SquareMatrix<Complex,4> 
00294 directProduct(const SquareMatrix<Complex,2>& left, 
00295               const SquareMatrix<Complex,2>& right);
00296 
00297 // Return conjugate of SquareMatrix.
00298 SquareMatrix<Complex,2> conj(const SquareMatrix<Complex,2>& m);
00299 
00300 // Return conjugate of SquareMatrix.
00301 SquareMatrix<Complex,4> conj(const SquareMatrix<Complex,4>& m);
00302 
00303 // Return adjoint of SquareMatrix.
00304 SquareMatrix<Complex,2> adjoint(const SquareMatrix<Complex,2>& m);
00305 
00306 // Return adjoint of SquareMatrix.
00307 SquareMatrix<Complex,4> adjoint(const SquareMatrix<Complex,4>& m);
00308 
00309 // Write SquareMatrix to output, uses Matrix to do the work.
00310 ostream& operator<<(ostream& os, const SquareMatrix<Complex,2>& m);
00311 ostream& operator<<(ostream& os, const SquareMatrix<Complex,4>& m);
00312 ostream& operator<<(ostream& os, const SquareMatrix<Float,2>& m);
00313 ostream& operator<<(ostream& os, const SquareMatrix<Float,4>& m);
00314 // </group>
00315 
00316 
00317 
00318 } //# NAMESPACE CASA - END
00319 
00320 #ifndef CASACORE_NO_AUTO_TEMPLATES
00321 #include <scimath/Mathematics/SquareMatrix.tcc>
00322 #endif //# CASACORE_NO_AUTO_TEMPLATES
00323 #endif