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

Subversion Repositories gelsvn

(root)/branches/cpp11-devel/src/LinAlg/LapackFunc.cpp – Rev 595 → 630