[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Help-gsl] LAPACK,CBLAS matrix inversion
From: |
James Bergstra |
Subject: |
Re: [Help-gsl] LAPACK,CBLAS matrix inversion |
Date: |
Mon, 7 Aug 2006 18:59:53 +0200 |
User-agent: |
Mutt/1.5.11+cvs20060403 |
On Fri, Aug 04, 2006 at 02:06:39PM +0100, Wei Cheng wrote:
> Hi,
>
> The cholesky decomposition and LU decomposition in GSL linear algerbra
> section, is only suitable for small matrix, for large matrix LAPACK is
> better as stated on GSL website. I am wondering if there is anyone who has
> experience in using cholesky and LU decomposion in LAPACK?? My matrix is 500
> by 500. or where to find relavent information? Someone suggested to use
> CBLAS. What is the difference between LAPACK and CBLAS? Which is better??
> Regards
My previous plan of making some sort of guide is way beyond my expertise... so
scratch that. If you like, here's a bit of code that uses dpotrf, which is
[obviously!] the way you compute a cholesky factor with LAPACK. It is also
good to know that ATLAS implements a few LAPACK functions, and dpotrf is one of
the lucky ones.
gsl_matrix * covarL = gsl_matrix_alloc(_N,_N);
gsl_matrix * covarL2 = gsl_matrix_alloc(_N,_N);
for (size_t i = 0; i < _N; ++i)
{
for (size_t j = 0; j <= i; ++j)
{
double cov = gsl_stats_covariance_m( GMAT_COL(mindata,i),
GMAT_COL(mindata,j),
mindata->size1, GMAT_AT(_mu,0,i), GMAT_AT(_mu,0,j));
gsl_matrix_set(covarL,i,j,cov);
gsl_matrix_set(covarL,j,i,cov);
}
}
//save a backup in case covariance is not full-rank
gsl_matrix_memcpy(covarL2, covarL);
int null_dim = clapack_dpotrf(CblasRowMajor, CblasLower, _N, covarL->data,
_N);
if (null_dim < 0)
{
assert(!"illegal argument to dpotrf...");
}
else if (null_dim > 0)
{
fprintf(stderr, "WARNING: covariance is low on rank (%i of %zu),
normalizing with lambda diagonal.\n", null_dim-1, _N);
for (size_t i = 0; i < _N; ++i)
{
mAT(covarL2, i, i) += lambda;
}
gsl_matrix_memcpy(covarL, covarL2);
null_dim = clapack_dpotrf(CblasRowMajor, CblasLower, _N, covarL->data,
_N);
if (null_dim > 0) assert(!"even our hack could not make cholesky
work:(");
}
--
http://www-etud.iro.umontreal.ca/~bergstrj