SUBROUTINE CDSVBK (
*
*  inputs
*
     :                  U, W, V, M, N, MP, NP, B, TMP,
*
*  outputs
*
     :                  X)
*
*  Module number:
*
*  Module name:
*
*  Keyphrase:
*  ----------
*  'back-substitution' of the singular value decomposition (SVD) procedure
*
*  Description:
*  ------------
*  adapted from Numerical Recipes by William Press et. al., p. 57, (1986)
*  subroutine SVBKSB with double precision
*
*  solve A * X = B for a vector X, where A is specified by the arrays U, W, V as
*  returned from a singular value decomposition routine.  Input quantities are
*  preserved after the call of this routine
*
*  FORTRAN name: CDSVBK.FOR
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       none
*  SDAS:
*       none
*  Others:
*       None
*
*  History:
*  --------
*  Version      Date        Author          Description
*     1       12-07-88     J.-C. HSU        coding
*-------------------------------------------------------------------------------
*
*== input:
*                                       --logical dimensions of the matrix U
        INTEGER                 M, N,
*                                       --physical dimensions of the matrix U
     :                          MP, NP
*                                       --component matrices
        DOUBLE PRECISION        U(MP,NP), W(NP), V(NP,NP),
*                                       --input vector
     :                          B(MP),
*                                       --working array
     :                          TMP(NP)
*
*== output:
*                                       --output vector
        DOUBLE PRECISION        X(NP)
*
*== local:
*                                       --array indices
        INTEGER                 I, J, JJ
        DOUBLE PRECISION        S
*------------------------------------------------------------------------------
*
*  calculate transpose(U) * B
*
        DO 20 J = 1, N
            S = 0.
*
*  nonzero result only if w(j) is nonzero
*
            IF (W(J) .NE. 0.) THEN
                DO 10 I = 1, M
                    S = S + U(I,J) * B(I)
   10           CONTINUE
*
*  this is the divide by w(j)
*
                S = S / W(J)
            ENDIF
*
            TMP(J) = S
   20   CONTINUE
*
*  matrix multiply by V to get answer
*
        DO 40 J = 1, N
            S = 0.
            DO 30 JJ = 1, N
                S = S + V(J,JJ) * TMP(JJ)
   30       CONTINUE
*
            X(J) = S
   40   CONTINUE
*
        RETURN
        END