typedef unsigned int    size_t;

struct gsl_block_struct
{
  size_t size;
  double *data;
};

typedef struct gsl_block_struct gsl_block;

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  double * data;
  gsl_block * block;
  int owner;
} gsl_matrix;

gsl_matrix *gsl_matrix_alloc (const size_t n1, const size_t n2);
void gsl_matrix_free (gsl_matrix * m);

inline
double
gsl_matrix_get(const gsl_matrix * m, const size_t i, const size_t j)
{
  return m->data[i * m->tda + j] ;
}

inline
void
gsl_matrix_set(gsl_matrix * m, const size_t i, const size_t j, const double x)
{
  m->data[i * m->tda + j] = x ;
}

int gsl_matrix_transpose_memcpy (gsl_matrix * dest, const gsl_matrix * src);

extern double cos  ( double );
extern double sin  ( double );
extern double fabs  ( double );

typedef struct
{
  size_t size;
  size_t stride;
  double *data;
  gsl_block *block;
  int owner;
}
gsl_vector;

gsl_vector *gsl_vector_alloc (const size_t n);
void gsl_vector_free (gsl_vector * v);

inline
void
gsl_vector_set (gsl_vector * v, const size_t i, double x)
{
  v->data[i * v->stride] = x;
}

void gsl_test (int status, const char *test_description, ...);

int  gsl_blas_dger (double alpha,
                    const gsl_vector * X,
                    const gsl_vector * Y,
                    gsl_matrix * A);
int gsl_linalg_LQ_decomp (gsl_matrix * A, gsl_vector * tau);
int gsl_linalg_LQ_unpack (const gsl_matrix * LQ, const gsl_vector * tau,
			  gsl_matrix * Q, gsl_matrix * L);

enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113};
typedef  enum CBLAS_TRANSPOSE   CBLAS_TRANSPOSE_t;

int  gsl_blas_dgemv (CBLAS_TRANSPOSE_t TransA,
                     double alpha,
                     const gsl_matrix * A,
                     const gsl_vector * X,
                     double beta,
                     gsl_vector * Y);
int gsl_linalg_LQ_update (gsl_matrix * Q, gsl_matrix * R,
			  const gsl_vector * v, gsl_vector * w);

int  gsl_blas_dgemm (CBLAS_TRANSPOSE_t TransA,
                     CBLAS_TRANSPOSE_t TransB,
                     double alpha,
                     const gsl_matrix * A,
                     const gsl_matrix * B,
                     double beta,
                     gsl_matrix * C);

int
check (double x, double actual, double eps)
{
  if (x == actual)
    {
      return 0;
    }
  else if (actual == 0)
    {
      return fabs(x) > eps;
    }
  else
    {
      return (fabs(x - actual)/fabs(actual)) > eps;
    }
}

gsl_matrix *
create_general_matrix(unsigned long size1, unsigned long size2)
{
  unsigned long i, j;
  gsl_matrix * m = gsl_matrix_alloc(size1, size2);
  for(i=0; i<size1; i++) {
    for(j=0; j<size2; j++) {
      gsl_matrix_set(m, i, j, 1.0/(i+j+1.0));
    }
  }
  return m;
}

gsl_matrix * m53;

int
test_LQ_update_dim(const gsl_matrix * m, double eps)
{
  int s = 0;
  unsigned long i,j, M = m->size1, N = m->size2;

  gsl_matrix * lq1  = gsl_matrix_alloc(N,M);
  gsl_matrix * lq2  = gsl_matrix_alloc(N,M);
  gsl_matrix * q1  = gsl_matrix_alloc(M,M);
  gsl_matrix * l1  = gsl_matrix_alloc(N,M);
  gsl_matrix * q2  = gsl_matrix_alloc(M,M);
  gsl_matrix * l2  = gsl_matrix_alloc(N,M);
  gsl_vector * d2 = gsl_vector_alloc( ( ( M ) < ( N ) ? ( M ) : ( N ) ));
  gsl_vector * u = gsl_vector_alloc(M);
  gsl_vector * v = gsl_vector_alloc(N);
  gsl_vector * w = gsl_vector_alloc(M);

  gsl_matrix_transpose_memcpy(lq1,m);
  gsl_matrix_transpose_memcpy(lq2,m);
  for(i=0; i<M; i++) gsl_vector_set(u, i, sin(i+1.0));
  for(i=0; i<N; i++) gsl_vector_set(v, i, cos(i+2.0) + sin(i*i+3.0));

  gsl_blas_dger(1.0, v, u, lq1);

  s += gsl_linalg_LQ_decomp(lq2, d2);
  s += gsl_linalg_LQ_unpack(lq2, d2, q2, l2);

  gsl_blas_dgemv(CblasNoTrans, 1.0, q2, u, 0.0, w);

  s += gsl_linalg_LQ_update(q2, l2, v, w);

  gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,l2,q2,0.0,lq2);

  for(i=0; i<N; i++) {
    for(j=0; j<M; j++) {
      double s1 = gsl_matrix_get(lq1, i, j);
      double s2 = gsl_matrix_get(lq2, i, j);

      int foo = check(s1, s2, eps);
      s += foo;
    }
  }

  gsl_vector_free(d2);
  gsl_vector_free(u);
  gsl_vector_free(v);
  gsl_vector_free(w);
  gsl_matrix_free(lq1);
  gsl_matrix_free(lq2);
  gsl_matrix_free(q1);
  gsl_matrix_free(l1);
  gsl_matrix_free(q2);
  gsl_matrix_free(l2);

  return s;
}

int test_LQ_update(void)
{
  int f;
  int s = 0;

  f = test_LQ_update_dim(m53, 2 * 512.0 *  2.2204460492503131e-16 );
  gsl_test(f, "  LQ_update m(5,3)");
  s += f;

  return s;
}

int main(void)
{
  m53 = create_general_matrix(5,3);
  gsl_test(test_LQ_update(),             "LQ Rank-1 Update");
}