WebSVN – gelsvn – Diff – /branches/imesh-and-hmesh/src/LinAlg/LapackFunc.h


#if !defined(LAPACKFUNC_H_HAA_AGUST_2001)
#define LAPACKFUNC_H_HAA_AGUST_2001
 
#if defined(_MSC_VER)
#pragma message("Note: including lib: clapack.lib\n")
#pragma comment(lib, "clapack.lib")
 
#endif
 
#include "Matrix.h"
#include "Vector.h"
 
namespace LinAlg
{
 
/*!
\file LapackFunc.h
\brief Interface to some of the LAPACK functionality.
 
These are functions which more or less directly interface with the
Lapack provided algorithms.
 
For indepth reference to the LAPACK functions see:
  LAPACK Users' Guide - 3rd Edition,
  by E. Anderson et al.,
  ISBN 0-89871-447-8,
  Published by SIAM,
 
This book is also available at: \URL{http://www.netlib.org/lapack/lug/lapack_lug.html}
 
The official LAPACK sites where from the source can be downloaded are:
  \URL{http://www.netlib.org/clapack/} and
  \URL{http://www.netlib.org/lapack/}
 
 
NB: When running this in MS Visual C++ it is usually required to set the 
multithread "\MD" compiler option. This is to ensure correct linkage to the 
precompiled library "clapack.lib" and/or "clapackDB.lib".
 
 
\author  Henrik Aanæs
\version Aug 2001
*/
 
/*!
\name Singular Value Decomposition SVD 
 
These functions perform the Singular Value Decomposition SVD of 
the MxN matrix A. The SVD is defined by:
 
  A=U*S*V^T
 
where:
- U is a M by M orthogonal matrix 
- V is a N by N orthogonal matrix 
- S is a M by N diaggonal matrix. The values in the diagonal are the singular values
 
 
\param  A the matrix to perform SVD on  
\return U will be resized if it is does not have the correct dimensions
\return V will be resized if it is does not have the correct dimensions
\return S will be resized if it is does not have the correct dimensions. 
\exception assert(info==0) for Lapack. Add a throw statement later.
\version  Aug 2001 
\author  Henrik Aanæs
*/  
//@{ 
///SVD of A, where the singular values are returned in a Vector.
void SVD(const CMatrix& A,CMatrix& U,CVector& s,CMatrix& V);
///SVD of A, where the singular values are returned in a 'diagonal' Matrix.
void SVD(const CMatrix& A,CMatrix& U,CMatrix& S,CMatrix& V);
///SVD of A, returning only the singular values in a Vector.
CVector SVD(const CMatrix& A);
//@}
 
 
/*!
\name Linear Equations
These functions solve the system of linear equations
 
  A*x=b 
 
for x, where:
- A is a N by N matrix 
- b is a N vector
- x is a N vector
 
There a speceilaized functions for symetric positive definite (SPD) 
matrices yeilding better performance. These are denote by SPD in 
there function name.
 
\param A the NxN square matrix
\param b the N vector
\return x will be resized if it is does not have the correct dimensions
\exception assert(info==0) for Lapack. Add a throw statement later.
\exception assert(A.Row()==A.Col()). Add a throw statement later.
\exception assert(A.Row()==b.Length()). Add a throw statement later.
\version  Aug 2001 
\author  Henrik Aanæs
*/  
//@{
///Solves Ax=b for x.
void LinearSolve(const CMatrix& A,const CVector&b,CVector& x);
///Solves Ax=b for x and returns x.
CVector LinearSolve(const CMatrix& A,const CVector&b);
///Solves Ax=b for x, where A is SPD.
void LinearSolveSPD(const CMatrix& A,const CVector&b,CVector& x);
///Solves Ax=b for x and returns x, where A is SPD.
CVector LinearSolveSPD(const CMatrix& A,const CVector&b);
//@}
 
void LinearSolveSym(const CMatrix& A,
										const CVector&b,
										CVector& x);
 
/**
\name Linear Least Squares
These functions solve the Linear Least Squares problem:
 
  min_x ||Ax-b||^2
 
for x, where:
- || || denotes the 2-norm
- A	is a M by N matrix. For a well formed M>=N and rank (A)=N. See below.
- b	is a M vector.
- x	is a N vector 
 
If the solution is not \em well \em formed the algorithm provided will find a
solution, x, which is not unique, but which sets the objective function
to 0. The reson being that the underlining algorithm works by SVD.
 
\param A the MxN matrix
\param b the M vector
\return x will be resized if it is does not have the correct dimensions
\exception assert(info==0) for Lapack. Add a throw statement later.
\exception assert(A.Rows()==b.Length());. Add a throw statement later.
\version  Aug 2001 
\author  Henrik Aanæs
*/  
//@{
///Solves the Linear Least Squares problem min_x ||Ax=b||^2 for x.
void LinearLSSolve(const CMatrix& A,const CVector&b,CVector& x);
///Solves the Linear Least Squares problem min_x ||Ax=b||^2 for x, and returnes x.
CVector LinearLSSolve(const CMatrix& A,const CVector&b);
//@}
 
/**
\name Matrix Inversion
These functions inverts the square matrix A. This matrix A must have
full rank. 
\param A square matrix
\return InvA the invers of A for one instance.
\exception assert(info==0) for Lapack. This wil among others happen if A is rank deficient. Add a throw statement later. 
\exception assert(A.Rows()==A.Cols()). Add a throw statement later.
\version  Aug 2001 
\author  Henrik Aanæs
*/  
//@{
///Invertes the square matrix A. That is here A is altered as opposed to the other Invert functions.
void Invert(CMatrix& A);
/// Returns the inverse of the square matrix A in InvA.
void Inverted(const CMatrix& A,CMatrix& InvA);
/// Returns the inverse of the square matrix A.
CMatrix Inverted(const CMatrix& A);
//@}
 
 
/**
\name QR Factorization
This function returns the QR factorization of A, such that Q*R=A where 
Q is a orthonormal matrix and R is an upper triangular matrix. However, 
in the case of A.Col()>A.Row(), the last A.Col-A.Row columns of Q are 
'carbage' and as such not part of a orthonormal matrix.
 \param A  the input matrix
\return Q an orthonormal matrix. (See above)
\return R an upper triangular matrix.
\exception assert(info==0) for Lapack. This wil among others happen if A is rank deficient. Add a throw statement later. 
\exception assert(A.Rows()>0 && A.Cols()>0). Add a throw statement later.
\version  Aug 2001 
\author  Henrik Aanæs
*/ 
//@{ 
void QRfact(const CMatrix& A,CMatrix& Q, CMatrix& R);
//@}
 
 
/**
\name RQ Factorization
This function returns the RQ factorization of A, such that R*Q=A where 
Q is a orthonormal matrix and R is an upper triangular matrix. However, 
in the case of A not beeing a square matrix, there might be some fuck up of Q.
 \param A  the input matrix
\return Q an orthonormal matrix. (See above)
\return R an upper triangular matrix.
\exception assert(info==0) for Lapack. This wil among others happen if A is rank deficient. Add a throw statement later. 
\exception assert(A.Rows()>0 && A.Cols()>0). Add a throw statement later.
\version  Aug 2001 
\author  Henrik Aanæs
*/ 
//@{ 
void RQfact(const CMatrix& A,CMatrix& R, CMatrix& Q);
//@}
 
}
 
#endif // !defined(LAPACKFUNC_H_HAA_AGUST_2001)
 

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

Subversion Repositories gelsvn

(root)/branches/imesh-and-hmesh/src/LinAlg/LapackFunc.h – Rev 30 → 37