C
C************************************************************************
C************************************************************************
C*                                                                     **
C*  BLACKBOX.FOR          23 MAY 1998                                  **
C*                                                                     **
C*  THE BLACK BOX -- PERFORMS SINGULAR VALUE DECOMPOSITION ON A MATRIX **
C*    WITH MISSING ENTRIES.  DEVELOPED BY KEITH T. POOLE, JULY-OCTOBER **
C*    1982.  FURTHER DEVELOPMENT IN NOVEMBER-DECEMBER 1992.            **
C*                                                                     **
C*                                                                     **
C************************************************************************
C************************************************************************
C
      DIMENSION XBIGONE(3539,60),KISSUE(100),KID(3539),W(152,127),
     C          KMISS(60,10),XDATA(3539,125),XT(3539,25),RSAVE(25,25),
     C          YHAT(3539),UUU(3539,25),VVV(3539,25),LID(3539),
     C          DC(152)
      COMMON /WORKSP/ RWKSP
      REAL RWKSP(61913)
      integer*2 itim1,itim2,itim3,itim4,jtim1,jtim2,jtim3,jtim4
      CHARACTER*120 FNAME,FMT,TITLE
      CHARACTER*10 CAND(20)
      OPEN(4,FILE='BLACKBOX.DAT',STATUS='OLD')
      OPEN(23,FILE='BLACK23.DAT')
      OPEN(24,FILE='BLACK24.DAT')
      OPEN(28,FILE='BLACK28.DAT')
C
C
  70  FORMAT(' ITERATION RECORD'
     C/3X,'DIM',5X,'ERROR',4X,'EXPLAINED',
     C6X,'PERCENT',2X,'CUM PERCENT',3X,'R-SQUARE',2X,'STD ERR EST')
  113 FORMAT(1X,I5,8F12.4)
  195 FORMAT(' ITERATION SUMMARY')
  196 FORMAT(' NUMBER OF DIMENSIONS='I4)
  197 FORMAT(5X,'PSI-TRANSPOSE*PSI')
  198 FORMAT(I6,20F7.3)
  199 FORMAT(1X,'CONSTRAINT CHECKS ON PSI AND W'/
     C       5X,'SUM OF COLUMNS OF PSI'/
     C       6X,35F10.4)
  200 FORMAT(I6,35F10.4)
  201 FORMAT(2I7,35F7.3)
  202 FORMAT(5X,'W-TRANSPOSE*W')
  225 FORMAT(1X,A10,I7,35F7.3)
  301 FORMAT(' NUMBER OF CASES',3I6)
  315 FORMAT(' R-SQUARE ',I8,F7.3)
  330 FORMAT(I5,7F10.4)
  337 FORMAT(' RANK CHECK OF PSI*W',I5)
 1000 format(22x,4(i2.2,'.'))
 1001 format(' ELAPSED TIME OF JOB  ',4(i2.2,'.'))
 1002 FORMAT(8I5)
 1003 FORMAT(A10)
 1082 format(A120)
 2222 FORMAT(1X,78('*'))
C
C
      CALL IWKIN(61913)
C
C
      call gettim(itim1,itim2,itim3,itim4)
C
C
      read(4,1082)FNAME 
      write(23,1082)FNAME
      write(*,1082)FNAME
      open(22,file=FNAME,STATUS='OLD')
      READ(4,'(A)')TITLE
      WRITE(23,'(A)')TITLE
      WRITE(*,'(A)')TITLE
      READ(4,1002)NS,NY,NMISS,KVMIN
      WRITE(23,1002)NS,NY,NMISS,KVMIN
      WRITE(*,1002)NS,NY,NMISS,KVMIN
      READ(4,'(A)')FMT
      WRITE(23,'(A)')FMT
      WRITE(*,'(A)')FMT
      DO 111 I=1,NY
      READ(4,1003)CAND(I)
      WRITE(23,1003)CAND(I)
      WRITE(*,1003)CAND(I)
      READ(4,1002)(KMISS(I,J),J=1,NMISS)
      WRITE(23,1002)(KMISS(I,J),J=1,NMISS)
      WRITE(*,1002)(KMISS(I,J),J=1,NMISS)
  111 CONTINUE
C      NY=14
C      NS=3
C
      I=0
      II=0
  500 READ(22,FMT,END=400)MID,(KISSUE(J),J=1,NY)
      I=I+1
      IMARK=0
      DO 1 J=1,NY
      DO 112 K=1,NMISS
      IF(KISSUE(J).EQ.KMISS(J,K))THEN 
         KISSUE(J)=-999
         IMARK=IMARK+1
         GO TO 1
      ENDIF
  112 CONTINUE
  1   CONTINUE
      IF((NY-IMARK).LT.KVMIN)GO TO 500
      II=II+1
      KID(II)=MID
      DO 11 J=1,NY
      XBIGONE(II,J)=FLOAT(KISSUE(J))
  11  CONTINUE
      GO TO 500
  400 WRITE(*,301)I,II
      NOBS=II
      WRITE(23,301)NOBS
C
C  SUBROUTINE BLACKB--GENERATES MULTI-DIMENSIONAL STARTING VALUES 
C    FOR LEGISLATORS BY PERFORMING A PSEUDO-SINGULAR VALUE 
C    DECOMPOSITION OF NY ROLL CALLS
C
C
      WRITE(23,2222)
      WRITE(23,2222)
      DO 999 KKK=1,NS
      WRITE(23,2222)
      WRITE(23,196)KKK
      WRITE(*,196)KKK
      call BLACKB(NOBS,NY,KKK,XBIGONE,XDATA,W,DC,SVSUM)
C  
C  CONSTRAINT CHECKS
C
      DO 233 K=1,KKK
      SUM=0.0
      DO 234 I=1,NOBS
      SUM=SUM+XDATA(I,K)
  234 CONTINUE
      VVV(1,K)=SUM/FLOAT(NOBS)
  233 CONTINUE
      WRITE(23,199)(VVV(1,K),K=1,KKK)
      DO 227 J=1,KKK
      DO 227 K=1,KKK
      SUM=0.0
      DO 228 I=1,NOBS
      SUM=SUM+XDATA(I,J)*XDATA(I,K)
  228 CONTINUE
  227 VVV(J,K)=SUM
      WRITE(23,197)
      DO 229 I=1,KKK
  229 WRITE(23,200)I,(VVV(I,J),J=1,KKK)
      WRITE(23,202)
      DO 230 J=1,KKK
      DO 230 K=1,KKK
      SUM=0.0
      DO 231 I=1,NY
      SUM=SUM+W(I,J+1)*W(I,K+1)
  231 CONTINUE
  230 VVV(J,K)=SUM
      DO 232 I=1,KKK
  232 WRITE(23,200)I,(VVV(I,J),J=1,KKK)
      ASUM=0.0
      BSUM=0.0
      CSUM=0.0
      DSUM=0.0
      ESUM=0.0
      SUME=0.0
      KK=0
      DO 223 J=1,NY
      AASUM=0.0
      BBSUM=0.0
      CCSUM=0.0
      DDSUM=0.0
      EESUM=0.0
      KJJ=0
      DO 224 I=1,NOBS
      SUM=0.0
      DO 726 K=1,KKK
      SUM=SUM+XDATA(I,K)*W(J,K+1)
  726 CONTINUE
      XT(I,J)=SUM
      AA=SUM+W(J,1)
      IF(XBIGONE(I,J).EQ.-999.0)GO TO 224
      BB=XBIGONE(I,J)
      SUME=SUME+(AA-BB)**2
      AASUM=AASUM+AA
      BBSUM=BBSUM+BB
      CCSUM=CCSUM+AA*AA
      DDSUM=DDSUM+BB*BB
      EESUM=EESUM+AA*BB
      KJJ=KJJ+1
  224 CONTINUE
      AAA=FLOAT(KJJ)*EESUM-AASUM*BBSUM
      BBB=FLOAT(KJJ)*CCSUM-AASUM*AASUM
      CCC=FLOAT(KJJ)*DDSUM-BBSUM*BBSUM
      RRR=0.0
      IF(ABS(BBB*CCC).GT.0.0)RRR=(AAA*AAA)/(BBB*CCC)
      W(J,KKK+2)=RRR
      LID(J)=KJJ
      ASUM=ASUM+AASUM
      BSUM=BSUM+BBSUM
      CSUM=CSUM+CCSUM
      DSUM=DSUM+DDSUM
      ESUM=ESUM+EESUM
      KK=KK+KJJ
  223 CONTINUE
      AA=FLOAT(KK)*ESUM-ASUM*BSUM
      BB=FLOAT(KK)*CSUM-ASUM*ASUM
      CC=FLOAT(KK)*DSUM-BSUM*BSUM
      RRR=(AA*AA)/(BB*CC)
      WRITE(*,315)KK,RRR
      WRITE(23,315)KK,RRR
C
      RSAVE(KKK,1)=SUME
      RSAVE(KKK,2)=RRR
      RSAVE(KKK,3)=SQRT(SUME/FLOAT(KK - KKK*(NOBS+NY) + NY))
C
C  *****patched in code*****
C
      XTOL=.001
      CALL LSVRR(NOBS,NY,XT,3539,21,XTOL,IRANK,YHAT,UUU,
     C           3539,VVV,3539)
      WRITE(23,337)IRANK
      DO 464 I=1,NS+3
      WRITE(23,330)I,YHAT(I)
  464 CONTINUE
C
C
C  WRITE OUT COLUMN PARAMETERS
C
      DO 222 I=1,NY
      WRITE(28,225)CAND(I),LID(I),(W(I,JJ),JJ=1,KKK+2)
  222 CONTINUE
C
C  WRITE OUT ROW PARAMETERS
C
      DO 742 I=1,NOBS
      WRITE(24,201)I,KID(I),(XDATA(I,J),J=1,KKK)
C      WRITE(*,201)I,KID(I), (UUU(I,J),J=1,KKK)
  742 CONTINUE
  999 CONTINUE
      WRITE(23,2222)
      WRITE(23,70)
      WRITE(*,70)
      DO 69 J=1,NS
      AA=(SVSUM-RSAVE(J,1))/SVSUM
      CCC=SVSUM-RSAVE(J,1)
      IF(J.EQ.1)BB=SVSUM
      IF(J.GT.1)BB=RSAVE(J-1,1)
      SUM=(BB-RSAVE(J,1))/SVSUM
      WRITE(23,113)J,RSAVE(J,1),CCC,SUM*100.0,AA*100.0,RSAVE(J,2),
     C              RSAVE(J,3)
  69  WRITE(* ,113)J,RSAVE(J,1),CCC,SUM*100.0,AA*100.0,RSAVE(J,2),
     C              RSAVE(J,3)
      WRITE(23,2222)
C
      WRITE(23,2222)
      call gettim(jtim1,jtim2,jtim3,jtim4)
      write(*,1000)itim1,itim2,itim3,itim4
      write(*,1001)jtim1,jtim2,jtim3,jtim4
      STOP
      END
C
C  ***********************************************************************
C    SUBROUTINE BLACKB
C      
C      PERFORMS THE FOLLOWING DECOMPOSITION:
C
C           X = [PW' + Jc'] + E
C         
C          WHERE X IS NPxNY, P IS NPxNF, W IS NYxNF, J IS A NP LENGTH
C            VECTORS OF 1'S, c IS A NY VECTOR OF CONSTANTS,
C            AND E IS NPxNY MATRIX OF ERROR WHICH IS MINIMIZED
C
C  ***********************************************************************
C  
C
      SUBROUTINE BLACKB(NP,NY,NF,XBIGONE,XDATA,W,DC,SVSUM)
      DIMENSION X(3539,152),XX(408),FAL1(152),FAL2(152),
     C          PSIX(3539,152),XS(3539,152),LL(152),W(152,127),
     C          SQU(152),R(152,152),D(152),WK(500),ROOTC(3539,152),
     C          CROOT(3539,152),DD(152),DC(152),CC(152),XT(3539,152),
     C          XSS(3539,152),MPOS(152),TSUM(408),DDD(152),
     C          XBIGONE(3539,60),XDATA(3539,125),DX(152),
     C          UUU(3539,152),VVV(3539,152),WSAVE(127,127),
     C          XSAVE(3539,60),YHAT(3539)
      COMMON /AREA1/ JR,LR
      COMMON /XSREG/ AREG,BREG
      COMMON /WORKSP/ RWKSP
      REAL RWKSP(61913)
  23  FORMAT(' CORRELATION MATRIX FOR PSIX')
  70  FORMAT(' ITERATION RECORD'/3X,'DIM',5X,'ERROR',5X,'EXPLAINED',
     C6X,'PERCENT',2X,'CUM PERCENT',3X,'R-SQUARE')
  86  FORMAT(1X,'SINGULAR VALUES OF ESTIMATED MATRICES'/
     C       1X,'FIRST COLUMN:  ORIGINAL MATRIX WITH FILLED IN MISSING E
     CNTRIES'/
     C       1X,'SECOND COLUMN: REPRODUCED MATRIX -- PSI*W + Jc'/
     C       1X,'THIRD COLUMN:  ORIGINAL MATRIX WITH FILLED IN MISSING E
     CNTRIES MINUS THE ORIGINAL COLUMN MEANS'/
     C       1X,'FOURTH COLUMN:  PSI*W')
  113 FORMAT(1X,I5,8F12.4)
  199 FORMAT(' CONSTRAINT CHECKS ON PSI AND W'/6X,35F7.3)
  200 FORMAT(I6,35F10.4)
  300 FORMAT(1X,'NUMBER OF ROWS               =',I6/
     C       1X,'NUMBER OF COLUMNS            =',I6/
     C       1X,'TOTAL NUMBER OF DATA ENTRIES =',I6/
     C       1X,'NUMBER MISSING ENTRIES       =',I6/
     C       1X,'PERCENT MISSING DATA         =',F15.5/
     C       1X,'SUM OF SQUARES GRAND MEAN    =',F15.5)
  301 FORMAT(2I5,F7.3)
  335 FORMAT(' CORRECTION ROW=',I3,1X,' # POSITIVE=',I5)
  336 FORMAT(1X,25I3/1X,25I3)
  337 FORMAT(' RANK CHECK',I5)
  338 FORMAT(I5,7F10.4)
  521 FORMAT(' CROSS PRODUCT MATRIX FOR PSIX')
 1113 FORMAT(I3,12F10.3)
 2222 FORMAT(1X,78('*'))
 7200 FORMAT(' ERROR PARAMETER EIGENVALUE ROUTINE=',I4)
 9011 FORMAT(I5,10X,3F10.4)
 9012 FORMAT(3I5,3F10.4)
      JR=23
      LR=JR
      NF1=NF+1
      DO 349 J=1,NY
      DC(J)=0.0
      CC(J)=0.0
  349 LL(J)=0.0
      KTOT=0
      SVSUM=0.0
      SWSUM=0.0
      WRITE(23,2222)
C
C
      DO 600 I=1,NP
      SUM=0.0
      SUMA=0.0
      KK=0
      DO 30 J=1,NY
      X(I,J)=XBIGONE(I,J)
      IF(X(I,J).NE.-999.0)THEN
         SUM=SUM+X(I,J)*X(I,J)
         SUMA=SUMA+X(I,J)
         KK=KK+1
      ENDIF
      XS(I,J)=X(I,J)
  30  CONTINUE
      KTOT=KTOT+KK
      SVSUM=SVSUM+SUM
      SWSUM=SWSUM+SUMA
      DO 31 J=1,NY
      XSS(I,J)=X(I,J)
      IF(X(I,J).EQ.-999.0)GO TO 31
      LL(J)=LL(J)+1
      DC(J)=DC(J)+X(I,J)
      CC(J)=CC(J)+X(I,J)*X(I,J)
  31  CONTINUE
  600 CONTINUE
C
C   COMPUTE NY COLUMN MEANS PLACE IN DC(NY), COMPUTE SUM OF SQUARES
C     AND PERCENT MISSING DATA
C
      DO 32 J=1,NY
      D(J)=CC(J)-(DC(J)*DC(J))/FLOAT(LL(J))
      DC(J)=DC(J)/FLOAT(LL(J))
C      WRITE(*,301)J,LL(J),DC(J)
C      WRITE(23,301)J,LL(J),DC(J)
  32  CONTINUE
      SVSUM=SVSUM-(SWSUM**2)/FLOAT(KTOT)
      LTOT=NP*NY
      PMISS=FLOAT(LTOT-KTOT)/FLOAT(LTOT)
      LTOT=LTOT-KTOT
      IF(NF.EQ.1)THEN
         WRITE(23,300)NP,NY,KTOT,LTOT,PMISS*100.0,SVSUM
         WRITE(* ,300)NP,NY,KTOT,LTOT,PMISS*100.0,SVSUM
         WRITE(23,2222)
         WRITE(*,2222)
      ENDIF   
C
C  COMPUTE CORRELATION NY BY NY CORRELATION MATRIX AND VECTOR OF SIGN
C    CHANGES FOR COLUMNS OF X/XS/XSS
C
      call CORR2(NP,NY,X,R,LL,MPOS,KS,KPOS,1)
C
C  MASTER LOOP FOR NF DIMENSIONS--ONE DIMENSION IS EXTRACTED AT A TIME.
C     THE RESIDUAL MATRIX IS COMPUTED AFTER EACH DIMENSION IS EXTRACTED
C     AND IT IS USED AS X IN THE FORMULA BELOW.  AT THE END OF THIS
C     PROCESS:
C
C              P'P = I, WITH P BEING NOBSxNF, W NYxNF, c NYx1
C
C           X = [PW' + Jc'] + E
C
      DO 9999 JJJ=1,NF
      KKK=JJJ
      DO 10 J=1,NY
      XX(J)=1.0
      XX(J+NY)=-DC(J)
      IF(JJJ.GT.1)XX(J+NY)=0.0
  10  CONTINUE
      XXK=0.0
      TXB=0.0
      KKT=0
C
C  AT FIRST ITERATION JJJ = 1 THE COLUMN MEANS ARE SUBTRACTED
C     FOR SUBSEQUENT ITERATIONS, COLUMN MEANS = 0.  NOTE USE OF
C     COLUMN SIGN CHANGES (LL(NY)) FROM CORR2
C
      DO 3 I=1,NP
      KK=0
      SUM=0.0
      DO 1 J=1,NY
      IF(X(I,J).EQ.-999.0)GO TO 1
      SUM=SUM+(X(I,J)*XX(J)+XX(J+NY))*FLOAT(LL(J))
      KK=KK+1
  1   CONTINUE
C
C  XT(NOBS,1) = (X - COL. MEAN)xSIGN CHANGE, SUMMED OVER ITH ROW--THIS
C     IS THE STARTING ESTIMATE OF THE FIRST COLUMN OF P (PSIX BELOW).
C     THE MEAN OF XT(NOBS,1) IS THEN SUBTRACTED OUT BELOW SO THAT THE
C     STARTING ESTIMATE SUMS TO ZERO.
C     XT(NOBS,2) = 1 FOR THE LINEAR REGRESSIONS BELOW---TO PICK UP THE
C       INTERCEPT VECTOR, c
C
      WXB=SUM/(FLOAT(KK))
C
      KKT=KKT+1
      TXB=TXB+WXB
      XT(I,1)=WXB
      XXK=XXK+XT(I,1)*XT(I,1)
  3   CONTINUE
      TXB=TXB/FLOAT(KKT)
      XXK=XXK-FLOAT(KKT)*TXB*TXB
      DO 40 I=1,NP
      XT(I,1)=(XT(I,1)-TXB)
      XT(I,2)=1.0
  40  CONTINUE
C
C  ALTERNATING LEAST SQUARES IS NOW PERFORMED TO FIRST ESTIMATE W AND c--
C    SUBROUTINE REG--ESTIMATES OF W AND c ARE RETURNED IN W(NY,NF+1)
C  W AND c ARE THEN USED TO ESTIMATE P IN SUBROUTINE REG2--
C  RESIDUALS ARE RETURNED IN X(NOBS,NY)
C
      DO 444 MM=1,4
      call REG(NP,1,NY,TSUM,W,XS,X,XT,1,1,KKK)
      call REG2(NP,1,NY,W,XS,X,XT,SQU,PXB,PXS,KKK)
C
C  PXB IS THE MEAN OF THE ESTIMATED COLUMN OF P
C  PXS IS THE VARIANCE OF THE ESTIMATED COLUMN OF P
C
      QXS=PXS
      XCOR=SQRT(XXK/PXS)
C      WRITE(JR,113)MM,XXK,QXS,PXS,XCOR,PXB,TXB
      DO 444 I=1,NP
C
C  THE ESTIMATED MEAN, PXB, IS SUBTRACTED OFF AND THE COLUMN IS
C    ADJUSTED SO IT HAS THE ORIGINAL SUM OF SQUARES WITH MEAN ZERO
C
      PSIX(I,JJJ)=(XT(I,1)-PXB)*XCOR
      XT(I,1)=(XT(I,1)-PXB)*XCOR
  444 CONTINUE
C
C  END OF A.L.S. LOOP
C
C  IF THIS IS THE LAST ITERATION, I.E. JJJ = NF, THEN PUT A VECTOR
C    OF ONES INTO THE NF+1 COLUMN OF THE ESTIMATED P MATRIX
C    AND RESET X(NOBS,NY)--WHICH STARTED AS THE ORIGINAL DATA AND
C    THEN BECAME THE RESIDUAL MATRIX AT EACH STAGE--BACK TO THE
C    ORIGINAL DATA WHICH WAS SAVED IN XSS(NOBS,NY)
C
C  IF THIS IS NOT THE LAST ITERATION, PUT THE RESIDUALS X(NOBS,NY) INTO
C    XS(NOBS,NY), AND CALL CORR2 TO GET THE SIGN CHANGES FOR THE COLUMNS
C    WHICH ARE USED ABOVE
C
      DO 461 I=1,NP
      IF(KKK.EQ.NF)PSIX(I,NF+1)=1.0
      DO 461 J=1,NY
      IF(KKK.EQ.NF)X(I,J)=XSS(I,J)
  461 XS(I,J)=X(I,J)
      IF(KKK.EQ.NF)GO TO 9999
      call CORR2(NP,NY,X,R,LL,MPOS,KS,KPOS,1)
 9999 CONTINUE
C
C  WRITE OUT CORRELATION AND CROSS PRODUCT MATRIX OF P
C
C      WRITE(JR,23)
C      call CORR2(NP,NF,PSIX,R,LL,MPOS,KS,KPOS,1)
C      WRITE(LR,521)
C      call PSIPRM(NP,NF1,PSIX)
C
C  PERFORM A.L.S. WITH FULL P MATRIX UNTIL CONVERGENCE
C
      DO 222 NN=1,5
      call REG(NP,NF,NY,TSUM,W,XS,X,PSIX,1,1,NF)
      call REG2(NP,NF,NY,W,XS,X,PSIX,SQU,PXB,PXS,NF)
C 
C  CENTER THE PSI MATRIX AT ZERO
C
      DO 42 K=1,NF
      SUM=0.0
      DO 41 I=1,NP
      SUM=SUM+PSIX(I,K)
  41  CONTINUE
      SUM=SUM/FLOAT(NP)
      DO 43 I=1,NP
      PSIX(I,K)=PSIX(I,K)-SUM
  43  CONTINUE
  42  CONTINUE
      AA=ABS(AREG-BREG)
      IF(AA.LT..01)GO TO 223
  222 CONTINUE
  223 NFPS=NF+1
C
C  COMPUTE PW' + Jc' AND PUT INTO X(NOBS,NY)
C  FOR MISSING DATA, INSERT ESTIMATED VALUES IN X INTO THE 
C    ORIGINAL DATA MATRIX USING XT(NOBS,NY) SO THAT SINGULAR VALUES
C    CAN BE ESTRACTED FROM ORIGINAL DATA
C
      DO 60 I=1,NP
      DO 601 K=1,NY
  601 XT(I,K)=0.0
      DO 61 K=1,NY
      WK(K)=0.0
      WK(K+NY)=0.0
      SUM=0.0
      DO 62 J=1,NFPS
  62  SUM=SUM+PSIX(I,J)*W(K,J)
  61  X(I,K)=SUM
      DO 59 JJ=1,NY
      IF(XS(I,JJ).NE.-999.0)XT(I,JJ)=XS(I,JJ)
      IF(XS(I,JJ).EQ.-999.0)XT(I,JJ)=X(I,JJ)
  59  CONTINUE
  60  CONTINUE
C
C  EXTRACT SINGULAR VALUES FROM ORIGINAL DATA MATRIX WITH ESTIMATES
C   INSERTED FOR MISSING DATA
C
      DO 71 I=1,NY
      DO 71 J=1,NY
      SUM=0.0
      DO 73 K=1,NP
  73  SUM=SUM+XT(K,I)*XT(K,J)
  71  ROOTC(I,J)=SUM
      CALL RS(3539,NY,ROOTC,D,1,CROOT,FAL1,FAL2,IERR)
      DO 729 I=1,NP
      DO 729 J=1,NY
      ROOTC(I,J)=XT(I,J)
      CROOT(I,J)=XT(I,J)-DC(J)
  729 CONTINUE
      XTOL=.001
      CALL LSVRR(NP,NY,ROOTC,3539,21,XTOL,IRANK,YHAT,UUU,
     C           3539,VVV,3539)
C      WRITE(23,337)IRANK
      DO 464 I=1,NY
C      WRITE(23,338)I,YHAT(I)
      D(I)=YHAT(I)
  464 CONTINUE
      XTOL=.001
      CALL LSVRR(NP,NY,CROOT,3539,21,XTOL,IRANK,YHAT,UUU,
     C           3539,VVV,3539)
C      WRITE(23,337)IRANK
      DO 465 I=1,NY
C      WRITE(23,338)I,YHAT(I)
      DX(I)=YHAT(I)
  465 CONTINUE
C      WRITE(JR,7200)IERR
C      WRITE(*,7200)IERR
      ESUM=0.0
      SUMM=0.0
      SUMM1=0.0
      SUMM2=0.0
      DO 83 I=1,NP
      DO 83 J=1,NY
      SUMM=SUMM+X(I,J)**2
      IF(XS(I,J).NE.-999.0)ESUM=ESUM+(X(I,J)-XS(I,J))**2
      IF(XS(I,J).NE.-999.0)SUMM1=SUMM1+X(I,J)**2
      IF(XS(I,J).EQ.-999.0)SUMM2=SUMM2+X(I,J)**2
C
C  STORE PW' + Jc' IN XSS(NOBS,NY) AND COMPUTE PW' AND STORE IN XT(NOBS,NY)
C
      XSS(I,J)=X(I,J)
  83  XT(I,J)=X(I,J)-W(J,NF+1)
C
C  PERFORM SINGULAR VALUE DECOMPOSITION OF PW' + Jc'
C
      DO 75 I=1,NY
      DO 75 J=1,NY
      SUM=0.0
      DO 76 K=1,NP
      SUM=SUM+XSS(K,I)*XSS(K,J)
  76  CONTINUE
  75  ROOTC(I,J)=SUM
      CALL RS(3539,NY,ROOTC,DDD,1,CROOT,FAL1,FAL2,IERR)
      DO 728 I=1,NP
      DO 728 J=1,NY
      ROOTC(I,J)=XSS(I,J)
  728 CONTINUE
      XTOL=.001
      CALL LSVRR(NP,NY,ROOTC,3539,21,XTOL,IRANK,YHAT,UUU,
     C           3539,VVV,3539)
C      WRITE(23,337)IRANK
      DO 463 I=1,NY
C      WRITE(23,338)I,YHAT(I)
      DDD(I)=YHAT(I)
  463 CONTINUE
C      WRITE(JR,7200)IERR
C      WRITE(*,7200)IERR
C
C  PERFORM SINGULAR VALUE DECOMPOSITION OF PW'
C
      DO 77 I=1,NY
      DO 77 J=1,NY
      SUM=0.0
      DO 78 K=1,NP
      SUM=SUM+XT(K,I)*XT(K,J)
  78  CONTINUE
  77  ROOTC(I,J)=SUM
      CALL RS(3539,NY,ROOTC,DD,1,CROOT,FAL1,FAL2,IERR)
      DO 727 I=1,NP
      DO 727 J=1,NY
      ROOTC(I,J)=XT(I,J)
  727 CONTINUE
      XTOL=.001
      CALL LSVRR(NP,NY,ROOTC,3539,21,XTOL,IRANK,YHAT,UUU,
     C           3539,VVV,3539)
C      WRITE(23,337)IRANK
      DO 462 I=1,NP
      DD(I)=YHAT(I)
      DO 462 JJ=1,NF
      PSIX(I,JJ)=UUU(I,JJ)*SQRT(YHAT(JJ))
  462 CONTINUE
      DO 460 I=1,NY
      DO 460 JJ=1,NF
      W(I,JJ)=VVV(I,JJ)*SQRT(YHAT(JJ))
  460 CONTINUE
C      WRITE(JR,7200)IERR
C      WRITE(*,7200)IERR
C  
C  CONSTRAINT CHECKS
C
      DO 233 K=1,NFPS
      SUM=0.0
      DO 234 I=1,NP
      SUM=SUM+PSIX(I,K)
  234 CONTINUE
      VVV(1,K)=SUM/FLOAT(NP)
  233 CONTINUE
C      WRITE(23,199)(VVV(1,K),K=1,NF)
      DO 227 J=1,NF
      DO 227 K=1,NF
      SUM=0.0
      DO 228 I=1,NP
      SUM=SUM+PSIX(I,J)*PSIX(I,K)
  228 CONTINUE
  227 VVV(J,K)=SUM
      DO 229 I=1,NF
C      WRITE(23,200)I,(VVV(I,J),J=1,NF)
  229 CONTINUE
      DO 230 J=1,NF
      DO 230 K=1,NF
      SUM=0.0
      DO 231 I=1,NY
      SUM=SUM+W(I,J)*W(I,K)
  231 CONTINUE
  230 VVV(J,K)=SUM
      DO 232 I=1,NF
C      WRITE(23,200)I,(VVV(I,J),J=1,NF)
  232 CONTINUE
C
C  SAVE PSI AND W HERE
C
      DO 266 J=1,NF
      DO 266 I=1,NP
      XSAVE(I,J)=PSIX(I,J)
  266 CONTINUE
      DO 265 J=1,NF
      DO 265 I=1,NY
      WSAVE(I,J+1)=W(I,J)
      WSAVE(I,1)=W(I,NF+1)
  265 CONTINUE        
C
C
C  WRITE OUT SINGULAR VALUES FOR ORIGINAL DATA (USING ESTIMATES FOR
C   MISSING ENTRIES), PW' + Jc', AND PW'
C
      WRITE(JR,86)
      DO 87 J=1,NY
      IF(J.GT.NF+3)GO TO 87
      WRITE(23,1113)J,D(J),DDD(J),DX(J),DD(J)
      WRITE(*,1113)J,D(J),DDD(J),DX(J),DD(J)
  87  CONTINUE
C
C
      WRITE(JR,2222)
C
C  RESTORE PSI AND W HERE
C
      DO 264 J=1,NF
      DO 264 I=1,NP
      XDATA(I,J)=XSAVE(I,J)
  264 CONTINUE
      DO 263 J=1,NF+1
      DO 263 I=1,NY
      W(I,J)=WSAVE(I,J)
  263 CONTINUE        
C
      RETURN
      END
C
C  **********************************************************************
C    SUBROUTINE CORR2---CALLED BY BLACKB.  COMPUTES CORRELATION MATRIX
C     FOR INPUT MATRIX X AND COMPUTES THE VECTOR OF COLUMN SIGN CHANGES
C     WHICH IS STORED IN VECTOR LL.
C  **********************************************************************
C
      SUBROUTINE CORR2(NP,NY,X,R,LL,MPOS,KS,KPOS,IPRNT)
      DIMENSION X(3539,152),R(152,152),SA(152,152),SB(152,152),
     CSD(152,152),LL(152),MPOS(152),SC(152,152)
      COMMON /AREA1/ JR,LR
  200 FORMAT(1X,50F7.3)
      DO 1 J=1,NY
      DO 1 JJ=1,NY
      SA(J,JJ)=0.0
      SB(J,JJ)=0.0
      SC(J,JJ)=0.0
  1   SD(J,JJ)=0.0
      DO 40 I=1,NP
      DO 31 J=1,NY
      DO 31 JJ=1,J
      IF(X(I,J).EQ.-999.0)GO TO 31
      IF(X(I,JJ).EQ.-999.0)GO TO 31
      SA(J,JJ)=SA(J,JJ)+X(I,J)
      IF(J.NE.JJ)SA(JJ,J)=SA(JJ,J)+X(I,JJ)
      SB(J,JJ)=SB(J,JJ)+X(I,J)*X(I,J)
      IF(J.NE.JJ)SB(JJ,J)=SB(JJ,J)+X(I,JJ)*X(I,JJ)
      SC(J,JJ)=SC(J,JJ)+X(I,J)*X(I,JJ)
      SC(JJ,J)=SC(J,JJ)
      SD(J,JJ)=SD(J,JJ)+1.0
  31  CONTINUE
  40  CONTINUE
      DO 32 J=1,NY
      DO 32 JJ=1,J
      AA=SD(J,JJ)*SC(J,JJ)-SA(J,JJ)*SA(JJ,J)
      BB=SD(J,JJ)*SB(J,JJ)-SA(J,JJ)*SA(J,JJ)
      CC=SD(J,JJ)*SB(JJ,J)-SA(JJ,J)*SA(JJ,J)
      BBCC=BB*CC
      IF(BBCC.LE..0)R(JJ,J)=0.0
      IF(BBCC.LE..0)GO TO 343
      R(JJ,J)=AA/SQRT(BB*CC)
  343 CONTINUE
  32  R(J,JJ)=R(JJ,J)
      IF(IPRNT.EQ.1)GO TO 62
      DO 34 J=1,NY
  34  WRITE(JR,200)(R(J,JJ),JJ=1,NY)
  62  BB=-99.0
C
C  FIND ROW WITH LARGEST TOTAL SUM OF ABSOLUTE VALUED CORRELATIONS
C
      KS=0
      DO 50 J=1,NY
      SUM=0.0
      DO 51 JJ=1,NY
  51  SUM=SUM+ABS(R(J,JJ))
      BB=AMAX1(SUM,BB)
      IF(BB.EQ.SUM)KS=J
  50  CONTINUE
      DO 52 J=1,NY
      IF(R(KS,J).LE.0.0)LL(J)=-1
  52  IF(R(KS,J).GT.0.0)LL(J)=1
C
C  ITERATIVELY CHANGE SIGNS OF COLUMNS TO MAXIMIZE NUMBER OF ENTRIES
C   IN THE CORRELATION MATRIX WHICH ARE POSITIVE.  THIS VECTOR IS THEN
C   RETURNED AS THE SIGN CHANGES FOR THE DATA MATRIX X.
C
      NYD2=(NY-1)/2
      KPOS=0
      DO 60 JK=1,NY
      DO 60 J=1,NY
      KK=0
      KSUM=0
      DO 61 JJ=1,NY
      AA=R(J,JJ)*FLOAT(LL(JJ))*FLOAT(LL(J))
      IF(JK.EQ.NY.AND.AA.GE.0.0)KPOS=KPOS+1
      IF(JK.EQ.NY.AND.AA.GE.0.0)KSUM=KSUM+1
  61  IF(AA.LT.0.0)KK=KK+1
      IF(KK.GT.NYD2)LL(J)=LL(J)*(-1)
      IF(JK.EQ.NY)MPOS(J)=KSUM
  60  IF(KK.GT.NYD2)KS=999
      RETURN
      END
C
C  ************************************************************************
C    SUBROUTINE REG---CALLED BY BLACKB.  PERFORMS REGRESSION TO ESTIMATE
C      W AND c.
C  ************************************************************************
C
      SUBROUTINE REG(NP,NF,NY,TSUM,W,XS,X,PSI,IPRNT,ILAST,KKK)
      DIMENSION XS(3539,152),PSI(3539,152),A(127,127),B(127,127),
     CW(152,127),C(152),TSUM(408),X(3539,152),LLL(152),FV1(127),
     CXT(3539,152),R(152,152),LL(152),RSUM(152),MPOS(152),WVEC(127),
     CZMAT(127,127),FV2(127)
      COMMON /AREA1/ JR,LR
      COMMON /XSREG/ AREG,BREG
  100 FORMAT(' PERFORMANCE INDEX MATRIX INVERSE=',2I5)
  101 FORMAT(' PERFORMANCE INDEX EIGENVALUE/VECTOR ROUTINE=',3I5)
  102 FORMAT(2I5,50F13.5)
  111 FORMAT(1X,8A1,13F10.4)
  200 FORMAT(3I5,15F10.4/15x,15f10.4/15x,15f10.4)
  250 FORMAT(' DIMENSION=',I3,1X,'TOTAL SSE REG1=',F15.4)
  300 FORMAT(' SSE AND FINAL WEIGHT MATRIX')
      NF1=NF+1
      DO 15 K=1,NY
      LLL(K)=0
      TSUM(K)=0.0
      C(K)=0.0
      DO 15 J=1,NF1
  15  W(K,J)=0.0
C
C  IF ILAST = 1 THEN ONLY PERFORM THE ESTIMATE OF W AND c FOR THE 
C   CURRENT COLUMN OF P.
C  IF ILAST = 0 THEN ESTIMATE W AND c ASSUMING NF=1, THEN ASSUMING
C   NF=2, UP TO NF=NF.  THIS ALLOWS THE COMPUTATION OF THE SS ACCOUNTED
C   FOR BY A 1 DIMENSION PW' + Jc' MODEL, A 2 DIMENSION PW' + Jc' MODEL,
C   AND SO ON.
C
      IF(ILAST.EQ.1)KKA=NF
      IF(ILAST.EQ.0)KKA=1
      DO 11 KK=KKA,NF
      DO 16 K=1,NY
      DO 16 J=1,NF+1
  16  W(K,J)=0.0
C
C  RUN REGRESSION TO ESTIMATE W AND c ONE COLUMN AT A TIME
C   BECAUSE OF MISSING DATA
C
      NF1=KK+1
      NF11=KK+2
      DO 7 K=1,NY
      LLL(K)=0
C
C  COMPUTE [P'P]-1
C
      DO 2 J=1,NF1
      DO 2 JJ=1,NF1
      SUM=0.0
      DO 1 I=1,NP
      IF(XS(I,K).EQ.-999.0)GO TO 1
      SUM=SUM+PSI(I,J)*PSI(I,JJ)
  1   CONTINUE
  2   A(J,JJ)=SUM
      call rs(127,nf1,a,wvec,1,ZMAT,fv1,fv2,ier)
C      IF(ILAST.EQ.0)THEN
C         WRITE(23,101)NP,NF1,IER
C         WRITE(23,102)NP,NF1,(WVEC(J),J=1,NF1)
C         WRITE(*,101)NP,NF1,IER
C         WRITE(*,102)NP,NF1,(WVEC(J),J=1,NF1)
C      ENDIF
      DO 60 I=1,NF1
      DO 60 KX=1,NF1
      SUM=0.0
      DO 61 J=1,NF1
      IF(ABS(WVEC(J)).GT..0001)THEN
          SUM=SUM+ZMAT(KX,J)*(1.0/WVEC(J))*ZMAT(I,J)
      ENDIF
  61  CONTINUE
  60  B(I,KX)=SUM
C
C  COMPUTE [P'P]-1*P'X = W'/c'
C
      DO 3 I=1,NP
      IF(XS(I,K).EQ.-999.0)GO TO 3
      DO 5 J=1,NF1
      SUM=0.0
      DO 4 JJ=1,NF1
  4   SUM=SUM+B(J,JJ)*PSI(I,JJ)
  5   C(J)=SUM
      DO 6 J=1,NF1
  6   W(K,J)=W(K,J)+C(J)*XS(I,K)
  3   CONTINUE
C
C  CALCULATE R-SQUARE AND SSE FOR REGRESSION
C
      ESUM=0.0
      DO 8 I=1,NP
      IF(XS(I,K).EQ.-999.0)GO TO 8
      LLL(K)=LLL(K)+1
      SUM=0.0
      DO 9 J=1,NF1
  9   SUM=SUM+PSI(I,J)*W(K,J)
      ESUM=ESUM+(SUM-XS(I,K))**2
      X(I,K)=SUM-XS(I,K)
      XT(I,1)=SUM
      XT(I,2)=XS(I,K)
  8   CONTINUE
      W(K,NF11)=ESUM
      TSUM(KK)=TSUM(KK)+ESUM
      call CORR2(NP,2,XT,R,LL,MPOS,KS,KPOS,1)
      RSUM(K)=R(1,2)*R(1,2)
  7   CONTINUE
C
C  END OF REGRESSION LOOP FOR CURRENT DIMENSION
C
      IF(ILAST.EQ.1)GO TO 90
C
C  IF ILAST = 0 THEN CALCULATE THE R-SQUARE FOR THE CURRENT
C   NUMBER OF DIMENSIONS
C
      DO 888 I=1,NP
      DO 88 K=1,NY
      SUM=0.0
      DO 89 J=1,NF1
  89  SUM=SUM+PSI(I,J)*W(K,J)
  88  XT(I,K)=SUM
  888 CONTINUE
      call RSQUR(NP,NY,RR,XT,XS)
      TSUM(KK+NY)=RR
  90  CONTINUE
      IF(IPRNT.EQ.1)GO TO 333
      WRITE(JR,300)
      DO 10 J=1,NY
  10  WRITE(LR,200)KK,J,LLL(J),W(J,NF11),RSUM(J),(W(J,JJ),JJ=1
     C,NF1)
  333 CONTINUE
      IF(ILAST.EQ.1)WRITE(JR,250)KKK,TSUM(KK)
      IF(ILAST.EQ.0)WRITE(JR,250)KK,TSUM(KK)
      AREG=TSUM(KK)
  11  CONTINUE
C
C  END OF DIMENSION ITERATION LOOP--ONLY IF ILAST = 0
C
      RETURN
      END
C
C  **************************************************************************
C    SUBROUTINE REG2---CALLED BY BLACKB.  PERFORMS REGRESSION TO ESTIMATE P.
C  **************************************************************************
C
      SUBROUTINE REG2(NP,NF,NY,W,XS,X,PSI,SQU,PXB,PXS,KKK)
      DIMENSION W(152,127),X(3539,152),XS(3539,152),PSI(3539,152),
     CA(152,127),V(152),SQU(152),Y(152)
      COMMON /AREA1/ JR,LR
      COMMON /XSREG/ AREG,BREG
      NF1=NF+1
      ESUM=0.0
      PXB=0.0
      PXS=0.0
      XNS=0.0
C
C  REGRESS ONE ROW AT A TIME BECAUSE OF MISSING DATA
C
      DO 1 I=1,NP
      KK=0
C
C  SET UP DEPENDENT VARIABLE-- Y(NY) = X[I,NY] - NY COL. MEANS
C    AND PUT W(NY,NF) INTO MATRIX A(NY,NF)
C
      DO 2 J=1,NY
      V(J)=0.0
      IF(XS(I,J).EQ.-999.0)GO TO 2
      KK=KK+1
      Y(KK)=XS(I,J)-W(J,NF1)
      DO 3 JJ=1,NF
  3   A(KK,JJ)=W(J,JJ)
  2   CONTINUE
C
C  CALL SUBROUTINE REGA TO PERFORM LEAST SQUARES
C
      NS=KK
      call REGA(NS,NF,A,Y,V)
C
C  STORE ESTIMATED ROW ENTRIES OF P AND COMPUTE SSE
C
      DO 5 K=1,NY
      SUM=0.0
      IF(XS(I,K).EQ.-999.0)GO TO 5
      DO 4 J=1,NF
      PSI(I,J)=V(J)
      SUM=SUM+PSI(I,J)*W(K,J)
  4   CONTINUE
      SUM=SUM+W(K,NF1)
      X(I,K)=SUM-XS(I,K)
      ESUM=ESUM+(SUM-XS(I,K))**2
  5   CONTINUE
      PXB=PXB+PSI(I,1)
      PXS=PXS+PSI(I,1)**2
      XNS=XNS+1.0
  1   CONTINUE
C
C  END OF REGRESSION LOOP
C
C  PXB IS THE SUM OF THE ESTIMATED COLUMN OF P AND
C  PXS IS THE VARIANCE OF THE ESTIMATED COLUMN OF P.
C  THESE ARE RETURNED TO BLACKB FOR ADJUSTMENT PURPOSES
C
      PXB=PXB/XNS
      PXS=PXS-XNS*PXB*PXB
      SQU(KKK)=ESUM
      WRITE(JR,350)KKK,ESUM
  350 FORMAT(' DIMENSION=',I3,1X,'TOTAL SSE REG2=',F15.4)
      BREG=ESUM
      RETURN
      END
C
C  **************************************************************************
C    SUBROUTINE REGA---CALLED BY REG2.  PERFORMS THE REGRESSION:
C
C         [W'W]-1*W'[X(I,NY) - c']
C
C      TO GET THE ROW ENTRIES OF P.
C  **************************************************************************
C
      SUBROUTINE REGA(NS,NF,A,Y,V)
      DIMENSION A(152,127),Y(152),B(127,127),C(127,127),V(152),
     CBB(127,152),ZMAT(127,127),FV1(127),FV2(127),WVEC(127)
      COMMON /AREA1/ JR,LR
  101 FORMAT(' PERFORMANCE INDEX EIGENVALUE/VECTOR ROUTINE=',3I5)
  102 FORMAT(2I5,30F10.4)
      DO 1 J=1,NF
      DO 1 JJ=1,NF
      SUM=0.0
      DO 2 I=1,NS
  2   SUM=SUM+A(I,J)*A(I,JJ)
  1   B(J,JJ)=SUM
      call rs(127,nf,B,wvec,1,ZMAT,fv1,fv2,ier)
C      WRITE(*,101)NS,NF,IER
C      WRITE(*,102)NS,NF,(WVEC(J),J=1,NF)
      DO 60 I=1,NF
      DO 60 K=1,NF
      SUM=0.0
      DO 61 J=1,NF
      IF(ABS(WVEC(J)).GT..0001)THEN
          SUM=SUM+ZMAT(K,J)*(1.0/WVEC(J))*ZMAT(I,J)
      ENDIF
  61  CONTINUE
  60  C(I,K)=SUM
      DO 3 I=1,NS
      DO 3 J=1,NF
      SUM=0.0
      DO 4 JJ=1,NF
  4   SUM=SUM+C(J,JJ)*A(I,JJ)
  3   BB(J,I)=SUM
      DO 5 JJ=1,NF
      SUM=0.0
      DO 6 J=1,NS
  6   SUM=SUM+BB(JJ,J)*Y(J)
  5   V(JJ)=SUM
      RETURN
      END
C
C  ***************************************************************************
C    SUBROUTINE RSQUR---CALLED BY REG.  COMPUTES THE SIMPLE PEARSON R-SQUARE
C     BETWEEN TWO INPUT MATRICES---A(NOBS,NY) AND B(NOBS,NY)
C  ***************************************************************************
C
      SUBROUTINE RSQUR(NP,NY,R,A,B)
      DIMENSION A(3539,152),B(3539,152)
      COMMON /AREA1/ JR,LR
      ASUM=0.0
      BSUM=0.0
      CSUM=0.0
      DSUM=0.0
      ESUM=0.0
      XNT=0.0
      DO 1 I=1,NP
      DO 1 J=1,NY
      IF(A(I,J).EQ.-999.0)GO TO 1
      IF(B(I,J).EQ.-999.0)GO TO 1
      ASUM=ASUM+A(I,J)
      BSUM=BSUM+B(I,J)
      CSUM=CSUM+A(I,J)*A(I,J)
      DSUM=DSUM+B(I,J)*B(I,J)
      ESUM=ESUM+A(I,J)*B(I,J)
      XNT=XNT+1.0
  1   CONTINUE
      AA=XNT*ESUM-ASUM*BSUM
      BB=XNT*CSUM-ASUM*ASUM
      CC=XNT*DSUM-BSUM*BSUM
      R=(AA*AA)/(BB*CC)
      RETURN
      END
C
C  ***************************************************************************
C    SUBROUTINE PSIPRM---CALLED BY BLACKB.  CALCULATES P'P AND WRITES IT OUT.
C  ***************************************************************************
C
      SUBROUTINE PSIPRM(NP,NF,PSI)
      DIMENSION PSI(3539,152),R(127,127)
      COMMON /AREA1/ JR,LR
  200 FORMAT(1X,50F10.4)
      DO 1 J=1,NF
      DO 1 K=1,NF
      SUM=0.0
      DO 2 I=1,NP
      SUM=SUM+PSI(I,J)*PSI(I,K)
  2   CONTINUE
  1   R(J,K)=SUM
      DO 3 I=1,NF
  3   WRITE(LR,200)(R(I,J),J=1,NF)
      RETURN
      END

C
C
C  ************************************************************************
C    SORT SUBROUTINE--SORTS A VECTOR 'A' OF REAL ELEMENTS INTO ASCENDING
C    ORDER.  'LA' IS THE NUMBER OF ELEMENTS TO BE SORTED AND 'IR' IS A
C    VECTOR OF INTEGERS THAT RECORDS THE PERMUTATIONS--USUALLY SET TO
C    1,2,3,4,...
C  ************************************************************************
C
C
      SUBROUTINE RSORT(A,LA,IR)
      DIMENSION A(LA),IU(21),IL(21),IR(LA)
      IF (LA.LE.0) RETURN
      M = 1
      I = 1
      J = LA
      R = .375
    5 IF (I.EQ.J) GO TO 45
      IF (R.GT..5898437) GO TO 10
      R = R+3.90625E-2
      GO TO 15 
   10 R = R-.21875 
   15 K = I 
C
C SELECT A CENTRAL ELEMENT OF THE  
C ARRAY AND SAVE IT IN LOCATION T  
C
      IJ = I+(J-I)*R 
      T = A(IJ)   
      IT = IR(IJ) 
C
C FIRST ELEMENT OF ARRAY IS GREATER
C THAN T, INTERCHANGE WITH T       
C
      IF (A(I).LE.T) GO TO 20  
      A(IJ) = A(I) 
      A(I) = T     
      T = A(IJ) 
      IR(IJ) = IR(I)  
      IR(I) = IT
      IT = IR(IJ) 
   20 L = J
C
C IF LAST ELEMENT OF ARRAY IS LESS THAN
C T, INTERCHANGE WITH T
C
      IF (A(J).GE.T) GO TO 30
      A(IJ) = A(J)
      A(J) = T
      T = A(IJ)
      IR(IJ) = IR(J)
      IR(J) = IT
      IT = IR(IJ)
C
C IF FIRST ELEMENT OF ARRAY IS GREATER
C THAN T, INTERCHANGE WITH T
C
      IF (A(I).LE.T) GO TO 30
      A(IJ) = A(I)
      A(I) = T
      T = A(IJ)
      IR(IJ) = IR(I)
      IR(I) = IT
      IT = IR(IJ)
      GO TO 30
   25 IF (A(L).EQ.A(K)) GO TO 30
      TT = A(L)
      A(L) = A(K)
      A(K) = TT
      ITT = IR(L)
      IR(L) = IR(K)
      IR(K) = ITT
C
C FIND AN ELEMENT IN THE SECOND HALF OF
C THE ARRAY WHICH IS SMALLER THAN T
C
   30 L = L-1
      IF (A(L).GT.T) GO TO 30
C
C FIND AN ELEMENT IN THE FIRST HALF OF
C THE ARRAY WHICH IS GREATER THAN T
C
   35 K = K+1
      IF (A(K).LT.T) GO TO 35
C
C INTERCHANGE THESE ELEMENTS
C
      IF (K.LE.L) GO TO 25
C
C SAVE UPPER AND LOWER SUBSCRIPTS OF
C THE ARRAY YET TO BE SORTED
C
      IF (L-I.LE.J-K) GO TO 40
      IL(M) = I
      IU(M) = L
      I = K
      M = M+1
      GO TO 50
   40 IL(M) = K
      IU(M) = J
      J = L
      M = M+1
      GO TO 50
C
C BEGIN AGAIN ON ANOTHER PORTION OF
C THE UNSORTED ARRAY
C
   45 M = M-1
      IF (M.EQ.0) RETURN
      I = IL(M)
      J = IU(M)
   50 IF (J-I.GE.11) GO TO 15
      IF (I.EQ.1) GO TO 5
      I = I-1
   55 I = I+1
      IF (I.EQ.J) GO TO 45
      T = A(I+1)
      IT = IR(I+1)
      IF (A(I).LE.T) GO TO 55
      K = I
   60 A(K+1) = A(K)
      IR(K+1) = IR(K)
      K = K-1
      IF (T.LT.A(K)) GO TO 60
      A(K+1) = T
      IR(K+1) = IT
      GO TO 55
      END
C
C
C  **************************************************************************
C    EIGENVECTOR/EIGENVALUE DECOMPOSITION SUBROUTINES FOR A SYMMETRIC MATRIX
C    SUBROUTINES ARE FROM EISPACK
C  **************************************************************************
C
      SUBROUTINE RS(NM,N,A,W,MATZ,Z,FV1,FV2,IERR)
C
      INTEGER N,NM,IERR,MATZ
      REAL A(NM,N),W(N),Z(NM,N),FV1(N),FV2(N)
C
C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
C     OF A REAL SYMMETRIC MATRIX.
C
C     ON INPUT
C
C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C        DIMENSION STATEMENT.
C
C        N  IS THE ORDER OF THE MATRIX  A.
C
C        A  CONTAINS THE REAL SYMMETRIC MATRIX.
C
C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
C
C     ON OUTPUT
C
C        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER.
C
C        Z  CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO.
C
C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT
C           AND TQL2.  THE NORMAL COMPLETION CODE IS ZERO.
C
C        FV1  AND  FV2  ARE TEMPORARY STORAGE ARRAYS.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
C
C     THIS VERSION DATED AUGUST 1983.
C
C     ------------------------------------------------------------------
C
      IF (N .LE. NM) GO TO 10
      IERR = 10 * N
      GO TO 50
C
   10 IF (MATZ .NE. 0) GO TO 20
C     .......... FIND EIGENVALUES ONLY ..........
      CALL  TRED1(NM,N,A,W,FV1,FV2)
      CALL  TQLRAT(N,W,FV2,IERR)
      GO TO 50
C     .......... FIND BOTH EIGENVALUES AND EIGENVECTORS ..........
   20 CALL  TRED2(NM,N,A,W,FV1,Z)
      CALL  TQL2(NM,N,W,FV1,Z,IERR)
   50 RETURN
      END
      SUBROUTINE TRED1(NM,N,A,D,E,E2)
C
      INTEGER I,J,K,L,N,II,NM,JP1
      REAL A(NM,N),D(N),E(N),E2(N)
      REAL F,G,H,SCALE
C
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED1,
C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C
C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX
C     TO A SYMMETRIC TRIDIAGONAL MATRIX USING
C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
C
C     ON INPUT
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.
C
C        N IS THE ORDER OF THE MATRIX.
C
C        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE
C          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
C
C     ON OUTPUT
C
C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS-
C          FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER
C          TRIANGLE.  THE FULL UPPER TRIANGLE OF A IS UNALTERED.
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX.
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
C
C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
C
C     THIS VERSION DATED AUGUST 1983.
C
C     ------------------------------------------------------------------
C
      DO 100 I = 1, N
         D(I) = A(N,I)
         A(N,I) = A(I,I)
  100 CONTINUE
C     .......... FOR I=N STEP -1 UNTIL 1 DO -- ..........
      DO 300 II = 1, N
         I = N + 1 - II
         L = I - 1
         H = 0.0E0
         SCALE = 0.0E0
         IF (L .LT. 1) GO TO 130
C     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
         DO 120 K = 1, L
  120    SCALE = SCALE + ABS(D(K))
C
         IF (SCALE .NE. 0.0E0) GO TO 140
C
         DO 125 J = 1, L
            D(J) = A(L,J)
            A(L,J) = A(I,J)
            A(I,J) = 0.0E0
  125    CONTINUE
C
  130    E(I) = 0.0E0
         E2(I) = 0.0E0
         GO TO 300
C
  140    DO 150 K = 1, L
            D(K) = D(K) / SCALE
            H = H + D(K) * D(K)
  150    CONTINUE
C
         E2(I) = SCALE * SCALE * H
         F = D(L)
         G = -SIGN(SQRT(H),F)
         E(I) = SCALE * G
         H = H - F * G
         D(L) = F - G
         IF (L .EQ. 1) GO TO 285
C     .......... FORM A*U ..........
         DO 170 J = 1, L
  170    E(J) = 0.0E0
C
         DO 240 J = 1, L
            F = D(J)
            G = E(J) + A(J,J) * F
            JP1 = J + 1
            IF (L .LT. JP1) GO TO 220
C
            DO 200 K = JP1, L
               G = G + A(K,J) * D(K)
               E(K) = E(K) + A(K,J) * F
  200       CONTINUE
C
  220       E(J) = G
  240    CONTINUE
C     .......... FORM P ..........
         F = 0.0E0
C
         DO 245 J = 1, L
            E(J) = E(J) / H
            F = F + E(J) * D(J)
  245    CONTINUE
C
         H = F / (H + H)
C     .......... FORM Q ..........
         DO 250 J = 1, L
  250    E(J) = E(J) - H * D(J)
C     .......... FORM REDUCED A ..........
         DO 280 J = 1, L
            F = D(J)
            G = E(J)
C
            DO 260 K = J, L
  260       A(K,J) = A(K,J) - F * E(K) - G * D(K)
C
  280    CONTINUE
C
  285    DO 290 J = 1, L
            F = D(J)
            D(J) = A(L,J)
            A(L,J) = A(I,J)
            A(I,J) = F * SCALE
  290    CONTINUE
C
  300 CONTINUE
C
      RETURN
      END
      SUBROUTINE TRED2(NM,N,A,D,E,Z)
C
      INTEGER I,J,K,L,N,II,NM,JP1
      REAL A(NM,N),D(N),E(N),Z(NM,N)
      REAL F,G,H,HH,SCALE
C
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2,
C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C
C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A
C     SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING
C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
C
C     ON INPUT
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.
C
C        N IS THE ORDER OF THE MATRIX.
C
C        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE
C          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
C
C     ON OUTPUT
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX.
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
C
C        Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX
C          PRODUCED IN THE REDUCTION.
C
C        A AND Z MAY COINCIDE.  IF DISTINCT, A IS UNALTERED.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
C
C     THIS VERSION DATED AUGUST 1983.
C
C     ------------------------------------------------------------------
C
      DO 100 I = 1, N
C
         DO 80 J = I, N
   80    Z(J,I) = A(J,I)
C
         D(I) = A(N,I)
  100 CONTINUE
C
      IF (N .EQ. 1) GO TO 510
C     .......... FOR I=N STEP -1 UNTIL 2 DO -- ..........
      DO 300 II = 2, N
         I = N + 2 - II
         L = I - 1
         H = 0.0E0
         SCALE = 0.0E0
         IF (L .LT. 2) GO TO 130
C     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
         DO 120 K = 1, L
  120    SCALE = SCALE + ABS(D(K))
C
         IF (SCALE .NE. 0.0E0) GO TO 140
  130    E(I) = D(L)
C
         DO 135 J = 1, L
            D(J) = Z(L,J)
            Z(I,J) = 0.0E0
            Z(J,I) = 0.0E0
  135    CONTINUE
C
         GO TO 290
C
  140    DO 150 K = 1, L
            D(K) = D(K) / SCALE
            H = H + D(K) * D(K)
  150    CONTINUE
C
         F = D(L)
         G = -SIGN(SQRT(H),F)
         E(I) = SCALE * G
         H = H - F * G
         D(L) = F - G
C     .......... FORM A*U ..........
         DO 170 J = 1, L
  170    E(J) = 0.0E0
C
         DO 240 J = 1, L
            F = D(J)
            Z(J,I) = F
            G = E(J) + Z(J,J) * F
            JP1 = J + 1
            IF (L .LT. JP1) GO TO 220
C
            DO 200 K = JP1, L
               G = G + Z(K,J) * D(K)
               E(K) = E(K) + Z(K,J) * F
  200       CONTINUE
C
  220       E(J) = G
  240    CONTINUE
C     .......... FORM P ..........
         F = 0.0E0
C
         DO 245 J = 1, L
            E(J) = E(J) / H
            F = F + E(J) * D(J)
  245    CONTINUE
C
         HH = F / (H + H)
C     .......... FORM Q ..........
         DO 250 J = 1, L
  250    E(J) = E(J) - HH * D(J)
C     .......... FORM REDUCED A ..........
         DO 280 J = 1, L
            F = D(J)
            G = E(J)
C
            DO 260 K = J, L
  260       Z(K,J) = Z(K,J) - F * E(K) - G * D(K)
C
            D(J) = Z(L,J)
            Z(I,J) = 0.0E0
  280    CONTINUE
C
  290    D(I) = H
  300 CONTINUE
C     .......... ACCUMULATION OF TRANSFORMATION MATRICES ..........
      DO 500 I = 2, N
         L = I - 1
         Z(N,L) = Z(L,L)
         Z(L,L) = 1.0E0
         H = D(I)
         IF (H .EQ. 0.0E0) GO TO 380
C
         DO 330 K = 1, L
  330    D(K) = Z(K,I) / H
C
         DO 360 J = 1, L
            G = 0.0E0
C
            DO 340 K = 1, L
  340       G = G + Z(K,I) * Z(K,J)
C
            DO 360 K = 1, L
               Z(K,J) = Z(K,J) - G * D(K)
  360    CONTINUE
C
  380    DO 400 K = 1, L
  400    Z(K,I) = 0.0E0
C
  500 CONTINUE
C
  510 DO 520 I = 1, N
         D(I) = Z(N,I)
         Z(N,I) = 0.0E0
  520 CONTINUE
C
      Z(N,N) = 1.0E0
      E(1) = 0.0E0
      RETURN
      END
      SUBROUTINE TQL2(NM,N,D,E,Z,IERR)
C
      INTEGER I,J,K,L,M,N,II,L1,L2,NM,MML,IERR
      REAL D(N),E(N),Z(NM,N)
      REAL C,C2,C3,DL1,EL1,F,G,H,P,R,S,S2,TST1,TST2,PYTHAG
C
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2,
C     NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND
C     WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
C
C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
C     OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD.
C     THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
C     BE FOUND IF  TRED2  HAS BEEN USED TO REDUCE THIS
C     FULL MATRIX TO TRIDIAGONAL FORM.
C
C     ON INPUT
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.
C
C        N IS THE ORDER OF THE MATRIX.
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
C
C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
C          REDUCTION BY  TRED2, IF PERFORMED.  IF THE EIGENVECTORS
C          OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
C          THE IDENTITY MATRIX.
C
C      ON OUTPUT
C
C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
C          UNORDERED FOR INDICES 1,2,...,IERR-1.
C
C        E HAS BEEN DESTROYED.
C
C        Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
C          TRIDIAGONAL (OR FULL) MATRIX.  IF AN ERROR EXIT IS MADE,
C          Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
C          EIGENVALUES.
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS.
C
C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
C
C     THIS VERSION DATED AUGUST 1983.
C
C     ------------------------------------------------------------------
C
      IERR = 0
      IF (N .EQ. 1) GO TO 1001
C
      DO 100 I = 2, N
  100 E(I-1) = E(I)
C
      F = 0.0E0
      TST1 = 0.0E0
      E(N) = 0.0E0
C
      DO 240 L = 1, N
         J = 0
         H = ABS(D(L)) + ABS(E(L))
         IF (TST1 .LT. H) TST1 = H
C     .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT ..........
         DO 110 M = L, N
            TST2 = TST1 + ABS(E(M))
            IF (TST2 .EQ. TST1) GO TO 120
C     .......... E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C                THROUGH THE BOTTOM OF THE LOOP ..........
  110    CONTINUE
C
  120    IF (M .EQ. L) GO TO 220
  130    IF (J .EQ. 30) GO TO 1000
         J = J + 1
C     .......... FORM SHIFT ..........
         L1 = L + 1
         L2 = L1 + 1
         G = D(L)
         P = (D(L1) - G) / (2.0E0 * E(L))
         R = PYTHAG(P,1.0E0)
         D(L) = E(L) / (P + SIGN(R,P))
         D(L1) = E(L) * (P + SIGN(R,P))
         DL1 = D(L1)
         H = G - D(L)
         IF (L2 .GT. N) GO TO 145
C
         DO 140 I = L2, N
  140    D(I) = D(I) - H
C
  145    F = F + H
C     .......... QL TRANSFORMATION ..........
         P = D(M)
         C = 1.0E0
         C2 = C
         EL1 = E(L1)
         S = 0.0E0
         MML = M - L
C     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
         DO 200 II = 1, MML
            C3 = C2
            C2 = C
            S2 = S
            I = M - II
            G = C * E(I)
            H = C * P
            R = PYTHAG(P,E(I))
            E(I+1) = S * R
            S = E(I) / R
            C = P / R
            P = C * D(I) - S * G
            D(I+1) = H + S * (C * G + S * D(I))
C     .......... FORM VECTOR ..........
            DO 180 K = 1, N
               H = Z(K,I+1)
               Z(K,I+1) = S * Z(K,I) + C * H
               Z(K,I) = C * Z(K,I) - S * H
  180       CONTINUE
C
  200    CONTINUE
C
         P = -S * S2 * C3 * EL1 * E(L) / DL1
         E(L) = S * P
         D(L) = C * P
         TST2 = TST1 + ABS(E(L))
         IF (TST2 .GT. TST1) GO TO 130
  220    D(L) = D(L) + F
  240 CONTINUE
C     .......... ORDER EIGENVALUES AND EIGENVECTORS ..........
      DO 300 II = 2, N
         I = II - 1
         K = I
         P = D(I)
C
         DO 260 J = II, N
            IF (D(J) .GE. P) GO TO 260
            K = J
            P = D(J)
  260    CONTINUE
C
         IF (K .EQ. I) GO TO 300
         D(K) = D(I)
         D(I) = P
C
         DO 280 J = 1, N
            P = Z(J,I)
            Z(J,I) = Z(J,K)
            Z(J,K) = P
  280    CONTINUE
C
  300 CONTINUE
C
      GO TO 1001
C     .......... SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS ..........
 1000 IERR = L
 1001 RETURN
      END
      SUBROUTINE TQLRAT(N,D,E2,IERR)
C
      INTEGER I,J,L,M,N,II,L1,MML,IERR
      REAL D(N),E2(N)
      REAL B,C,F,G,H,P,R,S,T,EPSLON,PYTHAG
C
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT,
C     ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH.
C
C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
C     TRIDIAGONAL MATRIX BY THE RATIONAL QL METHOD.
C
C     ON INPUT
C
C        N IS THE ORDER OF THE MATRIX.
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
C
C        E2 CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF THE
C          INPUT MATRIX IN ITS LAST N-1 POSITIONS.  E2(1) IS ARBITRARY.
C
C      ON OUTPUT
C
C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
C          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
C          THE SMALLEST EIGENVALUES.
C
C        E2 HAS BEEN DESTROYED.
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS.
C
C     CALLS PYTHAG FOR  SQRT(A*A + B*B) .
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
C
C     THIS VERSION DATED AUGUST 1983.
C
C     ------------------------------------------------------------------
C
      IERR = 0
      IF (N .EQ. 1) GO TO 1001
C
      DO 100 I = 2, N
  100 E2(I-1) = E2(I)
C
      F = 0.0E0
      T = 0.0E0
      E2(N) = 0.0E0
C
      DO 290 L = 1, N
         J = 0
         H = ABS(D(L)) + SQRT(E2(L))
         IF (T .GT. H) GO TO 105
         T = H
         B = EPSLON(T)
         C = B * B
C     .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ..........
  105    DO 110 M = L, N
            IF (E2(M) .LE. C) GO TO 120
C     .......... E2(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C                THROUGH THE BOTTOM OF THE LOOP ..........
  110    CONTINUE
C
  120    IF (M .EQ. L) GO TO 210
  130    IF (J .EQ. 30) GO TO 1000
         J = J + 1
C     .......... FORM SHIFT ..........
         L1 = L + 1
         S = SQRT(E2(L))
         G = D(L)
         P = (D(L1) - G) / (2.0E0 * S)
         R = PYTHAG(P,1.0E0)
         D(L) = S / (P + SIGN(R,P))
         H = G - D(L)
C
         DO 140 I = L1, N
  140    D(I) = D(I) - H
C
         F = F + H
C     .......... RATIONAL QL TRANSFORMATION ..........
         G = D(M)
         IF (G .EQ. 0.0E0) G = B
         H = G
         S = 0.0E0
         MML = M - L
C     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
         DO 200 II = 1, MML
            I = M - II
            P = G * H
            R = P + E2(I)
            E2(I+1) = S * R
            S = E2(I) / R
            D(I+1) = H + S * (H + D(I))
            G = D(I) - E2(I) / G
            IF (G .EQ. 0.0E0) G = B
            H = G * P / R
  200    CONTINUE
C
         E2(L) = S * G
         D(L) = H
C     .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST ..........
         IF (H .EQ. 0.0E0) GO TO 210
         IF (ABS(E2(L)) .LE. ABS(C/H)) GO TO 210
         E2(L) = H * E2(L)
         IF (E2(L) .NE. 0.0E0) GO TO 130
  210    P = D(L) + F
C     .......... ORDER EIGENVALUES ..........
         IF (L .EQ. 1) GO TO 250
C     .......... FOR I=L STEP -1 UNTIL 2 DO -- ..........
         DO 230 II = 2, L
            I = L + 2 - II
            IF (P .GE. D(I-1)) GO TO 270
            D(I) = D(I-1)
  230    CONTINUE
C
  250    I = 1
  270    D(I) = P
  290 CONTINUE
C
      GO TO 1001
C     .......... SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS ..........
 1000 IERR = L
 1001 RETURN
      END
      REAL FUNCTION PYTHAG(A,B)
      REAL A,B
C
C     FINDS SQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW
C
      REAL P,R,S,T,U
      P = AMAX1(ABS(A),ABS(B))
      IF (P .EQ. 0.0E0) GO TO 20
      R = (AMIN1(ABS(A),ABS(B))/P)**2
   10 CONTINUE
         T = 4.0E0 + R
         IF (T .EQ. 4.0E0) GO TO 20
         S = R/T
         U = 1.0E0 + 2.0E0*S
         P = U*P
         R = (S/U)**2 * R
      GO TO 10
   20 PYTHAG = P
      RETURN
      END
C
      REAL FUNCTION EPSLON (X)
      REAL X
      REAL A,B,C,EPS
      A = 4.0E0/3.0E0
  10  B = A - 1.0E0
      C = B + B + B
      EPS = ABS(C-1.0E0)
      IF(EPS .EQ. 0.0E0)GO TO 10
      EPSLON = EPS*ABS(X)
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  LSVRR/DLSVRR (SINGLE/DOUBLE PRECISION VERSION)
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    JANUARY 1, 1985
C
C  PURPOSE:    COMPUTE THE SINGULAR VALUE DECOMPOSITION OF A REAL
C              MATRIX.
C
C  USAGE:      CALL LSVRR (NRA, NCA, A, LDA, IPATH, TOL, IRANK,
C                          S, U, LDU, V, LDV)
C
C  ARGUMENTS:
C     NRA    - NUMBER OF ROWS OF A.  (INPUT)
C     NCA    - NUMBER OF COLUMNS OF A.  (INPUT)
C     A      - NRA BY NCA MATRIX WHOSE SINGULAR VALUE DECOMPOSITION
C              IS TO BE COMPUTED.  (INPUT)
C     LDA    - LEADING DIMENSION OF A EXACTLY AS SPECIFIED IN THE
C              DIMENSION STATEMENT OF THE CALLING PROGRAM.  (INPUT)
C     IPATH  - FLAG USED TO CONTROL THE COMPUTATION OF THE SINGULAR
C              VECTORS.  (INPUT)
C              IPATH HAS THE DECIMAL EXPANSION IJ SUCH THAT:
C                 I = 0 MEANS DO NOT COMPUTE THE LEFT SINGULAR VECTORS,
C                 I = 1 MEANS RETURN THE NCA LEFT SINGULAR VECTORS IN U,
C                 I = 2 MEANS RETURN ONLY THE MIN(NRA,NCA) LEFT SINGULAR
C                       VECTORS IN U,
C                 J = 0 MEANS DO NOT COMPUTE THE RIGHT SINGULAR VECTORS,
C                 J = 1 MEANS RETURN THE RIGHT SINGULAR VECTORS IN V.
C              FOR EXAMPLE, IPATH = 20 MEANS I = 2 AND J = 0.
C     TOL    - SCALAR CONTAINING THE TOLERANCE USED TO DETERMINE WHEN A
C              SINGULAR VALUE IS NEGLIGIBLE.  (INPUT)
C              IF TOL IS POSITIVE THEN A SINGULAR VALUE SI IS CONSIDERED
C                 NEGLIGIBLE IF SI .LE. TOL.
C              IF TOL IS NEGATIVE THEN A SINGULAR VALUE SI IS CONSIDERED
C                 NEGLIGIBLE IF SI .LE. ABS(TOL)*(INFINITY NORM OF A).
C                 IN THIS CASE ABS(TOL) SHOULD GENERALLY CONTAIN AN
C                 ESTIMATE OF THE LEVEL OF RELATIVE ERROR IN THE DATA.
C     IRANK  - SCALAR CONTAINING AN ESTIMATE OF THE RANK OF A.  (OUTPUT)
C     S      - VECTOR OF LENGTH MIN(NRA+1,NCA) CONTAINING THE SINGULAR
C              VALUES OF A IN DESCENDING ORDER OF MAGNITUDE IN THE FIRST
C              MIN(NRA,NCA) POSITIONS.  (OUTPUT)
C     U      - NRA BY NCU MATRIX CONTAINING THE LEFT SINGULAR VECTORS OF
C              A.  (OUTPUT)
C              NCU MUST BE EQUAL TO NRA IF I IS EQUAL TO 1.
C              NCU MUST BE EQUAL TO MIN(NRA,NCA) IF I IS EQUAL TO 2.
C              U WILL NOT BE REFERENCED IF I IS EQUAL TO ZERO.  IF NRA
C              IS LESS THAN OR EQUAL TO NCU, THEN U CAN SHARE THE SAME
C              STORAGE LOCATIONS AS A.  SEE REMARKS.
C     LDU    - LEADING DIMENSION OF U EXACTLY AS SPECIFIED IN THE
C              DIMENSION STATEMENT OF THE CALLING PROGRAM.  (INPUT)
C     V      - NCA BY NCA MATRIX CONTAINING THE RIGHT SINGULAR VECTORS
C              OF A.  (OUTPUT)
C              V WILL NOT BE REFERENCED IF J IS EQUAL TO ZERO.  V CAN
C              SHARE THE SAME STORAGE LOCATION AS A; HOWEVER U AND V
C              CANNOT BOTH COINCIDE WITH A SIMULTANEOUSLY.
C     LDV    - LEADING DIMENSION OF V EXACTLY AS SPECIFIED IN THE
C              DIMENSION STATEMENT OF THE CALLING PROGRAM.  (INPUT)
C
C  REMARKS:
C  1. AUTOMATIC WORKSPACE USAGE IS
C              LSVRR    NRA*NCA + NRA + NCA + MAX(NRA,NCA) - 1 UNITS, OR
C              DLSVRR   2*(NRA*NCA + NRA + NCA + MAX(NRA,NCA) - 1)
C                       UNITS.
C     WORKSPACE MAY BE EXPLICITLY PROVIDED, IF DESIRED, BY USE OF
C     L2VRR/DL2VRR.  THE REFERENCE IS
C              CALL L2VRR (NRA, NCA, A, LDA, IPATH, TOL, IRANK,
C                          S, U, LDU, V, LDV, ACOPY, WK)
C     THE ADDITIONAL ARGUMENTS ARE AS FOLLOWS:
C     ACOPY  - WORK VECTOR OF LENGTH NRA*NCA FOR THE MATRIX A.  IF
C              A IS NOT NEEDED THEN A AND ACOPY MAY SHARE THE SAME
C              STORAGE LOCATIONS.
C     WK     - WORK VECTOR OF LENGTH NRA + NCA + MAX(NRA,NCA) - 1.
C
C  2. INFORMATIONAL ERROR
C     TYPE CODE
C       4   1  CONVERGENCE CANNOT BE ACHIEVED FOR ALL THE SINGULAR
C              VALUES AND THEIR CORRESPONDING SINGULAR VECTORS.
C
C  3. WHEN NRA IS MUCH GREATER THAN NCA, IT MIGHT NOT BE REASONABLE TO
C     STORE THE WHOLE MATRIX U.  IN THIS CASE IPATH WITH I = 2 ALLOWS
C     A SINGULAR VALUE FACTORIZATION OF A TO BE COMPUTED IN WHICH ONLY
C     THE FIRST NCA COLUMNS OF U ARE COMPUTED, AND IN MANY APPLICATIONS
C     THOSE ARE ALL THAT ARE NEEDED.
C
C  KEYWORDS:   LEAST SQUARES; COMPLETE ORTHOGONAL DECOMPOSITION;
C              RANK ESTIMATION
C
C  GAMS:       D6
C
C  CHAPTERS:   MATH/LIBRARY LINEAR SYSTEMS
C              STAT/LIBRARY MATHEMATICAL SUPPORT
C
C  COPYRIGHT:  1985 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE LSVRR (NRA, NCA, A, LDA, IPATH, TOL, IRANK, S, U,
     &                  LDU, V, LDV)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    NRA, NCA, LDA, IPATH, IRANK, LDU, LDV
      REAL       TOL, A(LDA,*), S(*), U(LDU,*), V(LDV,*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    INDA, INDW
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                                  SPECIFICATIONS FOR COMMON /WORKSP/
      REAL       RWKSP(61913)
      REAL       RDWKSP(5000)
      DOUBLE PRECISION DWKSP(2500)
      COMPLEX    CWKSP(2500)
      COMPLEX    CZWKSP(2500)
      COMPLEX    *16 ZWKSP(1250)
      INTEGER    IWKSP(5000)
      LOGICAL    LWKSP(5000)
      EQUIVALENCE (DWKSP(1), RWKSP(1))
      EQUIVALENCE (CWKSP(1), RWKSP(1)), (ZWKSP(1), RWKSP(1))
      EQUIVALENCE (IWKSP(1), RWKSP(1)), (LWKSP(1), RWKSP(1))
      EQUIVALENCE (RDWKSP(1), RWKSP(1)), (CZWKSP(1), RWKSP(1))
      COMMON     /WORKSP/ RWKSP
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  MAX0
      INTRINSIC  MAX0
      INTEGER    MAX0
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   E1MES, E1POP, E1PSH, E1STI, L2VRR
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   I1KGT, N1RCD
      INTEGER    I1KGT, N1RCD
C
      CALL E1PSH ('LSVRR ')
C
      IF (NRA.LE.0 .OR. NCA.LE.0) THEN
         CALL E1STI (1, NRA)
         CALL E1STI (2, NCA)
         CALL E1MES (5, 1, 'BOTH THE NUMBER OF ROWS AND THE '//
     &               'NUMBER OF COLUMNS OF THE INPUT MATRIX HAVE TO '//
     &               'BE POSITIVE WHILE NRA = %(I1) AND NCA = %(I2) '//
     &               'ARE GIVEN.')
      ELSE
         INDA = I1KGT(NRA*NCA,3)
         INDW = I1KGT(NRA+NCA+MAX0(NCA,NRA)-1,3)
         IF (N1RCD(0) .NE. 0) THEN
            CALL E1MES (5, 2, ' ')
            CALL E1STI (1, NRA)
            CALL E1STI (2, NCA)
            CALL E1MES (5, 2, 'THE WORKSPACE IS BASED ON NRA AND NCA '//
     &                  'WHERE NRA = %(I1) AND NCA = %(I2) ARE '//
     &                  'GIVEN.')
         ELSE
            CALL L2VRR (NRA, NCA, A, LDA, IPATH, TOL, IRANK, S, U,
     &                  LDU, V, LDV, RDWKSP(INDA), RDWKSP(INDW))
         END IF
      END IF
C
      CALL E1POP ('LSVRR ')
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  L2VRR/DL2VRR (SINGLE/DOUBLE PRECISION VERSION)
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    JANUARY 1, 1985
C
C  PURPOSE:    COMPUTE THE SINGULAR VALUE DECOMPOSITION OF A REAL
C              MATRIX.
C
C  USAGE:      CALL L2VRR (NRA, NCA, A, LDA, IPATH, TOL, IRANK,
C                          S, U, LDU, V, LDV, WKA, WK)
C
C  ARGUMENTS:  SEE LSVRR/DLSVRR.
C
C  REMARKS:    SEE LSVRR/DLSVRR.
C
C  CHAPTER:    MATH/LIBRARY LINEAR SYSTEMS
C
C  COPYRIGHT:  1985 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE L2VRR (NRA, NCA, A, LDA, IPATH, TOL, IRANK, S, U,
     &                  LDU, V, LDV, WKA, WK)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    NRA, NCA, LDA, IPATH, IRANK, LDU, LDV
      REAL       TOL, A(LDA,*), S(*), U(*), V(*), WKA(*), WK(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, INDE, INDW, JOBU, JOBV
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  MOD
      INTRINSIC  MOD
      INTEGER    MOD
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   E1MES, E1POP, E1PSH, E1STI, SCOPY, L3VRR
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   N1RCD
      INTEGER    N1RCD
C
      CALL E1PSH ('L2VRR   ')
C
      IF (NRA.LE.0 .OR. NCA.LE.0) THEN
         CALL E1STI (1, NRA)
         CALL E1STI (2, NCA)
         CALL E1MES (5, 1, 'BOTH THE NUMBER OF ROWS AND THE '//
     &               'NUMBER OF COLUMNS OF THE INPUT MATRIX HAVE TO '//
     &               'BE POSITIVE WHILE NRA = %(I1) AND NCA = %(I2) '//
     &               'ARE GIVEN.')
      ELSE
         IF (NRA .GT. LDA) THEN
            CALL E1STI (1, NRA)
            CALL E1STI (2, LDA)
            CALL E1MES (5, 2, 'THE NUMBER OF ROWS OF A MUST BE '//
     &                  'LESS THAN OR EQUAL TO ITS LEADING DIMENSION '//
     &                  'WHILE NRA = %(I1) AND LDA = %(I2) ARE '//
     &                  'GIVEN.')
         END IF
C
         IF (LDU .LE. 0) THEN
            CALL E1STI (1, LDU)
            CALL E1MES (5, 3, 'THE LEADING DIMENSION OF U MUST '//
     &                  'BE GREATER THAN ZERO.  LDU = %(I1) IS '//
     &                  'GIVEN.')
         END IF
C
         IF (LDV .LE. 0) THEN
            CALL E1STI (1, LDV)
            CALL E1MES (5, 4, 'THE LEADING DIMENSION OF V MUST '//
     &                  'BE GREATER THAN ZERO.  LDV = %(I1) IS '//
     &                  'GIVEN.')
         END IF
C
         IF (N1RCD(0) .EQ. 0) THEN
            JOBU = MOD(IPATH,100)/10
            JOBV = MOD(IPATH,10)
            IF ((JOBU.NE.0.AND.JOBU.NE.1.AND.JOBU.NE.2) .OR.
     &          (JOBV.NE.0.AND.JOBV.NE.1)) THEN
               CALL E1STI (1, JOBU)
               CALL E1STI (2, JOBV)
               CALL E1MES (5, 7, 'ERROR IN COMPUTATION CONTROL FLAG. '//
     &                     'THE IJ DECIMAL EXPANSION OF IPATH IS I = '//
     &                     '%(I1) AND J = %(I2).  I MUST BE EITHER '//
     &                     '0, 1 OR 2 AND J MUST BE EITHER 0 OR 1.')
            ELSE
C                                  MAKE A COPY OF A IN WKA AND WORK
C                                  WITH WKA ONLY.
               DO 10  I=1, NCA
                  CALL SCOPY (NRA, A(1,I), 1, WKA((I-1)*NRA+1), 1)
   10          CONTINUE
C
               INDE = 1
               INDW = NCA + 1
               CALL L3VRR (NRA, NCA, WKA, NRA, IPATH, TOL, IRANK, S,
     &                     U, LDU, V, LDV, WK(INDE), WK(INDW))
            END IF
         END IF
      END IF
C
      CALL E1POP ('L2VRR   ')
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  L3VRR/DL3VRR (SINGLE/DOUBLE PRECISION VERSION)
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    JANUARY 1, 1985
C
C  PURPOSE:    COMPUTE THE SINGULAR VALUE DECOMPOSITION OF A REAL
C              MATRIX.
C
C  USAGE:      CALL L3VRR (NRA, NCA, A, LDA, IPATH, TOL, IRANK,
C                          S, U, LDU, V, LDV, E, WORK)
C
C  ARGUMENTS:  SEE LSVRR/DLSVRR
C
C  REMARKS:    SEE LSVRR/DLSVRR.
C
C  CHAPTER:    MATH/LIBRARY LINEAR SYSTEMS
C
C  COPYRIGHT:  1985 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE L3VRR (NRA, NCA, A, LDA, IPATH, TOL, IRANK, S, U,
     &                  LDU, V, LDV, E, WORK)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    NRA, NCA, LDA, IPATH, IRANK, LDU, LDV
      REAL       TOL, A(LDA,*), S(*), U(LDU,*), V(LDV,*), E(*), WORK(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, IEND, INFO, ISTART, ITER, J, JOBU, K, KASE, KK, L,
     &           LL, LLS, LM1, LP1, LS, LU, M, MAXIT, MINMN, MM, MM1,
     &           MP1, NCT, NCTP1, NCU, NRT, NRTP1
      REAL       ANORM, B, C, CS, EL, EMM1, F, G, SCALE, SHIFT, SL,
     &           SM, SMM1, SN, STOL, T, T1, TEST, ZTEST
      LOGICAL    WANTU, WANTV
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  ABS,AMAX1,MAX0,MIN0,MOD,SIGN,SQRT
      INTRINSIC  ABS, AMAX1, MAX0, MIN0, MOD, SIGN, SQRT
      INTEGER    MAX0, MIN0, MOD
      REAL       ABS, AMAX1, SIGN, SQRT
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   E1MES, E1POP, E1PSH, E1STI, SCOPY, SGEMV, SGER, SROT,
     &           SROTG, SSCAL, SSET, SSWAP
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   SASUM, SNRM2
      REAL       SASUM, SNRM2
C
      CALL E1PSH ('L3VRR   ')
C                                  IF TOL.LT.0, COMPUTE THE INFINITY
C                                  NORM OF A AND SINGULAR VALUE
C                                  TOLERANCE FOR RANK ESTIMATION
      IF (TOL .LT. 0.0E0) THEN
         ANORM = 0.0E0
         DO 10  I=1, NRA
            ANORM = AMAX1(ANORM,SASUM(NCA,A(I,1),LDA))
   10    CONTINUE
         STOL = ABS(TOL)*ANORM
      ELSE
         STOL = TOL
      END IF
C                                  SET THE MAXIMUM NUMBER OF
C                                  ITERATIONS.
      MAXIT = 30
C                                  DETERMINE WHAT IS TO BE COMPUTED.
      WANTU = .FALSE.
      WANTV = .FALSE.
      JOBU = MOD(IPATH,100)/10
      NCU = NRA
      IF (JOBU .EQ. 2) NCU = MIN0(NRA,NCA)
      IF (JOBU .NE. 0) WANTU = .TRUE.
      IF (MOD(IPATH,10) .NE. 0) WANTV = .TRUE.
C                                  REDUCE A TO BIDIAGONAL FORM, STORING
C                                  THE DIAGONAL ELEMENTS IN S AND THE
C                                  SUPER-DIAGONAL ELEMENTS IN E.
      INFO = 0
      NCT = MIN0(NRA-1,NCA)
      NRT = MAX0(0,MIN0(NCA-2,NRA))
      LU = MAX0(NCT,NRT)
      IF (LU .GE. 1) THEN
         DO 20  L=1, LU
            LP1 = L + 1
            IF (L .LE. NCT) THEN
C                                  COMPUTE THE TRANSFORMATION FOR THE
C                                  L-TH NCA AND PLACE THE L-TH
C                                  DIAGONAL IN S(L).
               S(L) = SNRM2(NRA-L+1,A(L,L),1)
               IF (S(L) .NE. 0.0E0) THEN
                  IF (A(L,L) .NE. 0.0E0) S(L) = SIGN(S(L),A(L,L))
                  CALL SSCAL (NRA-L+1, 1.0E0/S(L), A(L,L), 1)
                  A(L,L) = 1.0E0 + A(L,L)
               END IF
               S(L) = -S(L)
            END IF
C
            IF (NCA .GE. LP1) THEN
               IF (L.LE.NCT .AND. S(L).NE.0.0E0) THEN
C                                  APPLY THE TRANSFORMATION.
                  CALL SGEMV ('T', NRA-L+1, NCA-LP1+1, -1.0E0/A(L,L),
     &                        A(L,LP1), LDA, A(L,L), 1, 0.0E0,
     &                        WORK(NRA+1), 1)
                  CALL SGER (NRA-L+1, NCA-LP1+1, 1.0E0, A(L,L), 1,
     &                       WORK(NRA+1), 1, A(L,LP1), LDA)
               END IF
C                                  PLACE THE L-TH ROW OF A INTO E FOR
C                                  THE SUBSEQUENT CALCULATION OF THE
C                                  ROW TRANSFORMATION.
               CALL SCOPY (NCA-LP1+1, A(L,LP1), LDA, E(LP1), 1)
            END IF
C                                  PLACE THE TRANSFORMATION IN U FOR
C                                  SUBSEQUENT BACK MULTIPLICATION.
            IF (WANTU .AND. L.LE.NCT) CALL SCOPY (NRA-L+1, A(L,L), 1,
     &          U(L,L), 1)
            IF (L .LE. NRT) THEN
C                                  COMPUTE THE L-TH ROW TRANSFORMATION
C                                  AND PLACE THE L-TH SUPER-DIAGONAL
C                                  IN E(L).
               E(L) = SNRM2(NCA-L,E(LP1),1)
               IF (E(L) .NE. 0.0E0) THEN
                  IF (E(LP1) .NE. 0.0E0) E(L) = SIGN(E(L),E(LP1))
                  CALL SSCAL (NCA-L, 1.0E0/E(L), E(LP1), 1)
                  E(LP1) = 1.0E0 + E(LP1)
               END IF
               E(L) = -E(L)
               IF (LP1.LE.NRA .AND. E(L).NE.0.0E0) THEN
C                                  APPLY THE TRANSFORMATION.
                  CALL SGEMV ('N', NRA-L, NCA-LP1+1, 1.0E0,
     &                        A(LP1,LP1), LDA, E(LP1), 1, 0.0E0,
     &                        WORK(LP1), 1)
                  CALL SGER (NRA-L, NCA-LP1+1, -1.0E0/E(LP1),
     &                       WORK(LP1), 1, E(LP1), 1, A(LP1,LP1), LDA)
               END IF
C                                  PLACE THE TRANSFORMATION IN V FOR
C                                  SUBSEQUENT BACK MULTIPLICATION.
               IF (WANTV) CALL SCOPY (NCA-LP1+1, E(LP1), 1, V(LP1,L),
     &             1)
            END IF
   20    CONTINUE
      END IF
C                                  SET UP THE FINAL BIDIAGONAL MATRIX
C                                  OR ORDER M.
      M = MIN0(NCA,NRA+1)
      NCTP1 = NCT + 1
      NRTP1 = NRT + 1
      IF (NCT .LT. NCA) S(NCTP1) = A(NCTP1,NCTP1)
      IF (NRA .LT. M) S(M) = 0.0E0
      IF (NRTP1 .LT. M) E(NRTP1) = A(NRTP1,M)
      E(M) = 0.0E0
C                                  IF REQUIRED, GENERATE U.
      IF (WANTU) THEN
         IF (NCU .GE. NCTP1) THEN
            DO 30  J=NCTP1, NCU
               CALL SSET (NRA, 0.0E0, U(1,J), 1)
               U(J,J) = 1.0E0
   30       CONTINUE
         END IF
         IF (NCT .GE. 1) THEN
            DO 40  LL=1, NCT
               L = NCT - LL + 1
               IF (S(L) .NE. 0.0E0) THEN
                  LP1 = L + 1
                  IF (NCU .GE. LP1) THEN
                     CALL SGEMV ('T', NRA-L+1, NCU-LP1+1,
     &                           -1.0E0/U(L,L), U(L,LP1), LDU, U(L,L),
     &                           1, 0.0E0, WORK(NRA+1), 1)
                     CALL SGER (NRA-L+1, NCU-LP1+1, 1.0E0, U(L,L), 1,
     &                          WORK(NRA+1), 1, U(L,LP1), LDU)
                  END IF
                  CALL SSCAL (NRA-L+1, -1.0E0, U(L,L), 1)
                  U(L,L) = 1.0E0 + U(L,L)
                  LM1 = L - 1
                  IF (LM1 .GE. 1) CALL SSET (LM1, 0.0E0, U(1,L), 1)
               ELSE
                  CALL SSET (NRA, 0.0E0, U(1,L), 1)
                  U(L,L) = 1.0E0
               END IF
   40       CONTINUE
         END IF
      END IF
C                                  IF IT IS REQUIRED, GENERATE V.
      IF (WANTV) THEN
         DO 50  LL=1, NCA
            L = NCA - LL + 1
            LP1 = L + 1
            IF (L.LE.NRT .AND. E(L).NE.0.0E0) THEN
               CALL SGEMV ('T', NCA-L, NCA-LP1+1, -1.0E0/V(LP1,L),
     &                     V(LP1,LP1), LDV, V(LP1,L), 1, 0.0E0,
     &                     WORK(NRA+1), 1)
               CALL SGER (NCA-L, NCA-LP1+1, 1.0E0, V(LP1,L), 1,
     &                    WORK(NRA+1), 1, V(LP1,LP1), LDV)
            END IF
            CALL SSET (NCA, 0.0E0, V(1,L), 1)
            V(L,L) = 1.0E0
   50    CONTINUE
      END IF
C                                  MAIN ITERATION LOOP FOR THE SINGULAR
C                                  VALUES.
      MM = M
      ITER = 0
   60 CONTINUE
C                                  QUIT IF ALL THE SINGULAR VALUES HAVE
C                                  BEEN FOUND.
      IF (M.NE.0 .AND. ITER.LE.MAXIT) THEN
C                                  THIS SECTION OF THE PROGRAM INSPECTS
C                                  FOR NEGLIGIBLE ELEMENTS IN THE S AND
C                                  E ARRAYS. ON COMPLETION THE VARIABLES
C                                  KASE AND L ARE SET AS FOLLOWS.
C
C                                  KASE = 1     IF S(M) AND E(L-1) ARE
C                                  NEGLIGIBLE AND L.LT.M
C
C                                  KASE = 2     IF S(L) IS NEGLIGIBLE
C                                  AND L.LT.M
C
C                                  KASE = 3     IF E(L-1) IS NEGLIGIBLE,
C                                  L.LT.M, AND S(L), ..., S(M) ARE NOT
C                                  NEGLIGIBLE (QR STEP).
C
C                                  KASE = 4     IF E(M-1) IS NEGLIGIBLE
C                                  (CONVERGENCE).
         DO 70  LL=1, M
            L = M - LL
            IF (L .EQ. 0) GO TO 80
            TEST = ABS(S(L)) + ABS(S(L+1))
            ZTEST = TEST + ABS(E(L))
            IF (ZTEST .EQ. TEST) THEN
               E(L) = 0.0E0
               GO TO 80
            END IF
   70    CONTINUE
   80    CONTINUE
         IF (L .EQ. M-1) THEN
            KASE = 4
         ELSE
            LP1 = L + 1
            MP1 = M + 1
            DO 90  LLS=LP1, MP1
               LS = M - LLS + LP1
               IF (LS .EQ. L) GO TO 100
               TEST = 0.0E0
               IF (LS .NE. M) TEST = TEST + ABS(E(LS))
               IF (LS .NE. L+1) TEST = TEST + ABS(E(LS-1))
               ZTEST = TEST + ABS(S(LS))
               IF (ZTEST .EQ. TEST) THEN
                  S(LS) = 0.0E0
                  GO TO 100
               END IF
   90       CONTINUE
  100       CONTINUE
            IF (LS .EQ. L) THEN
               KASE = 3
            ELSE IF (LS .EQ. M) THEN
               KASE = 1
            ELSE
               KASE = 2
               L = LS
            END IF
         END IF
         L = L + 1
C                                  PERFORM THE TASK INDICATED BY KASE.
         IF (KASE .EQ. 1) THEN
C                                  DEFLATE NEGLIGIBLE S(M).
            MM1 = M - 1
            F = E(M-1)
            E(M-1) = 0.0E0
            DO 110  KK=L, MM1
               K = MM1 - KK + L
               T1 = S(K)
               CALL SROTG (T1, F, CS, SN)
               S(K) = T1
               IF (K .NE. L) THEN
                  F = -SN*E(K-1)
                  E(K-1) = CS*E(K-1)
               END IF
               IF (WANTV) CALL SROT (NCA, V(1,K), 1, V(1,M), 1, CS, SN)
  110       CONTINUE
C                                  SPLIT AT NEGLIGIBLE S(L).
         ELSE IF (KASE .EQ. 2) THEN
            F = E(L-1)
            E(L-1) = 0.0E0
            DO 120  K=L, M
               T1 = S(K)
               CALL SROTG (T1, F, CS, SN)
               S(K) = T1
               F = -SN*E(K)
               E(K) = CS*E(K)
               IF (WANTU) CALL SROT (NRA, U(1,K), 1, U(1,L-1), 1, CS,
     &             SN)
  120       CONTINUE
C                                  PERFORM ONE QR STEP.
         ELSE IF (KASE .EQ. 3) THEN
C                                  CALCULATE THE SHIFT.
            SCALE = AMAX1(ABS(S(M)),ABS(S(M-1)),ABS(E(M-1)),ABS(S(L)),
     &              ABS(E(L)))
            SM = S(M)/SCALE
            SMM1 = S(M-1)/SCALE
            EMM1 = E(M-1)/SCALE
            SL = S(L)/SCALE
            EL = E(L)/SCALE
            B = ((SMM1+SM)*(SMM1-SM)+EMM1**2)/2.0E0
            C = (SM*EMM1)**2
            SHIFT = 0.0E0
            IF (B.NE.0.0E0 .OR. C.NE.0.0E0) THEN
               SHIFT = SQRT(B**2+C)
               IF (B .LT. 0.0E0) SHIFT = -SHIFT
               SHIFT = C/(B+SHIFT)
            END IF
            F = (SL+SM)*(SL-SM) - SHIFT
            G = SL*EL
C                                  CHASE ZEROS.
            MM1 = M - 1
            DO 130  K=L, MM1
               CALL SROTG (F, G, CS, SN)
               IF (K .NE. L) E(K-1) = F
               F = CS*S(K) + SN*E(K)
               E(K) = CS*E(K) - SN*S(K)
               G = SN*S(K+1)
               S(K+1) = CS*S(K+1)
               IF (WANTV) CALL SROT (NCA, V(1,K), 1, V(1,K+1), 1, CS,
     &             SN)
               CALL SROTG (F, G, CS, SN)
               S(K) = F
               F = CS*E(K) + SN*S(K+1)
               S(K+1) = -SN*E(K) + CS*S(K+1)
               G = SN*E(K+1)
               E(K+1) = CS*E(K+1)
               IF (WANTU .AND. K.LT.NRA) CALL SROT (NRA, U(1,K), 1,
     &             U(1,K+1), 1, CS, SN)
  130       CONTINUE
            E(M-1) = F
            ITER = ITER + 1
C                                  CONVERGENCE.
         ELSE IF (KASE .EQ. 4) THEN
C                                  MAKE THE SINGULAR VALUE POSITIVE.
            IF (S(L) .LT. 0.0E0) THEN
               S(L) = -S(L)
               IF (WANTV) CALL SSCAL (NCA, -1.0E0, V(1,L), 1)
            END IF
C                                  ORDER THE SINGULAR VALUE.
  140       IF (L .EQ. MM) GO TO 150
            IF (S(L) .GE. S(L+1)) GO TO 150
            T = S(L)
            S(L) = S(L+1)
            S(L+1) = T
            IF (WANTV .AND. L.LT.NCA) CALL SSWAP (NCA, V(1,L), 1, V(1,
     &          L+1), 1)
            IF (WANTU .AND. L.LT.NRA) CALL SSWAP (NRA, U(1,L), 1, U(1,
     &          L+1), 1)
            L = L + 1
            GO TO 140
  150       CONTINUE
            ITER = 0
            M = M - 1
         END IF
         GO TO 60
      END IF
      IF (ITER .GT. MAXIT) INFO = M
      IF (INFO .EQ. 0) THEN
C                                  ESTIMATE THE RANK OF A
         LS = MIN0(NCA+1,NRA)
         MINMN = MIN0(NCA,NRA)
         IRANK = MINMN
         DO 160  I=1, MINMN
            IF (ABS(S(I)) .LE. STOL) THEN
               IRANK = I - 1
               GO TO 170
            END IF
  160    CONTINUE
  170    CONTINUE
      ELSE
         ISTART = INFO + 1
         IEND = MIN0(NRA,NCA)
         CALL E1STI (1, ISTART)
         CALL E1STI (2, IEND)
C         CALL E1MES (4, 1, 'CONVERGENCE CAN ONLY BE OBTAINED '//
C     &               'FOR THE %(I1), ... , %(I2) SINGULAR VALUES '//
C     &               'AND THEIR CORRESPONDING SINGULAR VECTORS.')
      END IF
      CALL E1POP ('L3VRR ')
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  E1MES
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    MARCH 2, 1984
C
C  PURPOSE:    SET AN ERROR STATE FOR THE CURRENT LEVEL IN THE STACK.
C              THE MESSAGE IS PRINTED IMMEDIATELY IF THE ERROR TYPE IS
C              5, 6, OR 7 AND THE PRINT ATTRIBUTE FOR THAT TYPE IS YES.
C
C  USAGE:      CALL E1MES(IERTYP,IERCOD,MSGPKD)
C
C  ARGUMENTS:
C     IERTYP - INTEGER SPECIFYING THE ERROR TYPE.  (INPUT)
C                IERTYP=1,  INFORMATIONAL/NOTE
C                IERTYP=2,  INFORMATIONAL/ALERT
C                IERTYP=3,  INFORMATIONAL/WARNING
C                IERTYP=4,  INFORMATIONAL/FATAL
C                IERTYP=5,  TERMINAL
C                IERTYP=6,  PROTRAN/WARNING
C                IERTYP=7,  PROTRAN/FATAL
C     IERCOD - INTEGER SPECIFYING THE ERROR CODE.  (INPUT)
C     MSGPKD - A CHARACTER STRING CONTAINING THE MESSAGE.
C              (INPUT)  WITHIN THE MESSAGE, ANY OF FOLLOWING MAY APPEAR
C                %(A1),%(A2),...,%(A9) FOR CHARACTER ARRAYS
C                %(C1),%(C2),...,%(C9) FOR COMPLEX NUMBERS
C                %(D1),%(D2),...,%(D9) FOR DOUBLE PRECISION NUMBERS
C                %(I1),%(I2),...,%(I9) FOR INTEGER NUMBERS
C                %(K1),%(K2),...,%(K9) FOR KEYWORDS
C                %(L1),%(L2),...,%(L9) FOR LITERALS (STRINGS)
C                %(R1),%(R2),...,%(R9) FOR REAL NUMBERS
C                %(Z1),%(Z2),...,%(Z9) FOR DOUBLE COMPLEX NUMBERS
C              THIS PROVIDES A WAY TO INSERT CHARACTER ARRAYS, STRINGS,
C              NUMBERS, AND KEYWORDS INTO THE MESSAGE.  SEE REMARKS
C              BELOW.
C
C  REMARKS:
C     THE NUMBER OF CHARACTERS IN THE MESSAGE AFTER THE INSERTION OF
C     THE CORRESPONDING STRINGS, ETC. SHOULD NOT EXCEED 255.  IF THE
C     LIMIT IS EXCEEDED, ONLY THE FIRST 255 CHARACTERS WILL BE USED.
C     THE APPROPRIATE STRINGS, ETC. NEED TO HAVE BEEN PREVIOUSLY STORED
C     IN COMMON VIA CALLS TO E1STA, E1STD, ETC.  LINE BREAKS MAY BE
C     SPECIFIED BY INSERTING THE TWO CHARACTERS '%/' INTO THE MESSAGE
C     AT THE DESIRED LOCATIONS.
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE E1MES (IERTYP, IERCOD, MSGPKD)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    IERTYP, IERCOD
      CHARACTER  MSGPKD*(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    ERTYP2, I, IER, IPLEN, ISUB, LAST, LEN2, LOC, M, MS,
     &           NLOC, NUM, PBEG
      CHARACTER  MSGTMP(255)
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      INTEGER    IFINIT, NFORMS
      CHARACTER  BLNK, DBB(3), FIND(4), FORMS(9), INREF(25), LPAR,
     &           NCHECK(3), PERCNT, RPAR
      SAVE       BLNK, DBB, FIND, FORMS, IFINIT, INREF, LPAR, NCHECK,
     &           NFORMS, PERCNT, RPAR
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                              SPECIFICATIONS FOR COMMON /ERCOM1/
      INTEGER    CALLVL, MAXLEV, MSGLEN, ERTYPE(51), ERCODE(51),
     &           PRINTB(7), STOPTB(7), PLEN, IFERR6, IFERR7,
     &           IALLOC(51), HDRFMT(7), TRACON(7)
      COMMON     /ERCOM1/ CALLVL, MAXLEV, MSGLEN, ERTYPE, ERCODE,
     &           PRINTB, STOPTB, PLEN, IFERR6, IFERR7, IALLOC, HDRFMT,
     &           TRACON
      SAVE       /ERCOM1/
C                              SPECIFICATIONS FOR COMMON /ERCOM2/
      CHARACTER  MSGSAV(255), PLIST(300), RNAME(51)*6
      COMMON     /ERCOM2/ MSGSAV, PLIST, RNAME
      SAVE       /ERCOM2/
C                              SPECIFICATIONS FOR COMMON /ERCOM3/
      DOUBLE PRECISION ERCKSM
      COMMON     /ERCOM3/ ERCKSM
      SAVE       /ERCOM3/
C                              SPECIFICATIONS FOR COMMON /ERCOM4/
      LOGICAL    ISUSER(51)
      COMMON     /ERCOM4/ ISUSER
      SAVE       /ERCOM4/
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  LEN,MIN0
      INTRINSIC  LEN, MIN0
      INTEGER    LEN, MIN0
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   C1TCI, E1INIT, E1PRT, E1UCS, M1VE, M1VECH
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   I1DX
      INTEGER    I1DX
C
      DATA FORMS/'A', 'C', 'D', 'I', 'K', 'L', 'R', 'S', 'Z'/,
     &     NFORMS/9/
      DATA PERCNT/'%'/, LPAR/'('/, RPAR/')'/, BLNK/' '/
      DATA INREF/' ', 'I', 'N', ' ', 'R', 'E', 'F', 'E', 'R',
     &     'E', 'N', 'C', 'E', ' ', 'T', 'O', ' ', 'K', 'E',
     &     'Y', 'W', 'O', 'R', 'D', ' '/
      DATA NCHECK/'N', '1', '*'/, DBB/'.', ' ', ' '/
      DATA FIND/'*', ' ', ' ', '*'/
      DATA IFINIT/0/
C                                  INITIALIZE ERROR TABLE IF NECESSARY
      IF (IFINIT .EQ. 0) THEN
         CALL E1INIT
         IFINIT = 1
      END IF
C                                  CHECK AND SET ERROR TYPE IF NECESSARY
      IF (IERTYP .NE. -1) THEN
         ERTYPE(CALLVL) = IERTYP
      ELSE IF (IERTYP.LT.-1 .OR. IERTYP.GT.7) THEN
         MSGLEN = 51
         CALL M1VECH ('.  ERROR FROM E1MES.  ILLEGAL ERROR TYPE'//
     &                ' SPECIFIED. ', MSGLEN, MSGSAV, MSGLEN)
         CALL E1PRT
         STOP
      END IF
C
      ERTYP2 = ERTYPE(CALLVL)
C                                  SET ERROR CODE IF NECESSARY
      IF (IERCOD .GT. -1) ERCODE(CALLVL) = IERCOD
      LEN2 = LEN(MSGPKD)
C
      IF (IERTYP.EQ.0 .OR. IERCOD.EQ.0) THEN
C                                  REMOVE THE ERROR STATE
         MSGLEN = 0
      ELSE IF (LEN2.EQ.0 .OR. (LEN2.EQ.1.AND.MSGPKD(1:1).EQ.BLNK)) THEN
         IF (ERTYP2 .EQ. 6) IFERR6 = 1
         IF (ERTYP2 .EQ. 7) IFERR7 = 1
C                                  UPDATE CHECKSUM PARAMETER ERCKSM
         CALL E1UCS
C                                  PRINT MESSAGE IF NECESSARY
         IF (ERTYP2.GE.5 .AND. PRINTB(ERTYP2).EQ.1) CALL E1PRT
      ELSE
C                                  FILL UP MSGSAV WITH EXPANDED MESSAGE
         LEN2 = MIN0(LEN2,255)
         DO 10  I=1, LEN2
            MSGTMP(I) = MSGPKD(I:I)
   10    CONTINUE
         MS = 0
         M = 0
C                                  CHECK PLIST FOR KEYWORD NAME
         NLOC = I1DX(PLIST,PLEN,NCHECK,3)
         IF (NLOC.GT.0 .AND. HDRFMT(ERTYP2).EQ.3) THEN
C                                  M1VE INREF INTO MSGSAV
            CALL M1VE (INREF, 1, 25, 25, MSGSAV, 1, 25, 25, IER)
C                                  GET LENGTH OF KEYWORD NAME
            CALL C1TCI (PLIST(NLOC+3), 3, IPLEN, IER)
            PBEG = NLOC + 3 + IER
C                                  M1VE KEYWORD NAME INTO MSGSAV
            CALL M1VE (PLIST, PBEG, PBEG+IPLEN-1, PLEN, MSGSAV, 26,
     &                 IPLEN+25, 255, IER)
C                                  UPDATE POINTER
            MS = IPLEN + 25
         END IF
C                                  INSERT DOT, BLANK, BLANK
         CALL M1VE (DBB, 1, 3, 3, MSGSAV, MS+1, MS+3, 255, IER)
         MS = MS + 3
C                                  LOOK AT NEXT CHARACTER
   20    M = M + 1
         ISUB = 0
         IF (M .GT. LEN2-4) THEN
            LAST = LEN2 - M + 1
            DO 30  I=1, LAST
   30       MSGSAV(MS+I) = MSGTMP(M+I-1)
            MSGLEN = MS + LAST
            GO TO 40
         ELSE IF (MSGTMP(M).EQ.PERCNT .AND. MSGTMP(M+1).EQ.LPAR .AND.
     &           MSGTMP(M+4).EQ.RPAR) THEN
            CALL C1TCI (MSGTMP(M+3), 1, NUM, IER)
            IF (IER.EQ.0 .AND. NUM.NE.0 .AND. I1DX(FORMS,NFORMS,
     &          MSGTMP(M+2),1).NE.0) THEN
C                                  LOCATE THE ITEM IN THE PARAMETER LIST
               CALL M1VE (MSGTMP(M+2), 1, 2, 2, FIND, 2, 3, 4, IER)
               LOC = I1DX(PLIST,PLEN,FIND,4)
               IF (LOC .GT. 0) THEN
C                                  SET IPLEN = LENGTH OF STRING
                  CALL C1TCI (PLIST(LOC+4), 4, IPLEN, IER)
                  PBEG = LOC + 4 + IER
C                                  ADJUST IPLEN IF IT IS TOO BIG
                  IPLEN = MIN0(IPLEN,255-MS)
C                                  M1VE STRING FROM PLIST INTO MSGSAV
                  CALL M1VE (PLIST, PBEG, PBEG+IPLEN-1, PLEN, MSGSAV,
     &                       MS+1, MS+IPLEN, 255, IER)
                  IF (IER.GE.0 .AND. IER.LT.IPLEN) THEN
C                                  UPDATE POINTERS
                     M = M + 4
                     MS = MS + IPLEN - IER
C                                  BAIL OUT IF NO MORE ROOM
                     IF (MS .GE. 255) THEN
                        MSGLEN = 255
                        GO TO 40
                     END IF
C                                  SET FLAG TO SHOW SUBSTITION WAS MADE
                     ISUB = 1
                  END IF
               END IF
            END IF
         END IF
         IF (ISUB .EQ. 0) THEN
            MS = MS + 1
            MSGSAV(MS) = MSGTMP(M)
         END IF
         GO TO 20
   40    ERTYP2 = ERTYPE(CALLVL)
         IF (ERTYP2 .EQ. 6) IFERR6 = 1
         IF (ERTYP2 .EQ. 7) IFERR7 = 1
C                                  UPDATE CHECKSUM PARAMETER ERCKSM
         CALL E1UCS
C                                  PRINT MESSAGE IF NECESSARY
         IF (ERTYP2.GE.5 .AND. PRINTB(ERTYP2).EQ.1) CALL E1PRT
      END IF
C                                  CLEAR PARAMETER LIST
      PLEN = 1
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  E1POP
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    MARCH 13, 1984
C
C  PURPOSE:    TO POP A SUBROUTINE NAME FROM THE ERROR CONTROL STACK.
C
C  USAGE:      CALL E1POP(NAME)
C
C  ARGUMENTS:
C     NAME   - A CHARACTER STRING OF LENGTH SIX SPECIFYING THE NAME
C              OF THE SUBROUTINE.  (INPUT)
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE E1POP (NAME)
C                                  SPECIFICATIONS FOR ARGUMENTS
      CHARACTER  NAME*(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    IERTYP, IR
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                              SPECIFICATIONS FOR COMMON /ERCOM1/
      INTEGER    CALLVL, MAXLEV, MSGLEN, ERTYPE(51), ERCODE(51),
     &           PRINTB(7), STOPTB(7), PLEN, IFERR6, IFERR7,
     &           IALLOC(51), HDRFMT(7), TRACON(7)
      COMMON     /ERCOM1/ CALLVL, MAXLEV, MSGLEN, ERTYPE, ERCODE,
     &           PRINTB, STOPTB, PLEN, IFERR6, IFERR7, IALLOC, HDRFMT,
     &           TRACON
      SAVE       /ERCOM1/
C                              SPECIFICATIONS FOR COMMON /ERCOM2/
      CHARACTER  MSGSAV(255), PLIST(300), RNAME(51)*6
      COMMON     /ERCOM2/ MSGSAV, PLIST, RNAME
      SAVE       /ERCOM2/
C                              SPECIFICATIONS FOR COMMON /ERCOM3/
      DOUBLE PRECISION ERCKSM
      COMMON     /ERCOM3/ ERCKSM
      SAVE       /ERCOM3/
C                              SPECIFICATIONS FOR COMMON /ERCOM4/
      LOGICAL    ISUSER(51)
      COMMON     /ERCOM4/ ISUSER
      SAVE       /ERCOM4/
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   E1MES, E1PRT, E1PSH, E1STI, E1STL, I1KRL
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   I1KST
      INTEGER    I1KST
C
      IF (CALLVL .LE. 1) THEN
         CALL E1PSH ('E1POP ')
         CALL E1STL (1, NAME)
         CALL E1MES (5, 1, 'ERROR CONDITION IN E1POP.  CANNOT POP '//
     &               'FROM %(L1) BECAUSE STACK IS EMPTY.')
         STOP
      ELSE IF (NAME .NE. RNAME(CALLVL)) THEN
         CALL E1STL (1, NAME)
         CALL E1STL (2, RNAME(CALLVL))
         CALL E1MES (5, 2, 'ERROR CONDITION IN E1POP.  %(L1) DOES '//
     &               'NOT MATCH THE NAME %(L2) IN THE STACK.')
         STOP
      ELSE
         IERTYP = ERTYPE(CALLVL)
         IF (IERTYP .NE. 0) THEN
C                                  M1VE ERROR TYPE AND ERROR CODE TO
C                                    PREVIOUS LEVEL FOR ERROR TYPES 2-7
            IF (IERTYP.GE.2 .AND. IERTYP.LE.7) THEN
               ERTYPE(CALLVL-1) = ERTYPE(CALLVL)
               ERCODE(CALLVL-1) = ERCODE(CALLVL)
            END IF
C                                  CHECK PRINT TABLE TO DETERMINE
C                                    WHETHER TO PRINT STORED MESSAGE
            IF (IERTYP .LE. 4) THEN
               IF (ISUSER(CALLVL-1) .AND. PRINTB(IERTYP).EQ.1)
     &             CALL E1PRT
            ELSE
               IF (PRINTB(IERTYP) .EQ. 1) CALL E1PRT
            END IF
C                                  CHECK STOP TABLE AND ERROR TYPE TO
C                                    DETERMINE WHETHER TO STOP
            IF (IERTYP .LE. 4) THEN
               IF (ISUSER(CALLVL-1) .AND. STOPTB(IERTYP).EQ.1) THEN
                  STOP
               END IF
            ELSE IF (IERTYP .EQ. 5) THEN
               IF (STOPTB(IERTYP) .EQ. 1) THEN
                  STOP
               END IF
            ELSE IF (HDRFMT(IERTYP) .EQ. 1) THEN
               IF (ISUSER(CALLVL-1)) THEN
                  IF (N1RGB(0) .NE. 0) THEN
                     STOP
                  END IF
               END IF
            END IF
         END IF
C                                  SET ERROR TYPE AND CODE
         IF (CALLVL .LT. MAXLEV) THEN
            ERTYPE(CALLVL+1) = -1
            ERCODE(CALLVL+1) = -1
         END IF
C                                  SET IR = AMOUNT OF WORKSPACE
C                                  ALLOCATED AT THIS LEVEL
         IR = I1KST(1) - IALLOC(CALLVL-1)
         IF (IR .GT. 0) THEN
C                                  RELEASE WORKSPACE
            CALL I1KRL (IR)
            IALLOC(CALLVL) = 0
         ELSE IF (IR .LT. 0) THEN
            CALL E1STI (1, CALLVL)
            CALL E1STI (2, IALLOC(CALLVL-1))
            CALL E1STI (3, I1KST(1))
            CALL E1MES (5, 3, 'ERROR CONDITION IN E1POP. '//
     &                  ' THE NUMBER OF WORKSPACE ALLOCATIONS AT '//
     &                  'LEVEL %(I1) IS %(I2).  HOWEVER, THE TOTAL '//
     &                  'NUMBER OF WORKSPACE ALLOCATIONS IS %(I3).')
            STOP
         END IF
C                                  DECREASE THE STACK POINTER BY ONE
         CALLVL = CALLVL - 1
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  E1PSH
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    MARCH 2, 1984
C
C  PURPOSE:    TO PUSH A SUBROUTINE NAME ONTO THE ERROR CONTROL STACK.
C
C  USAGE:      CALL E1PSH(NAME)
C
C  ARGUMENTS:
C     NAME   - A CHARACTER STRING OF LENGTH SIX SPECIFING THE NAME OF
C              THE SUBROUTINE.  (INPUT)
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE E1PSH (NAME)
C                                  SPECIFICATIONS FOR ARGUMENTS
      CHARACTER  NAME*(*)
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      INTEGER    IFINIT
      SAVE       IFINIT
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                              SPECIFICATIONS FOR COMMON /ERCOM1/
      INTEGER    CALLVL, MAXLEV, MSGLEN, ERTYPE(51), ERCODE(51),
     &           PRINTB(7), STOPTB(7), PLEN, IFERR6, IFERR7,
     &           IALLOC(51), HDRFMT(7), TRACON(7)
      COMMON     /ERCOM1/ CALLVL, MAXLEV, MSGLEN, ERTYPE, ERCODE,
     &           PRINTB, STOPTB, PLEN, IFERR6, IFERR7, IALLOC, HDRFMT,
     &           TRACON
      SAVE       /ERCOM1/
C                              SPECIFICATIONS FOR COMMON /ERCOM2/
      CHARACTER  MSGSAV(255), PLIST(300), RNAME(51)*6
      COMMON     /ERCOM2/ MSGSAV, PLIST, RNAME
      SAVE       /ERCOM2/
C                              SPECIFICATIONS FOR COMMON /ERCOM3/
      DOUBLE PRECISION ERCKSM
      COMMON     /ERCOM3/ ERCKSM
      SAVE       /ERCOM3/
C                              SPECIFICATIONS FOR COMMON /ERCOM4/
      LOGICAL    ISUSER(51)
      COMMON     /ERCOM4/ ISUSER
      SAVE       /ERCOM4/
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   E1INIT, E1MES, E1STI
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   I1KST
      INTEGER    I1KST
C
      DATA IFINIT/0/
C                                  INITIALIZE ERROR TABLE IF NECESSARY
      IF (IFINIT .EQ. 0) THEN
         CALL E1INIT
         IFINIT = 1
      END IF
      IF (CALLVL .GE. MAXLEV) THEN
         CALL E1STI (1, MAXLEV)
         CALL E1MES (5, 1, 'ERROR CONDITION IN E1PSH.  PUSH WOULD '//
     &               'CAUSE STACK LEVEL TO EXCEED %(I1). ')
         STOP
      ELSE
C                                  STORE ALLOCATION LEVEL
         IALLOC(CALLVL) = I1KST(1)
C                                  INCREMENT THE STACK POINTER BY ONE
         CALLVL = CALLVL + 1
C                                  PUT SUBROUTINE NAME INTO STACK
         RNAME(CALLVL) = NAME
C                                  SET ERROR TYPE AND ERROR CODE
         ERTYPE(CALLVL) = 0
         ERCODE(CALLVL) = 0
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  E1STI
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    MARCH 6, 1984
C
C  PURPOSE:    TO STORE AN INTEGER FOR SUBSEQUENT USE WITHIN AN ERROR
C              MESSAGE.
C
C  USAGE:      CALL E1STI(II, IVALUE)
C
C  ARGUMENTS:
C     II     - INTEGER SPECIFYING THE SUBSTITUTION INDEX.  II MUST BE
C              BETWEEN 1 AND 9.  (INPUT)
C     IVALUE - THE INTEGER TO BE STORED.  (INPUT)
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE E1STI (II, IVALUE)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    II, IVALUE
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    IBEG, IER, ILEN
      CHARACTER  ARRAY(14)
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      INTEGER    IFINIT
      CHARACTER  BLANK(1)
      SAVE       BLANK, IFINIT
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                              SPECIFICATIONS FOR COMMON /ERCOM1/
      INTEGER    CALLVL, MAXLEV, MSGLEN, ERTYPE(51), ERCODE(51),
     &           PRINTB(7), STOPTB(7), PLEN, IFERR6, IFERR7,
     &           IALLOC(51), HDRFMT(7), TRACON(7)
      COMMON     /ERCOM1/ CALLVL, MAXLEV, MSGLEN, ERTYPE, ERCODE,
     &           PRINTB, STOPTB, PLEN, IFERR6, IFERR7, IALLOC, HDRFMT,
     &           TRACON
      SAVE       /ERCOM1/
C                              SPECIFICATIONS FOR COMMON /ERCOM2/
      CHARACTER  MSGSAV(255), PLIST(300), RNAME(51)*6
      COMMON     /ERCOM2/ MSGSAV, PLIST, RNAME
      SAVE       /ERCOM2/
C                              SPECIFICATIONS FOR COMMON /ERCOM3/
      DOUBLE PRECISION ERCKSM
      COMMON     /ERCOM3/ ERCKSM
      SAVE       /ERCOM3/
C                              SPECIFICATIONS FOR COMMON /ERCOM4/
      LOGICAL    ISUSER(51)
      COMMON     /ERCOM4/ ISUSER
      SAVE       /ERCOM4/
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   C1TIC, E1INIT, E1INPL
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   I1ERIF
      INTEGER    I1ERIF
C
      DATA BLANK/' '/, IFINIT/0/
C                                  INITIALIZE IF NECESSARY
      IF (IFINIT .EQ. 0) THEN
         CALL E1INIT
         IFINIT = 1
      END IF
      CALL C1TIC (IVALUE, ARRAY, 14, IER)
      IBEG = I1ERIF(ARRAY,14,BLANK,1)
      IF (II.GE.1 .AND. II.LE.9 .AND. IER.EQ.0) THEN
         ILEN = 15 - IBEG
         CALL E1INPL ('I', II, ILEN, ARRAY(IBEG))
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  I1KGT
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    JANUARY 17, 1984
C
C  PURPOSE:    ALLOCATE NUMERICAL WORKSPACE.
C
C  USAGE:      I1KGT(NELMTS,ITYPE)
C
C  ARGUMENTS:
C     NELMTS - NUMBER OF ELEMENTS OF DATA TYPE ITYPE TO BE
C              ALLOCATED.  (INPUT)
C     ITYPE  - DATA TYPE OF ARRAY TO BE ALLOCATED.  (INPUT)
C                 1 - LOGICAL
C                 2 - INTEGER
C                 3 - REAL
C                 4 - DOUBLE PRECISION
C                 5 - COMPLEX
C                 6 - DOUBLE COMPLEX
C     I1KGT  - INTEGER FUNCTION.  (OUTPUT)  RETURNS THE INDEX OF THE
C              FIRST ELEMENT IN THE CURRENT ALLOCATION.
C
C  REMARKS:
C  1. ON RETURN, THE ARRAY WILL OCCUPY
C     WKSP(I1KGT), WKSP(I1KGT+1), ..., WKSP(I1KGT+NELMTS-1) WHERE
C     WKSP IS AN ARRAY OF DATA TYPE ITYPE EQUIVALENCED TO RWKSP.
C
C  2. IF I1KGT IS NEGATIVE, THE ABSOLUTE VALUE OF I1KGT IS THE
C     ADDITIONAL WORKSPACE NEEDED FOR THE CURRENT ALLOCATION.
C
C  3. THE ALLOCATOR RESERVES THE FIRST SIXTEEN INTEGER LOCATIONS OF
C     THE STACK FOR ITS OWN INTERNAL BOOKKEEPING.  THESE ARE INITIALIZED
C     BY THE FUNCTION IWKIN UPON THE FIRST CALL TO THE ALLOCATION
C     PACKAGE.
C
C  4. THE USE OF THE FIRST TEN INTEGER LOCATIONS IS AS FOLLOWS:
C      WKSP( 1) - LOUT    THE NUMBER OF CURRENT ALLOCATIONS
C      WKSP( 2) - LNOW    THE CURRENT ACTIVE LENGTH OF THE STACK
C      WKSP( 3) - LUSED   THE MAXIMUM VALUE OF WKSP(2) ACHIEVED
C                         THUS FAR
C      WKSP( 4) - LBND    THE LOWER BOUND OF PERMANENT STORAGE WHICH
C                         IS ONE NUMERIC STORAGE UNIT MORE THAN THE
C                         MAXIMUM ALLOWED LENGTH OF THE STACK.
C      WKSP( 5) - LMAX    THE MAXIMUM LENGTH OF THE STORAGE ARRAY
C      WKSP( 6) - LALC    THE TOTAL NUMBER OF ALLOCATIONS HANDLED BY
C                         I1KGT
C      WKSP( 7) - LNEED   THE NUMBER OF NUMERIC STORAGE UNITS BY WHICH
C                         THE ARRAY SIZE MUST BE INCREASED FOR ALL PAST
C                         ALLOCATIONS TO SUCCEED
C      WKSP( 8) - LBOOK   THE NUMBER OF NUMERIC STORAGE UNITS USED FOR
C                         BOOKKEEPING
C      WKSP( 9) - LCHAR   THE POINTER TO THE PORTION OF THE PERMANENT
C                         STACK WHICH CONTAINS THE BOOKKEEPING AND
C                         POINTERS FOR THE CHARACTER WORKSPACE
C                         ALLOCATION.
C      WKSP(10) - LLCHAR  THE LENGTH OF THE ARRAY BEGINNING AT LCHAR
C                         SET ASIDE FOR CHARACTER WORKSPACE BOOKKEEPING
C                         AND POINTERS.
C                 NOTE -  IF CHARACTER WORKSPACE IS NOT BEING USED,
C                         LCHAR AND LLCHAR CAN BE IGNORED.
C  5. THE NEXT SIX INTEGER LOCATIONS CONTAIN VALUES DESCRIBING THE
C     AMOUNT OF STORAGE ALLOCATED BY THE ALLOCATION SYSTEM TO THE
C     VARIOUS DATA TYPES.
C      WKSP(11) - NUMERIC STORAGE UNITS ALLOCATED TO LOGICAL
C      WKSP(12) - NUMERIC STORAGE UNITS ALLOCATED TO INTEGER
C      WKSP(13) - NUMERIC STORAGE UNITS ALLOCATED TO REAL
C      WKSP(14) - NUMERIC STORAGE UNITS ALLOCATED TO DOUBLE PRECISION
C      WKSP(15) - NUMERIC STORAGE UNITS ALLOCATED TO COMPLEX
C      WKSP(16) - NUMERIC STORAGE UNITS ALLOCATED TO DOUBLE COMPLEX
C
C  COPYRIGHT:  1984 BY IMSL, INC. ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      INTEGER FUNCTION I1KGT (NELMTS, ITYPE)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    NELMTS, ITYPE
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, IDUMAL, IGAP, ILEFT, IPA, IPA7, ISA, ISA7,
     &           ISIZE(6), JTYPE, LALC, LBND, LBOOK, LMAX, LNEED,
     &           LNEED1, LNOW, LOUT, LUSED
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      LOGICAL    FIRST
      SAVE       FIRST
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                              SPECIFICATIONS FOR COMMON /ERCOM8/
      INTEGER    PROLVL, XXLINE(10), XXPLEN(10), ICALOC(10), INALOC(10)
      COMMON     /ERCOM8/ PROLVL, XXLINE, XXPLEN, ICALOC, INALOC
      SAVE       /ERCOM8/
C                              SPECIFICATIONS FOR COMMON /ERCOM9/
      CHARACTER  XXPROC(10)*31
      COMMON     /ERCOM9/ XXPROC
      SAVE       /ERCOM9/
C                                  SPECIFICATIONS FOR COMMON /WORKSP/
      REAL       RWKSP(61913)
      REAL       RDWKSP(5000)
      DOUBLE PRECISION DWKSP(2500)
      COMPLEX    CWKSP(2500)
      COMPLEX    CZWKSP(2500)
      COMPLEX    *16 ZWKSP(1250)
      INTEGER    IWKSP(5000)
      LOGICAL    LWKSP(5000)
      EQUIVALENCE (DWKSP(1), RWKSP(1))
      EQUIVALENCE (CWKSP(1), RWKSP(1)), (ZWKSP(1), RWKSP(1))
      EQUIVALENCE (IWKSP(1), RWKSP(1)), (LWKSP(1), RWKSP(1))
      EQUIVALENCE (RDWKSP(1), RWKSP(1)), (CZWKSP(1), RWKSP(1))
      COMMON     /WORKSP/ RWKSP
C                                  SPECIFICATIONS FOR EQUIVALENCE
      EQUIVALENCE (LOUT, IWKSP(1))
      EQUIVALENCE (LNOW, IWKSP(2))
      EQUIVALENCE (LUSED, IWKSP(3))
      EQUIVALENCE (LBND, IWKSP(4))
      EQUIVALENCE (LMAX, IWKSP(5))
      EQUIVALENCE (LALC, IWKSP(6))
      EQUIVALENCE (LNEED, IWKSP(7))
      EQUIVALENCE (LBOOK, IWKSP(8))
      EQUIVALENCE (ISIZE(1), IWKSP(11))
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  IABS,MAX0,MOD
      INTRINSIC  IABS, MAX0, MOD
      INTEGER    IABS, MAX0, MOD
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   E1MES, E1POP, E1POS, E1PSH, E1STI, IWKIN
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   I1KQU
      INTEGER    I1KQU
C
      DATA FIRST/.TRUE./
C
      CALL E1PSH ('I1KGT ')
C
      IF (FIRST) THEN
C                                  INITIALIZE WORKSPACE IF NEEDED
         FIRST = .FALSE.
         CALL IWKIN (0)
      END IF
C                                  NUMBER OF ELEMENTS LESS THAN 0
      IF (NELMTS .LT. 0) THEN
         CALL E1STI (1, NELMTS)
         CALL E1MES (5, 2, 'NUMBER OF ELEMENTS IS NOT POSITIVE.%/'//
     &               'NELMTS = %(I1).')
         CALL E1POP ('I1KGT ')
         GO TO 9000
      END IF
C                                  ILLEGAL DATA TYPE REQUESTED
      IF (ITYPE.EQ.0 .OR. IABS(ITYPE).GE.7) THEN
         CALL E1MES (5, 3, 'ILLEGAL DATA TYPE REQUESTED.')
         CALL E1POP ('I1KGT ')
         GO TO 9000
      END IF
C                                  BOOKKEEPING OVERWRITTEN
      IF (LNOW.LT.LBOOK .OR. LNOW.GT.LUSED .OR. LUSED.GT.LMAX .OR.
     &    LNOW.GE.LBND .OR. LOUT.GT.LALC) THEN
         CALL E1MES (5, 4, 'ONE OR MORE OF THE FIRST EIGHT '//
     &               'BOOKKEEPING LOCATIONS IN IWKSP HAVE BEEN '//
     &               'OVERWRITTEN.')
         CALL E1POP ('I1KGT ')
         GO TO 9000
      END IF
C
      CALL E1POP ('I1KGT ')
C                                  DETERMINE NUMBER OF LOCATIONS STILL
C                                  AVAILABLE FOR DATA TYPE ITYPE
C                                  NOTE: I1KQU ALLOWS FOR 2 INTEGER
C                                        POINTERS WHICH MUST BE HANDLED
C                                        ARTIFICIALLY IF ILEFT = 0.
      ILEFT = I1KQU(IABS(ITYPE))
C
      IF (ITYPE .GT. 0) THEN
C                                  RELEASABLE STORAGE
         IF (ILEFT .GE. NELMTS) THEN
            I1KGT = (LNOW*ISIZE(2)-1)/ISIZE(ITYPE) + 2
            I = ((I1KGT-1+NELMTS)*ISIZE(ITYPE)-1)/ISIZE(2) + 3
C                                  IWKSP(I-1) CONTAINS THE DATA TYPE FOR
C                                  THIS ALLOCATION. IWKSP(I) CONTAINS
C                                  LNOW FOR THE PREVIOUS ALLOCATION.
            IWKSP(I-1) = ITYPE
            IWKSP(I) = LNOW
            LOUT = LOUT + 1
            LALC = LALC + 1
            LNOW = I
            LUSED = MAX0(LUSED,LNOW)
            LNEED = 0
         ELSE
C                                  RELEASABLE STORAGE WAS REQUESTED
C                                  BUT THE STACK WOULD OVERFLOW.
C                                  THEREFORE, ALLOCATE RELEASABLE
C                                  SPACE THROUGH THE END OF THE STACK
            IF (LNEED .EQ. 0) THEN
               IDUMAL = (LNOW*ISIZE(2)-1)/ISIZE(ITYPE) + 2
               I = ((IDUMAL-1+ILEFT)*ISIZE(ITYPE)-1)/ISIZE(2) + 3
C                                  ADVANCE COUNTERS AND STORE POINTERS
C                                  IF THERE IS ROOM TO DO SO
               IF (I .LT. LBND) THEN
C                                  IWKSP(I-1) CONTAINS THE DATA TYPE FOR
C                                  THIS ALLOCATION. IWKSP(I) CONTAINS
C                                  LNOW FOR THE PREVIOUS ALLOCATION.
                  IWKSP(I-1) = ITYPE
                  IWKSP(I) = LNOW
                  LOUT = LOUT + 1
                  LALC = LALC + 1
                  LNOW = I
                  LUSED = MAX0(LUSED,LNOW)
               END IF
            END IF
C                                  CALCULATE AMOUNT NEEDED TO ACCOMODATE
C                                  THIS ALLOCATION REQUEST
            LNEED1 = (NELMTS-ILEFT)*ISIZE(ITYPE)
            IF (ILEFT .EQ. 0) THEN
               IGAP = ISIZE(ITYPE) - MOD(LNOW+LNEED,ISIZE(ITYPE))
               IF (IGAP .EQ. ISIZE(ITYPE)) IGAP = 0
               LNEED1 = LNEED1 + 2*ISIZE(2) + IGAP
            END IF
C                                  MODIFY LNEED ACCORDING TO THE SIZE
C                                  OF THE BASE BEING USED (D.P. HERE)
            LNEED = LNEED + ((LNEED1+ISIZE(3)-1)/ISIZE(3))
C                                  SINCE CURRENT ALLOCATION IS ILLEGAL,
C                                  RETURN THE NEGATIVE OF THE ADDITIONAL
C                                  AMOUNT NEEDED TO MAKE IT LEGAL
            I1KGT = -LNEED
         END IF
      ELSE
C                                  PERMANENT STORAGE
         IF (ILEFT .GE. NELMTS) THEN
            JTYPE = -ITYPE
            I1KGT = (LBND*ISIZE(2)-1)/ISIZE(JTYPE) + 1 - NELMTS
            I = ((I1KGT-1)*ISIZE(JTYPE))/ISIZE(2) - 1
C                                  IWKSP(I) CONTAINS LBND FOR PREVIOUS
C                                  PERMANENT STORAGE ALLOCATION.
C                                  IWKSP(I+1) CONTAINS THE DATA TYPE FOR
C                                  THIS ALLOCATION.
            IWKSP(I) = LBND
            IWKSP(I+1) = JTYPE
            LALC = LALC + 1
            LBND = I
            LNEED = 0
         ELSE
C                                  PERMANENT STORAGE WAS REQUESTED
C                                  BUT THE STACK WOULD OVERFLOW,
C                                  THEREFORE, ALLOCATE RELEASABLE
C                                  SPACE THROUGH THE END OF THE STACK
            IF (LNEED .EQ. 0) THEN
               JTYPE = -ITYPE
               IDUMAL = (LNOW*ISIZE(2)-1)/ISIZE(JTYPE) + 2
               I = ((IDUMAL-1+ILEFT)*ISIZE(JTYPE)-1)/ISIZE(2) + 3
C                                  ADVANCE COUNTERS AND STORE POINTERS
C                                  IF THERE IS ROOM TO DO SO
               IF (I .LT. LBND) THEN
C                                  IWKSP(I-1) CONTAINS THE DATA TYPE FOR
C                                  THIS ALLOCATION. IWKSP(I) CONTAINS
C                                  LNOW FOR THE PREVIOUS ALLOCATION.
                  IWKSP(I-1) = JTYPE
                  IWKSP(I) = LNOW
                  LOUT = LOUT + 1
                  LALC = LALC + 1
                  LNOW = I
                  LUSED = MAX0(LUSED,LNOW)
               END IF
            END IF
C                                  CALCULATE AMOUNT NEEDED TO ACCOMODATE
C                                  THIS ALLOCATION REQUEST
            LNEED1 = (NELMTS-ILEFT)*ISIZE(-ITYPE)
            IF (ILEFT .EQ. 0) THEN
               IGAP = ISIZE(-ITYPE) - MOD(LNOW+LNEED,ISIZE(-ITYPE))
               IF (IGAP .EQ. ISIZE(-ITYPE)) IGAP = 0
               LNEED1 = LNEED1 + 2*ISIZE(2) + IGAP
            END IF
C                                  MODIFY LNEED ACCORDING TO THE SIZE
C                                  OF THE BASE BEING USED (D.P. HERE)
            LNEED = LNEED + ((LNEED1+ISIZE(3)-1)/ISIZE(3))
C                                  SINCE CURRENT ALLOCATION IS ILLEGAL,
C                                  RETURN THE NEGATIVE OF THE ADDITIONAL
C                                  AMOUNT NEEDED TO MAKE IT LEGAL
            I1KGT = -LNEED
         END IF
      END IF
C                                  STACK OVERFLOW - UNRECOVERABLE ERROR
 9000 IF (LNEED .GT. 0) THEN
         CALL E1POS (-5, IPA, ISA)
         CALL E1POS (5, 0, 0)
         CALL E1POS (-7, IPA7, ISA7)
         CALL E1POS (7, 0, 0)
         CALL E1PSH ('I1KGT ')
         CALL E1STI (1, LNEED+(LMAX/ISIZE(3)))
         IF (XXLINE(PROLVL).GE.1 .AND. XXLINE(PROLVL).LE.999) THEN
            CALL E1MES (7, 1, 'INSUFFICIENT WORKSPACE FOR CURRENT '//
     &                  'ALLOCATION(S).  CORRECT BY INSERTING THE '//
     &                  'FOLLOWING PROTRAN LINE: $OPTIONS;WORKSPACE=%'//
     &                  '(I1)')
         ELSE
            CALL E1MES (5, 5, 'INSUFFICIENT WORKSPACE FOR CURRENT '//
     &                  'ALLOCATION(S). CORRECT BY CALLING IWKIN '//
     &                  'FROM MAIN PROGRAM WITH THE THREE FOLLOWING '//
     &                  'STATEMENTS:  (REGARDLESS OF PRECISION)%/'//
     &                  '      COMMON /WORKSP/  RWKSP%/      REAL '//
     &                  'RWKSP(%(I1))%/      CALL IWKIN(%(I1))')
         END IF
         CALL E1POP ('I1KGT ')
         CALL E1POS (5, IPA, ISA)
         CALL E1POS (7, IPA7, ISA7)
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  N1RCD
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    MARCH 6, 1984
C
C  PURPOSE:    RETRIEVE AN ERROR CODE.
C
C  USAGE:      N1RCD(IOPT)
C
C  ARGUMENTS:
C     IOPT   - INTEGER SPECIFYING THE LEVEL.  (INPUT)
C              IF IOPT=0 THE ERROR CODE FOR THE CURRENT LEVEL IS
C              RETURNED.  IF IOPT=1 THE ERROR CODE FOR THE MOST
C              RECENTLY CALLED ROUTINE (LAST POP) IS RETURNED.
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      INTEGER FUNCTION N1RCD (IOPT)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    IOPT
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                              SPECIFICATIONS FOR COMMON /ERCOM1/
      INTEGER    CALLVL, MAXLEV, MSGLEN, ERTYPE(51), ERCODE(51),
     &           PRINTB(7), STOPTB(7), PLEN, IFERR6, IFERR7,
     &           IALLOC(51), HDRFMT(7), TRACON(7)
      COMMON     /ERCOM1/ CALLVL, MAXLEV, MSGLEN, ERTYPE, ERCODE,
     &           PRINTB, STOPTB, PLEN, IFERR6, IFERR7, IALLOC, HDRFMT,
     &           TRACON
      SAVE       /ERCOM1/
C                              SPECIFICATIONS FOR COMMON /ERCOM2/
      CHARACTER  MSGSAV(255), PLIST(300), RNAME(51)*6
      COMMON     /ERCOM2/ MSGSAV, PLIST, RNAME
      SAVE       /ERCOM2/
C                              SPECIFICATIONS FOR COMMON /ERCOM3/
      DOUBLE PRECISION ERCKSM
      COMMON     /ERCOM3/ ERCKSM
      SAVE       /ERCOM3/
C                              SPECIFICATIONS FOR COMMON /ERCOM4/
      LOGICAL    ISUSER(51)
      COMMON     /ERCOM4/ ISUSER
      SAVE       /ERCOM4/
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   E1PRT, M1VECH
C
      IF (IOPT.NE.0 .AND. IOPT.NE.1) THEN
         ERTYPE(CALLVL) = 5
         ERCODE(CALLVL) = 1
         MSGLEN = 47
         CALL M1VECH ('.  THE ARGUMENT PASSED TO N1RCD MUST BE 0 OR '//
     &                '1. ', MSGLEN, MSGSAV, MSGLEN)
         CALL E1PRT
         STOP
      ELSE
         N1RCD = ERCODE(CALLVL+IOPT)
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  SASUM (SINGLE PRECISION VERSION)
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    AUGUST 9, 1986
C
C  PURPOSE:    SUM THE ABSOLUTE VALUES OF THE COMPONENTS OF A
C              SINGLE PRECISION VECTOR.
C
C  USAGE:      SASUM(N, SX, INCX)
C
C  ARGUMENTS:
C     N      - LENGTH OF VECTORS X.  (INPUT)
C     SX     - REAL VECTOR OF LENGTH N*INCX.  (INPUT)
C     INCX   - DISPLACEMENT BETWEEN ELEMENTS OF SX.  (INPUT)
C              X(I) IS DEFINED TO BE SX(1+(I-1)*INCX). INCX MUST BE
C              GREATER THAN 0.
C     SASUM  - SINGLE PRECISION SUM FROM I=1 TO N OF ABS(X(I)).
C              (OUTPUT)
C              X(I) REFERS TO A SPECIFIC ELEMENT OF SX.
C
C  GAMS:       D1A
C
C  CHAPTERS:   MATH/LIBRARY BASIC MATRIX/VECTOR OPERATIONS
C              STAT/LIBRARY MATHEMATICAL SUPPORT
C
C  COPYRIGHT:  1986 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      REAL FUNCTION SASUM (N, SX, INCX)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    N, INCX
      REAL       SX(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, M, MP1, NINCX
C                                  SPECIFICATIONS FOR SPECIAL CASES
C     INTRINSIC  MOD
      INTRINSIC  MOD
      INTEGER    MOD
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  ABS
      INTRINSIC  ABS
      REAL       ABS
C
      SASUM = 0.0E0
      IF (N .GT. 0) THEN
         IF (INCX .NE. 1) THEN
C                                  CODE FOR INCREMENT NOT EQUAL TO 1
            NINCX = N*INCX
            DO 10  I=1, NINCX, INCX
               SASUM = SASUM + ABS(SX(I))
   10       CONTINUE
         ELSE
C                                  CODE FOR INCREMENT EQUAL TO 1
            M = MOD(N,6)
C                                  CLEAN-UP LOOP
            DO 30  I=1, M
               SASUM = SASUM + ABS(SX(I))
   30       CONTINUE
            MP1 = M + 1
            DO 40  I=MP1, N, 6
               SASUM = SASUM + ABS(SX(I)) + ABS(SX(I+1)) +
     &                 ABS(SX(I+2)) + ABS(SX(I+3)) + ABS(SX(I+4)) +
     &                 ABS(SX(I+5))
   40       CONTINUE
         END IF
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  SCOPY (SINGLE PRECISION VERSION)
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    AUGUST 9, 1986
C
C  PURPOSE:    COPY A VECTOR X TO A VECTOR Y, BOTH SINGLE PRECISION.
C
C  USAGE:      CALL SCOPY (N, SX, INCX, SY, INCY)
C
C  ARGUMENTS:
C     N      - LENGTH OF VECTORS X AND Y.  (INPUT)
C     SX     - REAL VECTOR OF LENGTH MAX(N*IABS(INCX),1).  (INPUT)
C     INCX   - DISPLACEMENT BETWEEN ELEMENTS OF SX.  (INPUT)
C              X(I) IS DEFINED TO BE.. SX(1+(I-1)*INCX) IF INCX .GE. 0
C              OR SX(1+(I-N)*INCX) IF INCX .LT. 0.
C     SY     - REAL VECTOR OF LENGTH MAX(N*IABS(INCY),1).  (OUTPUT)
C              SCOPY COPIES X(I) TO Y(I) FOR I=1,...,N. X(I) AND Y(I)
C              REFER TO SPECIFIC ELEMENTS OF SX AND SY, RESPECTIVELY.
C              SEE INCX AND INCY ARGUMENT DESCRIPTIONS.
C     INCY   - DISPLACEMENT BETWEEN ELEMENTS OF SY.  (INPUT)
C              Y(I) IS DEFINED TO BE.. SY(1+(I-1)*INCY) IF INCY .GE. 0
C              OR SY(1+(I-N)*INCY) IF INCY .LT. 0.
C
C  GAMS:       D1A
C
C  CHAPTERS:   MATH/LIBRARY BASIC MATRIX/VECTOR OPERATIONS
C              STAT/LIBRARY MATHEMATICAL SUPPORT
C
C  COPYRIGHT:  1986 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE SCOPY (N, SX, INCX, SY, INCY)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    N, INCX, INCY
      REAL       SX(*), SY(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, IX, IY, M, MP1
C                                  SPECIFICATIONS FOR SPECIAL CASES
C     INTRINSIC  MOD
      INTRINSIC  MOD
      INTEGER    MOD
C
      IF (N .GT. 0) THEN
         IF (INCX.NE.1 .OR. INCY.NE.1) THEN
C                                  CODE FOR UNEQUAL INCREMENTS
            IX = 1
            IY = 1
            IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
            IF (INCY .LT. 0) IY = (-N+1)*INCY + 1
            DO 10  I=1, N
               SY(IY) = SX(IX)
               IX = IX + INCX
               IY = IY + INCY
   10       CONTINUE
         ELSE
C                                  CODE FOR BOTH INCREMENTS EQUAL TO 1
            M = MOD(N,7)
C                                  CLEAN-UP LOOP
            DO 30  I=1, M
               SY(I) = SX(I)
   30       CONTINUE
            MP1 = M + 1
            DO 40  I=MP1, N, 7
               SY(I) = SX(I)
               SY(I+1) = SX(I+1)
               SY(I+2) = SX(I+2)
               SY(I+3) = SX(I+3)
               SY(I+4) = SX(I+4)
               SY(I+5) = SX(I+5)
               SY(I+6) = SX(I+6)
   40       CONTINUE
         END IF
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  SGEMV  (SINGLE PRECISION VERSION)
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    JULY 15, 1986
C
C  PURPOSE:    PERFORM THE MATRIX-VECTOR MULTIPLICATION
C              Y = ALPHA*A*X + BETA*Y OR Y = ALPHA*A'*X + BETA*Y,
C              ALL SINGLE PRECISION.
C
C  USAGE:      CALL SGEMV (TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y,
C                          INCY)
C
C  ARGUMENTS:
C     TRANS  - CHARACTER SPECIFING THE OPERATION TO BE PERFORMED.
C              (INPUT)
C                 TRANS               OPERATION
C              'N' OR 'N'      Y = ALPHA*A*X + BETA*Y
C              'T' OR 'T'      Y = ALPHA*A'*X + BETA*Y
C              'C' OR 'C'      Y = ALPHA*A'*X + BETA*Y
C     M      - NUMBER OF ROWS IN A.  (INPUT)
C     N      - NUMBER OF COLUMNS IN A.  (INPUT)
C     ALPHA  - SCALAR.  (INPUT)
C     A      - ARRAY OF SIZE M BY N.  (INPUT)
C     LDA    - LEADING DIMENSION OF A EXACTLY AS SPECIFIED IN THE
C              CALLING ROUTINE.  (INPUT)
C     X      - VECTOR OF LENGTH (N-1)*IABS(INCX)+1 WHEN TRANS IS
C              'N' OR 'N' AND OF LENGTH (M-1)*IABS(INCX)+1 OTHERWISE.
C              (INPUT)
C     INCX   - DISPLACEMENT BETWEEN ELEMENTS OF X.  (INPUT)
C     BETA   - SCALAR.  (INPUT)
C              WHEN BETA IS ZERO, Y IS NOT REFERENCED.
C     Y      - VECTOR OF LENGTH (N-1)*IABS(INCY)+1 WHEN TRANS IS
C              'M' OR 'M' AND OF LENGTH (M-1)*IABS(INCY)+1 OTHERWISE.
C              (INPUT/OUTPUT)
C     INCY   - DISPLACEMENT BETWEEN ELEMENTS OF Y.  (INPUT)
C
C  GAMS:       D1B
C
C  CHAPTER:    MATH/LIBRARY BASIC MATRIX/VECTOR OPERATIONS
C
C  COPYRIGHT:  1986 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE SGEMV (TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y,
     &                  INCY)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    M, N, LDA, INCX, INCY
      REAL       ALPHA, BETA, X(*), Y(*)
      CHARACTER  TRANS*1
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, IX, IY, KY, LENX, LENY
C                                  SPECIFICATIONS FOR SPECIAL CASES
      REAL       A(*)
      EXTERNAL   SAXPY, SDOT
      INTEGER    KX
      REAL       SDOT
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  IABS
      INTRINSIC  IABS
      INTEGER    IABS
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   SSCAL, SSET
C                                  QUICK RETURN IF POSSIBLE
      IF (M.EQ.0 .OR. N.EQ.0 .OR. (ALPHA.EQ.0.0) .AND. (BETA.EQ.1.0))
     &    GO TO 9000
C
      IF (TRANS.EQ.'N' .OR. TRANS.EQ.'N') THEN
         LENX = N
         LENY = M
      ELSE
         LENX = M
         LENY = N
      END IF
C
      IX = 1
      IY = 1
      IF (INCX .LT. 0) IX = (-LENX+1)*INCX + 1
      IF (INCY .LT. 0) IY = (-LENY+1)*INCY + 1
C
      IF (BETA .EQ. 1) THEN
      ELSE IF (INCY .EQ. 0) THEN
         IF (BETA .EQ. 0.0) THEN
            Y(1) = 0.0
         ELSE
            Y(1) = BETA**LENY*Y(1)
         END IF
      ELSE IF (BETA .EQ. 0.0) THEN
         CALL SSET (LENY, 0.0, Y, IABS(INCY))
      ELSE
         CALL SSCAL (LENY, BETA, Y, IABS(INCY))
      END IF
C
      IF (ALPHA .EQ. 0.0) GO TO 9000
C                                  NOT TRANSPOSE
      IF (TRANS.EQ.'N' .OR. TRANS.EQ.'N') THEN
         KX = IX
         DO 10  I=1, N
            CALL SAXPY (M, ALPHA*X(KX), A(LDA*(I-1)+1), 1, Y, INCY)
            KX = KX + INCX
   10    CONTINUE
      ELSE
C                                  TRANSPOSE
         KY = IY
         DO 20  I=1, N
            Y(KY) = Y(KY) + ALPHA*SDOT(M,A(LDA*(I-1)+1),1,X,INCX)
            KY = KY + INCY
   20    CONTINUE
      END IF
C
 9000 RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  SGER  (SINGLE PRECISION VERSION)
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    JULY 15, 1986
C
C  PURPOSE:    PERFORM THE RANK-ONE MATRIX UPDATE A = ALPHA*X*Y' + A,
C              ALL SINGLE PRECISION.
C
C  USAGE:      CALL SGER (M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
C
C  ARGUMENTS:
C     M      - NUMBER OF ROWS IN A.  (INPUT)
C     N      - NUMBER OF COLUMNS IN A.  (INPUT)
C     ALPHA  - REAL SCALAR.  (INPUT)
C     X      - REAL VECTOR OF LENGTH (M-1)*IABS(INCX)+1.  (INPUT)
C     INCX   - DISPLACEMENT BETWEEN ELEMENTS OF X.  (INPUT)
C     Y      - REAL VECTOR OF LENGTH (N-1)*IABS(INCY)+1.  (INPUT)
C     INCY   - DISPLACEMENT BETWEEN ELEMENTS OF Y.  (INPUT)
C     A      - ARRAY OF SIZE M BY N.  (INPUT/OUTPUT)
C              ON INPUT, A CONTAINS THE MATRIX TO BE UPDATED.
C              ON OUTPUT, A CONTAINS THE UPDATED MATRIX.
C     LDA    - LEADING DIMENSION OF A EXACTLY AS SPECIFIED IN THE
C              CALLING ROUTINE.  (INPUT)
C
C  GAMS:       D1B
C
C  CHAPTER:    MATH/LIBRARY BASIC MATRIX/VECTOR OPERATIONS
C
C  COPYRIGHT:  1986 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE SGER (M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    M, N, INCX, INCY, LDA
      REAL       ALPHA, X(*), Y(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    IY, J
C                                  SPECIFICATIONS FOR SPECIAL CASES
      REAL       A(*)
      INTEGER    I1X
      EXTERNAL   SAXPY
C                                  QUICK RETURN IF POSSIBLE
      IF (M.EQ.0 .OR. N.EQ.0 .OR. ALPHA.EQ.0.0) GO TO 9000
C
      IY = 1
      IF (INCY .LT. 0) IY = (-N+1)*INCY + 1
C
      I1X = 1
      DO 10  J=1, N
         CALL SAXPY (M, ALPHA*Y(IY), X, INCX, A(I1X), 1)
         IY = IY + INCY
         I1X = I1X + LDA
   10 CONTINUE
C
 9000 RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  SNRM2 (SINGLE PRECISION VERSION)
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    AUGUST 9, 1986
C
C  PURPOSE:    COMPUTE THE EUCLIDEAN LENGTH OR L2 NORM OF A
C              SINGLE-PRECISION VECTOR.
C
C  USAGE:      SNRM2(N, SX, INCX)
C
C  ARGUMENTS:
C     N      - LENGTH OF VECTOR X.  (INPUT)
C     SX     - REAL VECTOR OF LENGTH N*INCX.  (INPUT)
C     INCX   - DISPLACEMENT BETWEEN ELEMENTS OF SX.  (INPUT)
C              X(I) IS DEFINED TO BE SX(1+(I-1)*INCX). INCX MUST BE
C              GREATER THAN ZERO.
C     SNRM2  - SQUARE ROOT OF THE SUM FROM I=1 TO N OF X(I)**2.
C              (OUTPUT)
C              X(I) REFERS TO A SPECIFIC ELEMENT OF SX. SEE INCX
C              ARGUMENT DESCRIPTION.
C
C  GAMS:       D1A3B
C
C  CHAPTERS:   MATH/LIBRARY BASIC MATRIX/VECTOR OPERATIONS
C              STAT/LIBRARY MATHEMATICAL SUPPORT
C
C  COPYRIGHT:  1986 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      REAL FUNCTION SNRM2 (N, SX, INCX)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    N, INCX
      REAL       SX(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, J, NEXT, NN
      REAL       HITEST, SUM, XMAX
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      REAL       CUTHI, CUTLO, ONE, ZERO
      SAVE       CUTHI, CUTLO, ONE, ZERO
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  ABS,SQRT
      INTRINSIC  ABS, SQRT
      REAL       ABS, SQRT
C
      DATA ZERO/0.0E0/, ONE/1.0E0/
      DATA CUTLO/4.441E-16/, CUTHI/1.304E19/
C
      IF (N .GT. 0) GO TO 10
      SNRM2 = ZERO
      GO TO 140
C
   10 ASSIGN 30 TO NEXT
      SUM = ZERO
      NN = N*INCX
C                                  BEGIN MAIN LOOP
      I = 1
   20 GO TO NEXT, (30, 40, 70, 80)
   30 IF (ABS(SX(I)) .GT. CUTLO) GO TO 110
      ASSIGN 40 TO NEXT
      XMAX = ZERO
C                                  PHASE 1. SUM IS ZERO
   40 IF (SX(I) .EQ. ZERO) GO TO 130
      IF (ABS(SX(I)) .GT. CUTLO) GO TO 110
C                                  PREPARE FOR PHASE 2.
      ASSIGN 70 TO NEXT
      GO TO 60
C                                  PREPARE FOR PHASE 4.
   50 I = J
      ASSIGN 80 TO NEXT
      SUM = (SUM/SX(I))/SX(I)
   60 XMAX = ABS(SX(I))
      GO TO 90
C                                  PHASE 2. SUM IS SMALL. SCALE TO
C                                  AVOID DESTRUCTIVE UNDERFLOW.
   70 IF (ABS(SX(I)) .GT. CUTLO) GO TO 100
C                                  COMMON CODE FOR PHASES 2 AND 4. IN
C                                  PHASE 4 SUM IS LARGE. SCALE TO
C                                  AVOID OVERFLOW.
   80 IF (ABS(SX(I)) .LE. XMAX) GO TO 90
      SUM = ONE + SUM*(XMAX/SX(I))**2
      XMAX = ABS(SX(I))
      GO TO 130
C
   90 SUM = SUM + (SX(I)/XMAX)**2
      GO TO 130
C                                  PREPARE FOR PHASE 3.
  100 SUM = (SUM*XMAX)*XMAX
C                                  FOR REAL OR D.P. SET HITEST =
C                                  CUTHI/N FOR COMPLEX SET HITEST =
C                                  CUTHI/(2*N)
  110 HITEST = CUTHI/N
C                                  PHASE 3. SUM IS MID-RANGE. NO
C                                  SCALING.
      DO 120  J=I, NN, INCX
         IF (ABS(SX(J)) .GE. HITEST) GO TO 50
  120 SUM = SUM + SX(J)*SX(J)
      SNRM2 = SQRT(SUM)
      GO TO 140
C
  130 CONTINUE
      I = I + INCX
      IF (I .LE. NN) GO TO 20
C                                  END OF MAIN LOOP. COMPUTE SQUARE
C                                  ROOT AND ADJUST FOR SCALING.
      SNRM2 = XMAX*SQRT(SUM)
  140 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  SROT (SINGLE PRECISION VERSION)
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    AUGUST 9, 1986
C
C  PURPOSE:    APPLY A GIVENS PLANE ROTATION IN SINGLE PRECISION.
C
C  USAGE:      CALL SROT (N, SX, INCX, SY, INCY, C, S)
C
C  ARGUMENTS:
C     N      - LENGTH OF VECTORS X AND Y.  (INPUT)
C     SX     - REAL VECTOR OF LENGTH MAX(N*IABS(INCX),1).
C              (INPUT/OUTPUT)
C              SROT REPLACES X(I) WITH SC*X(I) + SS*Y(I) FOR I=1,...,N.
C              X(I) AND Y(I) REFER TO SPECIFIC ELEMENTS OF SX AND SY.
C     INCX   - DISPLACEMENT BETWEEN ELEMENTS OF SX.  (INPUT)
C              X(I) IS DEFINED TO BE
C                 SX(1+(I-1)*INCX) IF INCX.GE.0  OR
C                 SX(1+(I-N)*INCX) IF INCX.LT.0.
C     SY     - REAL VECTOR OF LENGTH MAX(N*IABS(INCY),1).
C              (INPUT/OUTPUT)
C              SROT REPLACES Y(I) WITH -SS*X(I) + SC*Y(I) FOR I=1,...,N.
C              X(I) AND Y(I) REFER TO SPECIFIC ELEMENTS OF SX AND SY.
C     INCY   - DISPLACEMENT BETWEEN ELEMENTS OF SY.  (INPUT)
C              Y(I) IS DEFINED TO BE
C                 SY(1+(I-1)*INCY) IF INCY.GE.0  OR
C                 SY(1+(I-N)*INCY) IF INCY.LT.0.
C     C      - REAL SCALAR CONTAINING ELEMENTS OF THE ROTATION MATRIX.
C              (INPUT)
C     S      - REAL SCALAR CONTAINING ELEMENTS OF THE ROTATION MATRIX.
C              (INPUT)
C
C  REMARK:
C                  ( SC SS )    (X(1) ... X(N))
C     SROT APPLIES (       ) TO (             )
C                  (-SS SC )    (Y(1) ... Y(N))
C
C  GAMS:       D1A8
C
C  CHAPTERS:   MATH/LIBRARY BASIC MATRIX/VECTOR OPERATIONS
C              STAT/LIBRARY MATHEMATICAL SUPPORT
C
C  COPYRIGHT:  1986 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE SROT (N, SX, INCX, SY, INCY, C, S)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    N, INCX, INCY
      REAL       C, S, SX(*), SY(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, IX, IY
      REAL       STEMP
C
      IF (N .GT. 0) THEN
         IF (INCX.NE.1 .OR. INCY.NE.1) THEN
C                                  CODE FOR UNEQUAL INCREMENTS OR EQUAL
C                                    INCREMENTS NOT EQUAL TO 1
            IX = 1
            IY = 1
            IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
            IF (INCY .LT. 0) IY = (-N+1)*INCY + 1
            DO 10  I=1, N
               STEMP = C*SX(IX) + S*SY(IY)
               SY(IY) = C*SY(IY) - S*SX(IX)
               SX(IX) = STEMP
               IX = IX + INCX
               IY = IY + INCY
   10       CONTINUE
         ELSE
C                                  CODE FOR BOTH INCREMENTS EQUAL TO 1
            DO 20  I=1, N
               STEMP = C*SX(I) + S*SY(I)
               SY(I) = C*SY(I) - S*SX(I)
               SX(I) = STEMP
   20       CONTINUE
         END IF
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  SROTG (SINGLE PRECISION VERSION)
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    AUGUST 9, 1986
C
C  PURPOSE:    CONSTRUCT A GIVENS PLANE ROTATION IN SINGLE PRECISION.
C
C  USAGE:      CALL SROTG (SA, SB, SC, SS)
C
C  ARGUMENTS:
C     SA     - FIRST ELEMENT OF VECTOR.  (INPUT/OUTPUT)
C              ON OUTPUT, R = (+/-)SQRT(SA**2 + SB**2) OVERWRITES SA.
C     SB     - SECOND ELEMENT OF VECTOR.  (INPUT/OUTPUT)
C              ON OUTPUT, Z OVERWRITES SB WHERE Z IS DEFINED TO BE
C                SS        IF ABS(SA) .GT. ABS(SB)
C                1.0/SC    IF ABS(SB) .GE. ABS(SA) AND SC .NE. 0.0
C                1.0       IF SC .EQ. 0.0.
C     SC     - REAL SCALAR CONTAINING ELEMENTS OF THE ROTATION MATRIX.
C              (OUTPUT)
C     SS     - REAL SCALAR CONTAINING ELEMENTS OF THE ROTATION MATRIX.
C              (OUTPUT)
C
C  REMARK:
C     SROTG CONSTRUCTS THE GIVENS ROTATION
C                         ( SC  SS )
C                     G = (        ) , SC**2 + SS**2 = 1
C                         (-SS  SC )
C     WHICH ZEROS THE SECOND ELEMENT OF (SA SB)(TRANSPOSE).
C
C  GAMS:       D1A8
C
C  CHAPTERS:   MATH/LIBRARY BASIC MATRIX/VECTOR OPERATIONS
C              STAT/LIBRARY MATHEMATICAL SUPPORT
C
C  COPYRIGHT:  1986 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE SROTG (SA, SB, SC, SS)
C                                  SPECIFICATIONS FOR ARGUMENTS
      REAL       SA, SB, SC, SS
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      REAL       R, U, V
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  ABS,SQRT
      INTRINSIC  ABS, SQRT
      REAL       ABS, SQRT
C
      IF (ABS(SA) .GT. ABS(SB)) THEN
C                                  HERE ABS(SA) .GT. ABS(SB)
         U = SA + SA
         V = SB/U
C                                  NOTE THAT U AND R HAVE THE SIGN OF
C                                    SA
         R = SQRT(.25+V**2)*U
C                                  NOTE THAT SC IS POSITIVE
         SC = SA/R
         SS = V*(SC+SC)
         SB = SS
         SA = R
      ELSE
C                                  HERE ABS(SA) .LE. ABS(SB)
         IF (SB .NE. 0.0) THEN
            U = SB + SB
            V = SA/U
C                                  NOTE THAT U AND R HAVE THE SIGN OF
C                                    SB (R IS IMMEDIATELY STORED IN SA)
            SA = SQRT(.25+V**2)*U
C                                  NOTE THAT SS IS POSITIVE
            SS = SB/SA
            SC = V*(SS+SS)
            IF (SC .NE. 0.0) THEN
               SB = 1.0/SC
            ELSE
               SB = 1.0
            END IF
C                                  HERE SA = SB = 0.
         ELSE
            SC = 1.0
            SS = 0.0
            SA = 0.0
            SB = 0.0
         END IF
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  SSCAL (SINGLE PRECISION VERSION)
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    AUGUST 9, 1986
C
C  PURPOSE:    MULTIPLY A VECTOR BY A SCALAR, Y = AY, BOTH SINGLE
C              PRECISION.
C
C  USAGE:      CALL SSCAL (N, SA, SX, INCX)
C
C  ARGUMENTS:
C     N      - LENGTH OF VECTOR X.  (INPUT)
C     SA     - REAL SCALAR.  (INPUT)
C     SX     - REAL VECTOR OF LENGTH N*INCX.  (INPUT/OUTPUT)
C              SSCAL REPLACES X(I) WITH SA*X(I) FOR I=1,...,N. X(I)
C              REFERS TO A SPECIFIC ELEMENT OF SX. SEE INCX ARGUMENT
C              DESCRIPTION.
C     INCX   - DISPLACEMENT BETWEEN ELEMENTS OF SX.  (INPUT)
C              X(I) IS DEFINED TO BE SX(1+(I-1)*INCX). INCX MUST BE
C              GREATER THAN ZERO.
C
C  GAMS:       D1A6
C
C  CHAPTERS:   MATH/LIBRARY BASIC MATRIX/VECTOR OPERATIONS
C              STAT/LIBRARY MATHEMATICAL SUPPORT
C
C  COPYRIGHT:  1986 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE SSCAL (N, SA, SX, INCX)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    N, INCX
      REAL       SA, SX(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, M, MP1, NS
C
      IF (N .GT. 0) THEN
         IF (INCX .NE. 1) THEN
C                                  CODE FOR INCREMENTS NOT EQUAL TO 1.
            NS = N*INCX
            DO 10  I=1, NS, INCX
               SX(I) = SA*SX(I)
   10       CONTINUE
         ELSE
C                                  CODE FOR INCREMENTS EQUAL TO 1.
C                                  CLEAN-UP LOOP SO REMAINING VECTOR
C                                  LENGTH IS A MULTIPLE OF 5.
            M = N - (N/5)*5
            DO 30  I=1, M
               SX(I) = SA*SX(I)
   30       CONTINUE
            MP1 = M + 1
            DO 40  I=MP1, N, 5
               SX(I) = SA*SX(I)
               SX(I+1) = SA*SX(I+1)
               SX(I+2) = SA*SX(I+2)
               SX(I+3) = SA*SX(I+3)
               SX(I+4) = SA*SX(I+4)
   40       CONTINUE
         END IF
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  SSET (SINGLE PRECISION VERSION)
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    AUGUST 9, 1986
C
C  PURPOSE:    SET THE COMPONENTS OF A VECTOR TO A SCALAR, ALL
C              SINGLE PRECISION.
C
C  USAGE:      CALL SSET (N, SA, SX, INCX)
C
C  ARGUMENTS:
C     N      - LENGTH OF VECTOR X.  (INPUT)
C     SA     - REAL SCALAR.  (INPUT)
C     SX     - REAL VECTOR OF LENGTH N*INCX.  (INPUT/OUTPUT)
C              SSET REPLACES X(I) WITH SA FOR I=1,...,N. X(I) REFERS TO
C              A SPECIFIC ELEMENT OF SX. SEE INCX ARGUMENT DESCRIPTION.
C     INCX   - DISPLACEMENT BETWEEN ELEMENTS OF SX.  (INPUT)
C              X(I) IS DEFINED TO BE SX(1+(I-1)*INCX). INCX MUST BE
C              GREATER THAN ZERO.
C
C  GAMS:       D1A1
C
C  CHAPTERS:   MATH/LIBRARY BASIC MATRIX/VECTOR OPERATIONS
C              STAT/LIBRARY MATHEMATICAL SUPPORT
C
C  COPYRIGHT:  1986 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE SSET (N, SA, SX, INCX)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    N, INCX
      REAL       SA, SX(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, M, MP1, NINCX
C                                  SPECIFICATIONS FOR SPECIAL CASES
C     INTRINSIC  MOD
      INTRINSIC  MOD
      INTEGER    MOD
C
      IF (N .GT. 0) THEN
         IF (INCX .NE. 1) THEN
C                                  CODE FOR INCREMENT NOT EQUAL TO 1
            NINCX = N*INCX
            DO 10  I=1, NINCX, INCX
               SX(I) = SA
   10       CONTINUE
         ELSE
C                                  CODE FOR INCREMENT EQUAL TO 1
            M = MOD(N,8)
C                                  CLEAN-UP LOOP
            DO 30  I=1, M
               SX(I) = SA
   30       CONTINUE
            MP1 = M + 1
            DO 40  I=MP1, N, 8
               SX(I) = SA
               SX(I+1) = SA
               SX(I+2) = SA
               SX(I+3) = SA
               SX(I+4) = SA
               SX(I+5) = SA
               SX(I+6) = SA
               SX(I+7) = SA
   40       CONTINUE
         END IF
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  SSWAP (SINGLE PRECISION VERSION)
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    AUGUST 9, 1986
C
C  PURPOSE:    INTERCHANGE VECTORS X AND Y, BOTH SINGLE PRECISION.
C
C  USAGE:      CALL SSWAP (N, SX, INCX, SY, INCY)
C
C  ARGUMENTS:
C     N      - LENGTH OF VECTORS X AND Y.  (INPUT)
C     SX     - REAL VECTOR OF LENGTH MAX(N*IABS(INCX),1).
C              (INPUT/OUTPUT)
C     INCX   - DISPLACEMENT BETWEEN ELEMENTS OF SX.  (INPUT)
C              X(I) IS DEFINED TO BE
C                 SX(1+(I-1)*INCX) IF INCX.GE.0  OR
C                 SX(1+(I-N)*INCX) IF INCX.LT.0.
C     SY     - REAL VECTOR OF LENGTH MAX(N*IABS(INCY),1).
C              (INPUT/OUTPUT)
C     INCY   - DISPLACEMENT BETWEEN ELEMENTS OF SY.  (INPUT)
C              Y(I) IS DEFINED TO BE
C                 SY(1+(I-1)*INCY) IF INCY.GE.0  OR
C                 SY(1+(I-N)*INCY) IF INCY.LT.0.
C
C  GAMS:       D1A5
C
C  CHAPTERS:   MATH/LIBRARY BASIC MATRIX/VECTOR OPERATIONS
C              STAT/LIBRARY MATHEMATICAL SUPPORT
C
C  COPYRIGHT:  1986 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE SSWAP (N, SX, INCX, SY, INCY)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    N, INCX, INCY
      REAL       SX(*), SY(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, IX, IY, M, MP1
      REAL       STEMP
C                                  SPECIFICATIONS FOR SPECIAL CASES
C     INTRINSIC  MOD
      INTRINSIC  MOD
      INTEGER    MOD
C
      IF (N .GT. 0) THEN
         IF (INCX.NE.1 .OR. INCY.NE.1) THEN
C                                  CODE FOR UNEQUAL INCREMENTS OR EQUAL
C                                    INCREMENTS NOT EQUAL TO 1
            IX = 1
            IY = 1
            IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
            IF (INCY .LT. 0) IY = (-N+1)*INCY + 1
            DO 10  I=1, N
               STEMP = SX(IX)
               SX(IX) = SY(IY)
               SY(IY) = STEMP
               IX = IX + INCX
               IY = IY + INCY
   10       CONTINUE
         ELSE
C                                  CODE FOR BOTH INCREMENTS EQUAL TO 1
            M = MOD(N,3)
C                                  CLEAN-UP LOOP
            DO 30  I=1, M
               STEMP = SX(I)
               SX(I) = SY(I)
               SY(I) = STEMP
   30       CONTINUE
            MP1 = M + 1
            DO 40  I=MP1, N, 3
               STEMP = SX(I)
               SX(I) = SY(I)
               SY(I) = STEMP
               STEMP = SX(I+1)
               SX(I+1) = SY(I+1)
               SY(I+1) = STEMP
               STEMP = SX(I+2)
               SX(I+2) = SY(I+2)
               SY(I+2) = STEMP
   40       CONTINUE
         END IF
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  C1TCI
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    AUGUST 13, 1984
C
C  PURPOSE:    CONVERT CHARACTER STRING INTO CORRESPONDING INTEGER
C              FORM.
C
C  USAGE:      CALL C1TCI (CHRSTR, SLEN, NUM, IER)
C
C  ARGUMENTS:
C   CHRSTR  - CHARACTER ARRAY THAT CONTAINS THE NUMBER DESCRIPTION.
C             (INPUT)
C   SLEN    - LENGTH OF THE CHARACTER ARRAY.  (INPUT)
C   NUM     - THE ANSWER.  (OUTPUT)
C   IER     - COMPLETION CODE.  (OUTPUT)  WHERE
C                IER =-2  INDICATES THAT THE NUMBER IS TOO LARGE TO
C                         BE CONVERTED;
C                IER =-1  INDICATES THAT SLEN <= 0;
C                IER = 0  INDICATES NORMAL COMPLETION;
C                IER > 0  INDICATES THAT THE INPUT STRING CONTAINS A
C                         NONNUMERIC CHARACTER.  IER IS THE INDEX OF
C                         THE FIRST NONNUMERIC CHARACTER IN CHRSTR.
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE C1TCI (CHRSTR, SLEN, NUM, IER)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    SLEN, NUM, IER
      CHARACTER  CHRSTR(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    COUNT, I, IMACH5, J, N, S, SIGN
      CHARACTER  ZERO
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      CHARACTER  BLANK, DIGIT*10, MINUS, PLUS
      SAVE       BLANK, DIGIT, MINUS, PLUS
C                                  SPECIFICATIONS FOR EQUIVALENCE
      EQUIVALENCE (DIGIT, ZERO)
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  INDEX
      INTRINSIC  INDEX
      INTEGER    INDEX
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   IMACH
      INTEGER    IMACH
C
      DATA DIGIT/'0123456789'/
      DATA BLANK/' '/, MINUS/'-'/, PLUS/'+'/
C
C                                  CHECK SLEN
      NUM = 0
      IF (SLEN .LE. 0) THEN
         IER = -1
         GO TO 50
      END IF
C                                  HANDLE LEADING BLANKS
      SIGN = 1
      I = 1
   10 IF (I .LE. SLEN) THEN
         IF (CHRSTR(I) .EQ. BLANK) THEN
            I = I + 1
            GO TO 10
         END IF
      ELSE
         IER = 1
         GO TO 50
      END IF
C                                  CHECK FOR SIGN, IF ANY
      S = I
      IF (CHRSTR(I) .EQ. MINUS) THEN
         SIGN = -1
         I = I + 1
      ELSE IF (CHRSTR(I) .EQ. PLUS) THEN
         I = I + 1
      END IF
   20 IF (I .LE. SLEN) THEN
         IF (CHRSTR(I) .EQ. BLANK) THEN
            I = I + 1
            GO TO 20
         END IF
      ELSE
         IER = S
         GO TO 50
      END IF
C                                  SKIP LEADING ZERO
      J = I
   30 IF (I .LE. SLEN) THEN
         IF (CHRSTR(I) .EQ. ZERO) THEN
            I = I + 1
            GO TO 30
         END IF
      ELSE
         IER = 0
         GO TO 50
      END IF
C                                  CHECK FIRST NONBLANK CHARACTER
      COUNT = 0
C                                  CHECK NUMERIC CHARACTERS
      IMACH5 = IMACH(5)
   40 N = INDEX(DIGIT,CHRSTR(I))
      IF (N .NE. 0) THEN
         COUNT = COUNT + 1
         IF (NUM .GT. ((IMACH5-N)+1)/10) THEN
            IER = -2
            GO TO 50
         ELSE
            NUM = NUM*10 - 1 + N
            I = I + 1
            IF (I .LE. SLEN) GO TO 40
         END IF
      END IF
C
      IF (COUNT .EQ. 0) THEN
         IF (I .GT. J) THEN
            IER = I
         ELSE
            IER = S
         END IF
      ELSE IF (I .GT. SLEN) THEN
         NUM = SIGN*NUM
         IER = 0
      ELSE
         NUM = SIGN*NUM
         IER = I
      END IF
C
   50 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  E1INIT
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    MARCH 13, 1984
C
C  PURPOSE:    INITIALIZATION.
C
C  USAGE:      CALL E1INIT
C
C  ARGUMENTS:  NONE
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE E1INIT
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    L
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      INTEGER    ISINIT
      SAVE       ISINIT
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                              SPECIFICATIONS FOR COMMON /ERCOM1/
      INTEGER    CALLVL, MAXLEV, MSGLEN, ERTYPE(51), ERCODE(51),
     &           PRINTB(7), STOPTB(7), PLEN, IFERR6, IFERR7,
     &           IALLOC(51), HDRFMT(7), TRACON(7)
      COMMON     /ERCOM1/ CALLVL, MAXLEV, MSGLEN, ERTYPE, ERCODE,
     &           PRINTB, STOPTB, PLEN, IFERR6, IFERR7, IALLOC, HDRFMT,
     &           TRACON
      SAVE       /ERCOM1/
C                              SPECIFICATIONS FOR COMMON /ERCOM2/
      CHARACTER  MSGSAV(255), PLIST(300), RNAME(51)*6
      COMMON     /ERCOM2/ MSGSAV, PLIST, RNAME
      SAVE       /ERCOM2/
C                              SPECIFICATIONS FOR COMMON /ERCOM3/
      DOUBLE PRECISION ERCKSM
      COMMON     /ERCOM3/ ERCKSM
      SAVE       /ERCOM3/
C                              SPECIFICATIONS FOR COMMON /ERCOM4/
      LOGICAL    ISUSER(51)
      COMMON     /ERCOM4/ ISUSER
      SAVE       /ERCOM4/
C                              SPECIFICATIONS FOR COMMON /ERCOM8/
      INTEGER    PROLVL, XXLINE(10), XXPLEN(10), ICALOC(10), INALOC(10)
      COMMON     /ERCOM8/ PROLVL, XXLINE, XXPLEN, ICALOC, INALOC
      SAVE       /ERCOM8/
C                              SPECIFICATIONS FOR COMMON /ERCOM9/
      CHARACTER  XXPROC(10)*31
      COMMON     /ERCOM9/ XXPROC
      SAVE       /ERCOM9/
C
      DATA ISINIT/0/
C
      IF (ISINIT .EQ. 0) THEN
C                                  INITIALIZE
         CALLVL = 1
         ERCODE(1) = 0
         ERTYPE(1) = 0
         IALLOC(1) = 0
         ISUSER(1) = .TRUE.
         IFERR6 = 0
         IFERR7 = 0
         PLEN = 1
         MAXLEV = 50
         DO 10  L=2, 51
            ERTYPE(L) = -1
            ERCODE(L) = -1
            IALLOC(L) = 0
            ISUSER(L) = .FALSE.
   10    CONTINUE
         DO 20  L=1, 7
            HDRFMT(L) = 1
            TRACON(L) = 1
   20    CONTINUE
         PROLVL = 1
         DO 30  L=1, 10
   30    ICALOC(L) = 0
         XXLINE(1) = 0
         XXPLEN(1) = 1
         XXPROC(1) = '?'
         RNAME(1) = 'USER'
         PRINTB(1) = 0
         PRINTB(2) = 0
         DO 40  L=3, 7
   40    PRINTB(L) = 1
         STOPTB(1) = 0
         STOPTB(2) = 0
         STOPTB(3) = 0
         STOPTB(4) = 1
         STOPTB(5) = 1
         STOPTB(6) = 0
         STOPTB(7) = 1
         ERCKSM = 0.0D0
C                                  SET FLAG TO INDICATE THAT
C                                    INITIALIZATION HAS OCCURRED
         ISINIT = 1
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  E1PRT
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    MARCH 14, 1984
C
C  PURPOSE:    TO PRINT AN ERROR MESSAGE.
C
C  USAGE:      CALL E1PRT
C
C  ARGUMENTS:  NONE
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE E1PRT
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    ALL, I, IBEG, IBLOC, IBLOC2, IEND, IER, IHDR, J,
     &           LERTYP, LOC, LOCM1, LOCX, MAXLOC, MAXTMP, MLOC, MOD,
     &           NCBEG, NLOC, NOUT
      CHARACTER  MSGTMP(70), STRING(10)
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      CHARACTER  ATLINE(9), BLANK(1), DBB(3), FROM(6), MSGTYP(8,7),
     &           PERSLA(2), QMARK, UNKNOW(8)
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                              SPECIFICATIONS FOR COMMON /ERCOM1/
      INTEGER    CALLVL, MAXLEV, MSGLEN, ERTYPE(51), ERCODE(51),
     &           PRINTB(7), STOPTB(7), PLEN, IFERR6, IFERR7,
     &           IALLOC(51), HDRFMT(7), TRACON(7)
      COMMON     /ERCOM1/ CALLVL, MAXLEV, MSGLEN, ERTYPE, ERCODE,
     &           PRINTB, STOPTB, PLEN, IFERR6, IFERR7, IALLOC, HDRFMT,
     &           TRACON
      SAVE       /ERCOM1/
C                              SPECIFICATIONS FOR COMMON /ERCOM2/
      CHARACTER  MSGSAV(255), PLIST(300), RNAME(51)*6
      COMMON     /ERCOM2/ MSGSAV, PLIST, RNAME
      SAVE       /ERCOM2/
C                              SPECIFICATIONS FOR COMMON /ERCOM3/
      DOUBLE PRECISION ERCKSM
      COMMON     /ERCOM3/ ERCKSM
      SAVE       /ERCOM3/
C                              SPECIFICATIONS FOR COMMON /ERCOM4/
      LOGICAL    ISUSER(51)
      COMMON     /ERCOM4/ ISUSER
      SAVE       /ERCOM4/
C                              SPECIFICATIONS FOR COMMON /ERCOM8/
      INTEGER    PROLVL, XXLINE(10), XXPLEN(10), ICALOC(10), INALOC(10)
      COMMON     /ERCOM8/ PROLVL, XXLINE, XXPLEN, ICALOC, INALOC
      SAVE       /ERCOM8/
C                              SPECIFICATIONS FOR COMMON /ERCOM9/
      CHARACTER  XXPROC(10)*31
      COMMON     /ERCOM9/ XXPROC
      SAVE       /ERCOM9/
      SAVE       ATLINE, BLANK, DBB, FROM, MSGTYP, PERSLA, QMARK,
     &           UNKNOW
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  MIN0
      INTRINSIC  MIN0
      INTEGER    MIN0
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   C1TIC, M1VE, UMACH
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   I1DX, I1ERIF
      INTEGER    I1DX, I1ERIF
C
      DATA MSGTYP/'N', 'O', 'T', 'E', ' ', ' ', ' ', ' ', 'A',
     &     'L', 'E', 'R', 'T', ' ', ' ', ' ', 'W', 'A', 'R',
     &     'N', 'I', 'N', 'G', ' ', 'F', 'A', 'T', 'A', 'L',
     &     ' ', ' ', ' ', 'T', 'E', 'R', 'M', 'I', 'N', 'A',
     &     'L', 'W', 'A', 'R', 'N', 'I', 'N', 'G', ' ', 'F',
     &     'A', 'T', 'A', 'L', ' ', ' ', ' '/
      DATA UNKNOW/'U', 'N', 'K', 'N', 'O', 'W', 'N', ' '/
      DATA ATLINE/' ', 'A', 'T', ' ', 'L', 'I', 'N', 'E', ' '/
      DATA BLANK/' '/, FROM/' ', 'F', 'R', 'O', 'M', ' '/
      DATA DBB/'.', ' ', ' '/, PERSLA/'%', '/'/
      DATA QMARK/'?'/
C
      IF (MSGLEN .LE. 0) RETURN
      CALL UMACH (2, NOUT)
      MAXTMP = 70
      MOD = 0
      LERTYP = ERTYPE(CALLVL)
      IHDR = HDRFMT(LERTYP)
      IF (IHDR .EQ. 3) THEN
         IF (XXPROC(PROLVL)(1:1).EQ.QMARK .AND. XXLINE(PROLVL).EQ.0)
     &       THEN
            IHDR = 1
         END IF
      END IF
      IEND = 0
      IF (IHDR.EQ.1 .AND. ERTYPE(CALLVL).LE.4) THEN
         MSGTMP(1) = BLANK(1)
         IEND = 1
C                                  CONVERT ERROR CODE INTO CHAR STRING
         CALL C1TIC (ERCODE(CALLVL), STRING, 10, IER)
C                                  LOCATE START OF NON-BLANK CHARACTERS
         IBEG = I1ERIF(STRING,10,BLANK,1)
C                                  M1VE IT TO MSGTMP
         CALL M1VE (STRING, IBEG, 10, 10, MSGTMP, IEND+1,
     &              IEND+11-IBEG, MAXTMP, IER)
         IEND = IEND + 11 - IBEG
      END IF
      IF (IHDR .NE. 2) THEN
         CALL M1VE (FROM, 1, 6, 6, MSGTMP, IEND+1, IEND+6, MAXTMP, IER)
         IEND = IEND + 6
      END IF
      IF (IHDR .EQ. 3) THEN
C                                  THIS IS A PROTRAN RUN TIME ERROR MSG.
C                                  RETRIEVE THE PROCEDURE NAME
         CALL M1VE (XXPROC(PROLVL), 1, XXPLEN(PROLVL), 31, MSGTMP,
     &              IEND+1, IEND+XXPLEN(PROLVL), MAXTMP, IER)
         MLOC = IEND + XXPLEN(PROLVL) + 1
         MSGTMP(MLOC) = BLANK(1)
         IEND = IEND + I1DX(MSGTMP(IEND+1),XXPLEN(PROLVL)+1,BLANK,1) -
     &          1
         IF (XXLINE(PROLVL) .GT. 0) THEN
C                                  INSERT ATLINE
            CALL M1VE (ATLINE, 1, 9, 9, MSGTMP, IEND+1, IEND+9,
     &                 MAXTMP, IER)
            IEND = IEND + 9
C                                  CONVERT PROTRAN GLOBAL LINE NUMBER
            CALL C1TIC (XXLINE(PROLVL), STRING, 10, IER)
C                                  LOCATE START OF NON-BLANK CHARACTERS
            IBEG = I1ERIF(STRING,10,BLANK,1)
C                                  M1VE GLOBAL LINE NUMBER TO MSGTMP
            CALL M1VE (STRING, IBEG, 10, 10, MSGTMP, IEND+1,
     &                 IEND+11-IBEG, MAXTMP, IER)
            IEND = IEND + 11 - IBEG
         END IF
      ELSE
C                                  THIS IS EITHER A LIBRARY ERROR MSG
C                                  OR A PROTRAN PREPROCESSOR ERROR MSG
         IF (IHDR .EQ. 1) THEN
C                                  THIS IS A LIBRARY ERROR MESSAGE.
C                                  RETRIEVE ROUTINE NAME
            CALL M1VE (RNAME(CALLVL), 1, 6, 6, MSGTMP, IEND+1, IEND+6,
     &                 MAXTMP, IER)
            MSGTMP(IEND+7) = BLANK(1)
            IEND = IEND + I1DX(MSGTMP(IEND+1),7,BLANK,1) - 1
         END IF
C                                  ADD DOT, BLANK, BLANK IF NEEDED
         IF (I1DX(MSGSAV,3,DBB,3) .NE. 1) THEN
            CALL M1VE (DBB, 1, 3, 3, MSGTMP, IEND+1, IEND+3, MAXTMP,
     &                 IER)
            IEND = IEND + 3
            MOD = 3
         END IF
      END IF
C                                  MSGTMP AND MSGSAV NOW CONTAIN THE
C                                   ERROR MESSAGE IN FINAL FORM.
      NCBEG = 59 - IEND - MOD
      ALL = 0
      IBLOC = I1DX(MSGSAV,MSGLEN,PERSLA,2)
      IF (IBLOC.NE.0 .AND. IBLOC.LT.NCBEG) THEN
         LOCM1 = IBLOC - 1
         LOC = IBLOC + 1
      ELSE IF (MSGLEN .LE. NCBEG) THEN
         LOCM1 = MSGLEN
         ALL = 1
      ELSE
         LOC = NCBEG
C                                  CHECK FOR APPROPRIATE PLACE TO SPLIT
   10    CONTINUE
         IF (MSGSAV(LOC) .NE. BLANK(1)) THEN
            LOC = LOC - 1
            IF (LOC .GT. 1) GO TO 10
            LOC = NCBEG + 1
         END IF
         LOCM1 = LOC - 1
      END IF
C                                  NO BLANKS FOUND IN FIRST NCBEG CHARS
      IF (LERTYP.GE.1 .AND. LERTYP.LE.7) THEN
         WRITE (NOUT,99995) (MSGTYP(I,LERTYP),I=1,8),
     &                     (MSGTMP(I),I=1,IEND), (MSGSAV(I),I=1,LOCM1)
      ELSE
         WRITE (NOUT,99995) (UNKNOW(I),I=1,8), (MSGTMP(I),I=1,IEND),
     &                     (MSGSAV(I),I=1,LOCM1)
      END IF
      IF (ALL .EQ. 0) THEN
C                                  PREPARE TO WRITE CONTINUATION OF
C                                    MESSAGE
C
C                                  FIND WHERE TO BREAK MESSAGE
C                                    LOC = NUMBER OF CHARACTERS OF
C                                          MESSAGE WRITTEN SO FAR
   20    LOCX = LOC + 64
         NLOC = LOC + 1
         IBLOC2 = IBLOC
         MAXLOC = MIN0(MSGLEN-LOC,64)
         IBLOC = I1DX(MSGSAV(NLOC),MAXLOC,PERSLA,2)
         IF (MSGSAV(NLOC).EQ.BLANK(1) .AND. IBLOC2.EQ.0) NLOC = NLOC +
     &       1
         IF (IBLOC .GT. 0) THEN
C                                  PAGE BREAK FOUND AT IBLOC
            LOCX = NLOC + IBLOC - 2
            WRITE (NOUT,99996) (MSGSAV(I),I=NLOC,LOCX)
            LOC = NLOC + IBLOC
            GO TO 20
C                                  DON'T BOTHER LOOKING FOR BLANK TO
C                                    BREAK AT IF LOCX .GE. MSGLEN
         ELSE IF (LOCX .LT. MSGLEN) THEN
C                                  CHECK FOR BLANK TO BREAK THE LINE
   30       CONTINUE
            IF (MSGSAV(LOCX) .EQ. BLANK(1)) THEN
C                                  BLANK FOUND AT LOCX
               WRITE (NOUT,99996) (MSGSAV(I),I=NLOC,LOCX)
               LOC = LOCX
               GO TO 20
            END IF
            LOCX = LOCX - 1
            IF (LOCX .GT. NLOC) GO TO 30
            LOCX = LOC + 64
C                                  NO BLANKS FOUND IN NEXT 64 CHARS
            WRITE (NOUT,99996) (MSGSAV(I),I=NLOC,LOCX)
            LOC = LOCX
            GO TO 20
         ELSE
C                                  ALL THE REST WILL FIT ON 1 LINE
            LOCX = MSGLEN
            WRITE (NOUT,99996) (MSGSAV(I),I=NLOC,LOCX)
         END IF
      END IF
C                                  SET LENGTH OF MSGSAV AND PLEN
C                                    TO SHOW THAT MESSAGE HAS
C                                    ALREADY BEEN PRINTED
 9000 MSGLEN = 0
      PLEN = 1
      IF (TRACON(LERTYP).EQ.1 .AND. CALLVL.GT.2) THEN
C                                  INITIATE TRACEBACK
         WRITE (NOUT,99997)
         DO 9005  J=CALLVL, 1, -1
            IF (J .GT. 1) THEN
               IF (ISUSER(J-1)) THEN
                  WRITE (NOUT,99998) RNAME(J), ERTYPE(J), ERCODE(J)
               ELSE
                  WRITE (NOUT,99999) RNAME(J), ERTYPE(J), ERCODE(J)
               END IF
            ELSE
               WRITE (NOUT,99998) RNAME(J), ERTYPE(J), ERCODE(J)
            END IF
 9005    CONTINUE
      END IF
C
      RETURN
99995 FORMAT (/, ' *** ', 8A1, ' ERROR', 59A1)
99996 FORMAT (' *** ', 9X, 64A1)
99997 FORMAT (14X, 'HERE IS A TRACEBACK OF SUBPROGRAM CALLS',
     &       ' IN REVERSE ORDER:', /, 14X, '      ROUTINE    ERROR ',
     &       'TYPE    ERROR CODE', /, 14X, '      -------    ',
     &       '----------    ----------')
99998 FORMAT (20X, A6, 5X, I6, 8X, I6)
99999 FORMAT (20X, A6, 5X, I6, 8X, I6, 4X, '(CALLED INTERNALLY)')
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  E1UCS
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    MARCH 8, 1984
C
C  PURPOSE:    TO UPDATE THE CHECKSUM NUMBER FOR ERROR MESSAGES.
C
C  USAGE:      CALL E1UCS
C
C  ARGUMENTS:  NONE
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE E1UCS
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, IBEG, IBEG2, IEND, ILOC, IPOS, JLOC, NCODE, NLEN
      DOUBLE PRECISION DNUM
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      DOUBLE PRECISION DMAX
      CHARACTER  BLANK(1), COMMA(1), EQUAL(1), LPAR(1)
      SAVE       BLANK, COMMA, DMAX, EQUAL, LPAR
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                              SPECIFICATIONS FOR COMMON /ERCOM1/
      INTEGER    CALLVL, MAXLEV, MSGLEN, ERTYPE(51), ERCODE(51),
     &           PRINTB(7), STOPTB(7), PLEN, IFERR6, IFERR7,
     &           IALLOC(51), HDRFMT(7), TRACON(7)
      COMMON     /ERCOM1/ CALLVL, MAXLEV, MSGLEN, ERTYPE, ERCODE,
     &           PRINTB, STOPTB, PLEN, IFERR6, IFERR7, IALLOC, HDRFMT,
     &           TRACON
      SAVE       /ERCOM1/
C                              SPECIFICATIONS FOR COMMON /ERCOM2/
      CHARACTER  MSGSAV(255), PLIST(300), RNAME(51)*6
      COMMON     /ERCOM2/ MSGSAV, PLIST, RNAME
      SAVE       /ERCOM2/
C                              SPECIFICATIONS FOR COMMON /ERCOM3/
      DOUBLE PRECISION ERCKSM
      COMMON     /ERCOM3/ ERCKSM
      SAVE       /ERCOM3/
C                              SPECIFICATIONS FOR COMMON /ERCOM4/
      LOGICAL    ISUSER(51)
      COMMON     /ERCOM4/ ISUSER
      SAVE       /ERCOM4/
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  DMOD
      INTRINSIC  DMOD
      DOUBLE PRECISION DMOD
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   S1ANUM
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   ICASE, I1X
      INTEGER    ICASE, I1X
C
      DATA BLANK(1)/' '/, COMMA(1)/','/, LPAR(1)/'('/
      DATA EQUAL(1)/'='/, DMAX/1.0D+9/
C
      IF (MSGLEN .GT. 1) THEN
         IPOS = 0
         IBEG2 = 1
   10    IBEG = IBEG2
         IEND = MSGLEN
C                                  LOOK FOR BLANK, COMMA, LEFT PAREN.,
C                                  OR EQUAL SIGN
         ILOC = I1X(MSGSAV(IBEG),IEND-IBEG+1,BLANK,1)
         JLOC = I1X(MSGSAV(IBEG),IEND-IBEG+1,COMMA,1)
         IF (ILOC.EQ.0 .OR. (JLOC.GT.0.AND.JLOC.LT.ILOC)) ILOC = JLOC
         JLOC = I1X(MSGSAV(IBEG),IEND-IBEG+1,LPAR,1)
         IF (ILOC.EQ.0 .OR. (JLOC.GT.0.AND.JLOC.LT.ILOC)) ILOC = JLOC
         JLOC = I1X(MSGSAV(IBEG),IEND-IBEG+1,EQUAL,1)
         IF (ILOC.EQ.0 .OR. (JLOC.GT.0.AND.JLOC.LT.ILOC)) ILOC = JLOC
         IF (ILOC .GE. 1) THEN
            CALL S1ANUM (MSGSAV(IBEG+ILOC), IEND-IBEG-ILOC+1, NCODE,
     &                   NLEN)
            IF (NCODE.EQ.2 .OR. NCODE.EQ.3) THEN
C                                  FLOATING POINT NUMBER FOUND.
C                                  SET POINTERS TO SKIP OVER IT
               IBEG2 = IBEG + ILOC + NLEN
               IF (IBEG2 .LE. MSGLEN) THEN
                  CALL S1ANUM (MSGSAV(IBEG2), IEND-IBEG2+1, NCODE,
     &                         NLEN)
                  IF ((MSGSAV(IBEG2).EQ.'+'.OR.MSGSAV(IBEG2).EQ.
     &                '-') .AND. NCODE.EQ.1) THEN
C                                  INTEGER IMMEDIATELY FOLLOWS A REAL AS
C                                  WITH SOME CDC NOS. LIKE 1.2345678+123
C                                  SET POINTERS TO SKIP OVER IT
                     IBEG2 = IBEG2 + NLEN
                  END IF
               END IF
            ELSE
               IBEG2 = IBEG + ILOC
            END IF
            IEND = IBEG + ILOC - 1
         END IF
C                                  UPDATE CKSUM USING PART OF MESSAGE
         DO 20  I=IBEG, IEND
            IPOS = IPOS + 1
            DNUM = ICASE(MSGSAV(I))
            ERCKSM = DMOD(ERCKSM+DNUM*IPOS,DMAX)
   20    CONTINUE
C                                  GO BACK FOR MORE IF NEEDED
         IF (IEND.LT.MSGLEN .AND. IBEG2.LT.MSGLEN) GO TO 10
C                                  UPDATE CKSUM USING ERROR TYPE
         DNUM = ERTYPE(CALLVL)
         ERCKSM = DMOD(ERCKSM+DNUM*(IPOS+1),DMAX)
C                                  UPDATE CKSUM USING ERROR CODE
         DNUM = ERCODE(CALLVL)
         ERCKSM = DMOD(ERCKSM+DNUM*(IPOS+2),DMAX)
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  M1VE
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    MARCH 5, 1984
C
C  PURPOSE:    MOVE A SUBSET OF ONE CHARACTER ARRAY TO ANOTHER.
C
C  USAGE:      CALL M1VE(INSTR, INBEG, INEND, INLEN, OUTSTR, OUTBEG,
C                         OUTEND, OUTLEN, IER)
C
C  ARGUMENTS:
C     INSTR  - SOURCE CHARACTER ARRAY.  (INPUT)
C     INBEG  - FIRST ELEMENT OF INSTR TO BE MOVED.  (INPUT)
C     INEND  - LAST ELEMENT OF INSTR TO BE MOVED.  (INPUT)
C              THE SOURCE SUBSET IS INSTR(INBEG),...,INSTR(INEND).
C     INLEN  - LENGTH OF INSTR.  (INPUT)
C     OUTSTR - DESTINATION CHARACTER ARRAY.  (OUTPUT)
C     IUTBEG - FIRST ELEMENT OF OUTSTR DESTINATION.  (INPUT)
C     IUTEND - LAST ELEMENT OF OUTSTR  DESTINATION.  (INPUT)
C              THE DESTINATION SUBSET IS OUTSRT(IUTBEG),...,
C              OUTSTR(IUTEND).
C     IUTLEN - LENGTH OF OUTSTR.  (INPUT)
C     IER    - COMPLETION CODE.  (OUTPUT)
C              IER = -2  INDICATES THAT THE INPUT PARAMETERS, INBEG,
C                        INEND, INLEN, IUTBEG, IUTEND ARE NOT
C                        CONSISTENT.  ONE OF THE CONDITIONS
C                        INBEG.GT.0, INEND.GE.INBEG, INLEN.GE.INEND,
C                        IUTBEG.GT.0, OR IUTEND.GE.IUTBEG IS NOT
C                        SATISFIED.
C              IER = -1  INDICATES THAT THE LENGTH OF OUTSTR IS
C                        INSUFFICIENT TO HOLD THE SUBSET OF INSTR.
C                        THAT IS, IUTLEN IS LESS THAN IUTEND.
C              IER =  0  INDICATES NORMAL COMPLETION
C              IER >  0  INDICATES THAT THE SPECIFIED SUBSET OF OUTSTR,
C                        OUTSTR(IUTBEG),...,OUTSTR(IUTEND) IS NOT LONG
C                        ENOUGH TO HOLD THE SUBSET INSTR(INBEG),...,
C                        INSTR(INEND) OF INSTR.  IER IS SET TO THE
C                        NUMBER OF CHARACTERS THAT WERE NOT MOVED.
C
C  REMARKS:
C  1. IF THE SUBSET OF OUTSTR IS LONGER THAN THE SUBSET OF INSTR,
C     TRAILING BLANKS ARE MOVED TO OUTSTR.
C  2. IF THE SUBSET OF INSTR IS LONGER THAN THE SUBSET OF OUTSTR,
C     THE SHORTER SUBSET IS MOVED TO OUTSTR AND IER IS SET TO THE NUMBER
C     OF CHARACTERS THAT WERE NOT MOVED TO OUTSTR.
C  3. IF THE LENGTH OF OUTSTR IS INSUFFICIENT TO HOLD THE SUBSET,
C     IER IS SET TO -2 AND NOTHING IS MOVED.
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE M1VE (INSTR, INBEG, INEND, INLEN, OUTSTR, IUTBEG,
     &                 IUTEND, IUTLEN, IER)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    INBEG, INEND, INLEN, IUTBEG, IUTEND, IUTLEN, IER
      CHARACTER  INSTR(*), OUTSTR(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    IUTLAS, KI, KO
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      CHARACTER  BLANK
      SAVE       BLANK
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  MIN0
      INTRINSIC  MIN0
      INTEGER    MIN0
C
      DATA BLANK/' '/
C                                  CHECK INBEG, INEND, INLEN, IUTBEG,
C                                  AND IUTEND
C
      IF (INBEG.LE.0 .OR. INEND.LT.INBEG .OR. INLEN.LT.INEND .OR.
     &    IUTBEG.LE.0 .OR. IUTEND.LT.IUTBEG) THEN
         IER = -2
         RETURN
      ELSE IF (IUTLEN .LT. IUTEND) THEN
         IER = -1
         RETURN
      END IF
C                                  DETERMINE LAST CHARACTER TO M1VE
      IUTLAS = IUTBEG + MIN0(INEND-INBEG,IUTEND-IUTBEG)
C                                  M1VE CHARACTERS
      KI = INBEG
      DO 10  KO=IUTBEG, IUTLAS
         OUTSTR(KO) = INSTR(KI)
         KI = KI + 1
   10 CONTINUE
C                                   SET IER TO NUMBER OF CHARACTERS THAT
C                                   WHERE NOT MOVED
      IER = KI - INEND - 1
C                                   APPEND BLANKS IF NECESSARY
      DO 20  KO=IUTLAS + 1, IUTEND
         OUTSTR(KO) = BLANK
   20 CONTINUE
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  M1VECH
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    DECEMBER 31, 1984
C
C  PURPOSE:    CHARACTER SUBSTRING ASSIGNMENT.
C
C  USAGE:      CALL M1VECH (STR1, LEN1, STR2, LEN2)
C
C  ARGUMENTS:
C     STR1   - SOURCE SUBSTRING.  (INPUT)
C              THE SOURCE SUBSTRING IS STR1(1:LEN1).
C     LEN1   - LENGTH OF STR1.  (INPUT)
C     STR2   - DESTINATION SUBSTRING.  (OUTPUT)
C              THE DESTINATION SUBSTRING IS STR2(1:LEN2).
C     LEN2   - LENGTH OF STR2.  (INPUT)
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE M1VECH (STR1, LEN1, STR2, LEN2)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    LEN1, LEN2
      CHARACTER  STR1*(*), STR2*(*)
C
      STR2(1:LEN2) = STR1(1:LEN1)
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  I1DX (SINGLE PRECISION VERSION)
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    SEPTEMBER 9, 1985
C
C  PURPOSE:    DETERMINE THE ARRAY SUBSCRIPT INDICATING THE STARTING
C              ELEMENT AT WHICH A KEY CHARACTER SEQUENCE BEGINS.
C              (CASE-INSENSITIVE VERSION)
C
C  USAGE:      I1DX(CHRSTR, I1LEN, KEY, KLEN)
C
C  ARGUMENTS:
C     CHRSTR - CHARACTER ARRAY TO BE SEARCHED.  (INPUT)
C     I1LEN  - LENGTH OF CHRSTR.  (INPUT)
C     KEY    - CHARACTER ARRAY THAT CONTAINS THE KEY SEQUENCE.  (INPUT)
C     KLEN   - LENGTH OF KEY.  (INPUT)
C     I1DX   - INTEGER FUNCTION.  (OUTPUT)
C
C  REMARKS:
C  1. RETURNS ZERO WHEN THERE IS NO MATCH.
C
C  2. RETURNS ZERO IF KLEN IS LONGER THAN ISLEN.
C
C  3. RETURNS ZERO WHEN ANY OF THE CHARACTER ARRAYS HAS A NEGATIVE OR
C     ZERO LENGTH.
C
C  GAMS:       N5C
C
C  CHAPTER:    MATH/LIBRARY UTILITIES
C
C  COPYRIGHT:  1985 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      INTEGER FUNCTION I1DX (CHRSTR, I1LEN, KEY, KLEN)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    I1LEN, KLEN
      CHARACTER  CHRSTR(*), KEY(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, II, J
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   ICASE, I1CSTR
      INTEGER    ICASE, I1CSTR
C
      I1DX = 0
      IF (KLEN.LE.0 .OR. I1LEN.LE.0) GO TO 9000
      IF (KLEN .GT. I1LEN) GO TO 9000
C
      I = 1
      II = I1LEN - KLEN + 1
   10 IF (I .LE. II) THEN
         IF (ICASE(CHRSTR(I)) .EQ. ICASE(KEY(1))) THEN
            IF (KLEN .NE. 1) THEN
               J = KLEN - 1
               IF (I1CSTR(CHRSTR(I+1),J,KEY(2),J) .EQ. 0) THEN
                  I1DX = I
                  GO TO 9000
               END IF
            ELSE
               I1DX = I
               GO TO 9000
            END IF
         END IF
         I = I + 1
         GO TO 10
      END IF
C
 9000 RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  I1KRL
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    AUGUST 9, 1983
C
C  PURPOSE:    DEALLOCATE THE LAST N ALLOCATIONS MADE IN THE WORKSPACE.
C              STACK BY I1KGT
C
C  USAGE:      CALL I1KRL(N)
C
C  ARGUMENTS:
C     N      - NUMBER OF ALLOCATIONS TO BE RELEASED TOP DOWN (INPUT)
C
C  COPYRIGHT:  1983 BY IMSL, INC.  ALL RIGHTS RESERVED
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE I1KRL (N)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    N
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, IN, LALC, LBND, LBOOK, LMAX, LNEED, LNOW, LOUT,
     &           LUSED, NDX, NEXT
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      LOGICAL    FIRST
      SAVE       FIRST
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                                  SPECIFICATIONS FOR COMMON /WORKSP/
      REAL       RWKSP(61913)
      REAL       RDWKSP(5000)
      DOUBLE PRECISION DWKSP(2500)
      COMPLEX    CWKSP(2500)
      COMPLEX    CZWKSP(2500)
      COMPLEX    *16 ZWKSP(1250)
      INTEGER    IWKSP(5000)
      LOGICAL    LWKSP(5000)
      EQUIVALENCE (DWKSP(1), RWKSP(1))
      EQUIVALENCE (CWKSP(1), RWKSP(1)), (ZWKSP(1), RWKSP(1))
      EQUIVALENCE (IWKSP(1), RWKSP(1)), (LWKSP(1), RWKSP(1))
      EQUIVALENCE (RDWKSP(1), RWKSP(1)), (CZWKSP(1), RWKSP(1))
      COMMON     /WORKSP/ RWKSP
C                                  SPECIFICATIONS FOR EQUIVALENCE
      EQUIVALENCE (LOUT, IWKSP(1))
      EQUIVALENCE (LNOW, IWKSP(2))
      EQUIVALENCE (LUSED, IWKSP(3))
      EQUIVALENCE (LBND, IWKSP(4))
      EQUIVALENCE (LMAX, IWKSP(5))
      EQUIVALENCE (LALC, IWKSP(6))
      EQUIVALENCE (LNEED, IWKSP(7))
      EQUIVALENCE (LBOOK, IWKSP(8))
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   E1MES, E1STI, IWKIN
C
      DATA FIRST/.TRUE./
C
      IF (FIRST) THEN
C                                  INITIALIZE WORKSPACE IF NEEDED
         FIRST = .FALSE.
         CALL IWKIN (0)
      END IF
C                                  CALLING I1KRL(0) WILL CONFIRM
C                                  INTEGRITY OF SYSTEM AND RETURN
      IF (N .LT. 0) THEN
         CALL E1MES (5, 10, 'ERROR FROM SUBROUTINE I1KRL:  ATTEMPT'//
     &               ' TO RELEASE A NEGATIVE NUMBER OF WORKSPACE'//
     &               ' ALLOCATIONS. ')
         GO TO 9000
      END IF
C                                  BOOKKEEPING OVERWRITTEN
      IF (LNOW.LT.LBOOK .OR. LNOW.GT.LUSED .OR. LUSED.GT.LMAX .OR.
     &    LNOW.GE.LBND .OR. LOUT.GT.LALC) THEN
         CALL E1MES (5, 11, 'ERROR FROM SUBROUTINE I1KRL:  ONE OR '//
     &               'MORE OF THE FIRST EIGHT BOOKKEEPING LOCATIONS '//
     &               'IN IWKSP HAVE BEEN OVERWRITTEN.  ')
         GO TO 9000
      END IF
C                                  CHECK ALL THE POINTERS IN THE
C                                  PERMANENT STORAGE AREA.  THEY MUST
C                                  BE MONOTONE INCREASING AND LESS THAN
C                                  OR EQUAL TO LMAX, AND THE INDEX OF
C                                  THE LAST POINTER MUST BE LMAX+1.
      NDX = LBND
      IF (NDX .NE. LMAX+1) THEN
         DO 10  I=1, LALC
            NEXT = IWKSP(NDX)
            IF (NEXT .EQ. LMAX+1) GO TO 20
C
            IF (NEXT.LE.NDX .OR. NEXT.GT.LMAX) THEN
               CALL E1MES (5, 12, 'ERROR FROM SUBROUTINE I1KRL:  '//
     &                     'A POINTER IN PERMANENT STORAGE HAS BEEN '//
     &                     ' OVERWRITTEN. ')
               GO TO 9000
            END IF
            NDX = NEXT
   10    CONTINUE
         CALL E1MES (5, 13, 'ERROR FROM SUBROUTINE I1KRL:  A '//
     &               'POINTER IN PERMANENT STORAGE HAS BEEN '//
     &               'OVERWRITTEN. ')
         GO TO 9000
      END IF
   20 IF (N .GT. 0) THEN
         DO 30  IN=1, N
            IF (LNOW .LE. LBOOK) THEN
               CALL E1MES (5, 14, 'ERROR FROM SUBROUTINE I1KRL:  '//
     &                     'ATTEMPT TO RELEASE A NONEXISTANT '//
     &                     'WORKSPACE  ALLOCATION. ')
               GO TO 9000
            ELSE IF (IWKSP(LNOW).LT.LBOOK .OR. IWKSP(LNOW).GE.LNOW-1)
     &              THEN
C                                  CHECK TO MAKE SURE THE BACK POINTERS
C                                  ARE MONOTONE.
               CALL E1STI (1, LNOW)
               CALL E1MES (5, 15, 'ERROR FROM SUBROUTINE I1KRL:  '//
     &                     'THE POINTER AT IWKSP(%(I1)) HAS BEEN '//
     &                     'OVERWRITTEN.  ')
               GO TO 9000
            ELSE
               LOUT = LOUT - 1
               LNOW = IWKSP(LNOW)
            END IF
   30    CONTINUE
      END IF
C
 9000 RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  I1KST
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    AUGUST 9, 1983
C
C  PURPOSE:    RETURN CONTROL INFORMATION ABOUT THE WORKSPACE STACK.
C
C  USAGE:      I1KST(NFACT)
C
C  ARGUMENTS:
C     NFACT  - INTEGER VALUE BETWEEN 1 AND 6 INCLUSIVE RETURNS THE
C                 FOLLOWING INFORMATION: (INPUT)
C                   NFACT = 1 - LOUT: NUMBER OF CURRENT ALLOCATIONS
C                               EXCLUDING PERMANENT STORAGE. AT THE
C                               END OF A RUN, THERE SHOULD BE NO
C                               ACTIVE ALLOCATIONS.
C                   NFACT = 2 - LNOW: CURRENT ACTIVE LENGTH
C                   NFACT = 3 - LTOTAL: TOTAL STORAGE USED THUS FAR
C                   NFACT = 4 - LMAX: MAXIMUM STORAGE ALLOWED
C                   NFACT = 5 - LALC: TOTAL NUMBER OF ALLOCATIONS MADE
C                               BY I1KGT THUS FAR
C                   NFACT = 6 - LNEED: NUMBER OF NUMERIC STORAGE UNITS
C                               BY WHICH THE STACK SIZE MUST BE
C                               INCREASED FOR ALL PAST ALLOCATIONS
C                               TO SUCCEED
C     I1KST  - INTEGER FUNCTION. (OUTPUT) RETURNS A WORKSPACE STACK
C              STATISTIC ACCORDING TO VALUE OF NFACT.
C
C  COPYRIGHT:  1983 BY IMSL, INC.  ALL RIGHTS RESERVED
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      INTEGER FUNCTION I1KST (NFACT)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    NFACT
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    ISTATS(7)
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      LOGICAL    FIRST
      SAVE       FIRST
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                                  SPECIFICATIONS FOR COMMON /WORKSP/
      REAL       RWKSP(61913)
      REAL       RDWKSP(5000)
      DOUBLE PRECISION DWKSP(2500)
      COMPLEX    CWKSP(2500)
      COMPLEX    CZWKSP(2500)
      COMPLEX    *16 ZWKSP(1250)
      INTEGER    IWKSP(5000)
      LOGICAL    LWKSP(5000)
      EQUIVALENCE (DWKSP(1), RWKSP(1))
      EQUIVALENCE (CWKSP(1), RWKSP(1)), (ZWKSP(1), RWKSP(1))
      EQUIVALENCE (IWKSP(1), RWKSP(1)), (LWKSP(1), RWKSP(1))
      EQUIVALENCE (RDWKSP(1), RWKSP(1)), (CZWKSP(1), RWKSP(1))
      COMMON     /WORKSP/ RWKSP
C                                  SPECIFICATIONS FOR EQUIVALENCE
      EQUIVALENCE (ISTATS(1), IWKSP(1))
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   E1MES, IWKIN
C
      DATA FIRST/.TRUE./
C
      IF (FIRST) THEN
C                                  INITIALIZE WORKSPACE IF NEEDED
         FIRST = .FALSE.
         CALL IWKIN (0)
      END IF
C
      IF (NFACT.LE.0 .OR. NFACT.GE.7) THEN
         CALL E1MES (5, 9, 'ERROR FROM SUBROUTINE I1KST:  ARGUMENT'//
     &               ' FOR I1KST MUST BE BETWEEN 1 AND 6 INCLUSIVE.')
      ELSE IF (NFACT .EQ. 1) THEN
C                                  LOUT
         I1KST = ISTATS(1)
      ELSE IF (NFACT .EQ. 2) THEN
C                                  LNOW + PERMANENT
         I1KST = ISTATS(2) + (ISTATS(5)-ISTATS(4)+1)
      ELSE IF (NFACT .EQ. 3) THEN
C                                  LUSED + PERMANENT
         I1KST = ISTATS(3) + (ISTATS(5)-ISTATS(4)+1)
      ELSE IF (NFACT .EQ. 4) THEN
C                                  LMAX
         I1KST = ISTATS(5)
      ELSE IF (NFACT .EQ. 5) THEN
C                                  LALC
         I1KST = ISTATS(6)
      ELSE IF (NFACT .EQ. 6) THEN
C                                  LNEED
         I1KST = ISTATS(7)
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  I1KQU
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    JANUARY 17, 1984
C
C  PURPOSE:    RETURN NUMBER OF ELEMENTS OF DATA TYPE ITYPE THAT
C              REMAIN TO BE ALLOCATED IN ONE REQUEST.
C
C  USAGE:      I1KQU(ITYPE)
C
C  ARGUMENTS:
C     ITYPE  - TYPE OF STORAGE TO BE CHECKED (INPUT)
C                 1 - LOGICAL
C                 2 - INTEGER
C                 3 - REAL
C                 4 - DOUBLE PRECISION
C                 5 - COMPLEX
C                 6 - DOUBLE COMPLEX
C     I1KQU  - INTEGER FUNCTION. (OUTPUT) RETURNS NUMBER OF ELEMENTS
C              OF DATA TYPE ITYPE REMAINING IN THE STACK.
C
C  COPYRIGHT:  1983 BY IMSL, INC.  ALL RIGHTS RESERVED
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      INTEGER FUNCTION I1KQU (ITYPE)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    ITYPE
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    ISIZE(6), LALC, LBND, LBOOK, LMAX, LNEED, LNOW, LOUT,
     &           LUSED
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      LOGICAL    FIRST
      SAVE       FIRST
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                                  SPECIFICATIONS FOR COMMON /WORKSP/
      REAL       RWKSP(61913)
      REAL       RDWKSP(5000)
      DOUBLE PRECISION DWKSP(2500)
      COMPLEX    CWKSP(2500)
      COMPLEX    CZWKSP(2500)
      COMPLEX    *16 ZWKSP(1250)
      INTEGER    IWKSP(5000)
      LOGICAL    LWKSP(5000)
      EQUIVALENCE (DWKSP(1), RWKSP(1))
      EQUIVALENCE (CWKSP(1), RWKSP(1)), (ZWKSP(1), RWKSP(1))
      EQUIVALENCE (IWKSP(1), RWKSP(1)), (LWKSP(1), RWKSP(1))
      EQUIVALENCE (RDWKSP(1), RWKSP(1)), (CZWKSP(1), RWKSP(1))
      COMMON     /WORKSP/ RWKSP
C                                  SPECIFICATIONS FOR EQUIVALENCE
      EQUIVALENCE (LOUT, IWKSP(1))
      EQUIVALENCE (LNOW, IWKSP(2))
      EQUIVALENCE (LUSED, IWKSP(3))
      EQUIVALENCE (LBND, IWKSP(4))
      EQUIVALENCE (LMAX, IWKSP(5))
      EQUIVALENCE (LALC, IWKSP(6))
      EQUIVALENCE (LNEED, IWKSP(7))
      EQUIVALENCE (LBOOK, IWKSP(8))
      EQUIVALENCE (ISIZE(1), IWKSP(11))
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  MAX0
      INTRINSIC  MAX0
      INTEGER    MAX0
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   E1MES, E1POP, E1PSH, IWKIN
C
      DATA FIRST/.TRUE./
C
      CALL E1PSH ('I1KQU ')
C
      IF (FIRST) THEN
C                                  INITIALIZE WORKSPACE IF NEEDED
         FIRST = .FALSE.
         CALL IWKIN (0)
      END IF
C                                  BOOKKEEPING OVERWRITTEN
      IF (LNOW.LT.LBOOK .OR. LNOW.GT.LUSED .OR. LUSED.GT.LMAX .OR.
     &    LNOW.GE.LBND .OR. LOUT.GT.LALC) THEN
         CALL E1MES (5, 7, 'ONE OR MORE OF THE FIRST EIGHT '//
     &               'BOOKKEEPING LOCATIONS IN IWKSP HAVE BEEN '//
     &               'OVERWRITTEN.')
      ELSE IF (ITYPE.LE.0 .OR. ITYPE.GE.7) THEN
C                                  ILLEGAL DATA TYPE REQUESTED
         CALL E1MES (5, 8, 'ILLEGAL DATA TYPE REQUESTED.')
      ELSE
C                                  THIS CALCULATION ALLOWS FOR THE
C                                  TWO POINTER LOCATIONS IN THE STACK
C                                  WHICH ARE ASSIGNED TO EACH ALLOCATION
         I1KQU = MAX0(((LBND-3)*ISIZE(2))/ISIZE(ITYPE)-(LNOW*ISIZE(2)-
     &           1)/ISIZE(ITYPE)-1,0)
      END IF
C
      CALL E1POP ('I1KQU ')
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  IWKIN (SINGLE PRECISION VERSION)
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    JANUARY 17, 1984
C
C  PURPOSE:    INITIALIZE BOOKKEEPING LOCATIONS DESCRIBING THE
C              WORKSPACE STACK.
C
C  USAGE:      CALL IWKIN (NSU)
C
C  ARGUMENT:
C     NSU    - NUMBER OF NUMERIC STORAGE UNITS TO WHICH THE WORKSPACE
C              STACK IS TO BE INITIALIZED
C
C  GAMS:       N4
C
C  CHAPTERS:   MATH/LIBRARY REFERENCE MATERIAL
C              STAT/LIBRARY REFERENCE MATERIAL
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE IWKIN (NSU)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    NSU
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    ISIZE(6), LALC, LBND, LBOOK, LMAX, LNEED, LNOW, LOUT,
     &           LUSED, MELMTS, MTYPE
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      LOGICAL    FIRST
      SAVE       FIRST
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                                  SPECIFICATIONS FOR COMMON /WORKSP/
      REAL       RWKSP(61913)
      REAL       RDWKSP(5000)
      DOUBLE PRECISION DWKSP(2500)
      COMPLEX    CWKSP(2500)
      COMPLEX    CZWKSP(2500)
      COMPLEX    *16 ZWKSP(1250)
      INTEGER    IWKSP(5000)
      LOGICAL    LWKSP(5000)
      EQUIVALENCE (DWKSP(1), RWKSP(1))
      EQUIVALENCE (CWKSP(1), RWKSP(1)), (ZWKSP(1), RWKSP(1))
      EQUIVALENCE (IWKSP(1), RWKSP(1)), (LWKSP(1), RWKSP(1))
      EQUIVALENCE (RDWKSP(1), RWKSP(1)), (CZWKSP(1), RWKSP(1))
      COMMON     /WORKSP/ RWKSP
C                                  SPECIFICATIONS FOR EQUIVALENCE
      EQUIVALENCE (LOUT, IWKSP(1))
      EQUIVALENCE (LNOW, IWKSP(2))
      EQUIVALENCE (LUSED, IWKSP(3))
      EQUIVALENCE (LBND, IWKSP(4))
      EQUIVALENCE (LMAX, IWKSP(5))
      EQUIVALENCE (LALC, IWKSP(6))
      EQUIVALENCE (LNEED, IWKSP(7))
      EQUIVALENCE (LBOOK, IWKSP(8))
      EQUIVALENCE (ISIZE(1), IWKSP(11))
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  MAX0
      INTRINSIC  MAX0
      INTEGER    MAX0
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   E1MES, E1STI
C
      DATA FIRST/.TRUE./
C
      IF (.NOT.FIRST) THEN
         IF (NSU .NE. 0) THEN
            CALL E1STI (1, LMAX)
            CALL E1MES (5, 100, 'ERROR FROM SUBROUTINE IWKIN:  '//
     &                  'WORKSPACE STACK HAS PREVIOUSLY BEEN '//
     &                  'INITIALIZED TO %(I1). CORRECT BY MAKING THE '//
     &                  'CALL TO IWKIN THE FIRST EXECUTABLE '//
     &                  'STATEMENT IN THE MAIN PROGRAM.  ')
C
            STOP
C
         ELSE
            RETURN
         END IF
      END IF
C
      IF (NSU .EQ. 0) THEN
C                                  IF NSU=0 USE DEFAULT SIZE 5000
         MELMTS = 5000
      ELSE
         MELMTS = NSU
      END IF
C                                  NUMBER OF ITEMS .LT. 0
      IF (MELMTS .LE. 0) THEN
         CALL E1STI (1, MELMTS)
         CALL E1MES (5, 1, 'ERROR FROM SUBROUTINE IWKIN:  NUMBER '//
     &               'OF NUMERIC STORAGE UNITS IS NOT POSITIVE. NSU '//
     &               '= %(I1) ')
      ELSE
C
         FIRST = .FALSE.
C                                  HERE TO INITIALIZE
C
C                                  SET DATA SIZES APPROPRIATE FOR A
C                                  STANDARD CONFORMING FORTRAN SYSTEM
C                                  USING THE FORTRAN
C                                  *NUMERIC STORAGE UNIT* AS THE
C                                  MEASURE OF SIZE.
C
C                                  TYPE IS REAL
         MTYPE = 3
C                                  LOGICAL
         ISIZE(1) = 1
C                                  INTEGER
         ISIZE(2) = 1
C                                  REAL
         ISIZE(3) = 1
C                                  DOUBLE PRECISION
         ISIZE(4) = 2
C                                  COMPLEX
         ISIZE(5) = 2
C                                  DOUBLE COMPLEX
         ISIZE(6) = 4
C                                  NUMBER OF WORDS USED FOR BOOKKEEPING
         LBOOK = 16
C                                  CURRENT ACTIVE LENGTH OF THE STACK
         LNOW = LBOOK
C                                  MAXIMUM VALUE OF LNOW ACHIEVED THUS
C                                  FAR
         LUSED = LBOOK
C                                  MAXIMUM LENGTH OF THE STORAGE ARRAY
         LMAX = MAX0(MELMTS,((LBOOK+2)*ISIZE(2)+ISIZE(3)-1)/ISIZE(3))
C                                  LOWER BOUND OF THE PERMANENT STORAGE
C                                  WHICH IS ONE WORD MORE THAN THE
C                                  MAXIMUM ALLOWED LENGTH OF THE STACK
         LBND = LMAX + 1
C                                  NUMBER OF CURRENT ALLOCATIONS
         LOUT = 0
C                                  TOTAL NUMBER OF ALLOCATIONS MADE
         LALC = 0
C                                  NUMBER OF WORDS BY WHICH THE ARRAY
C                                  SIZE MUST BE INCREASED FOR ALL PAST
C                                  ALLOCATIONS TO SUCCEED
         LNEED = 0
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  I1X (SINGLE PRECISION VERSION)
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    AUGUST 30, 1985
C
C  PURPOSE:    DETERMINE THE ARRAY SUBSCRIPT INDICATING THE STARTING
C              ELEMENT AT WHICH A KEY CHARACTER SEQUENCE BEGINS.
C              (CASE-SENSITIVE VERSION)
C
C  USAGE:      I1X(CHRSTR, I1LEN, KEY, KLEN)
C
C  ARGUMENTS:
C     CHRSTR - CHARACTER ARRAY TO BE SEARCHED.  (INPUT)
C     I1LEN  - LENGTH OF CHRSTR.  (INPUT)
C     KEY    - CHARACTER ARRAY THAT CONTAINS THE KEY SEQUENCE.  (INPUT)
C     KLEN   - LENGTH OF KEY.  (INPUT)
C     I1X    - INTEGER FUNCTION.  (OUTPUT)
C
C  REMARKS:
C  1. RETURNS ZERO WHEN THERE IS NO MATCH.
C
C  2. RETURNS ZERO IF KLEN IS LONGER THAN ISLEN.
C
C  3. RETURNS ZERO WHEN ANY OF THE CHARACTER ARRAYS HAS A NEGATIVE OR
C     ZERO LENGTH.
C
C  GAMS:       N5C
C
C  CHAPTER:    MATH/LIBRARY UTILITIES
C
C  COPYRIGHT:  1985 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      INTEGER FUNCTION I1X (CHRSTR, I1LEN, KEY, KLEN)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    I1LEN, KLEN
      CHARACTER  CHRSTR(*), KEY(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, II, J
C
      I1X = 0
      IF (KLEN.LE.0 .OR. I1LEN.LE.0) GO TO 9000
      IF (KLEN .GT. I1LEN) GO TO 9000
C
      I = 1
      II = I1LEN - KLEN + 1
   10 IF (I .LE. II) THEN
         IF (CHRSTR(I) .EQ. KEY(1)) THEN
            DO 20  J=2, KLEN
               IF (CHRSTR(I+J-1) .NE. KEY(J)) GO TO 30
   20       CONTINUE
            I1X = I
            GO TO 9000
   30       CONTINUE
         END IF
         I = I + 1
         GO TO 10
      END IF
C
 9000 RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  C1TIC
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    MARCH 9, 1984
C
C  PURPOSE:    CONVERT AN INTEGER TO ITS CORRESPONDING CHARACTER FORM.
C              (RIGHT JUSTIFIED)
C
C  USAGE:      CALL C1TIC(NUM, CHRSTR, SLEN, IER)
C
C  ARGUMENTS:
C     NUM    - INTEGER NUMBER.  (INPUT)
C     CHRSTR - CHARACTER ARRAY THAT RECEIVES THE RESULT.  (OUTPUT)
C     SLEN   - LENGTH OF THE CHARACTER ARRAY.  (INPUT)
C     IER    - COMPLETION CODE.  (OUTPUT) WHERE
C                 IER < 0  INDICATES THAT SLEN <= 0,
C                 IER = 0  INDICATES NORMAL COMPLETION,
C                 IER > 0  INDICATES THAT THE CHARACTER ARRAY IS TOO
C                       SMALL TO HOLD THE COMPLETE NUMBER.  IER
C                       INDICATES HOW MANY SIGNIFICANT DIGITS ARE
C                       BEING TRUNCATED.
C
C  REMARKS:
C  1. THE CHARACTER ARRAY IS FILLED IN A RIGHT JUSTIFIED MANNER.
C  2. LEADING ZEROS ARE REPLACED BY BLANKS.
C  3. SIGN IS INSERTED ONLY FOR NEGATIVE NUMBER.
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE C1TIC (NUM, CHRSTR, SLEN, IER)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    NUM, SLEN, IER
      CHARACTER  CHRSTR(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, J, K, L
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      CHARACTER  BLANK(1), DIGIT(10), MINUS(1)
      SAVE       BLANK, DIGIT, MINUS
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  IABS
      INTRINSIC  IABS
      INTEGER    IABS
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   M1VE
C
      DATA DIGIT/'0', '1', '2', '3', '4', '5', '6', '7', '8',
     &     '9'/
      DATA BLANK/' '/, MINUS/'-'/
C                                  CHECK SLEN
      IF (SLEN .LE. 0) THEN
         IER = -1
         RETURN
      END IF
C                                  THE NUMBER IS ZERO
      IF (NUM .EQ. 0) THEN
         CALL M1VE (BLANK, 1, 1, 1, CHRSTR, 1, SLEN-1, SLEN, I)
         CHRSTR(SLEN) = DIGIT(1)
         IER = 0
         RETURN
      END IF
C                                  CONVERT NUMBER DIGIT BY DIGIT TO
C                                  CHARACTER FORM
      J = SLEN
      K = IABS(NUM)
   10 IF (K.GT.0 .AND. J.GE.1) THEN
         L = K
         K = K/10
         L = L - K*10
         CHRSTR(J) = DIGIT(L+1)
         J = J - 1
         GO TO 10
      END IF
C
   20 IF (K .EQ. 0) THEN
         IF (NUM .LT. 0) THEN
            CALL M1VE (MINUS, 1, 1, 1, CHRSTR, J, J, SLEN, I)
            IF (I .NE. 0) THEN
               IER = 1
               RETURN
            END IF
            J = J - 1
         END IF
         IER = 0
         CALL M1VE (BLANK, 1, 1, 1, CHRSTR, 1, J, SLEN, I)
         RETURN
      END IF
C                                  DETERMINE THE NUMBER OF SIGNIFICANT
C                                  DIGITS BEING TRUNCATED
      I = 0
   30 IF (K .GT. 0) THEN
         K = K/10
         I = I + 1
         GO TO 30
      END IF
C
      IF (NUM .LT. 0) I = I + 1
      IER = I
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  N1RGB
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    MARCH 2, 1984
C
C  PURPOSE:    RETURN A POSITIVE NUMBER AS A FLAG TO INDICATED THAT A
C              STOP SHOULD OCCUR DUE TO ONE OR MORE GLOBAL ERRORS.
C
C  USAGE:      N1RGB(IDUMMY)
C
C  ARGUMENTS:
C     IDUMMY - INTEGER SCALAR DUMMY ARGUMENT.
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      INTEGER FUNCTION N1RGB (IDUMMY)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    IDUMMY
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                              SPECIFICATIONS FOR COMMON /ERCOM1/
      INTEGER    CALLVL, MAXLEV, MSGLEN, ERTYPE(51), ERCODE(51),
     &           PRINTB(7), STOPTB(7), PLEN, IFERR6, IFERR7,
     &           IALLOC(51), HDRFMT(7), TRACON(7)
      COMMON     /ERCOM1/ CALLVL, MAXLEV, MSGLEN, ERTYPE, ERCODE,
     &           PRINTB, STOPTB, PLEN, IFERR6, IFERR7, IALLOC, HDRFMT,
     &           TRACON
      SAVE       /ERCOM1/
C                              SPECIFICATIONS FOR COMMON /ERCOM2/
      CHARACTER  MSGSAV(255), PLIST(300), RNAME(51)*6
      COMMON     /ERCOM2/ MSGSAV, PLIST, RNAME
      SAVE       /ERCOM2/
C                              SPECIFICATIONS FOR COMMON /ERCOM3/
      DOUBLE PRECISION ERCKSM
      COMMON     /ERCOM3/ ERCKSM
      SAVE       /ERCOM3/
C                              SPECIFICATIONS FOR COMMON /ERCOM4/
      LOGICAL    ISUSER(51)
      COMMON     /ERCOM4/ ISUSER
      SAVE       /ERCOM4/
C                                  INITIALIZE FUNCTION
      N1RGB = 0
C                                  CHECK FOR GLOBAL ERROR TYPE 6
      IF (IFERR6 .GT. 0) THEN
         N1RGB = STOPTB(6)
         IFERR6 = 0
      END IF
C                                  CHECK FOR GLOBAL ERROR TYPE 7
      IF (IFERR7 .GT. 0) THEN
         N1RGB = STOPTB(7)
         IFERR7 = 0
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  SDOT (SINGLE PRECISION VERSION)
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    AUGUST 9, 1986
C
C  PURPOSE:    COMPUTE THE SINGLE-PRECISION DOT PRODUCT X*Y.
C
C  USAGE:      SDOT(N, SX, INCX, SY, INCY)
C
C  ARGUMENTS:
C     N      - LENGTH OF VECTORS X AND Y.  (INPUT)
C     SX     - REAL VECTOR OF LENGTH MAX(N*IABS(INCX),1).  (INPUT)
C     INCX   - DISPLACEMENT BETWEEN ELEMENTS OF SX.  (INPUT)
C              X(I) IS DEFINED TO BE.. SX(1+(I-1)*INCX) IF INCX .GE. 0
C              OR SX(1+(I-N)*INCX) IF INCX .LT. 0.
C     SY     - REAL VECTOR OF LENGTH MAX(N*IABS(INCY),1).  (INPUT)
C     INCY   - DISPLACEMENT BETWEEN ELEMENTS OF SY.  (INPUT)
C              Y(I) IS DEFINED TO BE.. SY(1+(I-1)*INCY) IF INCY .GE. 0
C              OR SY(1+(I-N)*INCY) IF INCY .LT. 0.
C     SDOT   - SUM FROM I=1 TO N OF X(I)*Y(I).  (OUTPUT)
C              X(I) AND Y(I) REFER TO SPECIFIC ELEMENTS OF SX AND SY,
C              RESPECTIVELY.  SEE INCX AND INCY ARGUMENT DESCRIPTIONS.
C
C  GAMS:       D1A4
C
C  CHAPTERS:   MATH/LIBRARY BASIC MATRIX/VECTOR OPERATIONS
C              STAT/LIBRARY MATHEMATICAL SUPPORT
C
C  COPYRIGHT:  1986 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      REAL FUNCTION SDOT (N, SX, INCX, SY, INCY)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    N, INCX, INCY
      REAL       SX(*), SY(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, IX, IY, M, MP1
C                                  SPECIFICATIONS FOR SPECIAL CASES
C     INTRINSIC  MOD
      INTRINSIC  MOD
      INTEGER    MOD
C
      SDOT = 0.0E0
      IF (N .GT. 0) THEN
         IF (INCX.NE.1 .OR. INCY.NE.1) THEN
C                                  CODE FOR UNEQUAL INCREMENTS
            IX = 1
            IY = 1
            IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
            IF (INCY .LT. 0) IY = (-N+1)*INCY + 1
            DO 10  I=1, N
               SDOT = SDOT + SX(IX)*SY(IY)
               IX = IX + INCX
               IY = IY + INCY
   10       CONTINUE
         ELSE
C                                  CODE FOR BOTH INCREMENTS EQUAL TO 1
            M = MOD(N,5)
C                                  CLEAN-UP LOOP SO REMAINING VECTOR
            DO 30  I=1, M
               SDOT = SDOT + SX(I)*SY(I)
   30       CONTINUE
            MP1 = M + 1
            DO 40  I=MP1, N, 5
               SDOT = SDOT + SX(I)*SY(I) + SX(I+1)*SY(I+1) +
     &                SX(I+2)*SY(I+2) + SX(I+3)*SY(I+3) +
     &                SX(I+4)*SY(I+4)
   40       CONTINUE
         END IF
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  E1POS
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    MARCH 2, 1984
C
C  PURPOSE:    SET OR RETRIEVE PRINT AND STOP ATTRIBUTES.
C
C  USAGE:      CALL E1POS(IERTYP,IPATT,ISATT)
C
C  ARGUMENTS:
C     IERTYP - INTEGER SPECIFYING THE ERROR TYPE FOR WHICH PRINT AND
C              STOP ATTRIBUTES ARE TO BE SET OR RETRIEVED.  (INPUT)  IF
C              IERTYP IS 0 THEN THE SETTINGS APPLY TO ALL ERROR TYPES.
C              IF IERTYP IS BETWEEN 1 AND 7, THEN THE SETTINGS ONLY
C              APPLY TO THAT SPECIFIED ERROR TYPE.  IF IERTYP IS
C              NEGATIVE THEN THE CURRENT PRINT AND STOP ATTRIBUTES WILL
C              BE RETURNED IN IPATT AND ISATT.
C     IPATT  - IF IERTYP IS POSITIVE, IPATT IS AN INTEGER SPECIFYING THE
C              DESIRED PRINT ATTRIBUTE AS FOLLOWS: -1 MEANS NO CHANGE,
C              0 MEANS NO, 1 MEANS YES, AND 2 MEANS ASSIGN THE DEFAULT
C              SETTING.  (INPUT)  IF IERTYP IS NEGATIVE, IPATT IS
C              RETURNED AS 1 IF PRINT IS YES OR 0 IF PRINT IS NO FOR
C              ERROR TYPE IABS(IERTYP).  (OUTPUT)
C     ISATT  - IF IERTYP IS POSITIVE, ISATT IS AN INTEGER SPECIFYING THE
C              DESIRED STOP ATTRIBUTE AS FOLLOWS: -1 MEANS NO CHANGE,
C              0 MEANS NO, 1 MEANS YES, AND 2 MEANS ASSIGN THE DEFAULT
C              SETTING.  (INPUT)  IF IERTYP IS NEGATIVE, ISATT IS
C              RETURNED AS 1 IF PRINT IS YES OR 0 IF PRINT IS NO FOR
C              ERROR TYPE IABS(IERTYP).  (OUTPUT)
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE E1POS (IERTYP, IPATT, ISATT)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    IERTYP, IPATT, ISATT
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, IER
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      INTEGER    DEFLTP(7), DEFLTS(7), IFINIT
      SAVE       DEFLTP, DEFLTS, IFINIT
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                              SPECIFICATIONS FOR COMMON /ERCOM1/
      INTEGER    CALLVL, MAXLEV, MSGLEN, ERTYPE(51), ERCODE(51),
     &           PRINTB(7), STOPTB(7), PLEN, IFERR6, IFERR7,
     &           IALLOC(51), HDRFMT(7), TRACON(7)
      COMMON     /ERCOM1/ CALLVL, MAXLEV, MSGLEN, ERTYPE, ERCODE,
     &           PRINTB, STOPTB, PLEN, IFERR6, IFERR7, IALLOC, HDRFMT,
     &           TRACON
      SAVE       /ERCOM1/
C                              SPECIFICATIONS FOR COMMON /ERCOM2/
      CHARACTER  MSGSAV(255), PLIST(300), RNAME(51)*6
      COMMON     /ERCOM2/ MSGSAV, PLIST, RNAME
      SAVE       /ERCOM2/
C                              SPECIFICATIONS FOR COMMON /ERCOM3/
      DOUBLE PRECISION ERCKSM
      COMMON     /ERCOM3/ ERCKSM
      SAVE       /ERCOM3/
C                              SPECIFICATIONS FOR COMMON /ERCOM4/
      LOGICAL    ISUSER(51)
      COMMON     /ERCOM4/ ISUSER
      SAVE       /ERCOM4/
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  IABS
      INTRINSIC  IABS
      INTEGER    IABS
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   E1INIT, E1MES, E1STI
C
      DATA IFINIT/0/
      DATA DEFLTP/0, 0, 1, 1, 1, 1, 1/, DEFLTS/0, 0, 0, 1, 1, 0, 1/
C                                  INITIALIZE ERROR TABLE IF NECESSARY
      IF (IFINIT .EQ. 0) THEN
         CALL E1INIT
         IFINIT = 1
      END IF
      IER = 0
      IF (IERTYP .GE. 0) THEN
         IF (IPATT.LT.-1 .OR. IPATT.GT.2) THEN
            CALL E1STI (1, IPATT)
            CALL E1MES (5, 1, 'INVALID VALUE SPECIFIED FOR PRINT '//
     &                  'TABLE ATTRIBUTE.  IPATT MUST BE -1, 0, 1, '//
     &                  'OR 2.  IPATT = %(I1)')
            IER = 1
         END IF
         IF (ISATT.LT.-1 .OR. ISATT.GT.2) THEN
            CALL E1STI (1, ISATT)
            CALL E1MES (5, 1, 'INVALID VALUE SPECIFIED FOR STOP '//
     &                  'TABLE ATTRIBUTE.  ISATT MUST BE -1, 0, 1, '//
     &                  'OR 2.  ISATT = %(I1)')
            IER = 1
         END IF
      END IF
      IF (IER .EQ. 0) THEN
         IF (IERTYP .EQ. 0) THEN
            IF (IPATT.EQ.0 .OR. IPATT.EQ.1) THEN
               DO 10  I=1, 7
   10          PRINTB(I) = IPATT
            ELSE IF (IPATT .EQ. 2) THEN
C                                  ASSIGN DEFAULT SETTINGS
               DO 20  I=1, 7
   20          PRINTB(I) = DEFLTP(I)
            END IF
            IF (ISATT.EQ.0 .OR. ISATT.EQ.1) THEN
               DO 30  I=1, 7
   30          STOPTB(I) = ISATT
            ELSE IF (ISATT .EQ. 2) THEN
C                                  ASSIGN DEFAULT SETTINGS
               DO 40  I=1, 7
   40          STOPTB(I) = DEFLTS(I)
            END IF
         ELSE IF (IERTYP.GE.1 .AND. IERTYP.LE.7) THEN
            IF (IPATT.EQ.0 .OR. IPATT.EQ.1) THEN
               PRINTB(IERTYP) = IPATT
            ELSE IF (IPATT .EQ. 2) THEN
C                                  ASSIGN DEFAULT SETTING
               PRINTB(IERTYP) = DEFLTP(IERTYP)
            END IF
            IF (ISATT.EQ.0 .OR. ISATT.EQ.1) THEN
               STOPTB(IERTYP) = ISATT
            ELSE IF (ISATT .EQ. 2) THEN
C                                  ASSIGN DEFAULT SETTING
               STOPTB(IERTYP) = DEFLTS(IERTYP)
            END IF
         ELSE IF (IERTYP.LE.-1 .AND. IERTYP.GE.-7) THEN
            I = IABS(IERTYP)
            IPATT = PRINTB(I)
            ISATT = STOPTB(I)
         END IF
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  E1STL
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    NOVEMBER 8, 1985
C
C  PURPOSE:    TO STORE A STRING FOR SUBSEQUENT USE WITHIN AN ERROR
C              MESSAGE.
C
C  USAGE:      CALL E1STL(IL,STRING)
C
C  ARGUMENTS:
C     IL     - INTEGER SPECIFYING THE SUBSTITUTION INDEX.  IL MUST BE
C              BETWEEN 1 AND 9.  (INPUT)
C     STRING - A CHARACTER STRING.  (INPUT)
C
C  COPYRIGHT:  1985 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE E1STL (IL, STRING)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    IL
      CHARACTER  STRING*(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, LEN2
      CHARACTER  STRGUP(255)
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      INTEGER    IFINIT
      SAVE       IFINIT
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                              SPECIFICATIONS FOR COMMON /ERCOM1/
      INTEGER    CALLVL, MAXLEV, MSGLEN, ERTYPE(51), ERCODE(51),
     &           PRINTB(7), STOPTB(7), PLEN, IFERR6, IFERR7,
     &           IALLOC(51), HDRFMT(7), TRACON(7)
      COMMON     /ERCOM1/ CALLVL, MAXLEV, MSGLEN, ERTYPE, ERCODE,
     &           PRINTB, STOPTB, PLEN, IFERR6, IFERR7, IALLOC, HDRFMT,
     &           TRACON
      SAVE       /ERCOM1/
C                              SPECIFICATIONS FOR COMMON /ERCOM2/
      CHARACTER  MSGSAV(255), PLIST(300), RNAME(51)*6
      COMMON     /ERCOM2/ MSGSAV, PLIST, RNAME
      SAVE       /ERCOM2/
C                              SPECIFICATIONS FOR COMMON /ERCOM3/
      DOUBLE PRECISION ERCKSM
      COMMON     /ERCOM3/ ERCKSM
      SAVE       /ERCOM3/
C                              SPECIFICATIONS FOR COMMON /ERCOM4/
      LOGICAL    ISUSER(51)
      COMMON     /ERCOM4/ ISUSER
      SAVE       /ERCOM4/
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  IABS,LEN,MIN0
      INTRINSIC  IABS, LEN, MIN0
      INTEGER    IABS, LEN, MIN0
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   E1INIT, E1INPL
C
      DATA IFINIT/0/
C                                  INITIALIZE IF NECESSARY
      IF (IFINIT .EQ. 0) THEN
         CALL E1INIT
         IFINIT = 1
      END IF
      LEN2 = LEN(STRING)
      LEN2 = MIN0(LEN2,255)
      DO 10  I=1, LEN2
         STRGUP(I) = STRING(I:I)
   10 CONTINUE
      IF (IABS(IL).GE.1 .AND. IABS(IL).LE.9) THEN
         CALL E1INPL ('L', IL, LEN2, STRGUP)
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  SAXPY (SINGLE PRECISION VERSION)
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    AUGUST 9, 1986
C
C  PURPOSE:    COMPUTE THE SCALAR TIMES A VECTOR PLUS A VECTOR,
C              Y = AX + Y, ALL SINGLE PRECISION.
C
C  USAGE:      CALL SAXPY (N, SA, SX, INCX, SY, INCY)
C
C  ARGUMENTS:
C     N      - LENGTH OF VECTORS X AND Y.  (INPUT)
C     SA     - REAL SCALAR.  (INPUT)
C     SX     - REAL VECTOR OF LENGTH MAX(N*IABS(INCX),1).  (INPUT)
C     INCX   - DISPLACEMENT BETWEEN ELEMENTS OF SX.  (INPUT)
C              X(I) IS DEFINED TO BE
C                 SX(1+(I-1)*INCX) IF INCX.GE.0  OR
C                 SX(1+(I-N)*INCX) IF INCX.LT.0.
C     SY     - REAL VECTOR OF LENGTH MAX(N*IABS(INCY),1).
C              (INPUT/OUTPUT)
C              SAXPY REPLACES Y(I) WITH SA*X(I) + Y(I) FOR I=1,...,N.
C              X(I) AND Y(I) REFER TO SPECIFIC ELEMENTS OF SX AND SY.
C     INCY   - DISPLACEMENT BETWEEN ELEMENTS OF SY.  (INPUT)
C              Y(I) IS DEFINED TO BE
C                 SY(1+(I-1)*INCY) IF INCY.GE.0  OR
C                 SY(1+(I-N)*INCY) IF INCY.LT.0.
C
C  GAMS:       D1A7
C
C  CHAPTERS:   MATH/LIBRARY BASIC MATRIX/VECTOR OPERATIONS
C              STAT/LIBRARY MATHEMATICAL SUPPORT
C
C  COPYRIGHT:  1986 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE SAXPY (N, SA, SX, INCX, SY, INCY)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    N, INCX, INCY
      REAL       SA, SX(*), SY(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, IX, IY, M, MP1
C                                  SPECIFICATIONS FOR SPECIAL CASES
C     INTRINSIC  MOD
      INTRINSIC  MOD
      INTEGER    MOD
C
      IF (N .GT. 0) THEN
         IF (SA .NE. 0.0) THEN
            IF (INCX.NE.1 .OR. INCY.NE.1) THEN
C                                  CODE FOR UNEQUAL INCREMENTS OR EQUAL
C                                  INCREMENTS NOT EQUAL TO 1
               IX = 1
               IY = 1
               IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
               IF (INCY .LT. 0) IY = (-N+1)*INCY + 1
               DO 10  I=1, N
                  SY(IY) = SY(IY) + SA*SX(IX)
                  IX = IX + INCX
                  IY = IY + INCY
   10          CONTINUE
            ELSE
C                                  CODE FOR BOTH INCREMENTS EQUAL TO 1
               M = MOD(N,4)
C                                  CLEAN-UP LOOP
               DO 30  I=1, M
                  SY(I) = SY(I) + SA*SX(I)
   30          CONTINUE
               MP1 = M + 1
               DO 40  I=MP1, N, 4
                  SY(I) = SY(I) + SA*SX(I)
                  SY(I+1) = SY(I+1) + SA*SX(I+1)
                  SY(I+2) = SY(I+2) + SA*SX(I+2)
                  SY(I+3) = SY(I+3) + SA*SX(I+3)
   40          CONTINUE
            END IF
         END IF
      END IF
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  I1ERIF
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    MARCH 13, 1984
C
C  PURPOSE:    RETURN THE POSITION OF THE FIRST ELEMENT OF A GIVEN
C              CHARACTER ARRAY WHICH IS NOT AN ELEMENT OF ANOTHER
C              CHARACTER ARRAY.
C
C  USAGE:      I1ERIF(STR1, LEN1, STR2, LEN2)
C
C  ARGUMENTS:
C     STR1   - CHARACTER ARRAY TO BE SEARCHED.  (INPUT)
C     LEN1   - LENGTH OF STR1.  (INPUT)
C     STR2   - CHARACTER ARRAY TO BE SEARCHED FOR.  (INPUT)
C     LEN2   - LENGTH OF STR2.  (INPUT)
C     I1ERIF - INTEGER FUNCTION.  (OUTPUT)
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      INTEGER FUNCTION I1ERIF (STR1, LEN1, STR2, LEN2)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    LEN1, LEN2
      CHARACTER  STR1(*), STR2(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   I1X
      INTEGER    I1X
C                              FIRST EXECUTABLE STATEMENT
      IF (LEN1.LE.0 .OR. LEN2.LE.0) THEN
         I1ERIF = 1
      ELSE
         DO 10  I=1, LEN1
            IF (I1X(STR2,LEN2,STR1(I),1) .EQ. 0) THEN
               I1ERIF = I
               RETURN
            END IF
   10    CONTINUE
         I1ERIF = 0
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  ICASE (SINGLE PRECISION VERSION)
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    SEPTEMBER 9, 1985
C
C  PURPOSE:    CONVERT FROM CHARACTER TO THE INTEGER ASCII VALUE WITHOUT
C              REGARD TO CASE.
C
C  USAGE:      ICASE(CH)
C
C  ARGUMENTS:
C     CH     - CHARACTER TO BE CONVERTED.  (INPUT)
C     ICASE  - INTEGER ASCII VALUE FOR CH WITHOUT REGARD TO THE CASE
C              OF CH.  (OUTPUT)
C              ICASE RETURNS THE SAME VALUE AS IMSL ROUTINE IACHAR FOR
C              ALL BUT LOWERCASE LETTERS.  FOR THESE, IT RETURNS THE
C              IACHAR VALUE FOR THE CORRESPONDING UPPERCASE LETTER.
C
C  GAMS:       N3
C
C  CHAPTER:    MATH/LIBRARY UTILITIES
C              STAT/LIBRARY UTILITIES
C
C  COPYRIGHT:  1986 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      INTEGER FUNCTION ICASE (CH)
C                                  SPECIFICATIONS FOR ARGUMENTS
      CHARACTER  CH
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   IACHAR
      INTEGER    IACHAR
C
      ICASE = IACHAR(CH)
      IF (ICASE.GE.97 .AND. ICASE.LE.122) ICASE = ICASE - 32
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  S1ANUM
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    MARCH 28, 1984
C
C  PURPOSE:    SCAN A TOKEN AND IDENTIFY IT AS FOLLOWS: INTEGER, REAL
C              NUMBER (SINGLE/DOUBLE), FORTRAN RELATIONAL OPERATOR,
C              FORTRAN LOGICAL OPERATOR, OR FORTRAN LOGICAL CONSTANT.
C
C  USAGE:      CALL S1ANUM(INSTR, SLEN, CODE, OLEN)
C
C  ARGUMENTS:
C     INSTR  - CHARACTER STRING TO BE SCANNED.  (INPUT)
C     SLEN   - LENGTH OF INSTR.  (INPUT)
C     CODE   - TOKEN CODE.  (OUTPUT)  WHERE
C                 CODE =  0  INDICATES AN UNKNOWN TOKEN,
C                 CODE =  1  INDICATES AN INTEGER NUMBER,
C                 CODE =  2  INDICATES A (SINGLE PRECISION) REAL NUMBER,
C                 CODE =  3  INDICATES A (DOUBLE PRECISION) REAL NUMBER,
C                 CODE =  4  INDICATES A LOGICAL CONSTANT (.TRUE. OR
C                               .FALSE.),
C                 CODE =  5  INDICATES THE RELATIONAL OPERATOR .EQ.,
C                 CODE =  6  INDICATES THE RELATIONAL OPERATOR .NE.,
C                 CODE =  7  INDICATES THE RELATIONAL OPERATOR .LT.,
C                 CODE =  8  INDICATES THE RELATIONAL OPERATOR .LE.,
C                 CODE =  9  INDICATES THE RELATIONAL OPERATOR .GT.,
C                 CODE = 10  INDICATES THE RELATIONAL OPERATOR .GE.,
C                 CODE = 11  INDICATES THE LOGICAL OPERATOR .AND.,
C                 CODE = 12  INDICATES THE LOGICAL OPERATOR .OR.,
C                 CODE = 13  INDICATES THE LOGICAL OPERATOR .EQV.,
C                 CODE = 14  INDICATES THE LOGICAL OPERATOR .NEQV.,
C                 CODE = 15  INDICATES THE LOGICAL OPERATOR .NOT..
C     OLEN   - LENGTH OF THE TOKEN AS COUNTED FROM THE FIRST CHARACTER
C              IN INSTR.  (OUTPUT)  OLEN RETURNS A ZERO FOR AN UNKNOWN
C              TOKEN (CODE = 0).
C
C  REMARKS:
C  1. BLANKS ARE CONSIDERED SIGNIFICANT.
C  2. LOWER AND UPPER CASE LETTERS ARE NOT SIGNIFICANT.
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE S1ANUM (INSTR, SLEN, CODE, OLEN)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    SLEN, CODE, OLEN
      CHARACTER  INSTR(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, IBEG, IIBEG, J
      LOGICAL    FLAG
      CHARACTER  CHRSTR(6)
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      INTEGER    TABPTR(16), TDCNST, TICNST, TOKEN(13), TRCNST, TZERR
      CHARACTER  DIGIT(10), LETTER(52), MINUS, PERIOD, PLUS, TABLE(38)
      SAVE       DIGIT, LETTER, MINUS, PERIOD, PLUS, TABLE, TABPTR,
     &           TDCNST, TICNST, TOKEN, TRCNST, TZERR
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   I1X, I1CSTR
      INTEGER    I1X, I1CSTR
C
      DATA TOKEN/5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 4, 4/
      DATA TABLE/'D', 'E', 'E', 'Q', 'N', 'E', 'L', 'T', 'L',
     &     'E', 'G', 'T', 'G', 'E', 'A', 'N', 'D', 'O', 'R',
     &     'E', 'Q', 'V', 'N', 'E', 'Q', 'V', 'N', 'O', 'T',
     &     'T', 'R', 'U', 'E', 'F', 'A', 'L', 'S', 'E'/
      DATA TABPTR/1, 2, 3, 5, 7, 9, 11, 13, 15, 18, 20, 23, 27, 30,
     &     34, 39/
      DATA DIGIT/'0', '1', '2', '3', '4', '5', '6', '7', '8',
     &     '9'/
      DATA LETTER/'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I',
     &     'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S',
     &     'T', 'U', 'V', 'W', 'X', 'Y', 'Z', 'A', 'B', 'C',
     &     'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M',
     &     'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W',
     &     'X', 'Y', 'Z'/
      DATA PERIOD/'.'/, PLUS/'+'/, MINUS/'-'/
      DATA TZERR/0/, TICNST/1/
      DATA TRCNST/2/, TDCNST/3/
C
      IF (SLEN .LE. 0) THEN
         CODE = 0
         OLEN = 0
         RETURN
      END IF
C                                  STATE 0 - ASSUME ERROR TOKEN
      IBEG = 1
      CODE = TZERR
C                                  CHECK SIGN
      IF (INSTR(IBEG).EQ.MINUS .OR. INSTR(IBEG).EQ.PLUS) THEN
         FLAG = .TRUE.
         IIBEG = IBEG
         IBEG = IBEG + 1
      ELSE
         FLAG = .FALSE.
      END IF
C                                  STATE 1 - ASSUME INTEGER CONSTANT
      IF (I1X(DIGIT,10,INSTR(IBEG),1) .NE. 0) THEN
         CODE = TICNST
         IIBEG = IBEG
         IBEG = IBEG + 1
C
   10    IF (IBEG .LE. SLEN) THEN
C
            IF (I1X(DIGIT,10,INSTR(IBEG),1) .NE. 0) THEN
               IIBEG = IBEG
               IBEG = IBEG + 1
               GO TO 10
C
            END IF
C
         ELSE
            GO TO 80
C
         END IF
C
         IF (INSTR(IBEG) .NE. PERIOD) GO TO 80
      END IF
C                                  STATE 2 - ASSUME REAL CONSTANT
      IF (CODE .EQ. TICNST) THEN
         CODE = TRCNST
         IIBEG = IBEG
         IBEG = IBEG + 1
         IF (IBEG .GT. SLEN) GO TO 80
      ELSE IF (INSTR(IBEG).EQ.PERIOD .AND. SLEN.GE.2) THEN
         IF (I1X(DIGIT,10,INSTR(IBEG+1),1) .NE. 0) THEN
            CODE = TRCNST
            IIBEG = IBEG + 1
            IBEG = IBEG + 2
            IF (IBEG .GT. SLEN) GO TO 80
         END IF
      END IF
C
      IF (I1X(DIGIT,10,INSTR(IBEG),1) .NE. 0) THEN
         CODE = TRCNST
         IIBEG = IBEG
         IBEG = IBEG + 1
C
   20    IF (IBEG .LE. SLEN) THEN
C
            IF (I1X(DIGIT,10,INSTR(IBEG),1) .NE. 0) THEN
               IIBEG = IBEG
               IBEG = IBEG + 1
               GO TO 20
C
            END IF
C
         ELSE
            GO TO 80
C
         END IF
C
      END IF
C
      IF (CODE .EQ. TZERR) THEN
         IF (INSTR(IBEG) .NE. PERIOD) GO TO 80
         IBEG = IBEG + 1
         IF (IBEG .GT. SLEN) GO TO 80
      END IF
C
      IF (I1X(LETTER,52,INSTR(IBEG),1) .EQ. 0) GO TO 80
      CHRSTR(1) = INSTR(IBEG)
C
      DO 30  I=2, 6
         IBEG = IBEG + 1
         IF (IBEG .GT. SLEN) GO TO 80
         IF (I1X(LETTER,52,INSTR(IBEG),1) .EQ. 0) GO TO 40
         CHRSTR(I) = INSTR(IBEG)
   30 CONTINUE
C
      GO TO 80
C
   40 CONTINUE
C
      DO 50  J=1, 15
         IF (I1CSTR(CHRSTR,I-1,TABLE(TABPTR(J)),TABPTR(J+1)-TABPTR(J))
     &        .EQ. 0) GO TO 60
   50 CONTINUE
C
      GO TO 80
C                                  STATE 4 - LOGICAL OPERATOR
   60 IF (J .GT. 2) THEN
C
         IF (CODE .EQ. TRCNST) THEN
C
            IF (INSTR(IBEG) .EQ. PERIOD) THEN
               CODE = TICNST
               IIBEG = IIBEG - 1
            END IF
C
            GO TO 80
C
         ELSE IF (INSTR(IBEG) .NE. PERIOD) THEN
            GO TO 80
C
         ELSE IF (FLAG) THEN
            GO TO 80
C
         ELSE
            CODE = TOKEN(J-2)
            IIBEG = IBEG
            GO TO 80
C
         END IF
C
      END IF
C                                  STATE 5 - DOUBLE PRECISION CONSTANT
      IF (CODE .NE. TRCNST) GO TO 80
      IF (INSTR(IBEG).EQ.MINUS .OR. INSTR(IBEG).EQ.PLUS) IBEG = IBEG +
     &    1
      IF (IBEG .GT. SLEN) GO TO 80
C
      IF (I1X(DIGIT,10,INSTR(IBEG),1) .EQ. 0) THEN
         GO TO 80
C
      ELSE
         IIBEG = IBEG
         IBEG = IBEG + 1
C
   70    IF (IBEG .LE. SLEN) THEN
C
            IF (I1X(DIGIT,10,INSTR(IBEG),1) .NE. 0) THEN
               IIBEG = IBEG
               IBEG = IBEG + 1
               GO TO 70
C
            END IF
C
         END IF
C
      END IF
C
      IF (J .EQ. 1) CODE = TDCNST
C
   80 CONTINUE
C
      IF (CODE .EQ. TZERR) THEN
         OLEN = 0
C
      ELSE
         OLEN = IIBEG
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  IMACH (SINGLE PRECISION VERSION)
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    MARCH 26, 1984
C
C  PURPOSE:    RETRIEVE INTEGER MACHINE CONSTANTS.
C
C  USAGE:      IMACH(N)
C
C  ARGUMENTS:
C     N      - INDEX OF DESIRED CONSTANT.  (INPUT)
C     IMACH  - MACHINE CONSTANT.  (OUTPUT)
C
C  REMARK:
C     FOLLOWING IS A DESCRIPTION OF THE ASSORTED INTEGER MACHINE
C     CONSTANTS.
C
C     WORDS
C
C        IMACH( 1) = NUMBER OF BITS PER INTEGER STORAGE UNIT.
C        IMACH( 2) = NUMBER OF CHARACTERS PER INTEGER STORAGE UNIT.
C
C     INTEGERS
C
C        ASSUME INTEGERS ARE REPRESENTED IN THE S-DIGIT, BASE-A FORM
C        SIGN ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
C        WHERE 0 .LE. X(I) .LT. A FOR I=0,...,S-1.  THEN
C
C        IMACH( 3) = A, THE BASE.
C        IMACH( 4) = S, NUMBER OF BASE-A DIGITS.
C        IMACH( 5) = A**S - 1, LARGEST MAGNITUDE.
C
C     FLOATING-POINT NUMBERS
C
C        ASSUME FLOATING-POINT NUMBERS ARE REPRESENTED IN THE T-DIGIT,
C        BASE-B FORM SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
C        WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T,
C        0 .LT. X(1), AND EMIN .LE. E .LE. EMAX.  THEN
C
C        IMACH( 6) = B, THE BASE.
C
C        SINGLE PRECISION
C
C           IMACH( 7) = T, NUMBER OF BASE-B DIGITS.
C           IMACH( 8) = EMIN, SMALLEST EXPONENT E.
C           IMACH( 9) = EMAX, LARGEST EXPONENT E.
C
C        DOUBLE PRECISION
C
C           IMACH(10) = T, NUMBER OF BASE-B DIGITS.
C           IMACH(11) = EMIN, SMALLEST EXPONENT E.
C           IMACH(12) = EMAX, LARGEST EXPONENT E.
C
C  GAMS:       R1
C
C  CHAPTERS:   MATH/LIBRARY REFERENCE MATERIAL
C              STAT/LIBRARY REFERENCE MATERIAL
C              SFUN/LIBRARY REFERENCE MATERIAL
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      INTEGER FUNCTION IMACH (N)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    N
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    NOUT
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      INTEGER    IMACHV(12)
      SAVE       IMACHV
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   UMACH
C                                  DEFINE CONSTANTS
      DATA IMACHV(1)/32/
      DATA IMACHV(2)/4/
      DATA IMACHV(3)/2/
      DATA IMACHV(4)/31/
      DATA IMACHV(5)/2147483647/
      DATA IMACHV(6)/2/
      DATA IMACHV(7)/24/
      DATA IMACHV(8)/-125/
      DATA IMACHV(9)/128/
      DATA IMACHV(10)/53/
      DATA IMACHV(11)/-1021/
      DATA IMACHV(12)/1024/
C
      IF (N.LT.1 .OR. N.GT.12) THEN
C                                  ERROR.  INVALID RANGE FOR N.
         CALL UMACH (2, NOUT)
         WRITE (NOUT,99999) N
99999    FORMAT (/, ' *** TERMINAL ERROR 5 FROM IMACH.  THE ARGUMENT',
     &          /, ' ***          MUST BE BETWEEN 1 AND 12 INCLUSIVE.'
     &          , /, ' ***          N = ', I6, '.', /)
         IMACH = 0
         STOP
C
      ELSE
         IMACH = IMACHV(N)
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  E1INPL
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    MARCH 2, 1984
C
C  PURPOSE:    TO STORE A CHARACTER STRING IN THE PARAMETER LIST PLIST
C              FOR USE BY THE ERROR MESSAGE HANDLER.
C
C  USAGE:      CALL E1INPL(FORM,NUM,SLEN,STRUP)
C
C  ARGUMENTS:
C     FORM   - A CHARACTER STRING OF LENGTH ONE TO BE INSERTED INTO
C              PLIST WHICH SPECIFIES THE FORM OF THE STRING.  (INPUT)
C              FOR EXAMPLE, 'L' FOR STRING, 'A' FOR CHARACTER ARRAY,
C              'I' FOR INTEGER, 'K' FOR KEYWORD (PROTRAN ONLY).  AN
C              ASTERISK IS INSERTED INTO PLIST PRECEDING FORM.
C     NUM    - INTEGER TO BE INSERTED AS A CHARACTER INTO PLIST
C              IMMEDIATELY FOLLOWING FORM.  (INPUT)  NUM MUST BE BETWEEN
C              1 AND 9.
C     SLEN   - THE NUMBER OF CHARACTERS IN STRUP.  (INPUT)  LEN MUST BE
C              LESS THAN OR EQUAL TO 255.  THE CHARACTER REPRESENTATION
C              OF SLEN IS INSERTED INTO PLIST AFTER NUM AND AN ASTERISK.
C     STRUP  - A CHARACTER STRING OF LENGTH LEN WHICH IS TO BE INSERTED
C              INTO PLIST.  (INPUT)  TRAILING BLANKS ARE IGNORED.
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE E1INPL (FORM, NUM, SLEN, STRUP)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    NUM, SLEN
      CHARACTER  FORM, STRUP(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    IER, L, LEN2, LENCK, LOC, NLEN, NNUM
      CHARACTER  STRNCH(3)
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      CHARACTER  BLANK, PRCNT(1), TEMP(4)
      SAVE       BLANK, PRCNT, TEMP
C                                  SPECIFICATIONS FOR SPECIAL CASES
C                              SPECIFICATIONS FOR COMMON /ERCOM1/
      INTEGER    CALLVL, MAXLEV, MSGLEN, ERTYPE(51), ERCODE(51),
     &           PRINTB(7), STOPTB(7), PLEN, IFERR6, IFERR7,
     &           IALLOC(51), HDRFMT(7), TRACON(7)
      COMMON     /ERCOM1/ CALLVL, MAXLEV, MSGLEN, ERTYPE, ERCODE,
     &           PRINTB, STOPTB, PLEN, IFERR6, IFERR7, IALLOC, HDRFMT,
     &           TRACON
      SAVE       /ERCOM1/
C                              SPECIFICATIONS FOR COMMON /ERCOM2/
      CHARACTER  MSGSAV(255), PLIST(300), RNAME(51)*6
      COMMON     /ERCOM2/ MSGSAV, PLIST, RNAME
      SAVE       /ERCOM2/
C                              SPECIFICATIONS FOR COMMON /ERCOM3/
      DOUBLE PRECISION ERCKSM
      COMMON     /ERCOM3/ ERCKSM
      SAVE       /ERCOM3/
C                              SPECIFICATIONS FOR COMMON /ERCOM4/
      LOGICAL    ISUSER(51)
      COMMON     /ERCOM4/ ISUSER
      SAVE       /ERCOM4/
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  IABS
      INTRINSIC  IABS
      INTEGER    IABS
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   C1TIC, M1VE
C
      DATA TEMP/'*', ' ', ' ', '*'/, PRCNT/'%'/, BLANK/' '/
C
      NNUM = IABS(NUM)
      LENCK = PLEN + SLEN + 8
      IF (NNUM.GE.1 .AND. NNUM.LE.9 .AND. LENCK.LE.300) THEN
         TEMP(2) = FORM
         CALL C1TIC (NNUM, TEMP(3), 1, IER)
         LOC = PLEN + 1
         IF (LOC .EQ. 2) LOC = 1
         CALL M1VE (TEMP, 1, 4, 4, PLIST(LOC), 1, 4, 262, IER)
         LOC = LOC + 4
         IF (NUM .LT. 0) THEN
            LEN2 = SLEN
         ELSE
            DO 10  L=1, SLEN
               LEN2 = SLEN - L + 1
               IF (STRUP(LEN2) .NE. BLANK) GO TO 20
   10       CONTINUE
            LEN2 = 1
   20       CONTINUE
         END IF
         NLEN = 1
         IF (LEN2 .GE. 10) NLEN = 2
         IF (LEN2 .GE. 100) NLEN = 3
         CALL C1TIC (LEN2, STRNCH, NLEN, IER)
         CALL M1VE (STRNCH, 1, NLEN, 3, PLIST(LOC), 1, NLEN, 262, IER)
         LOC = LOC + NLEN
         CALL M1VE (PRCNT, 1, 1, 1, PLIST(LOC), 1, 1, 262, IER)
         LOC = LOC + 1
         CALL M1VE (STRUP, 1, LEN2, LEN2, PLIST(LOC), 1, LEN2, 262,
     &              IER)
         PLEN = LOC + LEN2 - 1
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  UMACH (SINGLE PRECISION VERSION)
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    MARCH 21, 1984
C
C  PURPOSE:    SET OR RETRIEVE INPUT OR OUTPUT DEVICE UNIT NUMBERS.
C
C  USAGE:      CALL UMACH (N, NUNIT)
C
C  ARGUMENTS:
C     N      - INDEX OF DESIRED UNIT.  (INPUT)
C              THE VALUES OF N ARE DEFINED AS FOLLOWS:
C              N = 1, CORRESPONDS TO THE STANDARD INPUT UNIT.
C              N = 2, CORRESPONDS TO THE STANDARD OUTPUT UNIT.
C     NUNIT  - I/O UNIT.  (INPUT OR OUTPUT)
C              IF THE VALUE OF N IS NEGATIVE, THE UNIT CORRESPONDING
C              TO THE INDEX IS RESET TO THE VALUE GIVEN IN NUNIT.
C              OTHERWISE, THE VALUE CORRESPONDING TO THE INDEX IS
C              RETURNED IN NUNIT.
C
C  GAMS:       R1
C
C  CHAPTERS:   MATH/LIBRARY REFERENCE MATERIAL
C              STAT/LIBRARY REFERENCE MATERIAL
C              SFUN/LIBRARY REFERENCE MATERIAL
C
C  COPYRIGHT:  1984 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE UMACH (N, NUNIT)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    N, NUNIT
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    NN, NOUT
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      INTEGER    UNIT(2)
      SAVE       UNIT
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  IABS
      INTRINSIC  IABS
      INTEGER    IABS
C
      DATA UNIT(1)/5/
      DATA UNIT(2)/6/
C
      NN = IABS(N)
      IF (NN.NE.1 .AND. NN.NE.2) THEN
C                                  ERROR.  INVALID RANGE FOR N.
         NOUT = UNIT(2)
         WRITE (NOUT,99999) NN
99999    FORMAT (/, ' *** TERMINAL ERROR 5 FROM UMACH.  THE ABSOLUTE',
     &          /, ' ***          VALUE OF THE INDEX VARIABLE MUST BE'
     &          , /, ' ***          1 OR 2.  IABS(N) = ', I6,
     &          '.', /)
         STOP
C                                  CHECK FOR RESET OR RETRIEVAL
      ELSE IF (N .LT. 0) THEN
C                                  RESET
         UNIT(NN) = NUNIT
      ELSE
C                                  RETRIEVE
         NUNIT = UNIT(N)
      END IF
C
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  IACHAR (SINGLE PRECISION VERSION)
C
C  COMPUTER:   PCDSMS/SINGLE
C
C  REVISED:    SEPTEMBER 9, 1985
C
C  PURPOSE:    RETURN THE INTEGER ASCII VALUE OF A CHARACTER ARGUMENT.
C
C  USAGE:      IACHAR(CH)
C
C  ARGUMENTS:
C     CH     - CHARACTER ARGUMENT FOR WHICH THE INTEGER ASCII VALUE
C              IS DESIRED.  (INPUT)
C     IACHAR - INTEGER ASCII VALUE FOR CH.  (OUTPUT)
C              THE CHARACTER CH IS IN THE IACHAR-TH POSITION OF THE
C              ASCII COLLATING SEQUENCE.
C
C  GAMS:       N3
C
C  CHAPTER:    MATH/LIBRARY UTILITIES
C              STAT/LIBRARY UTILITIES
C
C  COPYRIGHT:  1986 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      INTEGER FUNCTION IACHAR (CH)
C                                  SPECIFICATIONS FOR ARGUMENTS
      CHARACTER  CH
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      IACHAR = ICHAR(CH)
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL NAME:  I1CSTR (SINGLE PRECISION VERSION)
C
C  COMPUTER:   DGC/SINGLE
C
C  REVISED:    SEPTEMBER 10, 1985
C
C  PURPOSE:    CASE INSENSITIVE COMPARISON OF TWO CHARACTER ARRAYS.
C
C  USAGE:      I1CSTR(STR1, LEN1, STR2, LEN2)
C
C  ARGUMENTS:
C     STR1   - FIRST CHARACTER ARRAY.  (INPUT)
C     LEN1   - LENGTH OF STR1.  (INPUT)
C     STR2   - SECOND CHARACTER ARRAY.  (INPUT)
C     LEN2   - LENGTH OF STR2.  (INPUT)
C     I1CSTR - INTEGER FUNCTION.  (OUTPUT) WHERE
C              I1CSTR = -1  IF STR1 .LT. STR2,
C              I1CSTR =  0  IF STR1 .EQ. STR2,
C              I1CSTR =  1  IF STR1 .GT. STR2.
C
C  REMARKS:
C  1. IF THE TWO ARRAYS, STR1 AND STR2,  ARE OF UNEQUAL LENGTH, THE
C     SHORTER ARRAY IS CONSIDERED AS IF IT WERE EXTENDED WITH BLANKS
C     TO THE LENGTH OF THE LONGER ARRAY.
C
C  2. IF ONE OR BOTH LENGTHS ARE ZERO OR NEGATIVE THE I1CSTR OUTPUT IS
C     BASED ON COMPARISON OF THE LENGTHS.
C
C  GAMS:       N5C
C
C  CHAPTER:    MATH/LIBRARY UTILITIES
C
C  COPYRIGHT:  1985 BY IMSL, INC.  ALL RIGHTS RESERVED.
C
C  WARRANTY:   IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN APPLIED
C              TO THIS CODE.  NO OTHER WARRANTY, EXPRESSED OR IMPLIED,
C              IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      INTEGER FUNCTION I1CSTR (STR1, LEN1, STR2, LEN2)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    LEN1, LEN2
      CHARACTER  STR1(LEN1), STR2(LEN2)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    IC1, IC2, ICB, IS, L, LENM
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  ISIGN,MIN0
      INTRINSIC  ISIGN, MIN0
      INTEGER    ISIGN, MIN0
C                                  SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL   ICASE
      INTEGER    ICASE
C
      IF (LEN1.GT.0 .AND. LEN2.GT.0) THEN
C                                  COMPARE FIRST LENM CHARACTERS
         LENM = MIN0(LEN1,LEN2)
         DO 10  L=1, LENM
            IC1 = ICASE(STR1(L))
            IC2 = ICASE(STR2(L))
            IF (IC1 .NE. IC2) THEN
               I1CSTR = ISIGN(1,IC1-IC2)
               RETURN
            END IF
   10    CONTINUE
      END IF
C                                  COMPARISON BASED ON LENGTH OR
C                                  TRAILING BLANKS
      IS = LEN1 - LEN2
      IF (IS .EQ. 0) THEN
         I1CSTR = 0
      ELSE
         IF (LEN1.LE.0 .OR. LEN2.LE.0) THEN
C                                  COMPARISON BASED ON LENGTH
            I1CSTR = ISIGN(1,IS)
         ELSE
C                                  COMPARISON BASED ON TRAILING BLANKS
C                                  TO EXTEND SHORTER ARRAY
            LENM = LENM + 1
            ICB = ICASE(' ')
            IF (IS .GT. 0) THEN
C                                  EXTEND STR2 WITH BLANKS
               DO 20  L=LENM, LEN1
                  IC1 = ICASE(STR1(L))
                  IF (IC1 .NE. ICB) THEN
                     I1CSTR = ISIGN(1,IC1-ICB)
                     RETURN
                  END IF
   20          CONTINUE
            ELSE
C                                  EXTEND STR1 WITH BLANKS
               DO 30  L=LENM, LEN2
                  IC2 = ICASE(STR2(L))
                  IF (ICB .NE. IC2) THEN
                     I1CSTR = ISIGN(1,ICB-IC2)
                     RETURN
                  END IF
   30          CONTINUE
            END IF
C
            I1CSTR = 0
         END IF
      END IF
C
      RETURN
      END
