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

[Steven Vancoillie]

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]);   
>         /*reading matrix from file*/
>       file = fopen("matrix.dat", "r");
>       if(file==NULL){
>                      printf("unable to read file\n");
>                      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
> address@hidden
> http://lists.gnu.org/mailman/listinfo/help-gsl
> 

grtz
Steven

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