[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Help-gsl] Best way to Calculating inv(A) . B

From: John D Lamb
Subject: Re: [Help-gsl] Best way to Calculating inv(A) . B
Date: Fri, 02 Apr 2010 10:01:06 +0100

On Fri, 2010-04-02 at 14:27 +1100, Srimal Jayawardena wrote: 
> One way I can think of is to use,
>  gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);
> as explained in
> I would need to find the solution vector x for each col vector b (of
> matrix B) and and concatenate the x vectors to obtain my inv(A).B
> matrix.
> Is there a better / more direct approach that is more numerically
> stable and accurate ?

Two possibilities occur to me. Both use gsl_linalg_LU_refine(). This
uses some form of iterative improvement (perhaps Jacobi — I haven’t
checked) and computes a vector of residuals r = Ax – b, which should be
zero if x solves Ax = b. Something like the following ought to give more
numerically accurate solutions.

double d;
gsl_vector* r;

for( d = 10; d > TOLERANCE; ){
  gsl_linalg_LU_refine( A, LU, p, b, x, r );
  d = gsl_blas_dnrm2( r );

This should improve the value of x. It would make sense to check that d
is actually decreasing because (i) iterative improvements are not
guaranteed to converge (though they typically do), (ii) the improvement
will be limited by the numerical accuracy of double-precision

The first possibility is to improve each estimate x of a column of
inv(A).B in turn. The second is to estimate inv(A) directly and use
iterative improvement to improve the estimates of the columns of inv(A).
I think I would use the first method.

This is not more direct, but it should typically be more numerically


reply via email to

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