Subversion Repositories gelsvn

Rev

Rev 379 | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 379 Rev 393
1
#if !defined(LAPACKFUNC_H_HAA_AGUST_2001)
1
#if !defined(LAPACKFUNC_H_HAA_AGUST_2001)
2
#define LAPACKFUNC_H_HAA_AGUST_2001
2
#define LAPACKFUNC_H_HAA_AGUST_2001
3
 
3
 
4
#if defined(_MSC_VER)
4
#if defined(_MSC_VER)
5
#pragma message("Note: including lib: clapack.lib and ignoring defaultlib : LIBC\n")
5
#pragma message("Note: including lib: clapack.lib and ignoring defaultlib : LIBC\n")
6
#pragma comment(lib, "clapack.lib")
6
#pragma comment(lib, "clapack.lib")
7
#pragma comment(linker, "/NODEFAULTLIB:LIBC.LIB")
7
#pragma comment(linker, "/NODEFAULTLIB:LIBC.LIB")
8
 
8
 
9
#endif
9
#endif
10
 
10
 
11
#include "Matrix.h"
11
#include "Matrix.h"
12
#include "Vector.h"
12
#include "Vector.h"
13
 
13
 
14
namespace LinAlg
14
namespace LinAlg
15
{
15
{
16
 
16
 
17
/*!
17
/*!
18
\file LapackFunc.h
18
\file LapackFunc.h
19
\brief Interface to some of the LAPACK functionality.
19
\brief Interface to some of the LAPACK functionality.
20
 
20
 
21
These are functions which more or less directly interface with the
21
These are functions which more or less directly interface with the
22
Lapack provided algorithms.
22
Lapack provided algorithms.
23
 
23
 
24
For indepth reference to the LAPACK functions see:
24
For indepth reference to the LAPACK functions see:
25
  LAPACK Users' Guide - 3rd Edition,
25
  LAPACK Users' Guide - 3rd Edition,
26
  by E. Anderson et al.,
26
  by E. Anderson et al.,
27
  ISBN 0-89871-447-8,
27
  ISBN 0-89871-447-8,
28
  Published by SIAM,
28
  Published by SIAM,
29
 
29
 
30
This book is also available at: \URL{http://www.netlib.org/lapack/lug/lapack_lug.html}
30
This book is also available at: \URL{http://www.netlib.org/lapack/lug/lapack_lug.html}
31
 
31
 
32
The official LAPACK sites where from the source can be downloaded are:
32
The official LAPACK sites where from the source can be downloaded are:
33
  \URL{http://www.netlib.org/clapack/} and
33
  \URL{http://www.netlib.org/clapack/} and
34
  \URL{http://www.netlib.org/lapack/}
34
  \URL{http://www.netlib.org/lapack/}
35
 
35
 
36
 
36
 
37
NB: When running this in MS Visual C++ it is usually required to set the 
37
NB: When running this in MS Visual C++ it is usually required to set the 
38
multithread "\MD" compiler option. This is to ensure correct linkage to the 
38
multithread "\MD" compiler option. This is to ensure correct linkage to the 
39
precompiled library "clapack.lib" and/or "clapackDB.lib".
39
precompiled library "clapack.lib" and/or "clapackDB.lib".
40
 
40
 
41
 
41
 
42
\author  Henrik Aanæs
42
\author  Henrik Aanæs
43
\version Aug 2001
43
\version Aug 2001
44
*/
44
*/
45
 
45
 
46
/*!
46
/*!
47
\name Singular Value Decomposition SVD 
47
\name Singular Value Decomposition SVD 
48
 
48
 
49
These functions perform the Singular Value Decomposition SVD of 
49
These functions perform the Singular Value Decomposition SVD of 
50
the MxN matrix A. The SVD is defined by:
50
the MxN matrix A. The SVD is defined by:
51
 
51
 
52
  A=U*S*V^T
52
  A=U*S*V^T
53
 
53
 
54
where:
54
where:
55
- U is a M by M orthogonal matrix 
55
- U is a M by M orthogonal matrix 
56
- V is a N by N orthogonal matrix 
56
- V is a N by N orthogonal matrix 
57
- S is a M by N diaggonal matrix. The values in the diagonal are the singular values
57
- S is a M by N diaggonal matrix. The values in the diagonal are the singular values
58
 
58
 
59
 
59
 
60
\param  A the matrix to perform SVD on  
60
\param  A the matrix to perform SVD on  
61
\return U will be resized if it is does not have the correct dimensions
61
\return U will be resized if it is does not have the correct dimensions
62
\return V will be resized if it is does not have the correct dimensions
62
\return V will be resized if it is does not have the correct dimensions
63
\return S will be resized if it is does not have the correct dimensions. 
63
\return S will be resized if it is does not have the correct dimensions. 
64
\exception assert(info==0) for Lapack. Add a throw statement later.
64
\exception assert(info==0) for Lapack. Add a throw statement later.
65
\version  Aug 2001 
65
\version  Aug 2001 
66
\author  Henrik Aanæs
66
\author  Henrik Aanæs
67
*/  
67
*/  
68
//@{ 
68
//@{ 
69
///SVD of A, where the singular values are returned in a Vector.
69
///SVD of A, where the singular values are returned in a Vector.
70
void SVD(const CMatrix& A,CMatrix& U,CVector& s,CMatrix& V);
70
void SVD(const CMatrix& A,CMatrix& U,CVector& s,CMatrix& V);
71
///SVD of A, where the singular values are returned in a 'diagonal' Matrix.
71
///SVD of A, where the singular values are returned in a 'diagonal' Matrix.
72
void SVD(const CMatrix& A,CMatrix& U,CMatrix& S,CMatrix& V);
72
void SVD(const CMatrix& A,CMatrix& U,CMatrix& S,CMatrix& V);
73
///SVD of A, returning only the singular values in a Vector.
73
///SVD of A, returning only the singular values in a Vector.
74
CVector SVD(const CMatrix& A);
74
CVector SVD(const CMatrix& A);
75
//@}
75
//@}
76
 
76
 
77
 
77
 
78
/*!
78
/*!
79
\name Linear Equations
79
\name Linear Equations
80
These functions solve the system of linear equations
80
These functions solve the system of linear equations
81
 
81
 
82
  A*x=b 
82
  A*x=b 
83
 
83
 
84
for x, where:
84
for x, where:
85
- A is a N by N matrix 
85
- A is a N by N matrix 
86
- b is a N vector
86
- b is a N vector
87
- x is a N vector
87
- x is a N vector
88
 
88
 
89
There a speceilaized functions for symetric positive definite (SPD) 
89
There a speceilaized functions for symetric positive definite (SPD) 
90
matrices yeilding better performance. These are denote by SPD in 
90
matrices yeilding better performance. These are denote by SPD in 
91
there function name.
91
there function name.
92
 
92
 
93
\param A the NxN square matrix
93
\param A the NxN square matrix
94
\param b the N vector
94
\param b the N vector
95
\return x will be resized if it is does not have the correct dimensions
95
\return x will be resized if it is does not have the correct dimensions
96
\exception assert(info==0) for Lapack. Add a throw statement later.
96
\exception assert(info==0) for Lapack. Add a throw statement later.
97
\exception assert(A.Row()==A.Col()). Add a throw statement later.
97
\exception assert(A.Row()==A.Col()). Add a throw statement later.
98
\exception assert(A.Row()==b.Length()). Add a throw statement later.
98
\exception assert(A.Row()==b.Length()). Add a throw statement later.
99
\version  Aug 2001 
99
\version  Aug 2001 
100
\author  Henrik Aanæs
100
\author  Henrik Aanæs
101
*/  
101
*/  
102
//@{
102
//@{
103
///Solves Ax=b for x.
103
///Solves Ax=b for x.
104
void LinearSolve(const CMatrix& A,const CVector&b,CVector& x);
104
void LinearSolve(const CMatrix& A,const CVector&b,CVector& x);
105
///Solves Ax=b for x and returns x.
105
///Solves Ax=b for x and returns x.
106
CVector LinearSolve(const CMatrix& A,const CVector&b);
106
CVector LinearSolve(const CMatrix& A,const CVector&b);
107
///Solves Ax=b for x, where A is SPD.
107
///Solves Ax=b for x, where A is SPD.
108
void LinearSolveSPD(const CMatrix& A,const CVector&b,CVector& x);
108
void LinearSolveSPD(const CMatrix& A,const CVector&b,CVector& x);
109
///Solves Ax=b for x and returns x, where A is SPD.
109
///Solves Ax=b for x and returns x, where A is SPD.
110
CVector LinearSolveSPD(const CMatrix& A,const CVector&b);
110
CVector LinearSolveSPD(const CMatrix& A,const CVector&b);
111
//@}
111
//@}
112
 
112
 
113
void LinearSolveSym(const CMatrix& A,
113
void LinearSolveSym(const CMatrix& A,
114
										const CVector&b,
114
										const CVector&b,
115
										CVector& x);
115
										CVector& x);
116
 
116
 
117
/**
117
/**
118
\name Linear Least Squares
118
\name Linear Least Squares
119
These functions solve the Linear Least Squares problem:
119
These functions solve the Linear Least Squares problem:
120
 
120
 
121
  min_x ||Ax-b||^2
121
  min_x ||Ax-b||^2
122
 
122
 
123
for x, where:
123
for x, where:
124
- || || denotes the 2-norm
124
- || || denotes the 2-norm
125
- A	is a M by N matrix. For a well formed M>=N and rank (A)=N. See below.
125
- A	is a M by N matrix. For a well formed M>=N and rank (A)=N. See below.
126
- b	is a M vector.
126
- b	is a M vector.
127
- x	is a N vector 
127
- x	is a N vector 
128
 
128
 
129
If the solution is not \em well \em formed the algorithm provided will find a
129
If the solution is not \em well \em formed the algorithm provided will find a
130
solution, x, which is not unique, but which sets the objective function
130
solution, x, which is not unique, but which sets the objective function
131
to 0. The reson being that the underlining algorithm works by SVD.
131
to 0. The reson being that the underlining algorithm works by SVD.
132
 
132
 
133
\param A the MxN matrix
133
\param A the MxN matrix
134
\param b the M vector
134
\param b the M vector
135
\return x will be resized if it is does not have the correct dimensions
135
\return x will be resized if it is does not have the correct dimensions
136
\exception assert(info==0) for Lapack. Add a throw statement later.
136
\exception assert(info==0) for Lapack. Add a throw statement later.
137
\exception assert(A.Rows()==b.Length());. Add a throw statement later.
137
\exception assert(A.Rows()==b.Length());. Add a throw statement later.
138
\version  Aug 2001 
138
\version  Aug 2001 
139
\author  Henrik Aanæs
139
\author  Henrik Aanæs
140
*/  
140
*/  
141
//@{
141
//@{
142
///Solves the Linear Least Squares problem min_x ||Ax=b||^2 for x.
142
///Solves the Linear Least Squares problem min_x ||Ax=b||^2 for x.
143
void LinearLSSolve(const CMatrix& A,const CVector&b,CVector& x);
143
void LinearLSSolve(const CMatrix& A,const CVector&b,CVector& x);
144
///Solves the Linear Least Squares problem min_x ||Ax=b||^2 for x, and returnes x.
144
///Solves the Linear Least Squares problem min_x ||Ax=b||^2 for x, and returnes x.
145
CVector LinearLSSolve(const CMatrix& A,const CVector&b);
145
CVector LinearLSSolve(const CMatrix& A,const CVector&b);
146
//@}
146
//@}
147
 
147
 
148
/**
148
/**
149
\name Matrix Inversion
149
\name Matrix Inversion
150
These functions inverts the square matrix A. This matrix A must have
150
These functions inverts the square matrix A. This matrix A must have
151
full rank. 
151
full rank. 
152
\param A square matrix
152
\param A square matrix
153
\return InvA the invers of A for one instance.
153
\return InvA the invers of A for one instance.
154
\exception assert(info==0) for Lapack. This wil among others happen if A is rank deficient. Add a throw statement later. 
154
\exception assert(info==0) for Lapack. This wil among others happen if A is rank deficient. Add a throw statement later. 
155
\exception assert(A.Rows()==A.Cols()). Add a throw statement later.
155
\exception assert(A.Rows()==A.Cols()). Add a throw statement later.
156
\version  Aug 2001 
156
\version  Aug 2001 
157
\author  Henrik Aanæs
157
\author  Henrik Aanæs
158
*/  
158
*/  
159
//@{
159
//@{
160
///Invertes the square matrix A. That is here A is altered as opposed to the other Invert functions.
160
///Invertes the square matrix A. That is here A is altered as opposed to the other Invert functions.
161
void Invert(CMatrix& A);
161
void Invert(CMatrix& A);
162
/// Returns the inverse of the square matrix A in InvA.
162
/// Returns the inverse of the square matrix A in InvA.
163
void Inverted(const CMatrix& A,CMatrix& InvA);
163
void Inverted(const CMatrix& A,CMatrix& InvA);
164
/// Returns the inverse of the square matrix A.
164
/// Returns the inverse of the square matrix A.
165
CMatrix Inverted(const CMatrix& A);
165
CMatrix Inverted(const CMatrix& A);
166
//@}
166
//@}
167
 
167
 
168
 
168
 
169
/**
169
/**
170
\name QR Factorization
170
\name QR Factorization
171
This function returns the QR factorization of A, such that Q*R=A where 
171
This function returns the QR factorization of A, such that Q*R=A where 
172
Q is a orthonormal matrix and R is an upper triangular matrix. However, 
172
Q is a orthonormal matrix and R is an upper triangular matrix. However, 
173
in the case of A.Col()>A.Row(), the last A.Col-A.Row columns of Q are 
173
in the case of A.Col()>A.Row(), the last A.Col-A.Row columns of Q are 
174
'carbage' and as such not part of a orthonormal matrix.
174
'carbage' and as such not part of a orthonormal matrix.
175
 \param A  the input matrix
175
 \param A  the input matrix
176
\return Q an orthonormal matrix. (See above)
176
\return Q an orthonormal matrix. (See above)
177
\return R an upper triangular matrix.
177
\return R an upper triangular matrix.
178
\exception assert(info==0) for Lapack. This wil among others happen if A is rank deficient. Add a throw statement later. 
178
\exception assert(info==0) for Lapack. This wil among others happen if A is rank deficient. Add a throw statement later. 
179
\exception assert(A.Rows()>0 && A.Cols()>0). Add a throw statement later.
179
\exception assert(A.Rows()>0 && A.Cols()>0). Add a throw statement later.
180
\version  Aug 2001 
180
\version  Aug 2001 
181
\author  Henrik Aanæs
181
\author  Henrik Aanæs
182
*/ 
182
*/ 
183
//@{ 
183
//@{ 
184
void QRfact(const CMatrix& A,CMatrix& Q, CMatrix& R);
184
void QRfact(const CMatrix& A,CMatrix& Q, CMatrix& R);
185
//@}
185
//@}
186
 
186
 
187
 
187
 
188
/**
188
/**
189
\name RQ Factorization
189
\name RQ Factorization
190
This function returns the RQ factorization of A, such that R*Q=A where 
190
This function returns the RQ factorization of A, such that R*Q=A where 
191
Q is a orthonormal matrix and R is an upper triangular matrix. However, 
191
Q is a orthonormal matrix and R is an upper triangular matrix. However, 
192
in the case of A not beeing a square matrix, there might be some fuck up of Q.
192
in the case of A not beeing a square matrix, there might be some fuck up of Q.
193
 \param A  the input matrix
193
 \param A  the input matrix
194
\return Q an orthonormal matrix. (See above)
194
\return Q an orthonormal matrix. (See above)
195
\return R an upper triangular matrix.
195
\return R an upper triangular matrix.
196
\exception assert(info==0) for Lapack. This wil among others happen if A is rank deficient. Add a throw statement later. 
196
\exception assert(info==0) for Lapack. This wil among others happen if A is rank deficient. Add a throw statement later. 
197
\exception assert(A.Rows()>0 && A.Cols()>0). Add a throw statement later.
197
\exception assert(A.Rows()>0 && A.Cols()>0). Add a throw statement later.
198
\version  Aug 2001 
198
\version  Aug 2001 
199
\author  Henrik Aanæs
199
\author  Henrik Aanæs
200
*/ 
200
*/ 
201
//@{ 
201
//@{ 
202
void RQfact(const CMatrix& A,CMatrix& R, CMatrix& Q);
202
void RQfact(const CMatrix& A,CMatrix& R, CMatrix& Q);
203
//@}
203
//@}
204
 
204
 
205
/**
205
/**
206
\name Find eigensolutions of a symmetric real matrix.
206
\name Find eigensolutions of a symmetric real matrix.
207
This function accepts a real symmetric matrix Q and a vector b.
207
This function accepts a real symmetric matrix Q and a vector b.
208
When the function returns, the eigenvalues of the matrix Q will be stored in b and
208
When the function returns, the eigenvalues of the matrix Q will be stored in b and
209
the eigenvectors form the columns of Q. This function is based on the
209
the eigenvectors form the columns of Q. This function is based on the
210
Lapack function dsyev, and returns its info code. A code of 0 indicates
210
Lapack function dsyev, and returns its info code. A code of 0 indicates
211
success, and a code < 0 indicates an error. Probably Q is not a real symmetric matrix.
211
success, and a code < 0 indicates an error. Probably Q is not a real symmetric matrix.
212
If the code is > 0 "the algorithm failed  to  converge;  code  off-diagonal  elements  
212
If the code is > 0 "the algorithm failed  to  converge;  code  off-diagonal  elements  
213
of an intermediate tridiagonal form did not converge to zero." Presumably this means that 
213
of an intermediate tridiagonal form did not converge to zero." Presumably this means that 
214
code contains the number of eigenvalues which are ok. 
214
code contains the number of eigenvalues which are ok. 
215
\author Andreas B¾rentzen.
215
\author Andreas B¾rentzen.
216
*/
216
*/
217
//@{
217
//@{
218
int EigenSolutionsSym(CMatrix& Q, CVector& b);
218
int EigenSolutionsSym(CMatrix& Q, CVector& b);
219
//@}
219
//@}
220
}
220
}
221
 
221
 
222
#endif // !defined(LAPACKFUNC_H_HAA_AGUST_2001)
222
#endif // !defined(LAPACKFUNC_H_HAA_AGUST_2001)
223
 
223