[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
## Re: [Help-gsl] problem with gsl_blas_dtrsv

**From**: |
Martin Jansche |

**Subject**: |
Re: [Help-gsl] problem with gsl_blas_dtrsv |

**Date**: |
Sun, 9 Sep 2007 16:08:37 -0400 |

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*
>