/* 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'.
*
* 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
void svdbak(int m, int n, double *q,double **u,double **v,double *x,double *b)
{
double *t,s;
int i,j;
t = (double *)calloc(n,sizeof(double));
for (i=0;i