Re: [Help-gsl] LAPACK,CBLAS matrix inversion

[James Bergstra]

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