casacore
MatrixMathLA.h
Go to the documentation of this file.
1 //# MatrixMath.h: The Casacore linear algebra functions
2 //# Copyright (C) 1994,1995,1996,1999,2000,2002
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 //# $Id$
27 
28 #ifndef SCIMATH_MATRIXMATHLA_H
29 #define SCIMATH_MATRIXMATHLA_H
30 
31 
32 #include <casacore/casa/aips.h>
33 #include <casacore/casa/Arrays/Vector.h>
34 #include <casacore/casa/Arrays/Matrix.h>
35 #include <casacore/casa/BasicSL/Complex.h>
36 
37 
38 namespace casacore { //# NAMESPACE CASACORE - BEGIN
39 
40 //<summary>
41 // Linear algebra functions on Vectors and Matrices.
42 // </summary>
43 //
44 // <reviewed reviewer="UNKNOWN" date="before2004/08/25" tests="tLinAlgebra">
45 //
46 // <linkfrom anchor="Linear Algebra" classes="Vector Matrix">
47 // <here>Linear Algebra</here> -- Linear algebra functions
48 // on Vectors and Matrices.
49 // </linkfrom>
50 //
51 //<group name="Linear Algebra">
52 
53 // Routines which calculate the inverse of a matrix. The inverse is very
54 // often the worst way to do a calculation. Nevertheless it is often
55 // convenient. The present implementation uses LU decomposition implemented
56 // by LAPACK. The determinate can be calculated "for free" as it is the
57 // product of the diagonal terms after decomposition. If the input matrix is
58 // singular, a matrix with no rows or columns is returned. <src>in</src>
59 // must be a square matrix.
60 // <note role="warning">This function will only work for complex types if
61 // Complex and DComplex map onto their FORTRAN counterparts.</note>
62 //# We could special case small matrices for efficiency.
63 //<group>
64 template<class T> void invert(Matrix<T> & out, T& determinate,
65  const Matrix<T> &in);
66 template<class T> Matrix<T> invert(const Matrix<T> &in);
67 template<class T> T determinate(const Matrix<T> &in);
68 //</group>
69 
70 // This function inverts a symmetric positive definite matrix. It is
71 // written in C++, so it should work with any data type for which
72 // operators +, -, *, /, =, and sqrt are defined. The function uses
73 // the Cholesky decomposition method to invert the matrix. Cholesky
74 // decomposition is about a factor of 2 better than LU decomposition
75 // where symmetry is ignored.
76 template<class T> void invertSymPosDef(Matrix<T> & out, T& determinate,
77  const Matrix<T> &in);
78 template<class T> Matrix<T> invertSymPosDef(const Matrix<T> &in);
79 
80 //</group>
81 
82 //# These are actually used by invertSymPosDef. They will not
83 //# normally be called by the end user.
84 
85 //# This function performs Cholesky decomposition.
86 //# A is a positive-definite symmetric matrix. Only the upper triangle of
87 //# A is needed on input. On output, the lower triangle of A contains the
88 //# Cholesky factor L. The diagonal elements of L are returned in vector
89 //# diag.
90 template<class T> void CholeskyDecomp(Matrix<T> &A, Vector<T> &diag);
91 
92 //# Solve linear equation A*x = b, where A positive-definite symmetric.
93 //# On input, A contains Cholesky factor L in its low triangle except the
94 //# diagonal elements which are in vector diag. On return x contains the
95 //# solution. b and x can be the same vector to save memory space.
96 template<class T> void CholeskySolve(Matrix<T> &A, Vector<T> &diag,
97  Vector<T> &b, Vector<T> &x);
98 
99 
100 //# These are the LAPACK routines actually used by invert. They will not
101 //# normally be called by the end user.
102 
103 #if !defined(NEED_FORTRAN_UNDERSCORES)
104 #define NEED_FORTRAN_UNDERSCORES 1
105 #endif
106 
107 #if NEED_FORTRAN_UNDERSCORES
108 #define sgetrf sgetrf_
109 #define dgetrf dgetrf_
110 #define cgetrf cgetrf_
111 #define zgetrf zgetrf_
112 #define sgetri sgetri_
113 #define dgetri dgetri_
114 #define cgetri cgetri_
115 #define zgetri zgetri_
116 #define sposv sposv_
117 #define dposv dposv_
118 #define cposv cposv_
119 #define zposv zposv_
120 #define spotri spotri_
121 #define dpotri dpotri_
122 #define cpotri cpotri_
123 #define zpotri zpotri_
124 #endif
125 
126 extern "C" {
127  void sgetrf(const int *m, const int *n, float *a, const int *lda,
128  int *ipiv, int *info);
129  void dgetrf(const int *m, const int *n, double *a, const int *lda,
130  int *ipiv, int *info);
131  void cgetrf(const int *m, const int *n, Complex *a, const int *lda,
132  int *ipiv, int *info);
133  void zgetrf(const int *m, const int *n, DComplex *a, const int *lda,
134  int *ipiv, int *info);
135  void sgetri(const int *m, float *a, const int *lda, const int *ipiv,
136  float *work, const int *lwork, int *info);
137  void dgetri(const int *m, double *a, const int *lda, const int *ipiv,
138  double *work, const int *lwork, int *info);
139  void cgetri(const int *m, Complex *a, const int *lda, const int *ipiv,
140  Complex *work, const int *lwork, int *info);
141  void zgetri(const int *m, DComplex *a, const int *lda, const int *ipiv,
142  DComplex *work, const int *lwork, int *info);
143 
144 
145  void sposv(const char *uplo, const int *n, const int* nrhs, float *a,
146  const int *lda, float *b, const int *ldb, int *info);
147  void dposv(const char *uplo, const int *n, const int* nrhs, double *a,
148  const int *lda, double *b, const int *ldb, int *info);
149  void cposv(const char *uplo, const int *n, const int* nrhs, Complex *a,
150  const int *lda, Complex *b, const int *ldb, int *info);
151  void zposv(const char *uplo, const int *n, const int* nrhs, DComplex *a,
152  const int *lda, DComplex *b, const int *ldb, int *info);
153 
154 
155  void spotri(const char *uplo, const int *n, float *a,
156  const int *lda, int *info);
157  void dpotri(const char *uplo, const int *n, double *a,
158  const int *lda, int *info);
159  void cpotri(const char *uplo, const int *n, Complex *a,
160  const int *lda, int *info);
161  void zpotri(const char *uplo, const int *n, DComplex *a,
162  const int *lda, int *info);
163 
164 }
165 
166 //# Overloaded versions of the above to make templating work more easily
167 inline void getrf(const int *m, const int *n, float *a, const int *lda,
168  int *ipiv, int *info)
169  { sgetrf(m, n, a, lda, ipiv, info); }
170 inline void getrf(const int *m, const int *n, double *a, const int *lda,
171  int *ipiv, int *info)
172  { dgetrf(m, n, a, lda, ipiv, info); }
173 inline void getrf(const int *m, const int *n, Complex *a, const int *lda,
174  int *ipiv, int *info)
175  { cgetrf(m, n, a, lda, ipiv, info); }
176 inline void getrf(const int *m, const int *n, DComplex *a, const int *lda,
177  int *ipiv, int *info)
178  { zgetrf(m, n, a, lda, ipiv, info); }
179 inline void getri(const int *m, float *a, const int *lda, const int *ipiv,
180  float *work, const int *lwork, int *info)
181  { sgetri(m, a, lda, ipiv, work, lwork, info); }
182 inline void getri(const int *m, double *a, const int *lda, const int *ipiv,
183  double *work, const int *lwork, int *info)
184  { dgetri(m, a, lda, ipiv, work, lwork, info); }
185 inline void getri(const int *m, Complex *a, const int *lda, const int *ipiv,
186  Complex *work, const int *lwork, int *info)
187  { cgetri(m, a, lda, ipiv, work, lwork, info); }
188 inline void getri(const int *m, DComplex *a, const int *lda, const int *ipiv,
189  DComplex *work, const int *lwork, int *info)
190  { zgetri(m, a, lda, ipiv, work, lwork, info); }
191 
192 inline void posv(const char *uplo, const int *n, const int* nrhs, float *a,
193  const int *lda, float *b, const int *ldb, int *info)
194  { sposv(uplo, n, nrhs, a, lda, b, ldb, info); }
195 inline void posv(const char *uplo, const int *n, const int* nrhs, double *a,
196  const int *lda, double *b, const int *ldb, int *info)
197  { dposv(uplo, n, nrhs, a, lda, b, ldb, info); }
198 inline void posv(const char *uplo, const int *n, const int* nrhs, Complex *a,
199  const int *lda, Complex *b, const int *ldb, int *info)
200  { cposv(uplo, n, nrhs, a, lda, b, ldb, info); }
201 inline void posv(const char *uplo, const int *n, const int* nrhs, DComplex *a,
202  const int *lda, DComplex *b, const int *ldb, int *info)
203  { zposv(uplo, n, nrhs, a, lda, b, ldb, info); }
204 
205 inline void potri(const char *uplo, const int *n, float *a,
206  const int *lda, int *info)
207  { spotri(uplo, n, a, lda, info); }
208 inline void potri(const char *uplo, const int *n, double *a,
209  const int *lda, int *info)
210  { dpotri(uplo, n, a, lda, info); }
211 inline void potri(const char *uplo, const int *n, Complex *a,
212  const int *lda, int *info)
213  { cpotri(uplo, n, a, lda, info); }
214 inline void potri(const char *uplo, const int *n, DComplex *a,
215  const int *lda, int *info)
216  { zpotri(uplo, n, a, lda, info); }
217 
218 
219 
220 } //# NAMESPACE CASACORE - END
221 
222 #ifndef CASACORE_NO_AUTO_TEMPLATES
223 #include <casacore/scimath/Mathematics/MatrixMathLA.tcc>
224 #endif //# CASACORE_NO_AUTO_TEMPLATES
225 #endif
this file contains all the compiler specific defines
Definition: mainpage.dox:28
void zgetri(const int *m, DComplex *a, const int *lda, const int *ipiv, DComplex *work, const int *lwork, int *info)
void cpotri(const char *uplo, const int *n, Complex *a, const int *lda, int *info)
void spotri(const char *uplo, const int *n, float *a, const int *lda, int *info)
void sposv(const char *uplo, const int *n, const int *nrhs, float *a, const int *lda, float *b, const int *ldb, int *info)
void zgetrf(const int *m, const int *n, DComplex *a, const int *lda, int *ipiv, int *info)
void zpotri(const char *uplo, const int *n, DComplex *a, const int *lda, int *info)
void zposv(const char *uplo, const int *n, const int *nrhs, DComplex *a, const int *lda, DComplex *b, const int *ldb, int *info)
void CholeskySolve(Matrix< T > &A, Vector< T > &diag, Vector< T > &b, Vector< T > &x)
void sgetri(const int *m, float *a, const int *lda, const int *ipiv, float *work, const int *lwork, int *info)
void dposv(const char *uplo, const int *n, const int *nrhs, double *a, const int *lda, double *b, const int *ldb, int *info)
void dgetrf(const int *m, const int *n, double *a, const int *lda, int *ipiv, int *info)
void CholeskyDecomp(Matrix< T > &A, Vector< T > &diag)
void getri(const int *m, float *a, const int *lda, const int *ipiv, float *work, const int *lwork, int *info)
Definition: MatrixMathLA.h:179
void cposv(const char *uplo, const int *n, const int *nrhs, Complex *a, const int *lda, Complex *b, const int *ldb, int *info)
void potri(const char *uplo, const int *n, float *a, const int *lda, int *info)
Definition: MatrixMathLA.h:205
void cgetrf(const int *m, const int *n, Complex *a, const int *lda, int *ipiv, int *info)
void dpotri(const char *uplo, const int *n, double *a, const int *lda, int *info)
void sgetrf(const int *m, const int *n, float *a, const int *lda, int *ipiv, int *info)
void cgetri(const int *m, Complex *a, const int *lda, const int *ipiv, Complex *work, const int *lwork, int *info)
void posv(const char *uplo, const int *n, const int *nrhs, float *a, const int *lda, float *b, const int *ldb, int *info)
Definition: MatrixMathLA.h:192
void getrf(const int *m, const int *n, float *a, const int *lda, int *ipiv, int *info)
Definition: MatrixMathLA.h:167
void dgetri(const int *m, double *a, const int *lda, const int *ipiv, double *work, const int *lwork, int *info)
Matrix< T > invertSymPosDef(const Matrix< T > &in)
void invert(Matrix< T > &out, T &determinate, const Matrix< T > &in)
Routines which calculate the inverse of a matrix.
void invertSymPosDef(Matrix< T > &out, T &determinate, const Matrix< T > &in)
This function inverts a symmetric positive definite matrix.