[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