Re: [Help-gsl] Please, run the two programms.

[Jo_Jo]

Brian Gough wrote:
> At Thu, 18 Oct 2007 12:07:13 +0200,
> andrzej buczynski wrote:
>   
>> Please, run the two programms.
>> In both programms vector v is not equal to vector Q*v, or Qtraqnsp * v.
>>     
>
> The function requires the input matrix to be decomposed with
> gsl_linalg_QR_decomp first.
>
>   
Dear Brian,

In  the Chapter 13 of the User Manual in pdf -format, about  Linear
Algebra,  page 130
on reads:
" .....                                                                
                     It may be useful if the
product b′ = QT b has already been computed using gsl_linalg_QR_QTvec."

And about gsl_linalg_QR_QTvec one reads:

"int gsl_linalg_QR_QTvec (const gsl matrix * QR, const gsl vector
*                 [Function]
          tau, gsl vector * v )
     This function applies the matrix QT encoded in the decomposition
(QR,tau) to the
     vector v, storing the result QT v in v. The matrix multiplication
is carried out directly
     using the encoding of the Householder vectors without needing to
form the full matrix
     QT ."

That means, that these vector v = QT * b  ?

And about function int gsl_linalg_QR_Qvec one read :
" Function  int gsl_linalg_QR_Qvec (const gsl matrix * QR, const gsl
vector *                [Function]
         tau, gsl vector * v )
     This function applies the matrix Q encoded in the decomposition
(QR,tau) to the
     vector v, storing the result Qv in v. The matrix multiplication is
carried out directly
     using the encoding of the Householder vectors without needing to
form the full matrix
     Q."

That means that these vector v = Q * b ?

These two descriptions of both functions should be :

1.)  in the first  description:

"Function  int gsl_linalg_QR_QTvec (const gsl matrix * QR, const gsl
vector *                 [Function]
          tau, gsl vector * b )
     This function applies the matrix QT encoded in the decomposition
(QR,tau) to the
     vector b, storing the result QT* b in b. The matrix multiplication
is carried out directly
     using the encoding of the Householder vectors without needing to
form the full matrix
     QT ."

2.) in the  second description:

" Function  int gsl_linalg_QR_Qvec (const gsl matrix * QR, const gsl
vector *                [Function]
         tau, gsl vector * b )
     This function applies the matrix Q encoded in the decomposition
(QR,tau) to the
     vector b, storing the result Q*b in b. The matrix multiplication is
carried out directly
     using the encoding of the Householder vectors without needing to
form the full matrix
     Q."


In my two programs :

- program  QRQvec.c :


//------program QRQvec.c
/*

  int gsl_linalg_QR_Qvec (const gsl matrix * QR,
                          const gsl vector * tau,
                          gsl vector * v )

* This function applies the matrix Q encoded in the decomposition
*  (QR,tau) to the vector v, storing the result Qv in v.
*
*  Should be :   v = Qv
* 
*/

#include <stdio.h>
#include <gsl/gsl_linalg.h>

int main (void)
{
  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,
                                          0.33, 0.44, -0.66, 0.12 };
       
       double b_data[] = { 1.0, 1.0, 1.0, 1.0 };

       int M = 4, N = 4 , i ;

       gsl_vector * tau = gsl_vector_alloc( GSL_MIN_INT ( M, N ) );
       gsl_matrix_view m = gsl_matrix_view_array (a_data, M, N);
       gsl_vector_set_zero ( tau ) ;
       gsl_vector_view v = gsl_vector_view_array (b_data, N);
 
 printf("-----vector v   ---------\n");
for (i = 0; i < M; i++) {
    printf ("%f ",gsl_vector_get (&v.vector, i ) );
    printf("\n"); };

printf("now, we are calling the function gsl_linalg_QR_Qvec()--\n");

       (void) gsl_linalg_QR_Qvec (&m.matrix, tau, &v.vector ) ;

printf(" and the-vector v should be Q *v , but is't ---------\n");

for (i = 0; i < M; i++) {
    printf ("%f ",gsl_vector_get (&v.vector, i ) );
   
    printf("\n"); };



 printf(" the-vector tau is not chenged too---------\n");
 for (i = 0; i < M; i++) {
    printf ("%f ",gsl_vector_get (tau, i ) );
    printf("\n"); };
 
 
     gsl_vector_free(tau);

       return 0;
  } ;
//---------- end of program QRQvec.c



- and in program QRQTvec.c  :


//-----program QRQTvec.c

/*
  int gsl_linalg_QR_QTvec (const gsl_matrix * QR,
                           const gsl_vector * tau,
                                  gsl_vector * v )
 
This function applies the matrix QT encoded in the decomposition
(QR,tau) to the vector v, storing the result QT*v in v.
The matrix multiplication is carried out directly
using the encoding of the Householder vectors without
needing to form the full matrix  QT .
   
*/

#include <stdio.h>
#include <gsl/gsl_linalg.h>

int main (void)
{
  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,
                                         0.33, 0.44, -0.66, 0.12 };
       
       double b_data[] = { 1.0, 1.0, 1.0, 1.0 };

       int M = 4, N = 4 , i ;

       gsl_vector * tau = gsl_vector_alloc( GSL_MIN_INT ( M, N ) );
       gsl_matrix_view m = gsl_matrix_view_array (a_data, M, N);
//       gsl_vector_set_zero ( tau ) ;
       gsl_vector_view v = gsl_vector_view_array (b_data, N);

       (void) gsl_linalg_QR_QTvec (&m.matrix, tau, &v.vector ) ;

   printf("-----vector v = QTrans*v---------\n");

for (i = 0; i < M; i++) {
    printf ("%f ",gsl_vector_get (&v.vector, i ) );
   
    printf("\n"); };

 printf("-----vector tau---------\n");
 for (i = 0; i < M; i++) {
    printf ("%f ",gsl_vector_get (tau, i ) );
    printf("\n"); };

 
 
     gsl_vector_free(tau);
       return 0;
     }
//---end of programm QRQTvec-c-------------------------

first, I  have initialize the vector v=b  with the values: 1.0, 1.0,
1.0, 1.0   , and
then I've tried to obtain only the vectors  Q*b  or QT*b respectively.

But in both cases the vector v = b is unchanged and  v = b = {  1.0,
1.0, 1.0, 1.0  } in output.

Please, run these two  programs.  For some reason I need  from the
linear system Ax = b
only the right hand  sites  of the systems:
          Rx = QT* b , or
          Rx = Q * b. 

 These two functions

    gsl_linalg_QR_QTvec() and gsl_linalg_QR_Qvec() don't provides the
result.

Best regards


Andrzej Buczynski
Moosburg, Germany