# include "data.h"

int svd (int argc, char **argv) {
  
  int i, Nx, Ny, status;
  float *in, *out, *A, *U, *W, *V;
  Vector *Vw;
  Buffer *Ma, *Mu, *Mv;

  if (argc != 6) goto usage;
  if (strcmp (argv[2], "=")) goto usage;

  if ((Ma = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
  if ((Mu = SelectBuffer (argv[3], ANYBUFFER, TRUE)) == NULL) return (FALSE);
  if ((Vw = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
  if ((Mv = SelectBuffer (argv[5], ANYBUFFER, TRUE)) == NULL) return (FALSE);

  Nx = Ma[0].header.Naxis[0];
  Ny = Ma[0].header.Naxis[1];

  /* U is Nx, Ny */
  gfits_free_matrix (&Mu[0].matrix);
  gfits_free_header (&Mu[0].header);
  Mu[0].bitpix = Ma[0].bitpix;
  Mu[0].unsign = Ma[0].unsign;
  Mu[0].bscale = Ma[0].bscale;
  Mu[0].bzero  = Ma[0].bzero;
  gfits_copy_header (&Ma[0].header, &Mu[0].header);
  gfits_create_matrix (&Mu[0].header, &Mu[0].matrix);
  
  /* V is Nx, Nx */
  gfits_free_matrix (&Mv[0].matrix);
  gfits_free_header (&Mv[0].header);
  Mv[0].bitpix = Ma[0].bitpix;
  Mv[0].unsign = Ma[0].unsign;
  Mv[0].bscale = Ma[0].bscale;
  Mv[0].bzero  = Ma[0].bzero;
  gfits_copy_header (&Ma[0].header, &Mv[0].header);
  gfits_modify (&Mv[0].header, "NAXIS2", "%d", 1, Nx);
  Mv[0].header.Naxis[1] = Nx;
  gfits_create_matrix (&Mv[0].header, &Mv[0].matrix);

  /* w is Nx */
  Vw[0].Nelements = Nx;
  REALLOCATE (Vw[0].elements, float, Nx);

  /* pointers to the various arrays */
  A = (float *) Ma[0].matrix.buffer;
  U = (float *) Mu[0].matrix.buffer;
  W = (float *) Vw[0].elements;
  V = (float *) Mv[0].matrix.buffer;

  /* copy A to U (svdcmp replaces A with U) */
  in  = A;
  out = U;
  for (i = 0; i < Nx*Ny; i++, in++, out++) *out = *in;
  /* use a bcopy instead? */

  status = svdcmp (U, W, V, Nx, Ny);
  if (!status) {
    gprint (GP_ERR, "error running svdcmp\n");
    return (FALSE);
  }
  return (TRUE);

 usage:
  gprint (GP_ERR, "USAGE: svd A = U w Vt\n");
  return (FALSE);
  
}

