Subversion Repositories gelsvn

Rev

Rev 379 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 379 Rev 595
Line -... Line 1...
-
 
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" {
Line 69... Line 75...
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,
Line 112... Line 118...
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
 
Line 192... Line 198...
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)
Line 234... Line 240...
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
{
Line 314... Line 321...
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
Line 363... Line 370...
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
Line 393... Line 400...
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
}