/* This is -*- C -*- */
/* vim: set sw=2: */
/* $Id: guppi-vector.c,v 1.4 2002/01/14 05:01:23 trow Exp $ */

/*
 * guppi-vector.c
 *
 * Copyright (C) 2001 The Free Software Foundation, Inc.
 *
 * Developed by Jon Trowbridge <trow@gnu.org>
 */

/*
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA.
 */

#include <config.h>
#include "guppi-vector.h"

#include <math.h>
#include <string.h>
#include "guppi-memory.h"


GuppiVector *
guppi_vector_new (gint n)
{
  GuppiVector *v;
  gint i;

  g_return_val_if_fail (n > 0, NULL);

  v = guppi_new0 (GuppiVector, 1);
  v->n = n;
  v->v = guppi_new0 (double, n);
  
  for (i=0; i<n; ++i)
    v->v[i] = 0.0;

  v->epsilon = 1e-8; /* default epsilon... */

  return v;
}

GuppiVector *
guppi_vector_new_basis (gint n, gint i)
{
  GuppiVector *v;

  g_return_val_if_fail (n > 0, NULL);
  g_return_val_if_fail (0 <= i && i < n, NULL);

  v = guppi_vector_new (n);
  v->v[i] = 1.0;

  return v;
}

GuppiVector *
guppi_vector_copy (GuppiVector *v)
{
  GuppiVector *w;

  if (v == NULL)
    return NULL;

  g_return_val_if_fail (v->n > 0, NULL);
  g_return_val_if_fail (v->v, NULL);

  w = guppi_new0 (GuppiVector, 1);
  w->n = v->n;
  w->v = guppi_new0 (double, v->n);
  memcpy (w->v, v->v, sizeof (double) * v->n);

  w->epsilon = v->epsilon;

  return w;
}

void
guppi_vector_free (GuppiVector *v)
{
  if (v) {
    guppi_free (v->v);
    guppi_free (v);
  }
}

double
guppi_vector_norm_sq (GuppiVector *v)
{
  double run;
  double *p;
  gint i;
  g_return_val_if_fail (v != NULL, 0);

  run = 0;
  p = v->v;
  for (i = 0; i < guppi_vector_dim (v); ++i) {
    run += (*p) * (*p);
    ++p;
  }

  return run;
}

double
guppi_vector_norm (GuppiVector *v)
{
  return sqrt (guppi_vector_norm_sq (v));
}

gboolean
guppi_vector_is_zero (GuppiVector *v)
{
  gint i;
  double *p;
  g_return_val_if_fail (v != NULL, TRUE);
  
  p = v->v;
  for (i = 0; i < guppi_vector_dim (v); ++i) {
    if (fabs (*p) > v->epsilon)
      return FALSE;
  }

  return TRUE;
}

gboolean
guppi_vector_same_dim (GuppiVector *v, GuppiVector *w)
{
  if (v == NULL || w == NULL)
    return v == w;

  return guppi_vector_dim (v) == guppi_vector_dim (w);
}

gboolean
guppi_vector_equal (GuppiVector *v, GuppiVector *w)
{
  gint i;
  double *p, *q;
  double eps;

  if ((v == NULL || w == NULL) && v != w)
    return FALSE;
  if (! guppi_vector_same_dim (v, w))
    return FALSE;

  p = v->v;
  q = w->v;
  eps = v->epsilon + w->epsilon;

  for (i = 0; i < guppi_vector_dim (v); ++i) {
    if (fabs ((*p) - (*q)) > eps)
      return FALSE;
    ++p;
    ++q;
  }

  return TRUE;
}

void
guppi_vector_scalar_multiply (GuppiVector *v, double c)
{
  gint i;
  double *p;
  g_return_if_fail (v != NULL);

  p = v->v;
  for (i = 0; i < guppi_vector_dim (v); ++i) {
    *p *= c;
  }
}

double
guppi_vector_dot_product (GuppiVector *v, GuppiVector *w)
{
  gint i;
  double *p, *q;
  double run;
  g_return_val_if_fail (v != NULL, 0);
  g_return_val_if_fail (w != NULL, 0);
  g_return_val_if_fail (guppi_vector_same_dim (v, w), 0);

  p = v->v;
  q = w->v;
  run = 0;

  for (i = 0; i < guppi_vector_dim (v); ++i) {
    run += (*p) * (*q);
    ++p;
    ++q;
  }

  return run;
}

void
guppi_vector_normalize (GuppiVector *v)
{
  double norm;

  g_return_if_fail (v != NULL);

  norm = guppi_vector_norm (v);

  if (norm > v->epsilon)
    guppi_vector_scalar_multiply (v, 1/norm);
}

/* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** */

void
guppi_vector_spew (GuppiVector *v)
{
  gint i;
  g_return_if_fail (v != NULL);

  g_print ("[ ");
  for (i = 0; i < guppi_vector_dim (v); ++i) {
    g_print ("%g ", guppi_vector_entry (v, i));
  }
  g_print ("]\n");
}