1 files changed, 109 insertions, 0 deletions
diff --git a/src/slalib/svdsol.f b/src/slalib/svdsol.f
new file mode 100644
index 0000000..46ef6bb
--- /dev/null
+++ b/src/slalib/svdsol.f
@@ -0,0 +1,109 @@
+      SUBROUTINE sla_SVDSOL (M, N, MP, NP, B, U, W, V, WORK, X)
+*+
+*     - - - - - - -
+*      S V D S O L
+*     - - - - - - -
+*
+*  From a given vector and the SVD of a matrix (as obtained from
+*  the SVD routine), obtain the solution vector (double precision)
+*
+*  This routine solves the equation:
+*
+*     A . x = b
+*
+*  where:
+*
+*     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 factorised (for example by
+*  the routine sla_SVD) into the following components:
+*
+*     A = U x W x VT
+*
+*  where:
+*
+*     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)
+*
+*  Notes:
+*
+*  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 factorisation 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.
+*
+*  Given:
+*     M,N    i         numbers of rows and columns in matrix A
+*     MP,NP  i         physical dimensions of array containing matrix A
+*     B      d(M)      known vector b
+*     U      d(MP,NP)  array containing MxN matrix U
+*     W      d(N)      NxN diagonal matrix W (diagonal elements only)
+*     V      d(NP,NP)  array containing NxN orthogonal matrix V
+*
+*  Returned:
+*     WORK   d(N)      workspace
+*     X      d(N)      unknown vector x
+*
+*  Reference:
+*     Numerical Recipes, section 2.9.
+*
+*  P.T.Wallace   Starlink   29 October 1993
+*
+*  Copyright (C) 1995 Rutherford Appleton Laboratory
+*-
+
+      IMPLICIT NONE
+
+      INTEGER M,N,MP,NP
+      DOUBLE PRECISION B(M),U(MP,NP),W(N),V(NP,NP),WORK(N),X(N)
+
+      INTEGER J,I,JJ
+      DOUBLE PRECISION S
+
+
+
+*  Calculate [diag(1/Wj)] . transpose(U) . b (or zero for zero Wj)
+      DO J=1,N
+         S=0D0
+         IF (W(J).NE.0D0) THEN
+            DO I=1,M
+               S=S+U(I,J)*B(I)
+            END DO
+            S=S/W(J)
+         END IF
+         WORK(J)=S
+      END DO
+
+*  Multiply by matrix V to get result
+      DO J=1,N
+         S=0D0
+         DO JJ=1,N
+            S=S+V(J,JJ)*WORK(JJ)
+         END DO
+         X(J)=S
+      END DO
+
+      END