#define NR 6
#define NC 6
/*
   static double a[6][6]={
  {0.,0.,1.,-25.,15.,0.},
  {0.,0.,1.,-25.,-15.,0.},
  {0.,0.,1.,25.,-15.,0.},
  {0.,0.,1.,25.,15.,0.},
  {0.,0.,0.,0.,0.,0.},
  {0.,0.,0.,0.,0.,0.},
};
*/
static double a[6][6]={
  {0.,0.,1.,-25.,15.,0.},
  {0.,0.,1.,-25.,-15.,0.},
  {0.,0.,1.,25.,-15.,0.},
  {0.,0.,1.,25.,15.,0.},
  {0.,0.,0.,0.,0.,0.},
  {0.,0.,0.,0.,0.,0.}
};


#include <stdio.h>
#include "nrutil.h"
#include "nr.h"

void svdcmp_cnv(int n_max, double *a, double *u, int m, int n, double w[], double *v)
{
  double **anr,**unr,**vnr,*wnr;
  int i,j;
  int mnr;

  mnr = m;
  if (m < n) mnr = n;
  anr=convert_matrix(a,1,m,1,n_max);
  unr=convert_matrix(u,1,mnr,1,n_max);
  vnr=convert_matrix(v,1,n,1,n_max);
  wnr=w-1;

  for(i=1;i<=m;i++) {
    for (j=1;j<=n;j++) {
      unr[i][j]=anr[i][j];
    }
  }
  if (m < n) {
    for(i=m+1;i<=n;i++) {
      for (j=1;j<=n;j++) {
	unr[i][j]=0.0;
      }
    }
  }
  svdcmp(unr,mnr,n,wnr,vnr);
  free_convert_matrix(anr,1,m,1,n_max);
  free_convert_matrix(unr,1,mnr,1,n_max);
  free_convert_matrix(vnr,1,n,1,n_max);
}
void round_off_cnv(double w[], int n, double lim)
{
  double wmax,wmin;
  int j;

  wmax=0.0;
  for(j=0;j<n;j++) if (w[j] > wmax) wmax=w[j];
  wmin=wmax*(double)lim;
  for(j=0;j<n;j++) if (w[j] < wmin) w[j]=0.0;
}

void svbksb_cnv(int n_max, double *u, double w[], double *v, int m, int n, double b[], double x[])
{
  double **unr,**vnr,*wnr,*bnr,*xnr;
  int mnr;

  mnr = m;
  if (m < n) mnr = n;

  unr=convert_matrix(u,1,mnr,1,n_max);
  vnr=convert_matrix(v,1,n,1,n_max);
  wnr=w-1;
  bnr=b-1;
  xnr=x-1;
  svbksb(unr,wnr,vnr,mnr,n,bnr,xnr);
  free_convert_matrix(unr,1,mnr,1,n_max);
  free_convert_matrix(vnr,1,n,1,n_max);
}

int main()
{

  int i,j;
  double u[NR][NC],v[NR][NC],w[NC],b[NR],x[NC];

  svdcmp_cnv(NC,&a[0][0],&u[0][0],4,NC,w,&v[0][0]);

  printf("singular value=");
  for(i=0; i<NC; i++)	printf(" %f ",w[i]);
  printf("\n");

  for(i=0; i<NC; i++)
    printf("u[%d]= %f %f %f %f %f %f\n",i,
	   u[0][i],u[1][i],u[2][i],u[3][i],u[4][i],u[5][i]);

  for(i=0; i<NC; i++)
    printf("v[%d]= %f %f %f %f %f %f\n",i,
	   v[i][0],v[i][1],v[i][2],v[i][3],v[i][4],v[i][5]);

  round_off_cnv(w,NC,1.0e-6);
    
  for(j=0;j<NR;j++) b[j]=(double)j;

  svbksb_cnv(NR,&u[0][0],w,&v[0][0],4,NC,b,x);
    
  printf("x=   ");
  for(i=0;i<NC;i++) printf(" %f",x[i]);
  printf("\n");
  return 1;
}


