Re: [Help-gsl] problem with gsl_blas_dtrsv

[Martin Jansche]

That's because you're using dTRsv -- the TR stands for "triangular".
You're only using the lower triangle of A, i.e. you're implicitly
inverting the triangular matrix

 {{2, 0, 0},
   {1, 2, 0},
   {1, 1, 1}}

whose inverse is

 {{1/2, 0, 0},
   {-1/4, 1/2, 0},
   {-1/4, -1/2, 1}}

When you multiply the latter with (1, 1, 1)', you get exactly the
solution you found using dtrsv().

If you want to solve a GEneral matrix problem, you have to use the
corresponding GE subroutine, e.g. dgesv (from LAPACK, not BLAS).  I
think the closest GSL analog of dgesv is gsl_linalg_LU_svx, and there
are other solvers too in the gsl_linalg package.

-- mj

On 9/8/07, Tomas Hudik <address@hidden> wrote:
> Hi all,
>
> I've been trying gsl_blas_dtrsv which solves equation: x<--A^{-1}x.
> (A is matrix and x is vector)
>
> Here is a part of my code:
>
>        double a[] = {2,1,1,
>                           1,2,1,
>                           1,1,1};
>
>        double x[]={1,1,1};
>
>
>        gsl_matrix_view A = gsl_matrix_view_array(m, 3, 3);
>        gsl_vector_view X = gsl_vector_view_array(x,3);
>
>        
> gsl_blas_dtrsv(CblasLower,CblasNoTrans,CblasNonUnit,&A.matrix,&X.vector);
>
>        cout << "\n\nResult: ";
>        for(int i=0;i<3;i++){
>              cout<<x[i]<<",";
>        }
> ---------------------------------------------------------------------------------------------------------------------------
>
> output is:
> Result: 0.5,0.25,0.25,
>
> But when I make it by hand, so:
> inversion matrix A^{-1} is:
> 1,0,-1
> 0,1,-1
> -1,-1,3
>
> so equation (matrix * vector) A^{-1}*x is:
> 1,0,-1
> 0,1,-1    *  (1,1,1) = (0,0,1) and NOT (0.5,0.25,0.25) as function dtrsv gave
> -1,-1,3
>
>
> Can you tell me why gsl_blas_dtrsv gives diffrent result?
>
> Thanks for your comments, Tomas
>
>
> _______________________________________________
> Help-gsl mailing list
> address@hidden
> http://lists.gnu.org/mailman/listinfo/help-gsl
>