Subversion Repositories gelsvn

Rev

Rev 39 | Rev 379 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 39 Rev 223
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\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
 
8
 
8
#endif
9
#endif
9
 
10
 
10
#include "Matrix.h"
11
#include "Matrix.h"
11
#include "Vector.h"
12
#include "Vector.h"
12
 
13
 
13
namespace LinAlg
14
namespace LinAlg
14
{
15
{
15
 
16
 
16
/*!
17
/*!
17
\file LapackFunc.h
18
\file LapackFunc.h
18
\brief Interface to some of the LAPACK functionality.
19
\brief Interface to some of the LAPACK functionality.
19
 
20
 
20
These are functions which more or less directly interface with the
21
These are functions which more or less directly interface with the
21
Lapack provided algorithms.
22
Lapack provided algorithms.
22
 
23
 
23
For indepth reference to the LAPACK functions see:
24
For indepth reference to the LAPACK functions see:
24
  LAPACK Users' Guide - 3rd Edition,
25
  LAPACK Users' Guide - 3rd Edition,
25
  by E. Anderson et al.,
26
  by E. Anderson et al.,
26
  ISBN 0-89871-447-8,
27
  ISBN 0-89871-447-8,
27
  Published by SIAM,
28
  Published by SIAM,
28
 
29
 
29
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}
30
 
31
 
31
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:
32
  \URL{http://www.netlib.org/clapack/} and
33
  \URL{http://www.netlib.org/clapack/} and
33
  \URL{http://www.netlib.org/lapack/}
34
  \URL{http://www.netlib.org/lapack/}
34
 
35
 
35
 
36
 
36
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 
37
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 
38
precompiled library "clapack.lib" and/or "clapackDB.lib".
39
precompiled library "clapack.lib" and/or "clapackDB.lib".
39
 
40
 
40
 
41
 
41
\author  Henrik Aanæs
42
\author  Henrik Aanæs
42
\version Aug 2001
43
\version Aug 2001
43
*/
44
*/
44
 
45
 
45
/*!
46
/*!
46
\name Singular Value Decomposition SVD 
47
\name Singular Value Decomposition SVD 
47
 
48
 
48
These functions perform the Singular Value Decomposition SVD of 
49
These functions perform the Singular Value Decomposition SVD of 
49
the MxN matrix A. The SVD is defined by:
50
the MxN matrix A. The SVD is defined by:
50
 
51
 
51
  A=U*S*V^T
52
  A=U*S*V^T
52
 
53
 
53
where:
54
where:
54
- U is a M by M orthogonal matrix 
55
- U is a M by M orthogonal matrix 
55
- V is a N by N orthogonal matrix 
56
- V is a N by N orthogonal matrix 
56
- 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
57
 
58
 
58
 
59
 
59
\param  A the matrix to perform SVD on  
60
\param  A the matrix to perform SVD on  
60
\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
61
\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
62
\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. 
63
\exception assert(info==0) for Lapack. Add a throw statement later.
64
\exception assert(info==0) for Lapack. Add a throw statement later.
64
\version  Aug 2001 
65
\version  Aug 2001 
65
\author  Henrik Aanæs
66
\author  Henrik Aanæs
66
*/  
67
*/  
67
//@{ 
68
//@{ 
68
///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.
69
void SVD(const CMatrix& A,CMatrix& U,CVector& s,CMatrix& V);
70
void SVD(const CMatrix& A,CMatrix& U,CVector& s,CMatrix& V);
70
///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.
71
void SVD(const CMatrix& A,CMatrix& U,CMatrix& S,CMatrix& V);
72
void SVD(const CMatrix& A,CMatrix& U,CMatrix& S,CMatrix& V);
72
///SVD of A, returning only the singular values in a Vector.
73
///SVD of A, returning only the singular values in a Vector.
73
CVector SVD(const CMatrix& A);
74
CVector SVD(const CMatrix& A);
74
//@}
75
//@}
75
 
76
 
76
 
77
 
77
/*!
78
/*!
78
\name Linear Equations
79
\name Linear Equations
79
These functions solve the system of linear equations
80
These functions solve the system of linear equations
80
 
81
 
81
  A*x=b 
82
  A*x=b 
82
 
83
 
83
for x, where:
84
for x, where:
84
- A is a N by N matrix 
85
- A is a N by N matrix 
85
- b is a N vector
86
- b is a N vector
86
- x is a N vector
87
- x is a N vector
87
 
88
 
88
There a speceilaized functions for symetric positive definite (SPD) 
89
There a speceilaized functions for symetric positive definite (SPD) 
89
matrices yeilding better performance. These are denote by SPD in 
90
matrices yeilding better performance. These are denote by SPD in 
90
there function name.
91
there function name.
91
 
92
 
92
\param A the NxN square matrix
93
\param A the NxN square matrix
93
\param b the N vector
94
\param b the N vector
94
\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
95
\exception assert(info==0) for Lapack. Add a throw statement later.
96
\exception assert(info==0) for Lapack. Add a throw statement later.
96
\exception assert(A.Row()==A.Col()). Add a throw statement later.
97
\exception assert(A.Row()==A.Col()). Add a throw statement later.
97
\exception assert(A.Row()==b.Length()). Add a throw statement later.
98
\exception assert(A.Row()==b.Length()). Add a throw statement later.
98
\version  Aug 2001 
99
\version  Aug 2001 
99
\author  Henrik Aanæs
100
\author  Henrik Aanæs
100
*/  
101
*/  
101
//@{
102
//@{
102
///Solves Ax=b for x.
103
///Solves Ax=b for x.
103
void LinearSolve(const CMatrix& A,const CVector&b,CVector& x);
104
void LinearSolve(const CMatrix& A,const CVector&b,CVector& x);
104
///Solves Ax=b for x and returns x.
105
///Solves Ax=b for x and returns x.
105
CVector LinearSolve(const CMatrix& A,const CVector&b);
106
CVector LinearSolve(const CMatrix& A,const CVector&b);
106
///Solves Ax=b for x, where A is SPD.
107
///Solves Ax=b for x, where A is SPD.
107
void LinearSolveSPD(const CMatrix& A,const CVector&b,CVector& x);
108
void LinearSolveSPD(const CMatrix& A,const CVector&b,CVector& x);
108
///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.
109
CVector LinearSolveSPD(const CMatrix& A,const CVector&b);
110
CVector LinearSolveSPD(const CMatrix& A,const CVector&b);
110
//@}
111
//@}
111
 
112
 
112
void LinearSolveSym(const CMatrix& A,
113
void LinearSolveSym(const CMatrix& A,
113
										const CVector&b,
114
										const CVector&b,
114
										CVector& x);
115
										CVector& x);
115
 
116
 
116
/**
117
/**
117
\name Linear Least Squares
118
\name Linear Least Squares
118
These functions solve the Linear Least Squares problem:
119
These functions solve the Linear Least Squares problem:
119
 
120
 
120
  min_x ||Ax-b||^2
121
  min_x ||Ax-b||^2
121
 
122
 
122
for x, where:
123
for x, where:
123
- || || denotes the 2-norm
124
- || || denotes the 2-norm
124
- 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.
125
- b	is a M vector.
126
- b	is a M vector.
126
- x	is a N vector 
127
- x	is a N vector 
127
 
128
 
128
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
129
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
130
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.
131
 
132
 
132
\param A the MxN matrix
133
\param A the MxN matrix
133
\param b the M vector
134
\param b the M vector
134
\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
135
\exception assert(info==0) for Lapack. Add a throw statement later.
136
\exception assert(info==0) for Lapack. Add a throw statement later.
136
\exception assert(A.Rows()==b.Length());. Add a throw statement later.
137
\exception assert(A.Rows()==b.Length());. Add a throw statement later.
137
\version  Aug 2001 
138
\version  Aug 2001 
138
\author  Henrik Aanæs
139
\author  Henrik Aanæs
139
*/  
140
*/  
140
//@{
141
//@{
141
///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.
142
void LinearLSSolve(const CMatrix& A,const CVector&b,CVector& x);
143
void LinearLSSolve(const CMatrix& A,const CVector&b,CVector& x);
143
///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.
144
CVector LinearLSSolve(const CMatrix& A,const CVector&b);
145
CVector LinearLSSolve(const CMatrix& A,const CVector&b);
145
//@}
146
//@}
146
 
147
 
147
/**
148
/**
148
\name Matrix Inversion
149
\name Matrix Inversion
149
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
150
full rank. 
151
full rank. 
151
\param A square matrix
152
\param A square matrix
152
\return InvA the invers of A for one instance.
153
\return InvA the invers of A for one instance.
153
\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. 
154
\exception assert(A.Rows()==A.Cols()). Add a throw statement later.
155
\exception assert(A.Rows()==A.Cols()). Add a throw statement later.
155
\version  Aug 2001 
156
\version  Aug 2001 
156
\author  Henrik Aanæs
157
\author  Henrik Aanæs
157
*/  
158
*/  
158
//@{
159
//@{
159
///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.
160
void Invert(CMatrix& A);
161
void Invert(CMatrix& A);
161
/// Returns the inverse of the square matrix A in InvA.
162
/// Returns the inverse of the square matrix A in InvA.
162
void Inverted(const CMatrix& A,CMatrix& InvA);
163
void Inverted(const CMatrix& A,CMatrix& InvA);
163
/// Returns the inverse of the square matrix A.
164
/// Returns the inverse of the square matrix A.
164
CMatrix Inverted(const CMatrix& A);
165
CMatrix Inverted(const CMatrix& A);
165
//@}
166
//@}
166
 
167
 
167
 
168
 
168
/**
169
/**
169
\name QR Factorization
170
\name QR Factorization
170
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 
171
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, 
172
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 
173
'carbage' and as such not part of a orthonormal matrix.
174
'carbage' and as such not part of a orthonormal matrix.
174
 \param A  the input matrix
175
 \param A  the input matrix
175
\return Q an orthonormal matrix. (See above)
176
\return Q an orthonormal matrix. (See above)
176
\return R an upper triangular matrix.
177
\return R an upper triangular matrix.
177
\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. 
178
\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.
179
\version  Aug 2001 
180
\version  Aug 2001 
180
\author  Henrik Aanæs
181
\author  Henrik Aanæs
181
*/ 
182
*/ 
182
//@{ 
183
//@{ 
183
void QRfact(const CMatrix& A,CMatrix& Q, CMatrix& R);
184
void QRfact(const CMatrix& A,CMatrix& Q, CMatrix& R);
184
//@}
185
//@}
185
 
186
 
186
 
187
 
187
/**
188
/**
188
\name RQ Factorization
189
\name RQ Factorization
189
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 
190
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, 
191
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.
192
 \param A  the input matrix
193
 \param A  the input matrix
193
\return Q an orthonormal matrix. (See above)
194
\return Q an orthonormal matrix. (See above)
194
\return R an upper triangular matrix.
195
\return R an upper triangular matrix.
195
\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. 
196
\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.
197
\version  Aug 2001 
198
\version  Aug 2001 
198
\author  Henrik Aanæs
199
\author  Henrik Aanæs
199
*/ 
200
*/ 
200
//@{ 
201
//@{ 
201
void RQfact(const CMatrix& A,CMatrix& R, CMatrix& Q);
202
void RQfact(const CMatrix& A,CMatrix& R, CMatrix& Q);
202
//@}
203
//@}
203
 
204
 
204
}
205
}
205
 
206
 
206
#endif // !defined(LAPACKFUNC_H_HAA_AGUST_2001)
207
#endif // !defined(LAPACKFUNC_H_HAA_AGUST_2001)
207
 
208