help-gsl
[Top][All Lists]

## Re: [Help-gsl] nan error in matrix-vector multiplication

 From: Steven Vancoillie Subject: Re: [Help-gsl] nan error in matrix-vector multiplication Date: Tue, 6 Jan 2009 23:14:21 +0100 User-agent: KMail/1.9.9

```On Tuesday 06 January 2009 22:58:08 Tim Rambo wrote:
> The problem is this: given the two angles from a unit vector in spherical
> coordinates, I create a cartesian vector which I then multiply by the
> identity matrix. The first two elements of the vector are correct but the
> final element of the vector is nan. Any help is appreciated.
>
> Thanks
>
>
>
> g_rotate(double *az, double *el){
>       FILE *file;
>
>       double d[9];
>       /*creating cartesian vector*/
>       double v[]=
> {sin((*el)*M_PI/180.0)*cos((*az)*M_PI/180.0),sin((*el)*M_PI/180.0)*sin((*az)*M_PI/180.0),
>  cos((*el)*M_PI/180.0)};
>       printf("%lf    %lf     %lf\n", v[0], v[1], v[2]);
>       file = fopen("matrix.dat", "r");
>       if(file==NULL){
>                      return -1;
>                      }
>

The following line looks quite horrible
>       fscanf(file, "%lf%lf%lf%lf%lf%lf%lf%lf%lf", &d[0], &d[1], &d[2], &d[3],
> &d[4], &d[5], &d[6], &d[7], &d[8]);
I think this is more readable:
for (i=0; i<9; ++i) fscanf (file, "%lf", d+i);

>   fclose(file);
>
>        printf("%lf    %lf     %lf\n %lf    %lf     %lf\n %lf    %lf     %lf
> ", d[0], d[1], d[2], d[3], d[4], d[5], d[6], d[7], d[8]);
>        /*creating linear system*/
>       gsl_matrix_view R = gsl_matrix_view_array(d, 3,3);
>
>       gsl_vector_view V = gsl_vector_view_array(v,3);
>
>
>       gsl_vector *CV= gsl_vector_alloc(3);
>      /*multiplying matrix and vector*/

looks like you are using CV uninitialized and adding R*V to it in this line:
>       gsl_blas_dgemv(CblasNoTrans, 1.0, &R.matrix, &V.vector, 1.0, CV);
so I think you should use instead:
gsl_blas_dgemv(CblasNoTrans, 1.0, &R.matrix, &V.vector, 0.0, CV);

>       printf("corrected vector in rectangular coordinates\n");
>       gsl_vector_fprintf(stdout, CV, "%g");
>        /*decomposing vecor into angles and storing the information*/
>        double xsq=gsl_vector_get(CV, 0)*gsl_vector_get(CV, 0);
>        double ysq=gsl_vector_get(CV, 1)*gsl_vector_get(CV, 1);
>        double zsq=gsl_vector_get(CV, 2)*gsl_vector_get(CV, 2);
>        *az=(float)(acos( gsl_vector_get(CV,
> 2)/sqrt(xsq+ysq+zsq))*180.00/M_PI);
>        *el=(float)(atan(gsl_vector_get(CV, 1)/gsl_vector_get(CV,
> 0))*180.00/M_PI);
>
>         /*garbage collection*/
>        gsl_vector_free(CV);
>
>       }
> _______________________________________________
> Help-gsl mailing list
> http://lists.gnu.org/mailman/listinfo/help-gsl
>

grtz
Steven
```

signature.asc
Description: This is a digitally signed message part.