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
      INTEGER LDX,N,P,LDU,LDV,JOB,INFO,I,J,JB
      PARAMETER(LDX=8,N=8,P=3,LDU=N,LDV=P,JOB=11)
      REAL POP(N),Y(N),X(LDX,P),S(P),E(P),U(LDU,LDU),V(LDV,P),W(N)
      REAL 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.994575,
     *   91.972266,
     *  105.710620,
     *  122.775046,
     *  131.669275,
     *  150.697361,
     *  179.323175,
     *  203.235298/
C      
      DO 1 I=1,8
C        Y(I)=(1935.0-1900*(I-1.))/1935.0
        Y(I)=1900.0+(I-1)*10
        X(I,1)=1.0
        X(I,2)=Y(I)
        X(I,3)=Y(I)**2
    1 CONTINUE
C
      CALL SSVDC(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(*,*) (S(I),I=1,P)
      DO 2 I=1,P
        C(I)=0.0
    2 CONTINUE
C
C       RELERR REFLECTS NUMBER OF ACCURATE DIGITS IN DATA
C        E.G. 6 DIGITS ==> RELERR=1.E-6, ETC.
C       MAKING RELERR LARGER INCREASES RESIDUALS
      RELERR=1.E-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.0
         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(*,*) (C(I),I=1,P)
C
C       EVALUATE MODEL (HORNER'S RULE) AND RESIDUALS AT YEAR =1900,...,1980
C
      RSQ=0.0
      DO 75 I=1,9
         YEAR=1900.0+(I-1)*10.0
         POP80=0.0
         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,3F10.2)') 
     *       ' FOR YEAR',IFIX(YEAR),' POP ESTM, MEAS, AND RESIDUAL ARE '
     *               ,POP80,POP(I),RI
         ELSE
            WRITE(*,'(A,I6,A,F10.3)') 
     *        ' FOR YEAR',IFIX(YEAR),' POP ESTMATE IS ',POP80
         ENDIF
   75 CONTINUE
C
      R=SQRT(RSQ)
      WRITE(*,*)'SQUARE ROOT OF RESIDUAL SUM OF SQUARES IS: ',R
C
      END