/* Program to perform Gram-Schmidt orthogonalization. */
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

typedef double real;

/* Prototypes */
static real **row_vectors	( const real *, int, int );
static real *vector_subtract_unit_component	( real *, const real *, int);
static real *vector_normalize		( real *, int);
static real inner_product	( const real *, const real *, int);
static real norm		( const real *, int);
static void gram_schmidt	( real **, int, int);
static void print_vector	( const char *, real *, int n);

/* Implementations */
static real **
row_vectors (
	     const real *matrix,
	     int rows,
	     int cols)
{
  int i,j;
  real **out;
  out = (real **)malloc(rows*sizeof(real *));  assert(out);
  for(i=0; i<rows; i++){
    out[i]=(real *)malloc(cols*sizeof(real)); assert(out[i]);
    for(j=0; j<cols; j++) out[i][j] = matrix[i*cols+j];
  }
  return out;
}

static real *
vector_subtract_unit_component(
			       real *out,
			       const real *in,
			       int dimen)
{
  int k;
  real c;
  assert(out); assert(in);
  assert(dimen>=0);
  c = inner_product(out, in, dimen);
  if(c)
    for( k=0; k<dimen; k++) out[k] -= c*in[k];
  return out;
}

static real *
vector_normalize(
		 real *inout,
		 int dimen)
{
  int k;  real c;
  assert(inout);
  assert(dimen>=0);
  c = norm(inout, dimen);  assert(c);
  for( k=0; k<dimen; k++) inout[k] *= 1/c;
  return inout;
}

static real
inner_product(
	      const real *in1,
	      const real *in2,
	      int dimen)
{
  int k;
  real out=0;
  assert(in1); assert(in2);
  assert(dimen>=0);
  for( k=0; k<dimen; k++) out += in1[k]*in2[k];
  return out;
}

static real
norm(
     const real *in,
     int dimen)
{
  return (real)sqrt(inner_product(in,in,dimen));
}

static void
gram_schmidt (
	      real **b,
	      int k,
	      int dimen)
{
  int i, j;
  real c;
  assert(b);			/* assume non-NULL input arrays */
  assert(k>=0);			/* assume non-negative subspace dimension */
  assert(dimen>=k);		/* assume adequate dimensionality */
  
  for(i=0; i<k; i++){
    for(j=0; j<i; j++)
      vector_subtract_unit_component( b[i], b[j], dimen);
    vector_normalize( b[i], dimen);
  }
  return;
}

static void
print_vector( 
	     const char *title,
	     real *x,
	     int n)
{
  int i;
  printf("%s (dim=%d): ", title, n);
  for(i=0; i<n; i++ ) printf("%9.4g ", x[i]);
  putchar('\n');
  return;
}


int
main(void)
{
  int dim=4, k=3, n;
  real mat[3*4]={1,1,1,1,  1,0,2,0,  1,0,0,3}, **vecs;

  vecs = row_vectors(mat, k, dim);
  gram_schmidt(vecs, k, dim);
  for(n=0; n<k; n++)
    print_vector("", vecs[n], dim);
  return 0;
}