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

Subversion Repositories gelsvn

(root)/branches/wireframe-experiments/src/LinAlg/LapackFunc.cpp – Rev 379 → 393