/* solve_linear_eqs.c */
/* Solving a set of simultaneous linear equations by SVD */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "nrutil.h"
#include "svd.h"

#define TOL 1.0e-10


main()
{
  int i, j, m, n;
  double wmax, wmin, **a, **u, *w, **v, *b, *x;

  /* Specification of a problem */
  m = n = 2;
  a = dmatrix(1, m, 1, n);
  u = dmatrix(1, m, 1, n);
  v = dmatrix(1, n, 1, n);
  w = dvector(1, n);
  b = dvector(1, m);
  x = dvector(1, n);

  a[1][1] = 3;
  a[1][2] = -4;
  a[2][1] = 2;
  a[2][2] = 5;
  b[1] = 1;
  b[2] = 16;

  /* copy a into u */
  for (i = 1; i <= m; i++) {
    for (j = 1; j <= n; j++) {
      u[i][j] = a[i][j];
    }
  }

  printf("Coefficients :\n");
  for (i = 1; i <= m; i++) {
    for (j = 1; j <= n; j++) {
      printf("a[%d][%d] = %lf\t", i, j, a[i][j]);
    }
    printf("\n");
  }
  for (i = 1; i <= m; i++) printf("b[%d] = %lf\t", i, b[i]);
  printf("\n");

  /* SVD the matrix u */
  dsvdcmp(u, m, n, w, v);

  /* Set the threshold for singular values allowed to be nonzero. */
  wmax = 0.0;

  for (j = 1; j <= n; j++)
    if (w[j] > wmax) wmax = w[j];

  wmin = wmax*TOL;

  for (j = 1; j <= n; j++)
    if (w[j] < wmin) w[j] = 0.0;

  /* Now we can backsubstitute. */
  dsvbksb(u, w, v, m, n, b, x);

  printf("\nAnswer : \n");
  for (i = 1; i <= n; i++) printf("x[%d] = %lf\t ", i, x[i]);
  printf("\n");

  free_dmatrix(a, 1, m, 1, n);
  free_dmatrix(u, 1, m, 1, n);
  free_dmatrix(v, 1, n, 1, n);
  free_dvector(w, 1, n);
  free_dvector(b, 1, m);
  free_dvector(x, 1, n);

  return;
}