[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), 
                    mindata->size1, GMAT_AT(_mu,0,i), GMAT_AT(_mu,0,j));

    //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, 

    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, 

        if (null_dim > 0) assert(!"even our hack could not make cholesky 


reply via email to

[Prev in Thread] Current Thread [Next in Thread]