Subversion Repositories gelsvn

Rev

Rev 39 | Rev 595 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

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