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
38namespace 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">
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>
64template<class T> void invert(Matrix<T> & out, T& determinate,
65 const Matrix<T> &in);
66template<class T> Matrix<T> invert(const Matrix<T> &in);
67template<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.
76template<class T> void invertSymPosDef(Matrix<T> & out, T& determinate,
77 const Matrix<T> &in);
78template<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.
90template<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.
96template<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
126extern "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
167inline 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); }
170inline 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); }
173inline 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); }
176inline 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); }
179inline 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); }
182inline 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); }
185inline 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); }
188inline 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
192inline 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); }
195inline 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); }
198inline 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); }
201inline 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
205inline void potri(const char *uplo, const int *n, float *a,
206 const int *lda, int *info)
207 { spotri(uplo, n, a, lda, info); }
208inline void potri(const char *uplo, const int *n, double *a,
209 const int *lda, int *info)
210 { dpotri(uplo, n, a, lda, info); }
211inline void potri(const char *uplo, const int *n, Complex *a,
212 const int *lda, int *info)
213 { cpotri(uplo, n, a, lda, info); }
214inline 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)
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.
Matrix< T > invertSymPosDef(const Matrix< T > &in)