slaSvdsol -


void slaSvdsol(m, n, mp, np, b, u, w, v, work, x)


int m
int n
int mp
int np
double *b
double *u
double *w
double *v
double *work
double *x


  From a given vector and the SVD of a matrix (as obtained from
  the slaSvd routine), obtain the solution vector.

  (double precision)

This routine solves the equation

     a . x = b


     a   is a given m (rows) x n (columns) matrix, where
     x   is the n-vector we wish to find
     b   is a given m-vector

  By means of the singular value decomposition method (SVD).  In
  this method, the matrix a is first factorized (for example by

the routine slaSvd) into the following components

     a = u x w x vt


     a   is the m (rows) x n (columns) matrix
     u   is an m x n column-orthogonal matrix
     w   is an n x n diagonal matrix with w(i,i).ge.0
     vt  is the transpose of an nxn orthogonal matrix

     Note that m and n, above, are the logical dimensions of the
     matrices and vectors concerned, which can be located in
     arrays of larger physical dimensions mp and np.

The solution is found from the expression

     x = v . [diag(1/wj)] . ( transpose(u) . b )


  1)  If matrix a is square, and if the diagonal matrix w is not
      adjusted, the method is equivalent to conventional solution
      of simultaneous equations.

  2)  If m>n, the result is a least-squares fit.

  3)  If the solution is poorly determined, this shows up in the
      SVD factorization as very small or zero wj values.  Where
      a wj value is small but non-zero it can be set to zero to
      avoid ill effects.  The present routine detects such zero
      wj values and produces a sensible solution, with highly
      correlated terms kept under control rather than being allowed
      to elope to infinity, and with meaningful values for the
      other terms.


     m,n    int            numbers of rows and columns in matrix a
     mp,np  int            physical dimensions of array containing matrix a
     b      double[m]      known vector b
     u      double[mp][np] array containing mxn matrix u
     w      double[n]      nxn diagonal matrix w (diagonal elements only)
     v      double[np][np] array containing nxn orthogonal matrix v


     work   double[n]      workspace
     x      double[n]      unknown vector x

  Note:  If the relative sizes of m, n, mp and np are inconsistent,
         the vector x is returned unaltered.  This condition should
         have been detected when the SVD was performed using slaSvd.


     Numerical Recipes, Section 2.9.

  P.T.Wallace   Starlink   21 December 1993