Subversion Repositories gelsvn

Rev

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

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