WebSVN – gelsvn – Diff – /trunk/src/LinAlg/LapackFunc.cpp


/* ----------------------------------------------------------------------- *
 * This file is part of GEL, http://www.imm.dtu.dk/GEL
 * Copyright (C) the authors and DTU Informatics
 * For license and list of authors, see ../../doc/intro.pdf
 * ----------------------------------------------------------------------- */
 
#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 [] BX;
	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);
    
    delete [] work;
    
	return info;
}
 
}
 

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

Subversion Repositories gelsvn

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