WebSVN – gelsvn – Diff – /trunk/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);
	
	int dsytrd_(char *uplo, int *n, double *a, int * lda, double *d__, double *e, double *tau, double *work, int *lwork, int *info);
	int dsteqr_(char *compz, int *n, double *d__, double *e, double *z__, int *ldz, double *work, int *info);
    int dsyev_(char *jobz, char *uplo, int *n, double *a, int *lda, double *w, 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;
		}
	}
 
}
 
int EigenSolutionsSym(CMatrix& Q, CVector& b)
{
	assert(Q.Rows() == Q.Cols());
	int n = Q.Rows();
	char jobz='V', uplo='U';
	int info, lwork=6*n;
	b.Resize(n);
	double* work = new double[lwork];
	
	dsyev_(&jobz, &uplo, &n, Q[0], &n, &b[0], work, &lwork, &info);
	return info;
}
 
}
 

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

Subversion Repositories gelsvn

(root)/trunk/src/LinAlg/LapackFunc.cpp – Rev 39 → 379