/*  svdbak.c -- Solves Ax = b where A has been decomposed by svd.
 *              Inputs to this routine are the svd of A, and the
 *              vectors: 'x' and 'b'
 *  (C) 2000, C. Bond. All rights reserved..
 *
 *  Solve least squares problem by the following means:
 *                                        T         T        T
 *  Compute svd of 'A', yielding A = U.q.V , where U . U  = V . V  = I
 *  and  'q' is a diagonal matrix containing the singular values.
 *
 *  Then,
 *                                    T
 *      A.x = b is equivalent to U.q.V .x = b
 *         T      T        T      T
 *  So,   U .U.q.V .x = q.V .x = U .b
 *  
 *  Since 'q' is a diagonal matrix, its inverse simply consists of the
 *  reciprocals of the elements along the diagonal.
 *                 T                  T
 *  We now have,  V .x = (diag(1/q)).U .b
 *
 *  and, finally, since 'V' is orthogonal,
 *                         T
 *      x = V.(diag(1/q)).U .b
 *
 * This last formulation is implemented below.
 */                               
#include <malloc.h>
void svdbak(int m, int n, double *q,double **u,double **v,double *x,double *b)

{
	double *t;
	int i,j;

	t = (double *)calloc(n,sizeof(double));
	for (i=0;i<n;i++) {    
		t[i] = 0.0;
		if (q[i]) {
			for (j=0;j<m;j++) {
				t[i] += (u[j][i] * b[j]);  /*   U'.b       */
			}
			t[i] /= q[i];                  /*  (1/q).U'.b  */
		}
	}
	for (i=0;i<n;i++) {
		x[i] = 0.0;
		for (j=0;j<n;j++) {
			x[i] += (v[i][j] * t[j]);      /* V.(1/q).U'.b */
		}
	}
	free(t);
}