C       CENSUS DATA ANALYSIS WITH THE SVD
C       FROM SECTION 8.1
C     From the book "Numerical Methods and Software"
C          by  D. Kahaner, C. Moler, S. Nash
C               Prentice Hall 1988
C
      INTEGER LDX,N,P,LDU,LDV,JOB,INFO,I,J,JB
      PARAMETER(LDX=8,N=8,P=3,LDU=N,LDV=P,JOB=11)
      DOUBLE PRECISION POP(N),Y(N),X(LDX,P),S(P),E(P),
     *U(LDU,LDU),V(LDV,P),W(N),C(P),YEAR,POP80,
     *TOL,RELERR,SUM,R,RSQ,RI
C
C       C CONTAINS COEFFICIENTS OF POLYNOMIAL C(1)*1+C(2)*T+C(3)*T*T
C         T=YEAR (1900 ETC.)
C
      DATA POP/
     *   75.994575D0,
     *   91.972266D0,
     *  105.710620D0,
     *  122.775046D0,
     *  131.669275D0,
     *  150.697361D0,
     *  179.323175D0,
     *  203.235298D0/
C      
      DO 1 I=1,8
C        Y(I)=(1935.0D0-1900D0*(I-1))/1935.0D0
        Y(I)=1900.0D0+(I-1)*10D0
        X(I,1)=1.0D0
        X(I,2)=Y(I)
        X(I,3)=Y(I)**2
    1 CONTINUE
C
      CALL DSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,W,JOB,INFO)
C       NOTE: X IS DESTROYED BY ABOVE CALL!!!
C
C       WRITE SINGULAR VALUES (DESCENDING ORDER)
      WRITE(*,*) 'SINGULAR VALUES ARE: '
      WRITE(*,'(2X,D16.8)') (S(I),I=1,P)
      DO 2 I=1,P
        C(I)=0.0D0
    2 CONTINUE
C
C       RELERR REFLECTS NUMBER OF ACCURATE DIGITS IN DATA
C        E.G. 6 DIGITS ==> RELERR=1.D-6, ETC.
C       MAKING RELERR LARGER INCREASES RESIDUALS
      RELERR=1.D-6
      TOL=RELERR*S(1)
C
C       MULTIPLY U-TRANS * POP, AND SOLVE FOR COEFFICIENTS C(I)
C
      DO 60 J=1,P
         IF(S(J).LE.TOL)GO TO 60
         SUM=0.0D0
         DO 40 I=1,N
            SUM=SUM+U(I,J)*POP(I)
   40    CONTINUE
         SUM=SUM/S(J)
         DO 50 I=1,P
            C(I)=C(I)+SUM*V(I,J)
   50    CONTINUE
   60 CONTINUE
C
      WRITE(*,*) 'COEFFICIENTS (ASSUMING DATA GOOD TO 6 DIGITS) ARE:'
      WRITE(*,'(2X,D16.8)') (C(I),I=1,P)
      WRITE(*,*)
C
C       EVALUATE MODEL (HORNER'S RULE) AND RESIDUALS AT YEAR =1900,...,1980
C
      RSQ=0.0D0
      DO 75 I=1,9
         YEAR=1900.0D0+(I-1)*10.0D0
         POP80=0.0D0
         DO 70 JB=1,P
            J=P+1-JB
            POP80=YEAR*POP80+C(J)
   70    CONTINUE
         IF(I.LT.9) THEN
            RI=POP(I)-POP80
            RSQ=RSQ+RI*RI
            WRITE(*,'(A,I6,A)') 
     *       ' FOR YEAR',IDINT(YEAR),
     *' POP ESTM, MEAS, AND RESIDUAL ARE'
            WRITE(*,'(3D16.8)')POP80,POP(I),RI
         ELSE
            WRITE(*,'(A,I6,A,D16.8)') 
     *        ' FOR YEAR',IDINT(YEAR),' POP ESTMATE IS ',POP80
         ENDIF
   75 CONTINUE
C
      R=SQRT(RSQ)
      WRITE(*,'(1X,A,D16.8)')
     *'SQUARE ROOT OF RESIDUAL SUM OF SQUARES IS: ',R
C
      WRITE(*,*)
      WRITE(*,*)'REFERENCE RESULTS (PARTIAL) FROM IBM PC/AT'
      WRITE(*,*)'SINGULAR VALUES ARE: '
      WRITE(*,*)'   0.10594723E+08'
      WRITE(*,*)'   0.64774566E+02'
      WRITE(*,*)'   0.34620247E-03'
      WRITE(*,*)'COEFFICIENTS (ASSUMING DATA GOOD TO 6 DIGITS) ARE:'
      WRITE(*,*)'  -0.16714353E-02'
      WRITE(*,*)'  -0.16169731E+01'
      WRITE(*,*)'   0.87095700E-03'
      WRITE(*,*)'      .  '
      WRITE(*,*)'      .  '
      WRITE(*,*)'      .  '    
      WRITE(*,*)'FOR YEAR  1980 POP ESTMATE IS    0.21289144E+03'
      WRITE(*,*)'SQUARE ROOT OF RESIDUAL SUM OF SQUARES IS:  ',
     *'  0.16096596E+02'
      END