* Last processed by NICE on 12-Jun-2000 15:53:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
*	PROGRAM TEST
*       INTEGER NMAX
*       PARAMETER (NMAX=10000)
*	REAL*8 X(NMAX),Y(NMAX),WEIGHT(NMAX),A(20), B(NMAX)
*	REAL*8 U(NMAX,20),V(20,20),W(20),C(20)
*	INTEGER I,MA,NP
*	REAL*8 CHISQ
*	LOGICAL ERROR
*	EXTERNAL FPOLY,FCHEB
**
*	NP = NMAX
*	DO I=1,NP
*	    X(I) = (2*I-(1+NP))/FLOAT(NP-1)	! Normalize range
**	    Y(I) = x(i)**2+12*X(I)-13
*	    Y(I) = SIN(X(I))**3+SIN(X(I))**2
*	    WEIGHT(I) = 1.
*	ENDDO
*	write(6,*) 'MA ?'
*	read (5,*) MA
*	CALL SVDFIT(X,Y,WEIGHT,NP,A,MA,U,V,W,NP,MA,CHISQ,FPOLY,B,C,ERROR)
*	WRITE(6,*) 'Chi square: ',CHISQ
*       DO I=1,MA
*	  WRITE(6,'(A,E25.15,A,I2.2,A)') 'LET D ',A(I),'*A^',I-1,'+D'
*       ENDDO
*	CALL SVDFIT(X,Y,WEIGHT,NP,A,MA,U,V,W,NP,MA,CHISQ,FCHEB,B,C,ERROR)
*	WRITE(6,*) 'Chi square: ',CHISQ
*	WRITE(6,*) (A(I),I=1,MA)
*	END
*<FF>
      SUBROUTINE SVDFIT(X,Y,WEIGHT,NDATA,A,MA,U,V,W,MP,NP,CHISQ,
     $FUNCS,B,AFUNC,ERROR)
C---------------------------------------------------------------------------
C Singular value decomposition
C Given a set of NDATA points X(I), Y(I) with individual weights
C WEIGHT(I), use Chi2 minimization to determine the MA coefficients A of the
C fitting function. Here we solve the fitting
C equations using singular value decomposition of the NDATA by MA matrix.
C Arrays U,V,W provide workspace on input. On output they define the singular
C value decomposition and can be used to obtain the covariance matrix. MP, NP
C are the physical dimensions of the matrices U, V, W, as indicated below. It
C is necessary that MP>=NDATA, NP>=MA. The programs returns values for the MA
C fit parameters and Chi2, CHISQ. The user supplies a subroutine
C FUNCS(X,AFUNC,MA) that returns the MA basis functions evaluated at x=X in the
C array AFUNC.
C---------------------------------------------------------------------------
      INTEGER NDATA,MA,MP,NP
      REAL*8 X(NDATA),Y(NDATA),WEIGHT(NDATA),B(NDATA)
      REAL*8 A(MA),V(NP,NP),U(MP,NP),W(NP),AFUNC(MA)
      EXTERNAL FUNCS
      LOGICAL ERROR
*
      REAL*8 TOL
      PARAMETER (TOL=1D-5)
*
      INTEGER I,J
      REAL*8 WMAX,THRESH,CHISQ,SUM
*
      ERROR = .FALSE.
      DO I=1,NDATA
         CALL FUNCS (X(I),AFUNC,MA)
         DO J=1,MA
            U(I,J) = AFUNC(J)*WEIGHT(I)
         ENDDO
         B(I) = Y(I)*WEIGHT(I)
      ENDDO
      CALL SVDCMP(U,NDATA,MA,MP,NP,W,V,ERROR)
      IF (ERROR) RETURN
* Edit out the (nearly) singular values
      WMAX = 0.
      DO J=1,MA
         IF (W(J).GT.WMAX) WMAX=W(J)
      ENDDO
      THRESH = TOL*WMAX
      DO J=1,MA
         IF (W(J).LT.THRESH) W(J) = 0
      ENDDO
* Solve the equations
      CALL SVBKSB(U,W,V,NDATA,MA,MP,NP,B,A,ERROR)
      IF (ERROR) RETURN
* Evaluate chi-square
      CHISQ = 0.
      DO I=1,NDATA
         CALL FUNCS(X(I),AFUNC,MA)
         SUM=0.
         DO J=1,MA
            SUM = SUM+A(J)*AFUNC(J)
         ENDDO
         CHISQ = CHISQ+(Y(I)-SUM)**2*WEIGHT(I)
      ENDDO
      END
*<FF>
      SUBROUTINE SVDCMP(A,M,N,MP,NP,W,V,ERROR)
C---------------------------------------------------------------------------
C Singular value decomposition
C Given a matrix A with logical dimensions M by N and physical dimensions MP
C by NP, this routine computes its singular value decomposition
C	A = U.W.Vt
C The matrix U replaces A on output. The diagonal matrix of singular values W
C is a output as a vector W. The matrix V (not the transpose Vt) is output as
C V. M must be greater or equal to N. If it is smaller, then A should be filled
C up to square with zero rows.
C Typed from "numerical recipes".  TF
C---------------------------------------------------------------------------
      INTEGER M,N,MP,NP
      REAL*8 A(MP,NP),W(NP),V(NP,NP)
      LOGICAL ERROR
*
      INTEGER NMAX
      PARAMETER (NMAX=21)                ! Maximum anticipated value of N
      REAL*8 RV1(NMAX)
      INTEGER I,J,K,L,NM,ITS,JJ
      REAL*8 SCALE,ANORM,C,S,X,Y,Z
      REAL*8 F,G,H
*
      IF (N.GT.NMAX) THEN
         WRITE(6,'(A)')
     $   'E-SVDCMP,  NMAX dimension too small. Will need to recompile.'
         ERROR = .TRUE.
         RETURN
      ELSEIF (M.LT.N) THEN
         WRITE(6,'(A)')
     $   'E-SVDCMP, You must augment A with extra zero rows.'
         ERROR = .TRUE.
         RETURN
      ENDIF
* Householder reduction to diagonal form
      G = 0.0
      SCALE = 0.0
      ANORM = 0.0
      DO I=1,N
         L=I+1
         RV1(I) = SCALE*G
         G = 0.0
         S = 0.0
         SCALE = 0.0
         IF (I.LE.M) THEN
            DO K=I,M
               SCALE = SCALE+ABS(A(K,I))
            ENDDO
            IF (SCALE.NE.0.0) THEN
               DO K=I,M
                  A(K,I) = A(K,I)/SCALE
                  S = S+A(K,I)*A(K,I)
               ENDDO
               F = A(I,I)
               G = -SIGN(SQRT(S),F)
               H = F*G-S
               A(I,I) = F-G
               IF(I.NE.N) THEN
                  DO J=L,N
                     S=0.0
                     DO K=I,M
                        S = S+A(K,I)*A(K,J)
                     ENDDO
                     F = S/H
                     DO K=I,M
                        A(K,J) = A(K,J)+F*A(K,I)
                     ENDDO
                  ENDDO
               ENDIF
               DO K=I,M
                  A(K,I) = SCALE*A(K,I)
               ENDDO
            ENDIF
         ENDIF
         W(I) = SCALE*G
         G = 0.0
         S = 0.0
         SCALE = 0.0
         IF ((I.LE.M).AND.(I.NE.N)) THEN
            DO K=L,N
               SCALE = SCALE+ABS(A(I,K))
            ENDDO
            IF (SCALE.NE.0) THEN
               DO K=L,N
                  A(I,K)=A(I,K)/SCALE
                  S=S+A(I,K)*A(I,K)
               ENDDO
               F=A(I,L)
               G=-SIGN(SQRT(S),F)
               H=F*G-S
               A(I,L)=F-G
               DO K=L,N
                  RV1(K)=A(I,K)/H
               ENDDO
               IF (I.NE.M) THEN
                  DO J=L,M
                     S=0.0
                     DO K=L,N
                        S=S+A(J,K)*A(I,K)
                     ENDDO
                     DO K=L,N
                        A(J,K)=A(J,K)+S*RV1(K)
                     ENDDO
                  ENDDO
               ENDIF
               DO K=L,N
                  A(I,K)=SCALE*A(I,K)
               ENDDO
            ENDIF
         ENDIF
         ANORM = MAX(ANORM,(ABS(W(I))+ABS(RV1(I))))
      ENDDO
* Accumulation of right hand transformation
      DO I=N,1,-1
         IF (I.LT.N) THEN
            IF (G.NE.0.0) THEN
               DO J=L,N                  ! Double division avoids possible under
                  V(J,I) = (A(I,J)/A(I,L))/G
               ENDDO
               DO J=L,N
                  S=0.0
                  DO K=L,N
                     S=S+A(I,K)*V(K,J)
                  ENDDO
                  DO K=L,N
                     V(K,J)=V(K,J)+S*V(K,I)
                  ENDDO
               ENDDO
            ENDIF
            DO J=L,N
               V(I,J) = 0.0
               V(J,I) = 0.0
            ENDDO
         ENDIF
         V(I,I) = 1.0
         G = RV1(I)
         L = I
      ENDDO
* Accumulation of left hand transformations.
      DO I=N,1,-1
         L = I+1
         G = W(I)
         IF (I.LT.N) THEN
            DO J=L,N
               A(I,J) = 0.0
            ENDDO
         ENDIF
         IF (G.NE.0.0) THEN
            G = 1.0/G
            IF (I.NE.N) THEN
               DO J=L,N
                  S = 0.0
                  DO K=L,M
                     S = S+A(K,I)*A(K,J)
                  ENDDO
                  F = (S/A(I,I))*G
                  DO K=I,M
                     A(K,J)=A(K,J)+F*A(K,I)
                  ENDDO
               ENDDO
            ENDIF
            DO J=I,M
               A(J,I) = A(J,I)*G
            ENDDO
         ELSE
            DO J=I,M
               A(J,I) = 0.0
            ENDDO
         ENDIF
         A(I,I) = A(I,I)+1.0
      ENDDO
* Diagonalization of the bidiagonal form
      DO K=N,1,-1                        ! Loop over singular values
         DO ITS=1,30                     ! Loop over allowed iterations
            DO L=K,1,-1                  ! Test for splitting:
               NM=L-1                    ! Note that RV1(1) is always zero
               IF ((ABS(RV1(L))+ANORM).EQ.ANORM) GOTO 2
               IF ((ABS(W(NM))+ANORM).EQ.ANORM) GOTO 1
            ENDDO
1           C=0.0                        ! Cancellation of RV1(L), if L>1
            S=1.0
            DO I=L,K
               F=S*RV1(I)
               IF ((ABS(F)+ANORM).NE.ANORM) THEN
                  G = W(I)
                  H = SQRT(F*F+G*G)
                  W(I) = H
                  H = 1.0/H
                  C = (G*H)
                  S = -F*H
                  DO J=1,M
                     Y = A(J,NM)
                     Z = A(J,I)
                     A(J,NM)=(Y*C)+(Z*S)
                     A(J,I)=-(Y*S)+(Z*C)
                  ENDDO
               ENDIF
            ENDDO
2           Z=W(K)
            IF (L.EQ.K) THEN             ! Convergence
               IF (Z.LT.0.0) THEN        ! Singular value is made non negative
                  W(K)=-Z
                  DO J=1,N
                     V(J,K)=-V(J,K)
                  ENDDO
               ENDIF
               GOTO 3
            ENDIF
            IF (ITS.EQ.30) THEN
               WRITE(6,'(A)')
     $         'E-SVDCMP, No convergence in 30 iterations.'
               ERROR = .TRUE.
               RETURN
            ENDIF
            X=W(L)                       ! Shift from bottom 2-by-2 minor
            NM = K-1
            Y = W(NM)
            G = RV1(NM)
            H = RV1(K)
            F = ((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y)
            G = SQRT(F*F+1.0)
            F = ((X-Z)*(X+Z)+H*((Y/(F+SIGN(G,F)))-H))/X
* Next QR transformation
            C=1.0
            S=1.0
            DO J=L,NM
               I=J+1
               G=RV1(I)
               Y=W(I)
               H=S*G
               G=C*G
               Z=SQRT(F*F+H*H)
               RV1(J) = Z
               C=F/Z
               S=H/Z
               F= (X*C)+(G*S)
               G=-(X*S)+(G*C)
               H=Y*S
               Y=Y*C
               DO JJ=1,N
                  X=V(JJ,J)
                  Z=V(JJ,I)
                  V(JJ,J)= (X*C)+(Z*S)
                  V(JJ,I)=-(X*S)+(Z*C)
               ENDDO
               Z=SQRT(F*F+H*H)
               W(J) = Z
               IF (Z.NE.0) THEN          ! Rotation can be arbitrary if Z=0
                  Z=1.0/Z
                  C=F*Z
                  S=H*Z
               ENDIF
               F= (C*G)+(S*Y)
               X=-(S*G)+(C*Y)
               DO JJ=1,M
                  Y=A(JJ,J)
                  Z=A(JJ,I)
                  A(JJ,J) = (Y*C)+(Z*S)
                  A(JJ,I) =-(Y*S)+(Z*C)
               ENDDO
            ENDDO
            RV1(L)=0.0
            RV1(K)=F
            W(K)=X
         ENDDO
3        CONTINUE
      ENDDO
      RETURN
      END
*<FF>
      SUBROUTINE SVBKSB(U,W,V,M,N,MP,NP,B,X,ERROR)
C---------------------------------------------------------------------------
C Singular value decomposition
C Solves A.X = B for a vector X, where A is specified by the arrays U,W,V as
C returned by SVDCMP. M and N are the logical dimensions of A, and will be
C equal for square matrices. MP and NP are the physical dimensions of A. B
C is the input right hand side. X is the output soolution vector. No input
C quantities are destroyed, so the routine may be called sequentially with
C different B's.
C Typed from "numerical recipes".  TF
C---------------------------------------------------------------------------
      INTEGER M,N,MP,NP
      REAL*8 U(MP,NP),W(NP),V(NP,NP),B(MP),X(NP)
      LOGICAL ERROR
*
      INTEGER NMAX
      PARAMETER (NMAX=21)
      REAL*8 TMP(NMAX),S
      INTEGER I,J,JJ
*
      IF (N.GT.NMAX) THEN
         WRITE(6,'(A)')
     $   'E-SVDCMP,  NMAX dimension too small. Will need to recompile.'
         ERROR = .TRUE.
         RETURN
      ENDIF
* Calculate UtB
      DO J=1,N
         S = 0.
         IF (W(J).NE.0.) THEN
            DO I=1,M
               S=S+U(I,J)*B(I)
            ENDDO
            S = S/W(J)
         ENDIF
         TMP(J) = S
      ENDDO
* Matrix multiply by V to get answer
      DO J=1,N
         S = 0.0
         DO JJ=1,N
            S = S+V(J,JJ)*TMP(JJ)
         ENDDO
         X(J) = S
      ENDDO
      END
*<FF>
      SUBROUTINE FPOLY(X,P,NP)
* Fitting routine for a polynomial of degree NP-1.
      INTEGER NP
      REAL*8 X,P(NP)
*
      INTEGER J
*
      P(1) = 1.
      DO J=2,NP
         P(J)=P(J-1)*X
      ENDDO
      END
*<FF>
      SUBROUTINE FCHEB(X,P,NP)
* Fitting routine for a Chebychev polynomial of degree NP-1.
* Compute all coefficients of the polynomials, from the usual recurrence
* Beware of first point for "true" Chebychev polynomials
      INTEGER NP
      REAL*8 X,P(NP)
*
      INTEGER J
*
      P(1) = 1.D0
      P(2) = X
      DO J=3,NP
         P(J)=2*P(J-1)*X-P(J-2)
      ENDDO
      P(1) = .5D0
      END
