[Help-gsl] stray #IND00

[le a]

Hi all,
   
  I am a newbie so please, any sort of help will be greatly appreciated.
   
  Here attached is a program which gives me a #IND00 value for the Q matrix. 
What is disturbing is the fact that if I comment the code related with the 
inversion of the R matrix (up in the code) the previous #IND00 values disapear 
(and Q gets the corect ones).
  I am using GNUWin32 with DevC++, gsl 1.8, sorry, I know this is not the 
proper list, but the sourceforge list is a little sluggish. I am trying to 
build a XBasic wrapper which will be used on Linux machines also.
   
  Am I doing something wrong?
   
  Thanks,
  Lee
   
       #include <stdio.h>
     #include <gsl/gsl_matrix.h>
     #include <gsl/gsl_permutation.h>
     
     int
     main (void)
     {
unsigned int i, j, k, jj, kk; 
double  Stack[1][5];
double E11, E22, G12, Niu12, Teta, Niu21, m, n, Thk;
double Ajk, Bjk, Djk;
  
double Pidouble = 3.1415926535897932384626433832795;
   gsl_permutation * perm = gsl_permutation_alloc (3);
 gsl_matrix * dat = gsl_matrix_alloc (4, 6);
   gsl_matrix * R = gsl_matrix_alloc (3, 3);
 gsl_matrix * Rinv = gsl_matrix_alloc (3, 3);
 gsl_matrix_set_zero(R);
 gsl_matrix_set(R, 0, 0, 1);
 gsl_matrix_set(R, 1, 1, 1);
 gsl_matrix_set(R, 2, 2, 2);
   gsl_matrix * Temp1 = gsl_matrix_alloc(3, 3);
 gsl_permutation * perm1 = gsl_permutation_alloc (3); 
  //' ************************************************************************
//' **Resolve Inverse Matrix R
//' ************************************************************************
 gsl_matrix * Temp = gsl_matrix_alloc(3, 3);
 gsl_matrix_memcpy(Temp, R);
 int s = 0;
 // Get determinant of matrix R
  gsl_linalg_LU_det(Temp, s);
 // Make LU decomposition of matrix R
  gsl_linalg_LU_decomp (Temp, perm, &s);
 // Invert the matrix R
  gsl_linalg_LU_invert (Temp, perm, Rinv);
   gsl_matrix * S = gsl_matrix_alloc (3, 3);
 gsl_matrix_set_zero(S);
 gsl_matrix * Q = gsl_matrix_alloc (3, 3);
 gsl_matrix * T = gsl_matrix_alloc (3, 3);
   FILE * f = fopen("laminate.lam", "r");
 // Fill the matrix dat
 int ret = gsl_matrix_fscanf(f, dat);
 fclose (f);
//' ************************************************************************
//' **Resolve Stack THK Matrix STACK
//' ************************************************************************
 Thk = 0;
 Stack[0][0] = Thk;
 for (i = 0; i < 4; i++)
     {
  Thk = Thk + gsl_matrix_get(dat, i, 5);
  Stack[0][i+1] = Thk;
 }
    // START
for (i = 0; i < 4; i++)
{
        printf("############\n");
  E11 = gsl_matrix_get(dat, i, 0);
  E22 = gsl_matrix_get(dat, i, 1);
  G12 = gsl_matrix_get(dat, i, 2);
  Niu12 = gsl_matrix_get(dat, i, 3);
  Teta = gsl_matrix_get(dat, i, 4);
    Niu21 = Niu12 * E22 / E11;
  m = cos (Pidouble * Teta / 180);
  n = sin (Pidouble * Teta / 180);
  //' ************************************************************************
//' **Resolve Matrix S
//' ************************************************************************
  gsl_matrix_set(S, 0, 0, 1/E11);
  gsl_matrix_set(S, 0, 1, -Niu21/E22);
  gsl_matrix_set(S, 1, 0, -Niu12/E11);
  gsl_matrix_set(S, 1, 1, 1/E22);
  gsl_matrix_set(S, 2, 2, 1/G12);
      j=0;
    for (k=0; k<3; k++)
         printf("%e   %e   %e\n", gsl_matrix_get(S, k, j), gsl_matrix_get(S, k, 
j+1), gsl_matrix_get(S, k, j+2));
         printf("############\n");
//' ************************************************************************
//' **Resolve Matrix Q
//' ************************************************************************
  gsl_matrix_memcpy(Temp1, S);
  int p = 0;
 // Get determinant of matrix S
  double get = gsl_linalg_LU_det(Temp1, p);
        printf("get = %e\n", get);
 // Make LU decomposition of matrix S
  gsl_linalg_LU_decomp (Temp1, perm1, &p);
 // Invert the matrix S
  gsl_linalg_LU_invert (Temp1, perm1, Q);
      j=0;
    for (k=0; k<3; k++)
         printf("%e   %e   %e\n", gsl_matrix_get(Q, k, j), gsl_matrix_get(Q, k, 
j+1), gsl_matrix_get(Q, k, j+2));
    getchar();
}
gsl_matrix_free (dat);
gsl_matrix_free (R);
gsl_matrix_free (Rinv);
gsl_matrix_free (Temp);
gsl_permutation_free(perm);
  gsl_matrix_free (Temp1);
gsl_permutation_free (perm1);
  
gsl_matrix_free (S);
gsl_matrix_free (Q);
gsl_matrix_free (T);
  return 0;
}

   
  laminate.lam
   
  24545e6     5861e6     22960e6  0.3  0      0.00051
24545e6     5861e6     22960e6   0.3     45     0.00051
24545e6     5861e6     22960e6   0.3     -45    0.00051
24545e6     5861e6     22960e6   0.3     90     0.00051
   

 Send instant messages to your online friends http://au.messenger.yahoo.com