WebSVN – gelsvn – Diff – /branches/vc2005beta-branch/src/LinAlg/LapackFunc.cpp


#include <iostream>
#include <vector>
#include "LapackFunc.h"
 
extern "C" {
	extern void dgetrf_(const int *m, const int *n, double *a, const int *lda, int *ipiv, int *info);
	extern void dgetri_(const int *n, double *a, const int *lda, int *ipiv,double *work, int *lwork, int *info);
	extern int dgesvd_(const char *jobu, const char *jobvt, const int *m, const int *n, double *a, const int *lda, double *s,
		double *u, const int *ldu, double *vt, const int *ldvt, double *work, const int *lwork, int *info);
	extern void dgelss_(const int *, const int *, const int *, double *, const int *, double *, const int *, double *, double *, int *, double *, const int *, int *);
	extern void dgesv_(const int *N, const int *NRHS, double *A, const int *LDA, int *IPIV, double *B, const int *LDB, int *INFO);
	extern void dposv_(const char *UPLO, const int *N, const int *NRHS, double *A, const int *LDA, double *B, const int *LDB, int *INFO);
	extern void dgeqrf_(const int *m, const int *n, double* A,const int* lda,double*tau,double *work, const int* lwork, int* info);
	extern void dorgqr_(const int *m, const int *n, const int *k, double* A, const int *lda,const double *tau,double* work, const int *lwork, int *info);
	extern void dgerqf_(const int *m, const int *n, double* A,const int* lda,double*tau,double *work, const int* lwork, int* info);
	extern void dorgrq_(const int *m, const int *n, const int *k, double* A, const int *lda,const double *tau,double* work, const int *lwork, int *info);
 
	int dsysv_(char *uplo, int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, double *work, int *lwork, int *info);
}
 
namespace{
	
	template <class T>
		inline T MIN(const T& a,const T& b)
	{
		return a<b?a:b;
	}
	
	template <class T>
		inline T MAX(const T& a,const T& b)
	{
		return a>b?a:b;
	}
	
}
 
namespace LinAlg
{
 
void SVD(const CMatrix& A,
		 CMatrix& U,
		 CVector& s,
		 CMatrix& V)
{
	//Notice, that the Lapack/Fortran CMatrix representation is the transposed version of the C/C++
	CMatrix a;
	A.Transposed(a);
	
	int nRows=A.Rows();
	int nCols=A.Cols();
	int info;
	int lda = MAX(nRows,1);
	int lwork= 6*MIN(nRows,nCols) + nRows + nCols;
	
	s.Resize(MIN(nRows,nCols));
	double *work = new double[lwork];
	
	char jobu='A';
	int ldu=nRows;
	U.Resize(nRows,nRows);
	
	char jobvt='A';
	int ldvt=nCols;
	V.Resize(nCols,nCols);
	
	
	dgesvd_ (&jobu, &jobvt, &nRows, &nCols, a[0], &lda, &(s[0]), U[0], &ldu, V[0], &ldvt,work, &lwork, &info);
	
	assert(info==0);
	
	delete work;
	U.Transpose();
}
 
void SVD(const CMatrix& A,
		 CMatrix& U,
		 CMatrix& S,
		 CMatrix& V)
{
	CVector s;
	SVD(A,U,s,V);
	S.Resize(A.Rows(),A.Cols());
	S=0;
	for(int cS=MIN(A.Rows(),A.Cols())-1;cS>=0;cS--)
	{
		S[cS][cS]=s[cS];
	}
}
 
CVector SVD(const CMatrix& A)
{
	CMatrix a;
	A.Transposed(a);
	
	int nRows=A.Rows();
	int nCols=A.Cols();
	int info;
	int lda = MAX(nRows,1);
	int lwork= 6*MIN(nRows,nCols) + nRows + nCols;
	
	CVector s(MIN(nRows,nCols));
	double *work = new double[lwork];
	
	char jobu='N';
	int ldu=1;
	
	char jobvt='N';
	int ldvt=1;
	
	dgesvd_ (&jobu, &jobvt, &nRows, &nCols, a[0], &lda, &(s[0]), NULL, &ldu, NULL, &ldvt,work, &lwork, &info);
	
	assert(info==0);
	
	delete work;
	
	return s;
}
 
 
void LinearSolve(const CMatrix& A,
				 const CVector&b,
				 CVector& x)
{
	assert(A.Rows()==b.Length());
	assert(A.Rows()==A.Cols());
	
	CMatrix a;
	A.Transposed(a);
	x=b;
	int nRows=A.Rows();
	int nrhs=1;	//only one CVector, change here to make b a CMatrix.
	int info;
	int *ipiv =new int[nRows];
	
	dgesv_(&nRows, &nrhs, a[0], &nRows, ipiv, &(x[0]), &nRows, &info);
	
	assert(info==0);
	delete [] ipiv;
}
 
CVector LinearSolve(const CMatrix& A,
				   const CVector&b)
{
	CVector x;
	LinearSolve(A,b,x);
	return x;
}
 
 
void LinearSolveSPD(const CMatrix& A,
					const CVector&b,
					CVector& x)
{
	assert(A.Rows()==b.Length());
	assert(A.Rows()==A.Cols());
	
	CMatrix a(A);
	x=b;
	
	char uplo='U';
	int nRows=A.Rows();
	int info;
	int nrhs=1;
	
	dposv_(&uplo, &nRows, &nrhs, a[0], &nRows, &(x[0]), &nRows, &info);
	
	assert(info==0);
}
 
void LinearSolveSym(const CMatrix& A,
										const CVector&b,
										CVector& x)
{
	assert(A.Rows()==b.Length());
	assert(A.Rows()==A.Cols());
	
	CMatrix a(A);
	x=b;
	
	char uplo='U';
	int nRows=A.Rows();
	int nCols=A.Cols();
	int info;
	int nrhs=1;
	
	int lwork= 10*6*MIN(nRows,nCols) + nRows + nCols;
	double *work = new double[lwork];
 
	std::vector<int> ipiv(A.Rows());
	dsysv_(&uplo, &nRows, &nrhs, a[0], &nRows, &ipiv[0], &(x[0]), &nRows, 
				 work, &lwork, &info);
	delete work;
	assert(info==0);
}
 
 
CVector LinearSolveSPD(const CMatrix& A,const CVector&b)
{
	CVector x;
	LinearSolveSPD(A,b,x);
	return x;
}
 
void LinearLSSolve(const CMatrix& A,
				   const CVector& b,
				   CVector& x)
{
	assert(A.Rows()==b.Length());
	
	int nRows=A.Rows();
	int nCols=A.Cols();
	CMatrix a;
	A.Transposed(a);
	
	int ldb=MAX(nRows,nCols);
	double* BX=new double[ldb];
	memcpy(BX,&(b[0]),nRows*sizeof(double));
	
	int nrhs=1;
	double* s=new double[MIN(nRows,nCols)];
	double rcond = -1.0; // using machine precision
	int info,rank;
	int lwork = 5*nRows*nCols + 1; // larger than necessary
	double *work = new double[lwork];
	
	dgelss_(&nRows, &nCols, &nrhs, a[0], &nRows, BX, &ldb, s, &rcond, &rank, work, &lwork, &info);
	
	x.Resize(nCols);
	memcpy(&(x[0]),BX,nCols*sizeof(double));
	
	
	assert(info==0);
	
	delete s;
	delete work;
}
 
CVector LinearLSSolve(const CMatrix& A,
					 const CVector& b)
{
	CVector x;
	LinearLSSolve(A,b,x);
	return x;
}
 
void Invert(CMatrix& A)
{
	assert(A.Rows()==A.Cols());
	int nRows=A.Rows();
	int info;
	int *ipiv = new int[nRows];
	//Perform th LU factorization of A
	dgetrf_ (&nRows, &nRows, A[0], &nRows, ipiv, &info);
	if (info != 0) { //An error occured
		delete [] ipiv;
		assert(info==0);	//info will be < 0 if A is rank deficient.
		return;
	}
 
	//Calculate the inverse
	int lwork = nRows * 5; 
	double *work = new double[lwork];
	dgetri_ (&nRows, A[0], &nRows, ipiv, work, &lwork, &info);
	delete [] ipiv;
	delete [] work;
	assert(info==0);
}
 
void Inverted(const CMatrix& A,
			  CMatrix& InvA)
{
	InvA=A;
	Invert(InvA);
}
 
CMatrix Inverted(const CMatrix& A)
{
	CMatrix InvA(A);
	Invert(InvA);
	return InvA;
}
 
void QRfact(const CMatrix& A,CMatrix& Q, CMatrix& R)
{
 
	
	int nRows=A.Rows();
	int nCols=A.Cols();
	assert(A.Rows()>0 && A.Cols()>0); //herby lda=nRows.
	int lwork= 6*nCols;
	int info;
	int k=MIN(nRows,nCols);
 
 
	A.Transposed(R);
 
	double *work = new double[lwork];
 
 
	CVector TAU(k);
 
	dgeqrf_(&nRows,&nCols,R[0],&nRows,&TAU[0],work,&lwork,&info);
	assert(info==0);
 
 
	Q=R;
 
	nCols=MIN(nCols,nRows);
 
	dorgqr_(&nRows,&nCols,&k,Q[0],&nRows,&TAU[0],work,&lwork,&info);
	assert(info==0);
 
 
	delete work;
 
	R.Transpose();
	Q.Transpose();
 
	//Set the zeros of R
	double * pRow;
	int ColStop;
	for(int cRow=1;cRow<R.Rows();cRow++)
	{
		pRow=R[cRow];
		ColStop=MIN(cRow,R.Cols());
		for(int cCol=0;cCol<ColStop;cCol++)
		{
			pRow[cCol]=0;
		}
	}
 
}
 
 
void RQfact(const CMatrix& A,CMatrix& R, CMatrix& Q)
{
	int nRows=A.Rows();
	int nCols=A.Cols();
	assert(A.Rows()>0 && A.Cols()>0); //herby lda=nRows.
	int lwork= 6*nCols;
	int info;
	int k=MIN(nRows,nCols);
 
 
	A.Transposed(R);
 
	double *work = new double[lwork];
 
 
	CVector TAU(k);
 
	dgerqf_(&nRows,&nCols,R[0],&nRows,&TAU[0],work,&lwork,&info);
	assert(info==0);
 
	Q=R;
 
	nCols=MAX(nCols,nRows);
 
	dorgrq_(&nRows,&nCols,&k,Q[0],&nRows,&TAU[0],work,&lwork,&info);
	assert(info==0);
 
	
	delete work;
 
	R.Transpose();
	Q.Transpose();
 
	//Set the zeros of R
	double * pRow;
	int ColStop;
	for(int cRow=1;cRow<R.Rows();cRow++)
	{
		pRow=R[cRow];
		ColStop=MIN(cRow,R.Cols());
		for(int cCol=0;cCol<ColStop;cCol++)
		{
			pRow[cCol]=0;
		}
	}
 
}
 
 
}
 

Rev 39	Rev 119
1	`#include <iostream>`	1	`#include <iostream>`
2	`#include <vector>`	2	`#include <vector>`
3	`#include "LapackFunc.h"`	3	`#include "LapackFunc.h"`
4		4
5	`extern "C" {`	5	`extern "C" {`
6	`extern void dgetrf_(const int m, const int n, double a, const int lda, int ipiv, int info);`	6	`extern void dgetrf_(const int m, const int n, double a, const int lda, int ipiv, int info);`
7	`extern void dgetri_(const int n, double a, const int lda, int ipiv,double work, int lwork, int *info);`	7	`extern void dgetri_(const int n, double a, const int lda, int ipiv,double work, int lwork, int *info);`
8	`extern int dgesvd_(const char jobu, const char jobvt, const int m, const int n, double a, const int lda, double *s,`	8	`extern int dgesvd_(const char jobu, const char jobvt, const int m, const int n, double a, const int lda, double *s,`
9	`double u, const int ldu, double vt, const int ldvt, double work, const int lwork, int *info);`	9	`double u, const int ldu, double vt, const int ldvt, double work, const int lwork, int *info);`
10	`extern void dgelss_(const int , const int , const int , double , const int , double , const int , double , double , int , double , const int , int *);`	10	`extern void dgelss_(const int , const int , const int , double , const int , double , const int , double , double , int , double , const int , int *);`
11	`extern void dgesv_(const int N, const int NRHS, double A, const int LDA, int IPIV, double B, const int LDB, int INFO);`	11	`extern void dgesv_(const int N, const int NRHS, double A, const int LDA, int IPIV, double B, const int LDB, int INFO);`
12	`extern void dposv_(const char UPLO, const int N, const int NRHS, double A, const int LDA, double B, const int LDB, int INFO);`	12	`extern void dposv_(const char UPLO, const int N, const int NRHS, double A, const int LDA, double B, const int LDB, int INFO);`
13	`extern void dgeqrf_(const int m, const int n, double* A,const int* lda,doubletau,double work, const int* lwork, int* info);`	13	`extern void dgeqrf_(const int m, const int n, double* A,const int* lda,doubletau,double work, const int* lwork, int* info);`
14	`extern void dorgqr_(const int m, const int n, const int k, double A, const int lda,const double tau,double* work, const int lwork, int info);`	14	`extern void dorgqr_(const int m, const int n, const int k, double A, const int lda,const double tau,double* work, const int lwork, int info);`
15	`extern void dgerqf_(const int m, const int n, double* A,const int* lda,doubletau,double work, const int* lwork, int* info);`	15	`extern void dgerqf_(const int m, const int n, double* A,const int* lda,doubletau,double work, const int* lwork, int* info);`
16	`extern void dorgrq_(const int m, const int n, const int k, double A, const int lda,const double tau,double* work, const int lwork, int info);`	16	`extern void dorgrq_(const int m, const int n, const int k, double A, const int lda,const double tau,double* work, const int lwork, int info);`
17		17
18	`int dsysv_(char uplo, int n, int nrhs, double a, int lda, int ipiv, double b, int ldb, double work, int lwork, int *info);`	18	`int dsysv_(char uplo, int n, int nrhs, double a, int lda, int ipiv, double b, int ldb, double work, int lwork, int *info);`
19	`}`	19	`}`
20		20
21	`namespace{`	21	`namespace{`
22		22
23	`template <class T>`	23	`template <class T>`
24	`inline T MIN(const T& a,const T& b)`	24	`inline T MIN(const T& a,const T& b)`
25	`{`	25	`{`
26	`return a<b?a:b;`	26	`return a<b?a:b;`
27	`}`	27	`}`
28		28
29	`template <class T>`	29	`template <class T>`
30	`inline T MAX(const T& a,const T& b)`	30	`inline T MAX(const T& a,const T& b)`
31	`{`	31	`{`
32	`return a>b?a:b;`	32	`return a>b?a:b;`
33	`}`	33	`}`
34		34
35	`}`	35	`}`
36		36
37	`namespace LinAlg`	37	`namespace LinAlg`
38	`{`	38	`{`
39		39
40	`void SVD(const CMatrix& A,`	40	`void SVD(const CMatrix& A,`
41	`CMatrix& U,`	41	`CMatrix& U,`
42	`CVector& s,`	42	`CVector& s,`
43	`CMatrix& V)`	43	`CMatrix& V)`
44	`{`	44	`{`
45	`//Notice, that the Lapack/Fortran CMatrix representation is the transposed version of the C/C++`	45	`//Notice, that the Lapack/Fortran CMatrix representation is the transposed version of the C/C++`
46	`CMatrix a;`	46	`CMatrix a;`
47	`A.Transposed(a);`	47	`A.Transposed(a);`
48		48
49	`int nRows=A.Rows();`	49	`int nRows=A.Rows();`
50	`int nCols=A.Cols();`	50	`int nCols=A.Cols();`
51	`int info;`	51	`int info;`
52	`int lda = MAX(nRows,1);`	52	`int lda = MAX(nRows,1);`
53	`int lwork= 6*MIN(nRows,nCols) + nRows + nCols;`	53	`int lwork= 6*MIN(nRows,nCols) + nRows + nCols;`
54		54
55	`s.Resize(MIN(nRows,nCols));`	55	`s.Resize(MIN(nRows,nCols));`
56	`double *work = new double[lwork];`	56	`double *work = new double[lwork];`
57		57
58	`char jobu='A';`	58	`char jobu='A';`
59	`int ldu=nRows;`	59	`int ldu=nRows;`
60	`U.Resize(nRows,nRows);`	60	`U.Resize(nRows,nRows);`
61		61
62	`char jobvt='A';`	62	`char jobvt='A';`
63	`int ldvt=nCols;`	63	`int ldvt=nCols;`
64	`V.Resize(nCols,nCols);`	64	`V.Resize(nCols,nCols);`
65		65
66		66
67	`dgesvd_ (&jobu, &jobvt, &nRows, &nCols, a[0], &lda, &(s[0]), U[0], &ldu, V[0], &ldvt,work, &lwork, &info);`	67	`dgesvd_ (&jobu, &jobvt, &nRows, &nCols, a[0], &lda, &(s[0]), U[0], &ldu, V[0], &ldvt,work, &lwork, &info);`
68		68
69	`assert(info==0);`	69	`assert(info==0);`
70		70
71	`delete work;`	71	`delete work;`
72	`U.Transpose();`	72	`U.Transpose();`
73	`}`	73	`}`
74		74
75	`void SVD(const CMatrix& A,`	75	`void SVD(const CMatrix& A,`
76	`CMatrix& U,`	76	`CMatrix& U,`
77	`CMatrix& S,`	77	`CMatrix& S,`
78	`CMatrix& V)`	78	`CMatrix& V)`
79	`{`	79	`{`
80	`CVector s;`	80	`CVector s;`
81	`SVD(A,U,s,V);`	81	`SVD(A,U,s,V);`
82	`S.Resize(A.Rows(),A.Cols());`	82	`S.Resize(A.Rows(),A.Cols());`
83	`S=0;`	83	`S=0;`
84	`for(int cS=MIN(A.Rows(),A.Cols())-1;cS>=0;cS--)`	84	`for(int cS=MIN(A.Rows(),A.Cols())-1;cS>=0;cS--)`
85	`{`	85	`{`
86	`S[cS][cS]=s[cS];`	86	`S[cS][cS]=s[cS];`
87	`}`	87	`}`
88	`}`	88	`}`
89		89
90	`CVector SVD(const CMatrix& A)`	90	`CVector SVD(const CMatrix& A)`
91	`{`	91	`{`
92	`CMatrix a;`	92	`CMatrix a;`
93	`A.Transposed(a);`	93	`A.Transposed(a);`
94		94
95	`int nRows=A.Rows();`	95	`int nRows=A.Rows();`
96	`int nCols=A.Cols();`	96	`int nCols=A.Cols();`
97	`int info;`	97	`int info;`
98	`int lda = MAX(nRows,1);`	98	`int lda = MAX(nRows,1);`
99	`int lwork= 6*MIN(nRows,nCols) + nRows + nCols;`	99	`int lwork= 6*MIN(nRows,nCols) + nRows + nCols;`
100		100
101	`CVector s(MIN(nRows,nCols));`	101	`CVector s(MIN(nRows,nCols));`
102	`double *work = new double[lwork];`	102	`double *work = new double[lwork];`
103		103
104	`char jobu='N';`	104	`char jobu='N';`
105	`int ldu=1;`	105	`int ldu=1;`
106		106
107	`char jobvt='N';`	107	`char jobvt='N';`
108	`int ldvt=1;`	108	`int ldvt=1;`
109		109
110	`dgesvd_ (&jobu, &jobvt, &nRows, &nCols, a[0], &lda, &(s[0]), NULL, &ldu, NULL, &ldvt,work, &lwork, &info);`	110	`dgesvd_ (&jobu, &jobvt, &nRows, &nCols, a[0], &lda, &(s[0]), NULL, &ldu, NULL, &ldvt,work, &lwork, &info);`
111		111
112	`assert(info==0);`	112	`assert(info==0);`
113		113
114	`delete work;`	114	`delete work;`
115		115
116	`return s;`	116	`return s;`
117	`}`	117	`}`
118		118
119		119
120	`void LinearSolve(const CMatrix& A,`	120	`void LinearSolve(const CMatrix& A,`
121	`const CVector&b,`	121	`const CVector&b,`
122	`CVector& x)`	122	`CVector& x)`
123	`{`	123	`{`
124	`assert(A.Rows()==b.Length());`	124	`assert(A.Rows()==b.Length());`
125	`assert(A.Rows()==A.Cols());`	125	`assert(A.Rows()==A.Cols());`
126		126
127	`CMatrix a;`	127	`CMatrix a;`
128	`A.Transposed(a);`	128	`A.Transposed(a);`
129	`x=b;`	129	`x=b;`
130	`int nRows=A.Rows();`	130	`int nRows=A.Rows();`
131	`int nrhs=1; //only one CVector, change here to make b a CMatrix.`	131	`int nrhs=1; //only one CVector, change here to make b a CMatrix.`
132	`int info;`	132	`int info;`
133	`int *ipiv =new int[nRows];`	133	`int *ipiv =new int[nRows];`
134		134
135	`dgesv_(&nRows, &nrhs, a[0], &nRows, ipiv, &(x[0]), &nRows, &info);`	135	`dgesv_(&nRows, &nrhs, a[0], &nRows, ipiv, &(x[0]), &nRows, &info);`
136		136
137	`assert(info==0);`	137	`assert(info==0);`
138	`delete [] ipiv;`	138	`delete [] ipiv;`
139	`}`	139	`}`
140		140
141	`CVector LinearSolve(const CMatrix& A,`	141	`CVector LinearSolve(const CMatrix& A,`
142	`const CVector&b)`	142	`const CVector&b)`
143	`{`	143	`{`
144	`CVector x;`	144	`CVector x;`
145	`LinearSolve(A,b,x);`	145	`LinearSolve(A,b,x);`
146	`return x;`	146	`return x;`
147	`}`	147	`}`
148		148
149		149
150	`void LinearSolveSPD(const CMatrix& A,`	150	`void LinearSolveSPD(const CMatrix& A,`
151	`const CVector&b,`	151	`const CVector&b,`
152	`CVector& x)`	152	`CVector& x)`
153	`{`	153	`{`
154	`assert(A.Rows()==b.Length());`	154	`assert(A.Rows()==b.Length());`
155	`assert(A.Rows()==A.Cols());`	155	`assert(A.Rows()==A.Cols());`
156		156
157	`CMatrix a(A);`	157	`CMatrix a(A);`
158	`x=b;`	158	`x=b;`
159		159
160	`char uplo='U';`	160	`char uplo='U';`
161	`int nRows=A.Rows();`	161	`int nRows=A.Rows();`
162	`int info;`	162	`int info;`
163	`int nrhs=1;`	163	`int nrhs=1;`
164		164
165	`dposv_(&uplo, &nRows, &nrhs, a[0], &nRows, &(x[0]), &nRows, &info);`	165	`dposv_(&uplo, &nRows, &nrhs, a[0], &nRows, &(x[0]), &nRows, &info);`
166		166
167	`assert(info==0);`	167	`assert(info==0);`
168	`}`	168	`}`
169		169
170	`void LinearSolveSym(const CMatrix& A,`	170	`void LinearSolveSym(const CMatrix& A,`
171	`const CVector&b,`	171	`const CVector&b,`
172	`CVector& x)`	172	`CVector& x)`
173	`{`	173	`{`
174	`assert(A.Rows()==b.Length());`	174	`assert(A.Rows()==b.Length());`
175	`assert(A.Rows()==A.Cols());`	175	`assert(A.Rows()==A.Cols());`
176		176
177	`CMatrix a(A);`	177	`CMatrix a(A);`
178	`x=b;`	178	`x=b;`
179		179
180	`char uplo='U';`	180	`char uplo='U';`
181	`int nRows=A.Rows();`	181	`int nRows=A.Rows();`
182	`int nCols=A.Cols();`	182	`int nCols=A.Cols();`
183	`int info;`	183	`int info;`
184	`int nrhs=1;`	184	`int nrhs=1;`
185		185
186	`int lwork= 106MIN(nRows,nCols) + nRows + nCols;`	186	`int lwork= 106MIN(nRows,nCols) + nRows + nCols;`
187	`double *work = new double[lwork];`	187	`double *work = new double[lwork];`
188		188
189	`std::vector<int> ipiv(A.Rows());`	189	`std::vector<int> ipiv(A.Rows());`
190	`dsysv_(&uplo, &nRows, &nrhs, a[0], &nRows, &ipiv[0], &(x[0]), &nRows,`	190	`dsysv_(&uplo, &nRows, &nrhs, a[0], &nRows, &ipiv[0], &(x[0]), &nRows,`
191	`work, &lwork, &info);`	191	`work, &lwork, &info);`
192	`delete work;`	192	`delete work;`
193	`assert(info==0);`	193	`assert(info==0);`
194	`}`	194	`}`
195		195
196		196
197	`CVector LinearSolveSPD(const CMatrix& A,const CVector&b)`	197	`CVector LinearSolveSPD(const CMatrix& A,const CVector&b)`
198	`{`	198	`{`
199	`CVector x;`	199	`CVector x;`
200	`LinearSolveSPD(A,b,x);`	200	`LinearSolveSPD(A,b,x);`
201	`return x;`	201	`return x;`
202	`}`	202	`}`
203		203
204	`void LinearLSSolve(const CMatrix& A,`	204	`void LinearLSSolve(const CMatrix& A,`
205	`const CVector& b,`	205	`const CVector& b,`
206	`CVector& x)`	206	`CVector& x)`
207	`{`	207	`{`
208	`assert(A.Rows()==b.Length());`	208	`assert(A.Rows()==b.Length());`
209		209
210	`int nRows=A.Rows();`	210	`int nRows=A.Rows();`
211	`int nCols=A.Cols();`	211	`int nCols=A.Cols();`
212	`CMatrix a;`	212	`CMatrix a;`
213	`A.Transposed(a);`	213	`A.Transposed(a);`
214		214
215	`int ldb=MAX(nRows,nCols);`	215	`int ldb=MAX(nRows,nCols);`
216	`double* BX=new double[ldb];`	216	`double* BX=new double[ldb];`
217	`memcpy(BX,&(b[0]),nRows*sizeof(double));`	217	`memcpy(BX,&(b[0]),nRows*sizeof(double));`
218		218
219	`int nrhs=1;`	219	`int nrhs=1;`
220	`double* s=new double[MIN(nRows,nCols)];`	220	`double* s=new double[MIN(nRows,nCols)];`
221	`double rcond = -1.0; // using machine precision`	221	`double rcond = -1.0; // using machine precision`
222	`int info,rank;`	222	`int info,rank;`
223	`int lwork = 5nRowsnCols + 1; // larger than necessary`	223	`int lwork = 5nRowsnCols + 1; // larger than necessary`
224	`double *work = new double[lwork];`	224	`double *work = new double[lwork];`
225		225
226	`dgelss_(&nRows, &nCols, &nrhs, a[0], &nRows, BX, &ldb, s, &rcond, &rank, work, &lwork, &info);`	226	`dgelss_(&nRows, &nCols, &nrhs, a[0], &nRows, BX, &ldb, s, &rcond, &rank, work, &lwork, &info);`
227		227
228	`x.Resize(nCols);`	228	`x.Resize(nCols);`
229	`memcpy(&(x[0]),BX,nCols*sizeof(double));`	229	`memcpy(&(x[0]),BX,nCols*sizeof(double));`
230		230
231		231
232	`assert(info==0);`	232	`assert(info==0);`
233		233
234	`delete s;`	234	`delete s;`
235	`delete work;`	235	`delete work;`
236	`}`	236	`}`
237		237
238	`CVector LinearLSSolve(const CMatrix& A,`	238	`CVector LinearLSSolve(const CMatrix& A,`
239	`const CVector& b)`	239	`const CVector& b)`
240	`{`	240	`{`
241	`CVector x;`	241	`CVector x;`
242	`LinearLSSolve(A,b,x);`	242	`LinearLSSolve(A,b,x);`
243	`return x;`	243	`return x;`
244	`}`	244	`}`
245		245
246	`void Invert(CMatrix& A)`	246	`void Invert(CMatrix& A)`
247	`{`	247	`{`
248	`assert(A.Rows()==A.Cols());`	248	`assert(A.Rows()==A.Cols());`
249	`int nRows=A.Rows();`	249	`int nRows=A.Rows();`
250	`int info;`	250	`int info;`
251	`int *ipiv = new int[nRows];`	251	`int *ipiv = new int[nRows];`
252	`//Perform th LU factorization of A`	252	`//Perform th LU factorization of A`
253	`dgetrf_ (&nRows, &nRows, A[0], &nRows, ipiv, &info);`	253	`dgetrf_ (&nRows, &nRows, A[0], &nRows, ipiv, &info);`
254	`if (info != 0) { //An error occured`	254	`if (info != 0) { //An error occured`
255	`delete [] ipiv;`	255	`delete [] ipiv;`
256	`assert(info==0); //info will be < 0 if A is rank deficient.`	256	`assert(info==0); //info will be < 0 if A is rank deficient.`
257	`return;`	257	`return;`
258	`}`	258	`}`
259		259
260	`//Calculate the inverse`	260	`//Calculate the inverse`
261	`int lwork = nRows * 5;`	261	`int lwork = nRows * 5;`
262	`double *work = new double[lwork];`	262	`double *work = new double[lwork];`
263	`dgetri_ (&nRows, A[0], &nRows, ipiv, work, &lwork, &info);`	263	`dgetri_ (&nRows, A[0], &nRows, ipiv, work, &lwork, &info);`
264	`delete [] ipiv;`	264	`delete [] ipiv;`
265	`delete [] work;`	265	`delete [] work;`
266	`assert(info==0);`	266	`assert(info==0);`
267	`}`	267	`}`
268		268
269	`void Inverted(const CMatrix& A,`	269	`void Inverted(const CMatrix& A,`
270	`CMatrix& InvA)`	270	`CMatrix& InvA)`
271	`{`	271	`{`
272	`InvA=A;`	272	`InvA=A;`
273	`Invert(InvA);`	273	`Invert(InvA);`
274	`}`	274	`}`
275		275
276	`CMatrix Inverted(const CMatrix& A)`	276	`CMatrix Inverted(const CMatrix& A)`
277	`{`	277	`{`
278	`CMatrix InvA(A);`	278	`CMatrix InvA(A);`
279	`Invert(InvA);`	279	`Invert(InvA);`
280	`return InvA;`	280	`return InvA;`
281	`}`	281	`}`
282		282
283	`void QRfact(const CMatrix& A,CMatrix& Q, CMatrix& R)`	283	`void QRfact(const CMatrix& A,CMatrix& Q, CMatrix& R)`
284	`{`	284	`{`
285		285
286		286
287	`int nRows=A.Rows();`	287	`int nRows=A.Rows();`
288	`int nCols=A.Cols();`	288	`int nCols=A.Cols();`
289	`assert(A.Rows()>0 && A.Cols()>0); //herby lda=nRows.`	289	`assert(A.Rows()>0 && A.Cols()>0); //herby lda=nRows.`
290	`int lwork= 6*nCols;`	290	`int lwork= 6*nCols;`
291	`int info;`	291	`int info;`
292	`int k=MIN(nRows,nCols);`	292	`int k=MIN(nRows,nCols);`
293		293
294		294
295	`A.Transposed(R);`	295	`A.Transposed(R);`
296		296
297	`double *work = new double[lwork];`	297	`double *work = new double[lwork];`
298		298
299		299
300	`CVector TAU(k);`	300	`CVector TAU(k);`
301		301
302	`dgeqrf_(&nRows,&nCols,R[0],&nRows,&TAU[0],work,&lwork,&info);`	302	`dgeqrf_(&nRows,&nCols,R[0],&nRows,&TAU[0],work,&lwork,&info);`
303	`assert(info==0);`	303	`assert(info==0);`
304		304
305		305
306	`Q=R;`	306	`Q=R;`
307		307
308	`nCols=MIN(nCols,nRows);`	308	`nCols=MIN(nCols,nRows);`
309		309
310	`dorgqr_(&nRows,&nCols,&k,Q[0],&nRows,&TAU[0],work,&lwork,&info);`	310	`dorgqr_(&nRows,&nCols,&k,Q[0],&nRows,&TAU[0],work,&lwork,&info);`
311	`assert(info==0);`	311	`assert(info==0);`
312		312
313		313
314	`delete work;`	314	`delete work;`
315		315
316	`R.Transpose();`	316	`R.Transpose();`
317	`Q.Transpose();`	317	`Q.Transpose();`
318		318
319	`//Set the zeros of R`	319	`//Set the zeros of R`
320	`double * pRow;`	320	`double * pRow;`
321	`int ColStop;`	321	`int ColStop;`
322	`for(int cRow=1;cRow<R.Rows();cRow++)`	322	`for(int cRow=1;cRow<R.Rows();cRow++)`
323	`{`	323	`{`
324	`pRow=R[cRow];`	324	`pRow=R[cRow];`
325	`ColStop=MIN(cRow,R.Cols());`	325	`ColStop=MIN(cRow,R.Cols());`
326	`for(int cCol=0;cCol<ColStop;cCol++)`	326	`for(int cCol=0;cCol<ColStop;cCol++)`
327	`{`	327	`{`
328	`pRow[cCol]=0;`	328	`pRow[cCol]=0;`
329	`}`	329	`}`
330	`}`	330	`}`
331		331
332	`}`	332	`}`
333		333
334		334
335	`void RQfact(const CMatrix& A,CMatrix& R, CMatrix& Q)`	335	`void RQfact(const CMatrix& A,CMatrix& R, CMatrix& Q)`
336	`{`	336	`{`
337	`int nRows=A.Rows();`	337	`int nRows=A.Rows();`
338	`int nCols=A.Cols();`	338	`int nCols=A.Cols();`
339	`assert(A.Rows()>0 && A.Cols()>0); //herby lda=nRows.`	339	`assert(A.Rows()>0 && A.Cols()>0); //herby lda=nRows.`
340	`int lwork= 6*nCols;`	340	`int lwork= 6*nCols;`
341	`int info;`	341	`int info;`
342	`int k=MIN(nRows,nCols);`	342	`int k=MIN(nRows,nCols);`
343		343
344		344
345	`A.Transposed(R);`	345	`A.Transposed(R);`
346		346
347	`double *work = new double[lwork];`	347	`double *work = new double[lwork];`
348		348
349		349
350	`CVector TAU(k);`	350	`CVector TAU(k);`
351		351
352	`dgerqf_(&nRows,&nCols,R[0],&nRows,&TAU[0],work,&lwork,&info);`	352	`dgerqf_(&nRows,&nCols,R[0],&nRows,&TAU[0],work,&lwork,&info);`
353	`assert(info==0);`	353	`assert(info==0);`
354		354
355	`Q=R;`	355	`Q=R;`
356		356
357	`nCols=MAX(nCols,nRows);`	357	`nCols=MAX(nCols,nRows);`
358		358
359	`dorgrq_(&nRows,&nCols,&k,Q[0],&nRows,&TAU[0],work,&lwork,&info);`	359	`dorgrq_(&nRows,&nCols,&k,Q[0],&nRows,&TAU[0],work,&lwork,&info);`
360	`assert(info==0);`	360	`assert(info==0);`
361		361
362		362
363	`delete work;`	363	`delete work;`
364		364
365	`R.Transpose();`	365	`R.Transpose();`
366	`Q.Transpose();`	366	`Q.Transpose();`
367		367
368	`//Set the zeros of R`	368	`//Set the zeros of R`
369	`double * pRow;`	369	`double * pRow;`
370	`int ColStop;`	370	`int ColStop;`
371	`for(int cRow=1;cRow<R.Rows();cRow++)`	371	`for(int cRow=1;cRow<R.Rows();cRow++)`
372	`{`	372	`{`
373	`pRow=R[cRow];`	373	`pRow=R[cRow];`
374	`ColStop=MIN(cRow,R.Cols());`	374	`ColStop=MIN(cRow,R.Cols());`
375	`for(int cCol=0;cCol<ColStop;cCol++)`	375	`for(int cCol=0;cCol<ColStop;cCol++)`
376	`{`	376	`{`
377	`pRow[cCol]=0;`	377	`pRow[cCol]=0;`
378	`}`	378	`}`
379	`}`	379	`}`
380		380
381	`}`	381	`}`
382		382
383		383
384	`}`	384	`}`
385		385

Subversion Repositories gelsvn

(root)/branches/vc2005beta-branch/src/LinAlg/LapackFunc.cpp – Rev 39 → 119