real function cpsec()
      implicit real(a-h,o-z)
      cpsec=0
      return
      end
      DOUBLE PRECISION FUNCTION DBLINF(NRA,NCA,A,LDA,X,Y)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,*),X(*),Y(*),XTRAN(100),M(100)
      DO 10 J=1,NCA
         SUM=0D0
         DO 20 I=1,NRA
            SUM=SUM+X(I)*A(I,J)
 20         CONTINUE
         M(J)=SUM   
 10   CONTINUE
      SUM=0D0
      DO 30 I=1,NCA
         SUM=SUM+M(I)*Y(I)
 30      CONTINUE
      DBLINF=SUM
      RETURN
      END
      DOUBLE PRECISION FUNCTION DCHIDF(CHSQ,DF)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL QSTAT
      NWH=1
      CALL CDFCHI(NWH,DCHIDF,CHSQ,DF,QSTAT,BOUND)
      RETURN
      END


      DOUBLE PRECISION FUNCTION DCHIIN(P,DF)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER STATUS
      CALL CDFCHI(2,P,DCHIIN,DF,STATUS,BOUND)
      RETURN
      END
      real*8 function dconst(name)
      implicit real*8(a-h,o-z)
      character*15 name
      dconst=3.1415926535898d0
      return
      end
      SUBROUTINE DCRGRG(N,A,LDA,B,LDB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(LDA,*), B(LDB,*)
      CALL DLACPY('M',N,N,A,LDA,B,LDB)
      RETURN
      END
      SUBROUTINE DEVLSF(N,A,LDA,EVAL)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(LDD1SCR=10000)
      COMMON/DSCRB/D1SCRB(LDD1SCR)
      PARAMETER(LDD2SCR=100)
      COMMON/DSCRB/D2SCRB(LDD2SCR,LDD2SCR)
      DIMENSION A(LDA,*), EVAL(*),info(100)
      IF(N.GT.LDD1SCR) GO TO 100
      IF(N.GT.LDD2SCR) GO TO 110
      CALL DCRGRG(N,A,LDA,D2SCRB,LDD2SCR)
      CALL DSYEV('N','U',N,D2SCRB,LDD2SCR,EVAL,D1SCRB,LDD1SCR,INFO)
      RETURN
 100  CALL ERR(N,'LDD1SCR','D1SCRBLK')
      RETURN
 110  CALL ERR(N,'LDD2SCR','D2SCRBLK')
      RETURN
      END
      SUBROUTINE DFULL1(N,A,LDA)
C  This routine converts a symmetric matrix in IMSL symmetric storeage mode
C  (upper triangular) to one in full storeage mode.
C  Inputs:
C  N   Order of symmetric matrix
C  A   Matrix in symmetric storeage mode
C  LDA Leading dimension of A
C  Output:
C  A   Matrix in full storeage mode
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(LDA,N)
      N1=N-1
      DO 10 I=1,N1
      DO 10 J=I+1,N
   10 A(J,I)=A(I,J)
      RETURN
      END
      SUBROUTINE DLFTDS(N,A,LDA,FAC,LDFAC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,*),FAC(LDFAC,*)
      CALL DCRGRG(N,A,LDA,FAC,LDFAC)
      CALL DPOTRF('U',N,FAC,LDFAC,INFO)
      RETURN
      END
      SUBROUTINE DLINDS(N,A,LDA,AINV,LDAINV)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,*),AINV(LDAINV,*)
      CALL DLFTDS(N,A,LDA,AINV,LDAINV)
      CALL DPOTRI('U',N,AINV,LDAINV,INFO)
      DO 30 I=1,N
         DO 40 J=1,I
            AINV(I,J)=AINV(J,I)
 40          CONTINUE
 30          CONTINUE
      RETURN
      END
      SUBROUTINE DLINRG(N,A,LDA,AINV,LDAINV)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      parameter(ldiscr=10000)
      common/iscrb/iscrb(ldiscr)
      parameter(ldd1scr=10000)
      common/dscrb/d1scrb(ldd1scr)
      DIMENSION A(LDA,*),AINV(LDAINV,*)
      if(n.gt.ldiscr) go to 100
      if(n.gt.ldd1scr) go to 150
      DO 10 I=1,N
        DO 20 J=1,N
           AINV(I,J)=A(I,J)
 20     CONTINUE
 10   CONTINUE
      CALL DGETRF(N,N,AINV,LDAINV,iscrb,INFO)
      CALL DGETRI(N,AINV,LDAINV,iscrb,d1scrb,ldd1scr,INFO)
      return
  100 call err(n,'ldiscr','iscrblk')
      return
  150 call err(n,'ldd1scr','d1scrblk')
      RETURN
      END
      SUBROUTINE DLINRT(N,A,LDA,IPATH,AINV,LDAINV)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(LDA,*),AINV(LDAINV,*)
      DO 10 I=1,N
         DO 20 J=1,N
            AINV(I,J)=A(I,J)
 20         CONTINUE
 10      CONTINUE   
      IF (IPATH.EQ.1) THEN
      CALL DTRTRI('L','N',N,AINV,LDAINV,INFO)
      ELSE
      CALL DTRTRI('U','N',N,AINV,LDAINV,INFO)
      ENDIF
      RETURN
      END
      SUBROUTINE DLSLDS(N,A,LDA,B,X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter(ldd2scr=100)
      common/dscrb/d2scrb(ldd2scr,ldd2scr)
      DIMENSION A(LDA,*),B(*),X(*)
      if(n.gt.ldd2scr) go to 100
      DO 10 I=1,N
         DO 20 J=1,N
            d2scrb(I,J)=A(I,J)
 20         CONTINUE
 10         CONTINUE
      DO 30 I=1,N
         X(I)=B(I)
 30      CONTINUE
      CALL DPOTRF('U',N,d2scrb,ldd2scr,INFO)
      CALL DPOTRS('U',N,1,d2scrb,ldd2scr,X,N,INFO)
      return
 100  call err(n,'ldd2scr','d2scrblk')
      RETURN
      END
      SUBROUTINE DLSLRG(N,A,LDA,B,IPATH,X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      parameter(ldiscr=10000)
      common/iscrb/iscrb(ldiscr)
      parameter(ldd2scr=100)
      common/dscrb/d2scrb(ldd2scr,ldd2scr)
      DIMENSION A(LDA,*),B(N),X(N)
      if(n.gt.ldiscr) go to 100
      if(n.gt.ldd2scr) go to 150
      DO 10 I=1,N
         X(I)=B(I)
 10      CONTINUE
      DO 20 I=1,N
         DO 30 J=1,N
            d2scrb(I,J)=A(I,J)
 30         CONTINUE
 20      CONTINUE
      CALL DGETRF(N,N,d2scrb,ldd2scr,iscrb,INFO)
      IF (IPATH.EQ.1) THEN
         CALL DGETRS('N',N,1,d2scrb,ldd2scr,iscrb,X,N,INFO)
      ELSE
         CALL DGETRS('T',N,1,d2scrb,ldd2scr,iscrb,X,N,INFO)
      ENDIF
      return
  100 call err(n,'ldiscr','iscrblk')
      return
  150 call err(n,'ldd2scr','d2scrblk')
      RETURN
      END
      SUBROUTINE DLSLTO(N,A,B,IPATH,X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      parameter(ldd1scr=10000)
      common/dscrb/d1scrb(ldd1scr)
      DIMENSION A(2*N-1),B(N),X(N),R(1000)
      if(n.gt.ldd1scr) go to 150
      DO 10 I=1,N
         X(I)=B(I)
 10      CONTINUE
      CALL TSLD(A,X,R,N)
      RETURN
  150 call err(2*n-2,'ldd1scr','d1scrblk')
      RETURN 
      END
      DOUBLE PRECISION FUNCTION DMACH(I)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IF (I.EQ.7) THEN
         DMACH=1.7976931348623D+308
      ENDIF
      IF (I.EQ.8) THEN
         DMACH=-1.7976931348623D+308
         ENDIF
      IF (I.EQ.1) THEN
         DMACH=0.22251D-307
      ENDIF
      IF (I.EQ.4) THEN
         DMACH=0.22204D-15		
      ENDIF
      RETURN
      END
      SUBROUTINE DMRRRR(NRA,NCA,A,LDA,NRB,NCB,B,LDB,
     +NRC,NCC,C,LDC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,*),B(LDB,*),C(LDC,*)
      DO 10 I=1,NRA
         DO 20 J=1,NCB
            SUM=0D0
            DO 30 K=1,NCA
               SUM=SUM+A(I,K)*B(K,J)
 30         CONTINUE
            C(I,J)=SUM
 20      CONTINUE
 10   CONTINUE
      RETURN
      END
      SUBROUTINE DMURRV(NRA,NCA,A,LDA,NX,X,IPATH,NY,Y)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(LDA,*),X(NX),Y(NY)
      IF (IPATH.EQ.1) THEN
         DO 10 I=1,NRA
            SUM=0D0
            DO 20 J=1,NCA
               SUM=SUM+A(I,J)*X(J)
 20         CONTINUE
            Y(I)=SUM
 10      CONTINUE
      ENDIF
      IF (IPATH.EQ.2) THEN
         DO 30 J=1,NCA
            SUM=0D0
            DO 40 I=1,NRA
               SUM=SUM+A(I,J)*X(I)
 40         CONTINUE
            Y(J)=SUM
 30      CONTINUE
      ENDIF
      RETURN
      END
      SUBROUTINE DMXTXF(M,N,A,LDA,nb,B,LDB)
      IMPLICIT double precision (A-H,O-Z)
      DIMENSION A(LDA,N),B(LDB,N)
      DO 20 J=1,N
      DO 20 I=j,N
      AA=0.0D0
      DO 10 L=1,M
   10 AA=AA+A(L,I)*A(L,J)
   20 B(j,i)=AA
      call dfull1(n,b,ldb)
      RETURN
      END
      SUBROUTINE DMXTyF(NRA,NCA,A,LDA,NRB,NCB,B,LDB,
     +NRC,NCC,C,LDC)
      IMPLICIT double precision (A-H,O-Z)
      DIMENSION A(LDA,NCA),B(LDB,NCB), C(LDc,NCC)
      DO 20 I=1,Nca
      DO 20 J=1,Ncb
      AA=0.0D0
      DO 10 L=1,nra
   10 AA=AA+A(L,I)*b(L,J)
   20 c(i,j)=AA
      RETURN
      END
      SUBROUTINE DMXyTF(NRA,NCA,A,LDA,NRB,NCB,B,LDB,
     +NRC,NCC,C,LDC)
      IMPLICIT double precision (A-H,O-Z)
      DIMENSION A(LDA,NCA),B(LDB,NCB), C(LDc,NCC)
      DO 20 I=1,Nra
      DO 20 J=1,Nrb
      AA=0.0D0
      DO 10 L=1,nca
   10 AA=AA+a(I,L)*b(J,L)
   20 c(i,j)=AA
      RETURN
      END
      DOUBLE PRECISION FUNCTION DNORDF(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      CALL CDFNOR(1,DNORDF,X,0D0,1D0,N,BOUND)
      RETURN
      END
      
      SUBROUTINE DQDAG(F,A,B,ERRABS,ERRREL,IRULE,RESULT,ERREST)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (LIMIT=500,LENW=3000)
      DIMENSION WORK(LENW),IWORK(LIMIT) 
      EXTERNAL F
      CALL DQAG(F,A,B,ERRABS,ERRREL,IRULE,RESULT,ABSERR,NEVAL,IER,
     +LIMIT,LENW,LAST,IWORK,WORK)
      RETURN
      END
       SUBROUTINE DRNCHI(NR,DF,R)
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DIMENSION R(nr)
       DO 10 I=1,NR
       R(I)=GENCHI(DF)
 10    CONTINUE
       RETURN
       END

        
      SUBROUTINE DRNMVN(NR,K,RSIG,LDRSIG,R,LDR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION RSIG(LDRSIG,*),R(LDR,*)
      call um0set(nr,k,r,ldr)
      do 20 ii=1,nr
      DO 10 J=1,K
      EPS=DRNNOF()
      DO 10 I=j,k
   10 r(ii,I)=r(ii,I)+rsig(j,i)*EPS
   20 continue
      RETURN
      END
      SUBROUTINE DRNNOA(NR,R)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION R(NR)
      AV=0D0
      SD=1D0
      DO 10 I=1,NR
         R(I)=GENNOR(AV,SD)
 10      CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION DRNNOF()
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DRNNOF=GENNOR(0D0,1D0)
      CALL ADVNST(7)
      RETURN
      END
      DOUBLE PRECISION FUNCTION DRNUNF()
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DRNUNF=RANF()
      CALL ADVNST(7)
      RETURN
      END
      DOUBLE PRECISION FUNCTION DSUM(N,DX,INCX)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION DX(*)
      IF (N.LE.0D0) DSUM=0D0
      SUM=0D0
      DO 10 I=1,N,INCX
         SUM=SUM+DX(I)
 10      CONTINUE
      DSUM=SUM
      RETURN
      END
      SUBROUTINE DSVRGN(N,RA,RB)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      parameter(ldiscr=10000)
      common/iscrb/iscrb(ldiscr)
      DIMENSION RA(N),RB(N)
      if(n.gt.ldiscr) go to 100
      DO 10 I=1,N
      RB(I)=RA(I)
 10   CONTINUE
      CALL DPSORT(RB,N,iscrb,2,IER)
      RETURN
 100  call err(n,'ldiscr','iscrblk')
      return      
      END

      SUBROUTINE DSVRGP(N,RA,RB,IPERM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION RA(N),RB(N),IPERM(N)
      DO 10 I=1,N
      RB(I)=RA(I)
 10   CONTINUE
      CALL DPSORT(RB,N,IPERM,2,IER)
      RETURN
      END
      SUBROUTINE DTRNRR(NRA,NCA,A,LDA,NRB,NCB,B,LDB)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,*),B(LDB,*)
      DO 10 I=1,NCA
         DO 20 J=1,NRA
            B(I,J)=A(J,I)
 20         CONTINUE
 10         CONTINUE
      RETURN
      END
      SUBROUTINE DWRRRN(TITLE,NRA,NCA,A,LDA,ITRING)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*(*) TITLE
      DIMENSION A(LDA,*)
      IF (ITRING.EQ.0) THEN
         WRITE(6,*) TITLE
         DO 10 I=1,NRA
            WRITE(6,999) (A(I,J),J=1,NCA)
 10         CONTINUE
      ENDIF
      IF (ITRING.EQ.1) THEN
         WRITE(6,991) TITLE
         DO 20 I=1,NRA
            WRITE(6,999) (A(I,J),J=I,NCA)
 20         CONTINUE
      ENDIF
      IF (ITRING.EQ.-1) THEN
         WRITE(6,992) TITLE
         DO 30 I=1,NRA
            WRITE(6,999) (A(I,J),J=1,I)
 30      CONTINUE
      ENDIF
      IF (ITRING.EQ.2) THEN
         WRITE(6,993) TITLE
         DO 40 I=1,NRA
            WRITE(6,999) (A(I,J),J=I+1,NCA)
 40         CONTINUE
      ENDIF
      IF (ITRING.EQ.-2) THEN
         WRITE(6,994) TITLE
         DO 50 I=1,NRA
            WRITE(6,999) (A(I,J),J=1,I-1)
 50         CONTINUE
      ENDIF      
 999  FORMAT(1X,15(1X,F8.5))
 991  FORMAT(1X,'THE UPPER TRIANGULAR PART OF ',A5,/,'
     +  :INCLUDING DIAGONAL')
 992  FORMAT(1X,'THE LOWER TRIANGULAR PART OF ',A5,/,'
     +  :INCLUDING DIAGONAL')
 993  FORMAT(1X,'THE UPPER TRIANGULAR PART OF ',A5,/,'
     +  :EXCLUDING DIAGONAL')
 994  FORMAT(1X,'THE LOWER TRIANGULAR PART OF ',A5,/,'
     +  :EXCLUDING DIAGONAL')
        RETURN  
      END
      SUBROUTINE DZBREN(F,ERRABS,ERRREL,A,B,MAXFN)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      EXTERNAL F
      CALL ZBRENT(F,A,B,ZERO,ERRABS,ERRREL,MAXFN,ITER1)
      B=ZERO
      MAXFN=ITER1
      RETURN
      END
      SUBROUTINE DZPORC(NDEG,COEFF,ROOT)
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 ROOT
      LOGICAL FAIL
      DIMENSION COEFF(101),ROOT(NDEG),DCOEFF(101),ZEROR(100),
     +ZEROI(100)
      EXTERNAL RPOLY
      DO 10 I=1,NDEG+1
         DCOEFF(I)=COEFF(NDEG+1-(I-1))
 10   CONTINUE
      CALL RPOLY(DCOEFF,NDEG,ZEROR,ZEROI,FAIL)
      DO 20 I=1,NDEG
         ROOT(I)=DCMPLX(ZEROR(I),ZEROI(I))
 20   CONTINUE
      RETURN
      END
      
      
         
      
      subroutine err(n,a1,a2)
      character*(*)a1,a2
      write(6,10)n,a1
   10 format(/,'***', i10, ' exceeds the constant ', a)
      write(6,20)a2,a1
   20 format('    Edit the file ',a, ' to change ',
     +a, ', recompile and link.')
      return
      end
      INTEGER FUNCTION IDYWK(IDATE,MONTH,IYEAR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      IDYWK=0
      RETURN
      END
      real function cpsec()
      implicit real(a-h,o-z)
      cpsec=0
      return
      end
      DOUBLE PRECISION FUNCTION DBLINF(NRA,NCA,A,LDA,X,Y)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,*),X(*),Y(*),XTRAN(100),M(100)
      DO 10 J=1,NCA
         SUM=0D0
         DO 20 I=1,NRA
            SUM=SUM+X(I)*A(I,J)
 20         CONTINUE
         M(J)=SUM   
 10   CONTINUE
      SUM=0D0
      DO 30 I=1,NCA
         SUM=SUM+M(I)*Y(I)
 30      CONTINUE
      DBLINF=SUM
      RETURN
      END
      DOUBLE PRECISION FUNCTION DCHIDF(CHSQ,DF)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL QSTAT
      NWH=1
      CALL CDFCHI(NWH,DCHIDF,CHSQ,DF,QSTAT,BOUND)
      RETURN
      END


      DOUBLE PRECISION FUNCTION DCHIIN(P,DF)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER STATUS
      CALL CDFCHI(2,P,DCHIIN,DF,STATUS,BOUND)
      RETURN
      END
      real*8 function dconst(name)
      implicit real*8(a-h,o-z)
      character*15 name
      dconst=3.1415926535898d0
      return
      end
      SUBROUTINE DCRGRG(N,A,LDA,B,LDB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(LDA,*), B(LDB,*)
      CALL DLACPY('M',N,N,A,LDA,B,LDB)
      RETURN
      END
      SUBROUTINE DEVLSF(N,A,LDA,EVAL)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(LDD1SCR=10000)
      COMMON/DSCRB/D1SCRB(LDD1SCR)
      PARAMETER(LDD2SCR=100)
      COMMON/DSCRB/D2SCRB(LDD2SCR,LDD2SCR)
      DIMENSION A(LDA,*), EVAL(*),info(100)
      IF(N.GT.LDD1SCR) GO TO 100
      IF(N.GT.LDD2SCR) GO TO 110
      CALL DCRGRG(N,A,LDA,D2SCRB,LDD2SCR)
      CALL DSYEV('N','U',N,D2SCRB,LDD2SCR,EVAL,D1SCRB,LDD1SCR,INFO)
      RETURN
 100  CALL ERR(N,'LDD1SCR','D1SCRBLK')
      RETURN
 110  CALL ERR(N,'LDD2SCR','D2SCRBLK')
      RETURN
      END
      SUBROUTINE DFULL1(N,A,LDA)
C  This routine converts a symmetric matrix in IMSL symmetric storeage mode
C  (upper triangular) to one in full storeage mode.
C  Inputs:
C  N   Order of symmetric matrix
C  A   Matrix in symmetric storeage mode
C  LDA Leading dimension of A
C  Output:
C  A   Matrix in full storeage mode
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(LDA,N)
      N1=N-1
      DO 10 I=1,N1
      DO 10 J=I+1,N
   10 A(J,I)=A(I,J)
      RETURN
      END
      SUBROUTINE DLFTDS(N,A,LDA,FAC,LDFAC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,*),FAC(LDFAC,*)
      CALL DCRGRG(N,A,LDA,FAC,LDFAC)
      CALL DPOTRF('U',N,FAC,LDFAC,INFO)
      RETURN
      END
      SUBROUTINE DLINDS(N,A,LDA,AINV,LDAINV)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,*),AINV(LDAINV,*)
      CALL DLFTDS(N,A,LDA,AINV,LDAINV)
      CALL DPOTRI('U',N,AINV,LDAINV,INFO)
      DO 30 I=1,N
         DO 40 J=1,I
            AINV(I,J)=AINV(J,I)
 40          CONTINUE
 30          CONTINUE
      RETURN
      END
      SUBROUTINE DLINRG(N,A,LDA,AINV,LDAINV)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      parameter(ldiscr=10000)
      common/iscrb/iscrb(ldiscr)
      parameter(ldd1scr=10000)
      common/dscrb/d1scrb(ldd1scr)
      DIMENSION A(LDA,*),AINV(LDAINV,*)
      if(n.gt.ldiscr) go to 100
      if(n.gt.ldd1scr) go to 150
      DO 10 I=1,N
        DO 20 J=1,N
           AINV(I,J)=A(I,J)
 20     CONTINUE
 10   CONTINUE
      CALL DGETRF(N,N,AINV,LDAINV,iscrb,INFO)
      CALL DGETRI(N,AINV,LDAINV,iscrb,d1scrb,ldd1scr,INFO)
      return
  100 call err(n,'ldiscr','iscrblk')
      return
  150 call err(n,'ldd1scr','d1scrblk')
      RETURN
      END
      SUBROUTINE DLINRT(N,A,LDA,IPATH,AINV,LDAINV)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(LDA,*),AINV(LDAINV,*)
      DO 10 I=1,N
         DO 20 J=1,N
            AINV(I,J)=A(I,J)
 20         CONTINUE
 10      CONTINUE   
      IF (IPATH.EQ.1) THEN
      CALL DTRTRI('L','N',N,AINV,LDAINV,INFO)
      ELSE
      CALL DTRTRI('U','N',N,AINV,LDAINV,INFO)
      ENDIF
      RETURN
      END
      SUBROUTINE DLSLDS(N,A,LDA,B,X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter(ldd2scr=100)
      common/dscrb/d2scrb(ldd2scr,ldd2scr)
      DIMENSION A(LDA,*),B(*),X(*)
      if(n.gt.ldd2scr) go to 100
      DO 10 I=1,N
         DO 20 J=1,N
            d2scrb(I,J)=A(I,J)
 20         CONTINUE
 10         CONTINUE
      DO 30 I=1,N
         X(I)=B(I)
 30      CONTINUE
      CALL DPOTRF('U',N,d2scrb,ldd2scr,INFO)
      CALL DPOTRS('U',N,1,d2scrb,ldd2scr,X,N,INFO)
      return
 100  call err(n,'ldd2scr','d2scrblk')
      RETURN
      END
      SUBROUTINE DLSLRG(N,A,LDA,B,IPATH,X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      parameter(ldiscr=10000)
      common/iscrb/iscrb(ldiscr)
      parameter(ldd2scr=100)
      common/dscrb/d2scrb(ldd2scr,ldd2scr)
      DIMENSION A(LDA,*),B(N),X(N)
      if(n.gt.ldiscr) go to 100
      if(n.gt.ldd2scr) go to 150
      DO 10 I=1,N
         X(I)=B(I)
 10      CONTINUE
      DO 20 I=1,N
         DO 30 J=1,N
            d2scrb(I,J)=A(I,J)
 30         CONTINUE
 20      CONTINUE
      CALL DGETRF(N,N,d2scrb,ldd2scr,iscrb,INFO)
      IF (IPATH.EQ.1) THEN
         CALL DGETRS('N',N,1,d2scrb,ldd2scr,iscrb,X,N,INFO)
      ELSE
         CALL DGETRS('T',N,1,d2scrb,ldd2scr,iscrb,X,N,INFO)
      ENDIF
      return
  100 call err(n,'ldiscr','iscrblk')
      return
  150 call err(n,'ldd2scr','d2scrblk')
      RETURN
      END
      SUBROUTINE DLSLTO(N,A,B,IPATH,X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      parameter(ldd1scr=10000)
      common/dscrb/d1scrb(ldd1scr)
      DIMENSION A(2*N-1),B(N),X(N),R(1000)
      if(n.gt.ldd1scr) go to 150
      DO 10 I=1,N
         X(I)=B(I)
 10      CONTINUE
      CALL TSLD(A,X,R,N)
      RETURN
  150 call err(2*n-2,'ldd1scr','d1scrblk')
      RETURN 
      END
      DOUBLE PRECISION FUNCTION DMACH(I)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IF (I.EQ.7) THEN
         DMACH=1.7976931348623D+308
      ENDIF
      IF (I.EQ.8) THEN
         DMACH=-1.7976931348623D+308
         ENDIF
      IF (I.EQ.1) THEN
         DMACH=0.22251D-307
      ENDIF
      IF (I.EQ.4) THEN
         DMACH=0.22204D-15		
      ENDIF
      RETURN
      END
      SUBROUTINE DMRRRR(NRA,NCA,A,LDA,NRB,NCB,B,LDB,
     +NRC,NCC,C,LDC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,*),B(LDB,*),C(LDC,*)
      DO 10 I=1,NRA
         DO 20 J=1,NCB
            SUM=0D0
            DO 30 K=1,NCA
               SUM=SUM+A(I,K)*B(K,J)
 30         CONTINUE
            C(I,J)=SUM
 20      CONTINUE
 10   CONTINUE
      RETURN
      END
      SUBROUTINE DMURRV(NRA,NCA,A,LDA,NX,X,IPATH,NY,Y)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(LDA,*),X(NX),Y(NY)
      IF (IPATH.EQ.1) THEN
         DO 10 I=1,NRA
            SUM=0D0
            DO 20 J=1,NCA
               SUM=SUM+A(I,J)*X(J)
 20         CONTINUE
            Y(I)=SUM
 10      CONTINUE
      ENDIF
      IF (IPATH.EQ.2) THEN
         DO 30 J=1,NCA
            SUM=0D0
            DO 40 I=1,NRA
               SUM=SUM+A(I,J)*X(I)
 40         CONTINUE
            Y(J)=SUM
 30      CONTINUE
      ENDIF
      RETURN
      END
      SUBROUTINE DMXTXF(M,N,A,LDA,nb,B,LDB)
      IMPLICIT double precision (A-H,O-Z)
      DIMENSION A(LDA,N),B(LDB,N)
      DO 20 J=1,N
      DO 20 I=j,N
      AA=0.0D0
      DO 10 L=1,M
   10 AA=AA+A(L,I)*A(L,J)
   20 B(j,i)=AA
      call dfull1(n,b,ldb)
      RETURN
      END
      SUBROUTINE DMXTyF(NRA,NCA,A,LDA,NRB,NCB,B,LDB,
     +NRC,NCC,C,LDC)
      IMPLICIT double precision (A-H,O-Z)
      DIMENSION A(LDA,NCA),B(LDB,NCB), C(LDc,NCC)
      DO 20 I=1,Nca
      DO 20 J=1,Ncb
      AA=0.0D0
      DO 10 L=1,nra
   10 AA=AA+A(L,I)*b(L,J)
   20 c(i,j)=AA
      RETURN
      END
      SUBROUTINE DMXyTF(NRA,NCA,A,LDA,NRB,NCB,B,LDB,
     +NRC,NCC,C,LDC)
      IMPLICIT double precision (A-H,O-Z)
      DIMENSION A(LDA,NCA),B(LDB,NCB), C(LDc,NCC)
      DO 20 I=1,Nra
      DO 20 J=1,Nrb
      AA=0.0D0
      DO 10 L=1,nca
   10 AA=AA+a(I,L)*b(J,L)
   20 c(i,j)=AA
      RETURN
      END
      DOUBLE PRECISION FUNCTION DNORDF(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      CALL CDFNOR(1,DNORDF,X,0D0,1D0,N,BOUND)
      RETURN
      END
      
      SUBROUTINE DQDAG(F,A,B,ERRABS,ERRREL,IRULE,RESULT,ERREST)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (LIMIT=500,LENW=3000)
      DIMENSION WORK(LENW),IWORK(LIMIT) 
      EXTERNAL F
      CALL DQAG(F,A,B,ERRABS,ERRREL,IRULE,RESULT,ABSERR,NEVAL,IER,
     +LIMIT,LENW,LAST,IWORK,WORK)
      RETURN
      END
       SUBROUTINE DRNCHI(NR,DF,R)
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DIMENSION R(nr)
       DO 10 I=1,NR
       R(I)=GENCHI(DF)
 10    CONTINUE
       RETURN
       END

        
      SUBROUTINE DRNMVN(NR,K,RSIG,LDRSIG,R,LDR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION RSIG(LDRSIG,*),R(LDR,*)
      call um0set(nr,k,r,ldr)
      do 20 ii=1,nr
      DO 10 J=1,K
      EPS=DRNNOF()
      DO 10 I=j,k
   10 r(ii,I)=r(ii,I)+rsig(j,i)*EPS
   20 continue
      RETURN
      END
      SUBROUTINE DRNNOA(NR,R)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION R(NR)
      AV=0D0
      SD=1D0
      DO 10 I=1,NR
         R(I)=GENNOR(AV,SD)
 10      CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION DRNNOF()
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DRNNOF=GENNOR(0D0,1D0)
      CALL ADVNST(7)
      RETURN
      END
      DOUBLE PRECISION FUNCTION DRNUNF()
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DRNUNF=RANF()
      CALL ADVNST(7)
      RETURN
      END
      DOUBLE PRECISION FUNCTION DSUM(N,DX,INCX)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION DX(*)
      IF (N.LE.0D0) DSUM=0D0
      SUM=0D0
      DO 10 I=1,N,INCX
         SUM=SUM+DX(I)
 10      CONTINUE
      DSUM=SUM
      RETURN
      END
      SUBROUTINE DSVRGN(N,RA,RB)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      parameter(ldiscr=10000)
      common/iscrb/iscrb(ldiscr)
      DIMENSION RA(N),RB(N)
      if(n.gt.ldiscr) go to 100
      DO 10 I=1,N
      RB(I)=RA(I)
 10   CONTINUE
      CALL DPSORT(RB,N,iscrb,2,IER)
      RETURN
 100  call err(n,'ldiscr','iscrblk')
      return      
      END

      SUBROUTINE DSVRGP(N,RA,RB,IPERM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION RA(N),RB(N),IPERM(N)
      DO 10 I=1,N
      RB(I)=RA(I)
 10   CONTINUE
      CALL DPSORT(RB,N,IPERM,2,IER)
      RETURN
      END
      SUBROUTINE DTRNRR(NRA,NCA,A,LDA,NRB,NCB,B,LDB)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,*),B(LDB,*)
      DO 10 I=1,NCA
         DO 20 J=1,NRA
            B(I,J)=A(J,I)
 20         CONTINUE
 10         CONTINUE
      RETURN
      END
      SUBROUTINE DWRRRN(TITLE,NRA,NCA,A,LDA,ITRING)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*(*) TITLE
      DIMENSION A(LDA,*)
      IF (ITRING.EQ.0) THEN
         WRITE(6,*) TITLE
         DO 10 I=1,NRA
            WRITE(6,999) (A(I,J),J=1,NCA)
 10         CONTINUE
      ENDIF
      IF (ITRING.EQ.1) THEN
         WRITE(6,991) TITLE
         DO 20 I=1,NRA
            WRITE(6,999) (A(I,J),J=I,NCA)
 20         CONTINUE
      ENDIF
      IF (ITRING.EQ.-1) THEN
         WRITE(6,992) TITLE
         DO 30 I=1,NRA
            WRITE(6,999) (A(I,J),J=1,I)
 30      CONTINUE
      ENDIF
      IF (ITRING.EQ.2) THEN
         WRITE(6,993) TITLE
         DO 40 I=1,NRA
            WRITE(6,999) (A(I,J),J=I+1,NCA)
 40         CONTINUE
      ENDIF
      IF (ITRING.EQ.-2) THEN
         WRITE(6,994) TITLE
         DO 50 I=1,NRA
            WRITE(6,999) (A(I,J),J=1,I-1)
 50         CONTINUE
      ENDIF      
 999  FORMAT(1X,15(1X,F8.5))
 991  FORMAT(1X,'THE UPPER TRIANGULAR PART OF ',A5,/,'
     +  :INCLUDING DIAGONAL')
 992  FORMAT(1X,'THE LOWER TRIANGULAR PART OF ',A5,/,'
     +  :INCLUDING DIAGONAL')
 993  FORMAT(1X,'THE UPPER TRIANGULAR PART OF ',A5,/,'
     +  :EXCLUDING DIAGONAL')
 994  FORMAT(1X,'THE LOWER TRIANGULAR PART OF ',A5,/,'
     +  :EXCLUDING DIAGONAL')
        RETURN  
      END
      SUBROUTINE DZBREN(F,ERRABS,ERRREL,A,B,MAXFN)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      EXTERNAL F
      CALL ZBRENT(F,A,B,ZERO,ERRABS,ERRREL,MAXFN,ITER1)
      B=ZERO
      MAXFN=ITER1
      RETURN
      END
      SUBROUTINE DZPORC(NDEG,COEFF,ROOT)
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 ROOT
      LOGICAL FAIL
      DIMENSION COEFF(101),ROOT(NDEG),DCOEFF(101),ZEROR(100),
     +ZEROI(100)
      EXTERNAL RPOLY
      DO 10 I=1,NDEG+1
         DCOEFF(I)=COEFF(NDEG+1-(I-1))
 10   CONTINUE
      CALL RPOLY(DCOEFF,NDEG,ZEROR,ZEROI,FAIL)
      DO 20 I=1,NDEG
         ROOT(I)=DCMPLX(ZEROR(I),ZEROI(I))
 20   CONTINUE
      RETURN
      END
      
      
         
      
      subroutine err(n,a1,a2)
      character*(*)a1,a2
      write(6,10)n,a1
   10 format(/,'***', i10, ' exceeds the constant ', a)
      write(6,20)a2,a1
   20 format('    Edit the file ',a, ' to change ',
     +a, ', recompile and link.')
      return
      end
      INTEGER FUNCTION IDYWK(IDATE,MONTH,IYEAR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      IDYWK=0
      RETURN
      END
      SUBROUTINE RNSET(ISEED)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CALL SETALL(123457,124579)
      RETURN
      END
      SUBROUTINE TDATE(IDATE,MONTH,IYEAR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      IDATE=0
      MONTH=0
      IYEAR=0
      RETURN
      END
      SUBROUTINE TIMDY(IHOUR,MINUTE,ISEC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      IHOUR=0
      MINUTE=0
      ISEC=0
      RETURN
      END
      SUBROUTINE RNSET(ISEED)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CALL SETALL(123457,124579)
      RETURN
      END
      SUBROUTINE TDATE(IDATE,MONTH,IYEAR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      IDATE=0
      MONTH=0
      IYEAR=0
      RETURN
      END
      SUBROUTINE TIMDY(IHOUR,MINUTE,ISEC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      IHOUR=0
      MINUTE=0
      ISEC=0
      RETURN
      END