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

int main ()
{
  /* a =

   0.918312   0.891045   0.575062   0.548165
   0.046109   0.638387   0.365323   0.680073
   0.816039   0.903395   0.227423   0.499313
   0.184944   0.199351   0.674672   0.625001

   [q,r,p] = qr(a)

   w =

   0.164588
   0.502226
   0.655178
   0.055821

   v =

   0.97662
   0.27795
   0.39237
   0.19906
     
  */

  double q_data[] = {  -0.621218,   0.166224,  -0.417979,  -0.641678,
                       -0.445071,   0.045123,   0.884341,  -0.133478,
                       -0.629828,  -0.394970,  -0.200523,   0.638048,
                       -0.138983,   0.902404,  -0.054994,   0.404138 };

  double r_data[] = {   -1.43435,  -0.75684,  -1.13066,  -1.04456,
                        0.00000,   0.63107,  -0.00069,   0.48859,
                        0.00000,   0.00000,  -0.51686,   0.23780,
                        0.00000,   0.00000,   0.00000,   0.12865 };

  double p_data[] = {   0,   0,   1,   0,      /* p = 1,2,0,3 */
                        1,   0,   0,   0,
                        0,   1,   0,   0,
                        0,   0,   0,   1 } ;

  double w_data[] = {   0.164588,
                        0.502226,
                        0.655178,
                        0.055821 } ;

  double v_data[] = {   0.97662,   
                        0.27795,   
                        0.39237,   
                        0.19906 } ;

  gsl_matrix_view q = gsl_matrix_view_array (q_data, 4, 4);
  gsl_matrix_view r = gsl_matrix_view_array (r_data, 4, 4);
  gsl_vector_view w = gsl_vector_view_array (w_data, 4);
  gsl_vector_view v = gsl_vector_view_array (v_data, 4);

  gsl_permutation * p = gsl_permutation_alloc(4);
  p->data[0] = 1; p->data[1] = 2; p->data[2] = 0; p->data[3] = 3;

  gsl_matrix * z = gsl_matrix_alloc(4,4);

  gsl_linalg_QRPT_update (&q.matrix, &r.matrix, p, &w.vector, &v.vector);

  printf("q2="); gsl_matrix_fprintf(stdout, &q.matrix, "%g");
  printf("r2="); gsl_matrix_fprintf(stdout, &r.matrix, "%g");

  /* q*(r+w*v'*p)

     ans =

     0.799757   0.446193   0.597557   0.482787
     0.783298   0.569891   0.555280   0.783854
     0.792831   0.071342   0.427551   0.420130
     0.315217   0.838238   0.592062   0.707928

  */

  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
                 1.0, &q.matrix, &r.matrix,
                 0.0, z);

  printf("z="); gsl_matrix_fprintf(stdout, z, "%g");

}