help-gsl
[Top][All Lists]

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

 From: Tim Rambo Subject: [Help-gsl] re:Re: nan error in matrix-vector multiplication Date: Thu, 08 Jan 2009 13:57:55 -0500

```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;
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
Subject: Re: [Help-gsl] nan error in matrix-vector multiplication
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;
>       /*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, v, v);
>       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, &d, &d, &d,
> &d, &d, &d, &d, &d);
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, d, d, d, d, d, d, d, d);
>        /*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);
>
>       }

```