[Help-gsl] Runtime error, Segmentation fault.

[Jo_Jo]

Hi help-gsl,

Please run these program. Runtime error, Segmentation fault. Why?

Regards
Andrzej Buczynski         Fri, 11/09/07


#include <stdio.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_math.h>

int main (void)
{
// first attempt with these data set    
/*
  double a_data[] = {-0.18, -0.60,  0.57, -0.96,
                      0.41, -0.24,  0.99,  0.58,
                     -0.14, -0.30, -0.97, -0.66,
                      0.51,  0.13,  0.19,  0.85 };
*/                       
 // second attempt with these data set
/*  double a_data[] = {-1.18, -9.60,  0.57, -0.96,
                      0.61, -0.24,  0.99,  8.58,
                     -0.14, -7.30, -0.97, -0.66,
                      0.51,  0.13,  7.19,  9.85 };
*/
 
// third attempt with these data set

 double a_data[] = {- 6.18, -9.60,  0.57, -0.01,
                      9.61, -0.24, -6.99,  8.58,
                     -3.14, -0.01, -0.97, -0.66,
                      0.51,  0.13, -7.19, -9.85 };
 
 
       
       double temp;
       int M =4, N =4, i, j ;
       int * signum ;
       gsl_matrix_view A = gsl_matrix_view_array (a_data, M, N);
       gsl_vector *norm = gsl_vector_alloc (M);
       gsl_vector *tau  = gsl_vector_alloc (GSL_MIN(M,N) );
       gsl_permutation * p = gsl_permutation_alloc (N);
       gsl_matrix *QRPT = gsl_matrix_alloc ( M, N);

printf("\n  ----------------------------------\n ");


// initialization with zeroes       
       gsl_vector_set_zero ( tau ) ;
       gsl_vector_set_zero ( norm ) ;
        
// copy A -> QRPT   
    
   for (i = 0; i < M; i++)
       for (j = 0; j < N; j++){
        temp = gsl_matrix_get (&A.matrix, i, j);
        gsl_matrix_set(QRPT, i, j, temp) ; }
 
printf("  ----------------------------------\n ");  
       
       
/* int gsl_linalg_QRPT_decomp (gsl matrix * A, 
                               gsl vector * tau,          
                               gsl permutation * p,
                               int * signum, 
                               gsl vector * norm )
This function factorizes the M -by-N matrix A into the QRPT
decomposition A = QRPT. On output the diagonal and upper triangular 
part of the input matrix contain the matrix R. 

The permutation matrix P is stored in the permutation p. 

The sign of the permutation is given by signum. It has the value (−1)n,
where n is the number of interchanges in the permutation. 

The vector tau and the columns of the lower triangular part of the
matrix A contain the Householder coefficients and vectors which
encode the orthogonal matrix Q. 
 
The vector tau must be of length k = min(M, N ).
                                                                                
                  T
The matrix Q is related to these components by, Q = Qk ...Q2 Q1 
where Qi = I −τi * vi *viT and vi is the Householder vector 
vi = (0, ..., 1, A(i + 1, i), A(i + 2, i), ..., A(m, i)). This
is the same storage scheme as used by lapack. 

The vector norm is a workspace of length N used for column pivoting.

The algorithm used to perform the decomposition is Householder QR 
with column pivoting 
(Golub & Van Loan, Matrix Computations, Algorithm 5.4.1).
*/      
    
 
printf("\n The matrix Q is related to the components by\n");
printf("\n  Q = Qk ...Q2 Q1 where Qi = I − taui* vi* viT\n");
printf("\n  and vi is the Householder vector\n");
printf("\n vi=(0, ..., 1, A(i + 1, i), A(i + 2, i), ..., A(m, i))\n"); 
printf("\n \n ");
printf("\n but \n ");

printf("\n in output: runtime error Segmentation fault \n\n\n");

 ( void) gsl_linalg_QRPT_decomp(QRPT,tau, p, signum, norm ) ;    
   
       printf ("vector tau = \n");
       gsl_vector_fprintf (stdout, tau, "%f");
       
       gsl_permutation_free (p);
       gsl_vector_free (tau);
       gsl_vector_free (norm);
       gsl_matrix_free (QRPT);
 
 
// runtime error  Segmentation fault
       
       
       return 0;
     }