[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*
>* *
>* http://www.gnu.org/software/gsl/manual/html_node/Linear-Algebra-Examples.html*
>* *
>* 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
calculations.
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
accurate.
--
JDL