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

[Tim Rambo]

Thanks, that fixed it, but it has exposed another problem. this one is just for 
solving a linear system and it spits out a a matrix of nans. This code reads 
three pairs (ideal and rotated) of spherical coordinate vectors from a file 
converts them to cartesian, and calculates the matrix that rotates between the 
ideal coordinate system and the rotated coordinate system.
 
 
void solar_correct(){    
  
  FILE *data_file;
  double az1, az2, az3, el1, el2, el3, caz1, caz2, caz3, cel1, cel2, cel3;
  /*reading data from file*/
  data_file=fopen("vectors.dat", "r");
  if(data_file==NULL){
                 printf("unable to open file v.dat\n");
                 exit(1); 
                 }
   fscanf(data_file,"%lf %lf %lf %lf", &az1, &el1, &caz1, &cel1);
   
 
   fscanf(data_file,"%lf %lf %lf %lf", &az2, &el2, &caz2, &cel2);
   
 
   fscanf(data_file,"%lf %lf %lf %lf", &az3, &el3, &caz3, &cel3);
   
  fclose(data_file);
  
     /*creating cartesian vectors*/
     double V1=sin((double)M_PI/((double)180)*el1)*cos(M_PI/((double)180)*az1);
     double V2=sin((double)M_PI/((double)180)*el1)*sin(M_PI/((double)180)*az1);
     double V3=cos((double)M_PI/((double)180)*el1);
     
     double W1=sin((double)M_PI/((double)180)*el2)*cos(M_PI/((double)180)*az2);
     double W2=sin((double)M_PI/((double)180)*el2)*sin(M_PI/((double)180)*az2);
     double W3=cos((double)M_PI/((double)180)*el2); 
      
     double X1=sin((double)M_PI/((double)180)*el3)*cos(M_PI/((double)180)*az3);
     double X2=sin((double)M_PI/((double)180)*el3)*sin(M_PI*az3/((double)180));
     double X3=cos((double)M_PI/((double)180)*el3);
     
     double 
CV1=sin((double)M_PI/((double)180)*caz1)*cos(M_PI/((double)180)*cel1);
     double 
CV2=sin((double)M_PI/((double)180)*caz1)*sin(M_PI/((double)180)*cel1);
     double CV3=cos((double)M_PI/((double)180)*caz1);
       
     double 
CW1=sin((double)M_PI/((double)180)*caz2)*cos(M_PI/((double)180)*cel2);
     double 
CW2=sin((double)M_PI/((double)180)*caz2)*sin(M_PI/((double)180)*cel2);
     double CW3=cos((double)M_PI/((double)180)*caz2);
      
     double 
CX1=sin((double)M_PI/((double)180)*caz3)*cos(M_PI/((double)180)*cel3);
     double 
CX2=sin((double)M_PI/((double)180)*caz3)*sin(M_PI/((double)180)*cel3);
     double CX3=cos((double)M_PI/((double)180)*caz3);
      /*forming linear system*/
       double a_data[] = { V1, V2, V3,  0,  0,  0,  0,  0,  0,
                            0,  0,  0, V1, V2, V3,  0,  0,  0,
                            0,  0,  0,  0,  0,  0, V1, V2, V3,
                           W1, W2, W3,  0,  0,  0,  0,  0,  0,
                            0,  0,  0, W1, W2, W3,  0,  0,  0,
                            0,  0,  0,  0,  0,  0, W1, W2, W3,
                           X1, X2, X3,  0,  0,  0,  0,  0,  0,
                            0,  0,  0, X1, X2, X3,  0,  0,  0, 
                            0,  0,  0,  0,  0,  0, X1, X2, X3};
     
       double b_data[] = { CV1, CV2, CV3, CW1, CW2, CW3, CX1, CX2, CX3 };
     
       gsl_matrix_view m 
         = gsl_matrix_view_array (a_data, 9, 9);
     
       gsl_vector_view b
         = gsl_vector_view_array (b_data, 9);
     
       gsl_vector *x = gsl_vector_alloc (9);
      /*solving linear system*/  
       int s;
     
       gsl_permutation * p = gsl_permutation_alloc (9);
     
       gsl_linalg_LU_decomp (&m.matrix, p, &s);
     
       gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);
       gsl_vector_fprintf(stdout, x, "%g");
      /*writing rotation matrix into file for permanent storage*/
       data_file =fopen("matrix.dat", "w");
       if(data_file==NULL){
                           printf("file matrix.dat could not be accessed or 
written to\n");
                           exit(1);
                           }
       gsl_vector_fprintf (data_file, x, "%lf ");
       fclose(data_file);
     
       /*garbage collection*/       
       gsl_permutation_free (p);
       gsl_vector_free (x);
     
     }
 
 
 
Message: 5
Date: Tue, 6 Jan 2009 23:14:21 +0100
From: Steven Vancoillie <address@hidden>
Subject: Re: [Help-gsl] nan error in matrix-vector multiplication
To: address@hidden 
Message-ID: <address@hidden>
Content-Type: text/plain; charset="iso-8859-1"
 
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);
>           
>       }