From a given vector and the SVD of a matrix (as obtained from the slaSvd routine), obtain the solution vector. (double precision)
a . x = b
a is a given m (rows) x n (columns) matrix, where m.ge.n 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
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.
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