
      PROGRAM DRIVER
C
C     TITLE: TEST OF SUBROUTINE 'LAPLAC'
C     ----------------------------------
C
C     THIS IS AN INTERACTIVE DRIVER PROGRAM FOR TESTING THE SUBROUTINE
C     'LAPLAC'.  THE SUBROUTINE SOLVES THE INTERIOR DIRICHLET PROBLEM
C     FOR LAPLACE'S EQUATION ON A SIMPLY-CONNECTED REGION D WITH A
C     SMOOTH BOUNDARY S.  THIS DRIVER PROGRAM IS SET UP TO SOLVE
C     PROBLEMS WITH KNOWN SOLUTIONS, GIVEN IN 'BDYFCN'; BUT IT CAN BE
C     USED WITH NO SIGNIFICANT CHANGE FOR THE CASE WHERE THE SOLUTION
C     IS UNKNOWN.
C
C     THE PROGRAM WILL QUIZ THE USER FOR NEEDED INFORMATION ABOUT THE
C     PROBLEM.  THE INFORMATION ASKED FOR IS AS FOLLOWS:
C
C     NUMSUR,A,B,C,ALPHA:  THESE VARIABLES DEFINE THE SURFACE S.  THEY
C          WERE USED EARLIER IN CALCULATING THE GALERKIN COEFFICIENTS
C          THAT ARE NEEDED BY 'LAPLAC'.  AS BEFORE, THE SURFACE IS
C          DEFINED BY SUBROUTINE 'SURFAC'.
C     FNAME:   THIS IS THE NAME OF THE FILE CONTAINING THE GALERKIN
C          COEFFICIENTS, CALCULATED PREVIOUSLY WITH SUBROUTINE 'GNRT'.
C          THIS FILE MUST BE THE FORMATFREE OUTPUT FILE FROM 'GNRT'.
C     FILOUT:   THIS IS THE NUMBER OF THE OUTPUT FILE, AND IT DEPENDS
C          ON WHETHER THE OUTPUT IS TO THE TERMINAL (FILOUT=TMLOUT)
C          OR TO A SEPARATE FILE (FILOUT=FILALT).  IN THE LATTER CASE,
C          THE USER WILL BE ASKED TO SUPPLY A NAME FOR THE OUTPUT FILE.
C          THE ACTUAL UNIT NUMBERS USED ARE GIVEN IN THE DATA
C          STATEMENT BELOW.
C     NUMPTS:   THE NUMBER OF POINTS IN THE REGION D AT WHICH THE
C          SOLUTION U (GIVEN IN USOLN) OF LAPLACE'S EQUATION IS BEING
C          SOUGHT.  THERE IS AN UPPER LIMIT ON 'NUMPTS'.  IT IS CALLED
C          'MAXPTS', AND IT IS GIVEN BELOW IN THE PARAMETER STATEMENT.
C     POINTS:   THIS ARRAY CONTAINS THE POINTS IN D AT WHICH THE
C          SOLUTION U IS TO BE FOUND.
C     NUMBF:   THE INDEX OF THE BOUNDARY FUNCTION, GIVEN IN 'BDYFCN'.
C     AUXPAR:   AN EXTRA PARAMETER THAT MAY BE USED IN DEFINING
C          THE BOUNDARY FUNCTION, IF DESIRED.
C     DEGREE:   THE DEGREE OF THE SPHERICAL HARMONIC POLYNOMIAL TO BE
C          USED IN SOLVING THE DOUBLE LAYER INTEGRAL EQUATION, IN
C          SUBROUTINE 'LAPLAC'.  THERE IS AN UPPER LIMIT FOR 'DEGREE'.
C          IT IS CALLED 'MAXDEG', AND IT IS GIVEN BELOW IN THE
C          PARAMETER STATEMENT.
C     NINTEG:   THE INTEGRATION PARAMETER USED IN DEFINING THE
C          GAUSSIAN QUADRATURE OF THE DOUBLE LAYER POTENTIAL, IN
C          SUBROUTINE 'LAPLAC'.  THERE IS AN UPPER LIMIT FOR 'NINTEG'.
C          IT IS CALLED 'MAXNIN', AND IT IS GIVEN BELOW IN THE
C          PARAMETER STATEMENT.  THERE ARE TWO GAUSSIAN QUADRATURE
C          SUBROUTINES, 'ZEROLG' AND 'GAUSS2' BEING USED HERE.  THE FIRST
C          COMPUTES THE GAUSS-LEGENDRE NODES AND WEIGHTS OF ORDER .LE. 20,
C          AND THE SECOND DOES SO FOR ORDERS .GE. 32, IN POWERS OF 2.
C          THUS IF 'NINTEG' IS CHOSEN LARGER THAN 20, IT SHOULD BE A
C          POWER OF 2.
C
C     THE ARRAY 'WORK' GIVES THE WORKING STORAGE NEEDED BY 'LAPLAC'.
C     THE DIMENSION OF 'WORK' IS GIVEN IN THE PARAMETER STATEMENT
C     BELOW, BASED ON DEFAULT UPPER LIMITS FOR VARIOUS INPUT VARIABLES.
C     THERE IS ALSO A CHECK AFTER THE INPUT OF ALL VARIABLES THAT
C     THE DIMENSION OF 'WORK' IS ADEQUATE.
C
C     MACHINE DEPENDENT CONSTANTS ARE SET THROUGH THE FUNCTION
C     'D1MACH', WHICH IS ONE OF THE SUBPROGRAMS GIVEN IN THE
C     'LAPLAC' PACKAGE.  THE ROUTINE 'D1MACH' SHOULD BE RESET
C     APPROPRIATELY FOR THE USER'S COMPUTER.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER FILEIN, FILALT, FILOUT, DEGREE, TML_IN, TMLOUT
      CHARACTER FNAME*50, T, ANAME*50
      PARAMETER(MAXDEG=9, MAXNIN=64 ,MAXPTS=20)
      PARAMETER(IWORK=(MAXDEG+1)**4+(MAXDEG+1)**2+2*MAXDEG
     *     +MAXNIN*(4+(MAXDEG*(MAXDEG+3))/2))
      DIMENSION USOLN(MAXPTS), POINTS(3,MAXPTS), WORK(IWORK)
      COMMON/BDYPAR/NUMBF, AUXPAR
      COMMON/BLKWRK/WORK
      COMMON/SURPAR/A, B, C, ALPHA, NUMSUR
      EXTERNAL SURFAC, BDYFCN
C
C     ******* INPUT/OUTPUT **********************************
C     THESE ARE THE I/O FILE NUMBERS, AND THEY MAY NEED
C     TO BE CHANGED FOR USE ON YOUR COMPUTER.
      DATA FILEIN/8/, TML_IN/5/, TMLOUT/6/, FILALT/9/
C     *******************************************************
C
C     INPUT THE SURFACE PARAMETERS.
      PRINT *, ' GIVE THE SURFACE PARAMETERS (NUMSUR,A,B,C,ALPHA).'
      READ *, NUMSUR,A,B,C,ALPHA
C
C     GET INFORMATION ABOUT INPUT AND OUTPUT FILES.
      PRINT *, ' GIVE THE NAME FOR YOUR FILE OF GALERKIN COEFFICIENTS.'
      READ(*,'(A)') FNAME
      L=LEN(FNAME)
      OPEN(FILEIN,FILE=FNAME(1:L),FORM='UNFORMATTED')
      PRINT *, ' DO YOU WANT THE ANSWERS WRITTEN TO THE TERMINAL (Y/N)?'
      READ(*,'(A1)') T
      IF((T .EQ. 'Y') .OR. (T .EQ. 'y')) THEN
          FILOUT=TMLOUT
      ELSE
          PRINT *,' GIVE THE NAME YOU WISH TO USE FOR YOUR OUTPUT FILE.'
          READ(*,'(A)') ANAME
          FILOUT=FILALT
          L=LEN(ANAME)
          OPEN(FILALT,FILE=ANAME(1:L))
      END IF
C
C     PRINT THE SURFACE AND GALERKIN COEFFICIENT PARAMETERS.
      REWIND FILEIN
      READ(FILEIN) NDEG,NINNER,NOUTER
      WRITE(FILOUT,1000) NUMSUR,A,B,C,ALPHA,NINNER,NOUTER
1000  FORMAT(//,' SURFACE PARAMETERS: SURFACE INDEX=',I1,/,' (A,B,C,)=',
     *     1PD17.10,2D20.10,/,' AUXILIARY SURFACE PARAMETER ALPHA = ',
     *     D17.10,/,' GALERKIN COEFFICIENT PARAMETERS: NINNER,NOUTER=',
     *     I3,I6,//)
C
C     READ THE POINTS AT WHICH THE SOLUTION IS TO BE EVALUATED.
      PRINT *, ' AT HOW MANY POINTS INSIDE THE SURFACE IS THE SOLUTION'
      PRINT *, ' TO BE EVALUATED?'
      READ *, NUMPTS
      PRINT *, ' GIVE THE POINTS (X,Y,Z).'
      READ *, ((POINTS(I,J),I=1,3),J=1,NUMPTS)
C
C     READ THE PARTICULAR PROBLEM PARAMETERS.  THE PROGRAM WILL LOOP
C     BACK TO THIS POINT AS OFTEN AS DESIRED.
10    PRINT *, '  GIVE THE INDEX OF THE BOUNDARY FUNCTION, NUMBF, AND'
      PRINT *, ' THE BOUNDARY PARAMETER AUXPAR.  TO STOP THE PROGRAM,'
      PRINT *, ' LET NUMBF=0.'
      READ *, NUMBF,AUXPAR
      IF(NUMBF .EQ. 0) THEN
          CLOSE(FILEIN)
          IF(FILOUT .EQ. FILALT) THEN
              CLOSE(FILALT)
          END IF
          STOP
      END IF
      PRINT *, ' GIVE THE DEGREE OF THE APPROXIMATION AND THE'
      PRINT *, ' INTEGRATION PARAMETER.  ALSO, THE'
      PRINT *, ' INTEGRATION PARAMETER SHOULD BE A POWER OF TWO'
      PRINT *, ' AND LESS THAT OR EQUAL TO 256.'
      READ *, DEGREE,NINTEG
C
C     CHECK ON THE NEEDED DIMENSION FOR 'WORK'.
      ICHECK=(DEGREE+1)**4+(DEGREE+1)**2+2*DEGREE
     *     +NINTEG*(4+(DEGREE*(DEGREE+3))/2)
      IF(ICHECK .GT. IWORK) THEN
C         THERE IS INSUFFICIENT SPACE IN THE ARRAY 'WORK'.
          PRINT 1001, ICHECK
1001      FORMAT(' THE DIMENSION OF THE ARRAY WORK IS NOT SUFFICIENT.',
     *         /,' THE DIMENSION SHOULD BE AT LEAST ',I5)
          STOP
      END IF
C
      CALL LAPLAC(BDYFCN,SURFAC,DEGREE,NINTEG,FILEIN,NUMPTS,POINTS,
     *     USOLN,0,WORK)
C
C     PRINT RESULTS.
      WRITE(FILOUT,1002) NUMBF,AUXPAR,DEGREE,NINTEG
1002  FORMAT(///,' BOUNDARY FUNCTION INDEX=',I2,/,' AUXILIARY BOUNDARY',
     *     ' FUNCTION PARAMETER=',1PD17.10,/,' DEGREE OF ',
     *     'APPROXIMATION=',I2,/,' INTEGRATION PARAMETER=',I3,/)
      WRITE(FILOUT,1003)
      IB=(DEGREE+1)**4+2*DEGREE+NINTEG*(4+(DEGREE*(DEGREE+3))/2)
1003  FORMAT(' LAPLACE EXPANSION COEFFICIENTS FOR DENSITY FUNCTION',
     *     //,7X,'INDICES',11X,'COEFFICIENT',//,3X,'I',4X,'N',4X,
     *     'M',3X,'L')
      DO 15 I=1,(DEGREE+1)**2
        CALL INDX(I,DEGREE,N,M,L)
15      WRITE(FILOUT,1004) I,N,M,L,WORK(IB+I)
1004  FORMAT(I4,I5,I5,I4,1PD21.10)
      WRITE(FILOUT,1005)
1005  FORMAT(//,' POTENTIAL SOLUTION U(X,Y,Z) AT DESIRED POINTS',//,
     *     2X,'I',7X,'X',9X,'Y',9X,'Z',12X,'U',18X,'ERROR')
      DO 20 I=1,NUMPTS
        TRUEU=BDYFCN(POINTS(1,I))
        ERROR=TRUEU-USOLN(I)
20      WRITE(FILOUT,1006) I,(POINTS(J,I),J=1,3),USOLN(I),ERROR
1006  FORMAT(I3,3F10.4,1PD21.10,D14.3)
      GO TO 10
      END

      FUNCTION BDYFCN(Q)
C     ---------------
C
C     THIS DEFINES THE BOUNDARY VALUE FUNCTION OF THE UNKNOWN
C     SOLUTION OF LAPLACE'S EQUATION.  IN THIS PARTICULAR INSTANCE,
C     IT ALSO GIVES THE TRUE SOLUTION OF LAPLACE'S EQUATION, FOR
C     DEMONSTRATION PURPOSES.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/BDYPAR/NUMBF,A
      DIMENSION Q(3)
      PARAMETER(ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,FOUR=4.0D0,
     *     FIVE=5.0D0)
C
      X=Q(1)
      Y=Q(2)
      Z=Q(3)
      GO TO (10,20,30,40,50),NUMBF
10    BDYFCN=ONE
      RETURN
20    BDYFCN=X
      RETURN
30    BDYFCN=X*X+Y*Y-TWO*Z*Z
      RETURN
40    BDYFCN=EXP(A*X)*COS(A*Y)+EXP(A*Z)*SIN(A*X)
      RETURN
50    BDYFCN=ONE/SQRT((X-FIVE)**2+(Y-FOUR)**2+(Z-THREE)**2)
      RETURN
      END

      SUBROUTINE SURFAC(Q,QCROSS,ST,CT,SP,CP,IFLAG)
C     -----------------
C
C     THIS CALCULATES A POINT Q ON THE SURFACE S, BASED ON S BEING THE
C     CONTINUOUS ONE-TO-ONE IMAGE OF THE UNIT SPHERE U.  LET Q=(X,Y,Z)
C     DENOTE THE DESIRED POINT ON S.  THEN IT IS TO CORRESPOND TO
C     THE POINT
C          UQ = (UX,UY,UZ)
C     ON U, WITH
C          UX = SIN(THETA)*COS(PHI) = ST*CP
C          UY = SIN(THETA)*SIN(PHI) = ST*SP
C          UZ = COS(THETA)          = CT
C
C     INPUT:
C     ST,CT,SP,CP   THESE ARE THE GIVEN TRIGONOMETRIC FUNCTION VALUES
C           FOR THE POINT UQ.
C     IFLAG = 0   COMPUTE Q.
C           = 1   COMPUTE Q AND THE CROSS-PRODUCT
C                 QCROSS = DPHI(Q) X DTHETA(Q)/SIN(THETA)
C           THE NOTATION DPHI(Q) MEANS THE DERIVATIVE OF Q REGARDED AS
C           AS A FUNCTION OF PHI, AND SIMILARLY FOR DTHETA(Q).  THE QUANTITY
C           QCROSS IS USED IN COMPUTING THE NORMAL TO S AT Q; AND ITS
C           MAGNITUDE IS THE JACOBIAN IN THE CHANGE OF THE SURFACE
C           DIFFERENTIAL FOR CHANGING AN INTEGRAL OVER S TO AN INTEGRAL
C           OVER U.  THE DIVISION BY SIN(THETA) REMOVES A SINGULARITY
C           IN THE MAPPING DUE TO THE USE OF SPHERICAL COORDINATES
C           IN REPRESENTING POINTS OF U.
C     OUTPUT:
C     Q       THE POINT ON S CORRESPONDING TO UQ ON THE SPHERE U.
C     QCROSS  THE CROSS-PRODUCT REFERRED TO ABOVE.
C
C     THIS PARTICULAR SUBROUTINE ASSUMES THAT THE SURFACE S IS GIVEN BY
C            Q = (X,Y,Z)
C              = R*(A*UX,B*UY,C*UZ)
C     WITH R = R(UQ) A POSITIVE SMOOTH CONTINUOUS FUNCTION ON U.  THE
C     CONSTANTS A,B,C ARE ASSUMED POSITIVE; AND FOR R=1, THE SURFACE
C     WILL BE AN ELLIPSOID WITH SEMI-AXES A,B,C.  THIS GENERAL FORM
C     IS EQUIVALENT TO ASSUMING THAT THE REGION ENCLOSED BY S IS
C     STARLIKE WITH RESPECT TO THE ORIGIN.  IT IS ALSO EQUIVALENT TO
C     THE STANDARD PARAMETRIC WAY OF DEFINING A SURFACE S IN TERMS
C     OF VARIABLES PHI AND THETA VARYING OVER  (O,2*PI) AND (0,PI),
C     RESPECTIVELY.
C
C     IT IS ASSUMED THAT R CAN BE EXTENDED AS A SMOOTH FUNCTION TO A
C     NEIGHBORHOOD OF U, SO THAT THE FIRST DERIVATIVES OF R(UX,UY,UZ)
C     CAN BE COMPUTED AT ALL POINTS OF U.  THE USER MUST SUPPLY R AND
C     ITS DERIVATIVES WITH RESPECT TO UX,UY, AND UZ, DENOTED BY DXR,
C     DYR, AND DZR, RESPECTIVELY.  EXAMPLES ARE GIVEN BELOW.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION Q(3),QCROSS(3)
      COMMON/SURPAR/A,B,C,ALPHA,NUMSUR
      COMMON/MACHCN/U100,U100SQ
C     U100 IS 100 TIMES THE MACHINE UNIT ROUND.
C
C     COMPUTE UQ.
      UX = ST*CP
      UY = ST*SP
      UZ = CT
C
C*******************************************************************
C     THIS BEGINS THE SECTION THAT THE USER MUST SPECIFY.
C
      GO TO (10,20,30), NUMSUR
C
C     SURFACE #1.
C     S IS AN ELLIPSOID.
10    R = 1.0D0
      IF(IFLAG .EQ. 0) GO TO 100
      DXR = 0.0D0
      DYR = 0.0D0
      DZR = 0.0D0
      GO TO 100
C
C     SURFACE #2.
C     S IS A PEANUT-SHAPED REGION, BASED ON THE OVALS OF CASSINI.
C     THE PARAMETER ALPHA IS A NON-ZERO POSITIVE CONSTANT.  AS IT
C     APPROACHES ZERO, THE SURFACE S BECOMES MORE CONSTRICTED IN
C     THE XY-PLANE.
20    C2T = 2.0D0*UZ*UZ - 1.0D0
      TEMP = SQRT(ALPHA + C2T*C2T)
      R = SQRT(C2T + TEMP)
      IF(IFLAG .EQ. 0) GO TO 100
      DXR = 0.0D0
      DYR = 0.0D0
      DZR=2.0D0*UZ*(1.0D0+C2T/TEMP)/R
      GO TO 100
C
C     SURFACE #3.
C     S IS A HEART-SHAPED REGION, INCREASINGLY SO AS THE POSITIVE
C     CONSTANT ALPHA INCREASES.
30    TEMP=1.0D0+ALPHA*(UZ+1.0D0)**2
      R=2.0D0-1.0D0/TEMP
      IF(IFLAG .EQ. 0) GO TO 100
      DXR=0.0D0
      DYR=0.0D0
      DZR=2.0D0*ALPHA*(UZ+1.0D0)/TEMP**2
      GO TO 100
C
C     THIS ENDS THE SECTION THAT THE USER MUST SUPPLY.
C*******************************************************************
C
C     COMPUTE THE POINT Q.
100   Q(1) = A*R*UX
      Q(2) = B*R*UY
      Q(3) = C*R*UZ
      IF(IFLAG .EQ. 0) RETURN
C
C     COMPUTE CROSS-PRODUCT.
C
      IF(ST .LE. U100) GO TO 200
C     BEGIN BY COMPUTING DERIVATIVES OF Q WITH RESPECT TO PHI AND THETA.
      DPHIR = -UY*DXR + UX*DYR
      DTHER = CT*(CP*DXR+SP*DYR) - ST*DZR
      DPHIQX = A*ST*(-SP*R + CP*DPHIR)
      DTHEQX = A*CP*(CT*R + ST*DTHER)
      DPHIQY = B*ST*(CP*R + SP*DPHIR)
      DTHEQY = B*SP*(CT*R + ST*DTHER)
      DPHIQZ = C*CT*DPHIR
      DTHEQZ = C*(-ST*R + CT*DTHER)
      QCROSS(1) = (DPHIQY*DTHEQZ-DTHEQY*DPHIQZ)/ST
      QCROSS(2) = (DPHIQZ*DTHEQX-DTHEQZ*DPHIQX)/ST
      QCROSS(3) = (DPHIQX*DTHEQY-DTHEQX*DPHIQY)/ST
      RETURN
C
C     COMPUTE CROSS-PRODUCT FOR Q AT A POLE.
200   QCROSS(1) = B*C*R*DXR
      QCROSS(2) = A*C*R*DYR
      QCROSS(3) = -A*B*UZ*R*R
      RETURN
      END

      SUBROUTINE LAPLAC(BDYFCN,SURFAC,DEGREE,NINTEG,FILEIN,NUMPTS,
     *     POINTS,USOLN,OPREPT,WORK)
C     -----------------
C
C     THIS CALCULATES THE SOLUTION OF THE INTERIOR DIRICHLET
C     PROBLEM FOR LAPLACE'S EQUATION ON A REGION D IN THREE
C     DIMENSIONS.  THE REGION D IS ASSUMED TO BE SIMPLY-CONNECTED
C     WITH A SMOOTH BOUNDARY S.  THE METHOD OF SOLUTION IS TO
C     REPRESENT THE UNKNOWN HARMONIC FUNCTION U(X,Y,Z) AS A DOUBLE
C     LAYER POTENTIAL.  THIS IS FOUND BY NUMERICALLY SOLVING FOR THE
C     DOUBLE LAYER DENSITY FUNCTION RHO(X,Y,Z) ON S BY MEANS OF
C     GALERKIN'S METHOD WITH SPHERICAL HARMONICS AS BASIS FUNCTIONS.
C     THE METHOD IS DESCRIBED IN THE ARTICLE
C          K. ATKINSON, 'THE NUMERICAL SOLUTION OF LAPLACE'S EQUATION
C          IN THREE DIMENSIONS', SIAM J. NUM. ANAL. 19(1982),263-274.
C     THE DOUBLE LAYER POTENTIAL IS THEN EVALUATED BY USING GAUSSIAN
C     QUADRATURE OVER THE SPHERE, AS DESCRIBED IN
C          K. ATKINSON, 'NUMERICAL INTEGRATION ON THE SPHERE', J.
C          AUSTRALIAN MATH. SOC., (SERIES B) 23(1982),332-347.
C
C     INPUT:
C     BDYFCN - THE BOUNDARY CONDITION FOR THE SOLUTION U.  IT MUST BE
C          GIVEN AS A FUNCTION BDYFCN(Q), WITH Q A VECTOR OF LENGTH 3.
C     SURFAC - THIS IS THE NAME OF THE SUBROUTINE THAT DEFINES THE
C          SURFACE S.  SEE THE EXAMPLE SUBROUTINE 'SURFAC' FOR DETAILS
C          OF THE DEFINITION AND THE CALLING SEQUENCE.
C     DEGREE - THE UPPER LIMIT ON THE DEGREE OF THE SPHERICAL HARMONICS
C          BEING USED IN THE APPROXIMATION OF RHO.
C     NINTEG - THE INTEGRATION PARAMETER BEING USED TO DETERMINE THE
C          GAUSSIAN QUADRATURE FORMULA FOR SPHERICAL INTEGRATION.  AT
C           PRESENT, IF 'NINTEG' IS GREATER THAN 20, THEN IT SHOULD BE A
C          POWER OF 2.  THIS IS DUE TO THE SUBROUTINE 'GAUSS2' BEING
C          USED TO PRODUCE THE GAUSS-LEGENDRE NODES AND WEIGHTS IN SUCH
C          A CASE.  THE TOTAL NUMBER OF INTEGRATION NODES ON THE SURFACE
C          S WILL BE 2*NINTEG**2.
C     FILEIN - THE NUMBER OF THE INPUT FILE CONTAINING THE GALERKIN
C          COEFFICIENTS, AS GENERATED BY THE SUBROUTINE 'GNRT'.  THIS
C          FILE MUST BE THE FORMATFREE FILE FROM 'GNRT'.
C     NUMPTS - THE NUMBER OF POINTS IN D AT WHICH THE SOLUTION U IS TO
C          BE EVALUATED.
C     POINTS - THE ARRAY CONTAINING THE POINTS FOR THE EVALUATION OF
C          U(X,Y,Z).  THE COORDINATES (X,Y,Z) OF POINT #J ARE STORED
C          IN POINTS(1,J), POINTS(2,J), POINTS(3,J).
C     OPREPT - OPREPT=0 MEANS THAT THIS IS THE FIRST CALL ON 'LAPLAC'
C          WITH THIS SET OF PROBLEM PARAMETERS.
C          OPREPT=1 MEANS 'LAPLAC' HAS BEEN CALLED PREVIOUSLY, AND
C          IT IS DESIRED TO EVALUATE U AT A NEW SET OF POINTS IN D.
C          ONLY 'NUMPTS' AND 'POINTS' SHOULD HAVE BEEN CHANGED SINCE
C          THE PREVIOUS CALL ON 'LAPLAC'.
C     WORK - A TEMPORARY ARRAY OF WORKING STORAGE.  ITS DIMENSION
C          SHOULD BE AT LEAST
C               (DEGREE+1)**4 + (DEGREE+1)**2 + 2*DEGREE
C                  +NINTEG*(4+(DEGREE*(DEGREE+3))/2)
C
C     OUTPUT:
C     USOLN - A ONE DIMENSIONAL ARRAY OF LENGTH AT LEAST 'NUMPTS'.
C          THE ELEMENT USOLN(I) CONTAINS THE APPROXIMATION TO U AT
C          POINT #I OF 'POINTS'.
C     GALERKIN COEFFICIENTS OF RHO - IF DESIRED, THESE ARE STORED IN
C          'WORK', BEGINNING AT POSITION
C               (DEGREE+1)**4 + 2*DEGREE
C                    + NINTEG*(4+(DEGREE*(DEGREE+3))/2)
C          THE SPHERICAL HARMONICS WITH WHICH THESE COEFFICIENTS ARE
C          ASSOCIATED CAN BE DETERMINED BY USING SUBROUTINE 'INDX' TO
C          CONVERT THE INDEX I OF THE COEFFICIENT #I.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER DEGREE,FILEIN,OPREPT,DEGSQ,DEG2
      DIMENSION POINTS(3,*),USOLN(*),WORK(*)
      COMMON/MACHCN/U100,U100SQ
C     THE VARIABLE U100 IS 100 TIMES THE MACHINE UNIT ROUND.
      SAVE
      EXTERNAL BDYFCN,SURFAC
C
C     SET MACHINE DEPENDENT CONSTANT.
      U100=100*D1MACH(4)
      U100SQ=U100**2
C
      IF(OPREPT .EQ. 1) GO TO 10
      DEGSQ=(DEGREE+1)**2
      DEG2=2*DEGREE+1
      LGDORD=1+(DEGREE*(DEGREE+3))/2
C
C     SUBDIVIDE 'WORK' INTO SMALLER ONE AND TWO DIMENSIONAL ARRAYS FOR
C     USE IN SUBROUTINE 'LPLC2'.
      I1=1
      I2=I1+DEGSQ**2
      I3=I2+NINTEG
      I4=I3+NINTEG
      I5=I4+NINTEG
      I6=I5+DEGREE
      I7=I6+DEGREE
      I8=I7+LGDORD*NINTEG
C
10    CALL LPLC2(BDYFCN,SURFAC,DEGREE,NINTEG,FILEIN,NUMPTS,POINTS,USOLN,
     *     OPREPT,WORK(I1),WORK(I2),WORK(I3),WORK(I4),WORK(I5),
     *     WORK(I6),WORK(I7),WORK(I8),DEGSQ,DEG2,LGDORD)
      RETURN
      END

      SUBROUTINE LPLC2(BDYFCN,SURFAC,DEGREE,NINTEG,FILEIN,NUMPTS,
     *     POINTS,USOLN,OPREPT,KERMAT,CSTH,SNTH,W,CSPH,SNPH,VALUES,
     *     RHS,DEGSQ,DEG2,LGDORD)
C     ----------------
C
C     THE ACTUAL SOLUTION OF LAPLACE'S EQUATION TAKES PLACE IN THIS ROUTINE.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION KERMAT,KERNEL
      INTEGER DEGREE,DEGSQ,DEG2,OPREPT,FILEIN
      DIMENSION KERMAT(DEGSQ,DEGSQ),CSTH(NINTEG),SNTH(NINTEG),
     *     W(NINTEG),CSPH(DEGREE),SNPH(DEGREE),VALUES(LGDORD,NINTEG),
     *     RHS(DEGSQ),POINTS(3,*),USOLN(*),Q(3),QCROSS(3)
      PARAMETER(ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0,FOUR=4.0D0)
      SAVE
C
      IF(OPREPT .EQ. 1) GO TO 200
C
C     MOVE POINTER OF INPUT FILE TO CORRECT SET OF TABLES.
      REWIND FILEIN
      DO 10 K=1,DEGREE-1
        KT=(K+1)**4
        READ(FILEIN)NDEG,NINNER,NOUTER
10      READ(FILEIN)(DUMMY,I=1,KT)
C
C     READ MATRIX OF GALERKIN COEFFICIENTS.
      READ(FILEIN)NDEG,NINNER,NOUTER
      READ(FILEIN)((KERMAT(I,J),J=1,DEGSQ),I=1,DEGSQ)
C
C     INITIALIZE.
      PI=FOUR*ATAN(ONE)
      TWOPI=TWO*PI
      HPHI=PI/NINTEG
      DO 20 I=1,DEGSQ
20      KERMAT(I,I)=KERMAT(I,I)+TWOPI
C
C     INITIALIZE INTEGRATION NODES AND THE ASSOCIATED TRIG FUNCTIONS.
      IF(NINTEG .LE. 20) THEN
          CALL ZEROLG(NINTEG,CSTH,W,-ONE,ONE)
      ELSE
          CALL GAUSS2(-ONE,ONE,NINTEG,CSTH,W)
      END IF
      DO 30 I=1,NINTEG
        THETA=ACOS(CSTH(I))
        SNTH(I)=SIN(THETA)
30      CALL LEGEND(CSTH(I),DEGREE,VALUES(1,I))
C
C     BEGIN SECTION TO CREATE RHS.
C
C     CREATE RIGHT-HAND SIDES OF GALERKIN SYSTEM.
      DO 70 I=1,DEGSQ
70      RHS(I)=ZERO
      DO 100 J=1,2*NINTEG
        PHI=J*HPHI
        CP=COS(PHI)
        SP=SIN(PHI)
        SNPH(1)=SP
        CSPH(1)=CP
        DO 80 M=2,DEGREE
          CSPH(M)=CP*CSPH(M-1)-SP*SNPH(M-1)
80        SNPH(M)=SP*CSPH(M-1)+CP*SNPH(M-1)
        DO 100 I=1,NINTEG
          CALL SURFAC(Q,DUMMY,SNTH(I),CSTH(I),SP,CP,0)
          BF=BDYFCN(Q)
          DO 90 N=1,DEGREE+1
90          RHS(N)=RHS(N)+W(I)*VALUES(N,I)*BF
          IB=DEGREE+1
          IBNM=DEGREE
          DO 100 M=1,DEGREE
          DO 100 N=M,DEGREE
            IB=IB+1
            IBNM=IBNM+2
            RHS(IBNM)=RHS(IBNM)+W(I)*BF*CSPH(M)*VALUES(IB,I)
100         RHS(IBNM+1)=RHS(IBNM+1)+W(I)*BF*SNPH(M)*VALUES(IB,I)
      DO 110 J=1,DEGSQ
110     RHS(J)=HPHI*RHS(J)
C
C     END SECTION ON CREATING RHS.
C
C     *********************************************************
C
C     SOLVE LINEAR SYSTEM FOR LAPLACE COEFFICIENTS
      CALL LINSYS(KERMAT,RHS,DEGSQ,DEGSQ,IER)
C     *********************************************************
C
C     INITIALIZE ALL VALUES OF U TO ZERO.
200   DO 210 L=1,NUMPTS
210     USOLN(L)=ZERO
C
C     BEGIN LOOP TO CREATE VALUES OF USOLN AT POINTS (X,Y,Z)
C     GIVEN IN THE ARRAY 'POINTS'.
C
      DO 250 J=1,2*NINTEG
        PHI=J*HPHI
        SP=SIN(PHI)
        CP=COS(PHI)
        SNPH(1)=SP
        CSPH(1)=CP
        DO 220 M=2,DEGREE
          SNPH(M)=SP*CSPH(M-1)+CP*SNPH(M-1)
220       CSPH(M)=CP*CSPH(M-1)-SP*SNPH(M-1)
        DO 250 I=1,NINTEG
C         CALCULATE RHO AT (THETA(I),PHI).
          RHO=ZERO
          DO 230 N=1,DEGREE+1
230         RHO=RHO+RHS(N)*VALUES(N,I)
          IB=DEGREE+1
          IBNM=DEGREE
          DO 240 M=1,DEGREE
          DO 240 N=M,DEGREE
            IB=IB+1
            IBNM=IBNM+2
240         RHO=RHO+VALUES(IB,I)*(CSPH(M)*RHS(IBNM)
     *         +SNPH(M)*RHS(IBNM+1))
C         CALCULATE VALUE OF U AT EACH (X,Y,Z) IN POINTS.
          DO 250 L=1,NUMPTS
            CALL SURFAC(Q,QCROSS,SNTH(I),CSTH(I),SP,CP,1)
            FM=KERNEL(POINTS(1,L),Q,QCROSS)
            USOLN(L)=USOLN(L)+W(I)*RHO*FM
250         CONTINUE
      DO 260 L=1,NUMPTS
260     USOLN(L)=HPHI*USOLN(L)
      RETURN
      END

      FUNCTION KERNEL(P,Q,QCROSS)
C     ---------------
C
C     EVALUATE THE NORMAL DERIVATIVE OF 1/ABS(P-Q) WITH RESPECT TO THE
C     POINT Q ON THE SURFACE S.  ALSO INCLUDE THE EFFECT OF THE CHANGE
C     OF SURFACE DIFFERENTIAL IN TRANSFORMING TO AN INTEGRAL OVER THE
C     UNIT SPHERE U.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION KERNEL
      DIMENSION P(3),Q(3),QCROSS(3)
      COMMON/MACHCN/U100,U100SQ
C     U100 IS 100 TIMES THE MACHINE UNIT ROUND.
      PARAMETER(ZERO=0.0D0)
C
      RS=(P(1)-Q(1))**2+(P(2)-Q(2))**2+(P(3)-Q(3))**2
      IF(RS .LE. U100SQ) THEN
          KERNEL=ZERO
      ELSE
          DENOM=RS*SQRT(RS)
          TEMP=ZERO
          DO 10 J=1,3
10          TEMP=TEMP+(P(J)-Q(J))*QCROSS(J)
          KERNEL=TEMP/DENOM
      END IF
      RETURN
      END

      SUBROUTINE LEGEND(X,NMAX,VALUES)
C     -----------------
C
C        THIS ROUTINE WILL EVALUATE THE NORMALIZED ASSOCIATED LEGENDRE
C     FUNCTIONS OF THE FIRST KIND AT X,
C                                 J
C                         Y(N,J)=P (X)
C                                 N
C
C     THE DEFINITION IS TAKEN FROM CHAPTER 8 OF ABRAMOWITZ AND STEGUN,
C     HANDBOOK OF MATHEMATICAL FUNCTIONS. THESE WILL BE EVALUATED AND
C     STORED IN THE VECTOR 'VALUES' FOR 0 .LE. J .LE. N .LE. NMAX, WITH
C     NMAX .GE. 1.  THE VALUE Y(N,J) IS STORED IN JB+N-J, WITH
C                 JB=(2+(2*NMAX+3)*J-J**2)/2
C     JB IS THE INDEX FOR Y(J,J).  THE NORMALIZATION IS DONE SO THAT
C     THE SUBSEQUENT CONSTRUCTION OF SPHERICAL HARMONICS WILL RESULT
C     IN FUNCTIONS OF NORM 1 OVER THE UNIT SPHERE.
C        THE EVALUATION IS DONE USING THE TRIPLE RECURSION RELATION
C     (8.5.3) ON PAGE 334 OF THE ABOVE REFERENCE. THIS IS PRECEDED BY
C     A SPECIAL LOOP TO SET Y(N,N) AND Y(N,N-1).  FOLLOWING THAT, THE
C     FUNCTIONS ARE NORMALIZED APPROPRIATELY.
C
C     NOTE: THERE IS AN UPPER LIMIT ON NMAX, SPECIFIED BY THE
C     VARIABLE 'MAXDEG' IN THE PARAMETER STATEMENT BELOW.  TO
C     INCREASE THE UPPER LIMIT, MERELY INCREASE 'MAXDEG'.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(MAXDEG=20,LFORD=2*MAXDEG)
      DIMENSION VALUES(*),FACT(0:LFORD)
      PARAMETER(ONE=1.0D0,TWO=2.0D0,FOUR=4.0D0)
C
C     SET INITIAL VALUES Y(N,N) AND Y(N,N-1), 0 .LE. N .LE. NMAX.
      RT=SQRT(ONE-X*X)
      VALUES(1)=ONE
      VALUES(2)=X
      IF(NMAX .EQ. 1) THEN
          VALUES(3)=-RT
          GO TO 30
      END IF
C
      JB=1
      DO 10 J=1,NMAX-1
        JBN=JB+NMAX-J+2
        VALUES(JBN)=-(2*J-1)*RT*VALUES(JB)
        VALUES(JBN+1)=-(2*J+1)*RT*VALUES(JB+1)
10      JB=JBN
      VALUES(JB+2)=-(2*NMAX-1)*RT*VALUES(JB)
C
C     PRODUCE THE REMAINDER OF THE TABLE VALUES Y(N,J).
      JB=1
      DO 20 J=0,NMAX-2
        JB=JB+2
        DO 20 N=J+2,NMAX
          VALUES(JB)=((2*N-1)*X*VALUES(JB-1)-(N+J-1)*VALUES(JB-2))/(N-J)
20        JB=JB+1
C
C     NORMALIZE THE LEGENDRE FUNCTIONS.
30    FACT(0)=ONE
      DO 40 J=1,2*NMAX
40      FACT(J)=J*FACT(J-1)
      IB=0
      PI=4*ATAN(ONE)
      TWOPI=TWO*PI
      FOURPI=FOUR*PI
      DO 50 N=0,NMAX
50      VALUES(N+1)=VALUES(N+1)*SQRT((2*N+1)/FOURPI)
      JB=NMAX+1
      DO 60 M=1,NMAX
      DO 60 N=M,NMAX
        JB=JB+1
        COEFF=SQRT(((2*N+1)/TWOPI)*FACT(N-M)/FACT(N+M))
60      VALUES(JB)=COEFF*VALUES(JB)
      RETURN
      END

      SUBROUTINE INDX(I,NDEG,N,M,L)
C     ---------------
C
C     THIS PROGRAM TAKES THE INDEX I OF A SPHERICAL HARMONIC,
C              1 .LE. I .LE. (NDEG+1)**2,
C     AND IT FINDS THE INDICES (N,M,L) TO WHICH I CORRESPONDS.
C     SEE FUNCTION 'INDX' FOR THE DEFINITION OF (N,M,L).
C
      KOL(MM)=NDEG*(MM-1)-(MM*(MM-3))/2-1
C
      IF(I .LE. NDEG+1) THEN
          N=I-1
          M=0
          L=0
          RETURN
      END IF
C
      J=I-NDEG-2
      IF(J .EQ. 2*(J/2)) THEN
          L=0
      ELSE
          L=1
      END IF
      J=(J-L)/2
      DO 40 M=1,NDEG
        IF(J .LT. KOL(M+1)) THEN
            N=J-KOL(M)+M
            RETURN
        END IF
40      CONTINUE
      END

      SUBROUTINE ZEROLG(N,ZZ,WW,A,B)
C     -----------------
C
C     PRODUCE THE NODES AND WEIGHTS FOR GAUSS-LEGENDRE QUADRATURE
C     ON (A,B), OF ORDER N. THE NUMBER OF NODES N IS RESTRICTED TO
C     BE .LE. 20.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION Z(109),W(109),ZZ(N),WW(N)
      DATA ONE/1.0D0/,TWO/2.0D0/
      DATA Z(1),Z(2),Z(3),Z(4),Z(5),Z(6),Z(7),Z(8),Z(9),Z(10)/
     *.577350269189626D0,.774596669241483D0,0.0D0,.861136311594053D0,
     *.339981043584856D0,.906179845938664D0,.538469310105683D0,0.0D0,
     *.932469514203152D0,.661209386466265D0/
      DATA Z(11),Z(12),Z(13),Z(14),Z(15),Z(16),Z(17),Z(18),Z(19),Z(20)/
     *.238619186083197D0,.949107912342759D0,.741531185599394D0,
     *.405845151377397D0,0.0D0,.960289856497536D0,.796666477413627D0,
     *.525532409916329D0,.183434642495650D0,.968160239507626D0/
      DATA Z(21),Z(22),Z(23),Z(24),Z(25),Z(26),Z(27),Z(28),Z(29)/
     *.836031107326636D0,.613371432700590D0,.324253423403809D0,
     *0.0D0,.973906528517172D0,.865063366688985D0,.679409568299024D0,
     *.433395394129247D0,.148874338981631D0/
      DATA Z(30),Z(31),Z(32),Z(33),Z(34),Z(35),Z(36),Z(37),Z(38)/
     *.978228658146057D0,.887062599768095D0,.730152005574049D0,
     *.519096129206812D0,.269543155952345D0,0.0D0,
     *.981560634246719D0,.904117256370475D0,.769902674194305D0/
      DATA Z(39),Z(40),Z(41),Z(42),Z(43),Z(44),Z(45),Z(46),Z(47)/
     *.587317954286617D0,.367831498998180D0,.125233408511469D0,
     *.984183054718588D0,.917598399222978D0,.801578090733310D0,
     *.642349339440340D0,.448492751036447D0,.230458315955135D0/
      DATA Z(48),Z(49),Z(50),Z(51),Z(52),Z(53),Z(54),Z(55),Z(56)/
     *0.0D0,.986283808696812D0,.928434883663574D0,
     *.827201315069765D0,.687292904811685D0,.515248636358154D0,
     *.319112368927890D0,.108054948707344D0,.987992518020485D0/
      DATA Z(57),Z(58),Z(59),Z(60),Z(61),Z(62),Z(63),Z(64),Z(65)/
     *.937273392400706D0,.848206583410427D0,.724417731360170D0,
     *.570972172608539D0,.394151347077563D0,.201194093997435D0,
     *0.0D0,.989400934991650D0,.944575023073233D0/
      DATA Z(66),Z(67),Z(68),Z(69),Z(70),Z(71),Z(72),Z(73),Z(74)/
     *.865631202387832D0,.755404408355003D0,.617876244402644D0,
     *.458016777657227D0,.281603550779259D0,.0950125098376374D0,
     *.990575475314417D0,.950675521768768D0,.880239153726986D0/
      DATA Z(75),Z(76),Z(77),Z(78),Z(79),Z(80),Z(81),Z(82),Z(83)/
     *.781514003896801D0,.657671159216691D0,.512690537086477D0,
     *.351231763453876D0,.178484181495848D0,0.0D0,
     *.991565168420931D0,.955823949571398D0,.892602466497556D0/
      DATA Z(84),Z(85),Z(86),Z(87),Z(88),Z(89),Z(90),Z(91),Z(92)/
     *.803704958972523D0,.691687043060353D0,.559770831073948D0,
     *.411751161462843D0,.251886225691506D0,.0847750130417353D0,
     *.992406843843584D0,.960208152134830D0,.903155903614818D0/
      DATA Z(93),Z(94),Z(95),Z(96),Z(97),Z(98),Z(99),Z(100),Z(101)/
     *.822714656537143D0,.720966177335229D0,.600545304661681D0,
     *.464570741375961D0,.316564099963630D0,.160358645640225D0,
     *0.0D0,.993128599185095D0,.963971927277914D0/
      DATA Z(102),Z(103),Z(104),Z(105),Z(106),Z(107),Z(108),Z(109)/
     *.912234428251326D0,.839116971822219D0,.746331906460151D0,
     *.636053680726515D0,.510867001950827D0,.373706088715420D0,
     *.227785851141645D0,.0765265211334973D0/
      DATA W(1),W(2),W(3),W(4),W(5),W(6),W(7),W(8),W(9),W(10)/
     *1.0D0,.555555555555556D0,.888888888888889D0,.347854845137454D0,
     *.652145154862546D0,.236926885056189D0,.478628670499366D0,
     *.568888888888889D0,.171324492379170D0,.360761573048139D0/
      DATA W(11),W(12),W(13),W(14),W(15),W(16),W(17),W(18),W(19),W(20)/
     *.467913934572691D0,.129484966168870D0,.279705391489277D0,
     *.381830050505119D0,.417959183673469D0,.101228536290376D0,
     *.222381034453374D0,.313706645877887D0,.362683783378362D0,
     *.0812743883615744D0/
      DATA W(21),W(22),W(23),W(24),W(25),W(26),W(27),W(28),W(29)/
     *.180648160694857D0,.260610696402935D0,.312347077040003D0,
     *.330239355001260D0,.0666713443086881D0,.149451349150581D0,
     *.219086362515982D0,.269266719309996D0,.295524224714753D0/
      DATA W(30),W(31),W(32),W(33),W(34),W(35),W(36),W(37),W(38)/
     *.0556685671161737D0,.125580369464905D0,.186290210927734D0,
     *.233193764591990D0,.262804544510247D0,.272925086777901D0,
     *.0471753363865118D0,.106939325995318D0,.160078328543346D0/
      DATA W(39),W(40),W(41),W(42),W(43),W(44),W(45),W(46),W(47)/
     *.203167426723066D0,.233492536538355D0,.249147045813403D0,
     *.0404840047653159D0,.0921214998377284D0,.138873510219787D0,
     *.178145980761946D0,.207816047536889D0,.226283180262897D0/
      DATA W(48),W(49),W(50),W(51),W(52),W(53),W(54),W(55),W(56)/
     *.232551553230874D0,.0351194603317519D0,.0801580871597602D0,
     *.121518570687903D0,.157203167158194D0,.185538397477938D0,
     *.205198463721296D0,.215263853463158D0,.0307532419961173D0/
      DATA W(57),W(58),W(59),W(60),W(61),W(62),W(63),W(64),W(65)/
     *.0703660474881081D0,.107159220467172D0,.139570677926154D0,
     *.166269205816994D0,.186161000015562D0,.198431485327112D0,
     *.202578241925561D0,.0271524594117541D0,.0622535239386478D0/
      DATA W(66),W(67),W(68),W(69),W(70),W(71),W(72),W(73),W(74)/
     *.0951585116824928D0,.124628971255534D0,.149595988816577D0,
     *.169156519395003D0,.182603415044924D0,.189450610455068D0,
     *.0241483028685479D0,.0554595293739872D0,.0850361483171792D0/
      DATA W(75),W(76),W(77),W(78),W(79),W(80),W(81),W(82),W(83)/
     *.111883847193404D0,.135136368468525D0,.154045761076810D0,
     *.168004102156450D0,.176562705366993D0,.179446470356207D0,
     *.0216160135264833D0,.0497145488949698D0,.0764257302548891D0/
      DATA W(84),W(85),W(86),W(87),W(88),W(89),W(90),W(91),W(92)/
     *.100942044106287D0,.122555206711478D0,.140642914670651D0,
     *.154684675126265D0,.164276483745833D0,.169142382963144D0,
     *.0194617882297265D0,.0448142267656996D0,.0690445427376412D0/
      DATA W(93),W(94),W(95),W(96),W(97),W(98),W(99),W(100),W(101)/
     *.0914900216224500D0,.111566645547334D0,.128753962539336D0,
     *.142606702173607D0,.152766042065860D0,.158968843393954D0,
     *.161054449848784D0,.0176140071391521D0,.0406014298003869D0/
      DATA W(102),W(103),W(104),W(105),W(106),W(107),W(108),W(109)/
     *.0626720483341091D0,.0832767415767047D0,.101930119817240D0,
     *.118194531961518D0,.131688638449177D0,.142096109318382D0,
     *.149172986472604D0,.152753387130726D0/
C
C     CALCULATE THE WEIGHTS AND NODES.
C
      SCALE=(B-A)/TWO
      IF((N/2)*2 .EQ. N) THEN
          M=N/2
          IBASE=M*M
      ELSE
          IBASE=(N*N-1)/4
          M=(N-1)/2
          ZZ(M+1)=(A+B)/TWO
          WW(M+1)=W(IBASE+M)*SCALE
      END IF
      DO 100 I=1,M
        T=Z(IBASE+I-1)
        ZZ(I)=(A*(ONE+T)+(ONE-T)*B)/TWO
        ZZ(N+1-I)=(A*(ONE-T)+(ONE+T)*B)/TWO
        WW(I)=W(IBASE+I-1)*SCALE
100     WW(N+1-I)=WW(I)
      RETURN
      END

      SUBROUTINE GAUSS2(A,B,N,TT,WW)
C     -----------------
C
C   THE GAUSSIAN QUADRATURE NODES AND WEIGHTS ARE CALCULATED FOR THE FORMULA
C
C             B                            I=N
C         INT   F(T)*DT   APPROX. =      SUM   WW(I)*F(TT(I))
C             A                            I=1
C
C   THE NODES WILL BE STORED IN TT IN INCREASING ORDER,
C         A .LT. TT(1) .LT.  ...  .LT. TT(N) .LT. B.
C   THESE NODES ARE THE ZEROS OF THE LEGENDRE POLYNOMIAL ON (A,B) OF
C   DEGREE N. THE CORRESPONDING WEIGHT FOR TT(I) WILL BE STORED IN WW(I).
C   ***WARNING***   THE INTEGER N MUST BE A POWER OF TWO,
C         2 .LE. N .LE. 256.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION WW(1),TT(1)
      DIMENSION W(255),T(255)
      DATA T(1),T(2),T(3),T(4),T(5),T(6),T(7),T(8),T(9),T(10),T(11),
     *     T(12),T(13),T(14),T(15)/.577350269189626D0,
     *  .861136311594053D0,.339981043584856D0,.960289856497536D0,
     *  .796666477413627D0,.525532409916329D0,.183434642495650D0,
     *  .989400934991650D0,.944575023073233D0,.865631202387832D0,
     *  .755404408355003D0,.617876244402644D0,.458016777657227D0,
     *  .281603550779259D0,.950125098376374D-1/
      DATA T(16),T(17),T(18),T(19),T(20),T(21),T(22),T(23),T(24),T(25),
     *     T(26),T(27),T(28),T(29),T(30),T(31)/.997263861849482D0,
     *  .985611511545268D0,.964762255587506D0,.934906075937740D0,
     *  .896321155766052D0,.849367613732570D0,.794483795967942D0,
     *  .732182118740290D0,.663044266930215D0,.587715757240762D0,
     *  .506899908932229D0,.421351276130635D0,.331868602282128D0,
     *  .239287362252137D0,.144471961582796D0,.483076656877383D-1/
      DATA T(32),T(33),T(34),T(35),T(36),T(37),T(38),T(39),T(40),T(41),
     *     T(42),T(43),T(44),T(45),T(46),T(47)/.999305041735772D0,
     *  .996340116771955D0,.991013371476744D0,.983336253884626D0,
     *  .973326827789911D0,.961008799652054D0,.946411374858403D0,
     *  .929569172131940D0,.910522137078503D0,.889315445995114D0,
     *  .865999398154093D0,.840629296252580D0,.813265315122798D0,
     *  .783972358943341D0,.752819907260532D0,.719881850171611D0/
      DATA T(48),T(49),T(50),T(51),T(52),T(53),T(54),T(55),T(56),T(57),
     *     T(58),T(59),T(60),T(61),T(62),T(63)/.685236313054233D0,
     *  .648965471254657D0,.611155355172393D0,.571895646202634D0,
     *  .531279464019895D0,.489403145707053D0,.446366017253464D0,
     *  .402270157963992D0,.357220158337668D0,.311322871990211D0,
     *  .264687162208767D0,.217423643740007D0,.169644420423993D0,
     *  .121462819296121D0,.729931217877990D-1,.243502926634244D-1/
      DATA T(64),T(65),T(66),T(67),T(68),T(69),T(70),T(71),T(72),T(73),
     *     T(74),T(75),T(76),T(77),T(78),T(79)/.999824887947132D0,
     *  .999077459977376D0,.997733248625514D0,.995792758534981D0,
     *  .993257112900213D0,.990127818491734D0,.986406742724586D0,
     *  .982096108435719D0,.977198491463907D0,.971716818747137D0,
     *  .965654366431965D0,.959014757853700D0,.951801961341264D0,
     *  .944020287830220D0,.935674388277916D0,.926769250878948D0/
      DATA T(80),T(81),T(82),T(83),T(84),T(85),T(86),T(87),T(88),T(89),
     *     T(90),T(91),T(92),T(93),T(94),T(95),T(96),T(97)/
     *  .917310198080961D0,.907302883401757D0,.896753288049158D0,
     *  .885667717345397D0,.874052796958032D0,.861915468939548D0,
     *  .849262987577969D0,.836102915060907D0,.822443116955644D0,
     *  .808291757507914D0,.793657294762193D0,.778548475506412D0,
     *  .762974330044095D0,.746944166797062D0,.730467566741909D0,
     *  .713554377683587D0,.696214708369514D0,.678458922447719D0/
      DATA T(98),T(99),T(100),T(101),T(102),T(103),T(104),T(105),T(106),
     *     T(107),T(108),T(109),T(110),T(111),T(112),T(113),T(114)/
     *  .660297632272646D0,.641741692562308D0,.622802193910585D0,
     *  .603490456158549D0,.583818021628763D0,.563796648226618D0,
     *  .543438302412810D0,.522755152051175D0,.501759559136144D0,
     *  .480464072404172D0,.458881419833552D0,.437024501037104D0,
     *  .414906379552275D0,.392540275033267D0,.369939555349859D0,
     *  .347117728597636D0,.324088435024413D0/
      DATA T(115),T(116),T(117),T(118),T(119),T(120),T(121),T(122),
     *     T(123),T(124),T(125),T(126),T(127)/.300865438877677D0,
     *  .277462620177904D0,.253893966422694D0,.230173564226660D0,
     *  .206315590902079D0,.182334305985337D0,.158244042714225D0,
     *  .134059199461188D0,.109794231127644D0,
     *  .854636405045155D-1,.610819696041396D-1,.366637909687335D-1,
     *  .122236989606158D-1/
      DATA T(128),T(129),T(130),T(131),T(132),T(133),T(134),T(135),
     *     T(136),T(137),T(138),T(139),T(140),T(141),T(142)/
     *  .999956050018992D0,.999768437409263D0,.999430937466261D0,
     *  .998943525843409D0,.998306266473006D0,.997519252756721D0,
     *  .996582602023382D0,.995496454481096D0,.994260972922410D0,
     *  .992876342608822D0,.991342771207583D0,.989660488745065D0,
     *  .987829747564861D0,.985850822286126D0,.983724009760315D0/
      DATA T(143),T(144),T(145),T(146),T(147),T(148),T(149),T(150),
     *     T(151),T(152),T(153),T(154),T(155),T(156),T(157)/
     *  .981449629025464D0,.979028021257622D0,.976459549719234D0,
     *  .973744599704370D0,.970883578480743D0,.967876915228489D0,
     *  .964725060975706D0,.961428488530732D0,.957987692411178D0,
     *  .954403188769716D0,.950675515316628D0,.946805231239127D0,
     *  .942792917117462D0,.938639174837814D0, .934344627502003D0/
      DATA T(158),T(159),T(160),T(161),T(162),T(163),T(164),T(165),
     *     T(166),T(167),T(168),T(169),T(170),T(171),T(172)/
     *  .929909919334006D0,.925335715583316D0,.920622702425146D0,
     *  .915771586857490D0,.910783096595065D0,.905657979960145D0,
     *  .900397005770304D0,.895000963223085D0,.889470661777611D0,
     *  .883806931033158D0,.878010620604707D0,.872082599995488D0,
     *  .866023758466555D0,.859835004903376D0,.853517267679503D0/
      DATA T(173),T(174),T(175),T(176),T(177),T(178),T(179),T(180),
     *     T(181),T(182),T(183),T(184),T(185),T(186),T(187)/
     *  .847071494517296D0,.840498652345763D0,.833799727155505D0,
     *  .826975723850813D0,.820027666098917D0,.812956596176432D0,
     *  .805763574812999D0,.798449681032171D0,.791016011989546D0,
     *  .783463682808184D0,.775793826411326D0,.768007593352446D0,
     *  .760106151642655D0,.752090686575492D0,.743962400549112D0/
      DATA T(188),T(189),T(190),T(191),T(192),T(193),T(194),T(195),
     *     T(196),T(197),T(198),T(199),T(200),T(201),T(202)/
     *  .735722512885918D0,.727372259649652D0,.718912893459971D0,
     *  .710345683304543D0,.701671914348685D0,.692892887742577D0,
     *  .684009920426076D0,.675024344931163D0,.665937509182049D0,
     *  .656750776292973D0,.647465524363725D0,.638083146272911D0,
     *  .628605049469015D0,.619032655759261D0,.609367401096334D0/
      DATA T(203),T(204),T(205),T(206),T(207),T(208),T(209),T(210),
     *     T(211),T(212),T(213),T(214),T(215),T(216),T(217)/
     *  .599610735362968D0,.589764122154454D0,.579829038559083D0,
     *  .569806974936569D0,.559699434694481D0,.549507934062719D0,
     *  .539234001866059D0,.528879179294822D0,.518445019673674D0,
     *  .507933088228616D0,.497344961852181D0,.486682228866890D0,
     *  .475946488786983D0,.465139352078479D0,.454262439917590D0/
      DATA T(218),T(219),T(220),T(221),T(222),T(223),T(224),T(225),
     *     T(226),T(227),T(228),T(229),T(230),T(231),T(232)/
     *  .443317383947527D0,.432305826033741D0,.421229418017624D0,
     *  .410089821468717D0,.398888707435459D0,.387627756194516D0,
     *  .376308656998716D0,.364933107823654D0,.353502815112970D0,
     *  .342019493522372D0,.330484865662417D0,.318900661840106D0,
     *  .307268619799319D0,.295590484460136D0,.283868007657082D0/
      DATA T(233),T(234),T(235),T(236),T(237),T(238),T(239),T(240),
     *     T(241),T(242),T(243),T(244),T(245),T(246),T(247)/
     *  .272102947876337D0,.260297069991943D0,.248452145001057D0,
     *  .236569949758284D0,.224652266709132D0,.212700883622626D0,
     *  .200717593323127D0,.188704193421389D0,.176662486044902D0,
     *  .164594277567554D0,.152501378338656D0,.140385602411376D0,
     *  .128248767270607D0,.116092693560333D0,.103919204810509D0/
      DATA   T(248),T(249),T(250),T(251),T(252),T(253),T(254),T(255)/
     *  .917301271635196D-1,.795272891002330D-1,.673125211657164D-1,
     *  .550876556946340D-1,.428545265363791D-1,.306149687799790D-1,
     *  .183708184788137D-1,.612391237518953D-2/
      DATA W(1),W(2),W(3),W(4),W(5),W(6),W(7),W(8),W(9),W(10),W(11),
     *     W(12),W(13),W(14),W(15)/1.0D0,.347854845137454D0,
     *  .652145154862546D0,.101228536290376D0,.222381034453374D0,
     *  .313706645877887D0,.362683783378362D0,.271524594117541D-1,
     *  .622535239386479D-1,.951585116824928D-1,.124628971255534D0,
     *  .149595988816577D0,.169156519395003D0,.182603415044924D0,
     *  .189450610455068D0/
      DATA W(16),W(17),W(18),W(19),W(20),W(21),W(22),W(23),W(24),W(25),
     *     W(26),W(27),W(28),W(29),W(30),W(31)/.701861000947010D-2,
     *  .162743947309057D-1,.253920653092621D-1,.342738629130214D-1,
     *  .428358980222267D-1,.509980592623762D-1,.586840934785355D-1,
     *  .658222227763618D-1,.723457941088485D-1,.781938957870703D-1,
     *  .833119242269468D-1,.876520930044038D-1,.911738786957639D-1,
     *  .938443990808046D-1,.956387200792749D-1,.965400885147278D-1/
      DATA W(32),W(33),W(34),W(35),W(36),W(37),W(38),W(39),W(40),W(41),
     *     W(42),W(43),W(44),W(45),W(46),W(47)/.178328072169643D-2,
     *  .414703326056247D-2,.650445796897836D-2,.884675982636395D-2,
     *  .111681394601311D-1,.134630478967186D-1,.157260304760247D-1,
     *  .179517157756973D-1,.201348231535302D-1,.222701738083833D-1,
     *  .243527025687109D-1,.263774697150547D-1,.283396726142595D-1,
     *  .302346570724025D-1,.320579283548516D-1,.338051618371416D-1/
      DATA W(48),W(49),W(50),W(51),W(52),W(53),W(54),W(55),W(56),W(57),
     *     W(58),W(59),W(60),W(61),W(62),W(63)/ .354722132568824D-1,
     *  .370551285402400D-1,.385501531786156D-1,.399537411327203D-1,
     *  .412625632426235D-1,.424735151236536D-1,.435837245293235D-1,
     *  .445905581637566D-1,.454916279274181D-1, .462847965813144D-1,
     *  .469681828162100D-1,.475401657148303D-1,.479993885964583D-1,
     *  .483447622348030D-1,.485754674415034D-1,.486909570091397D-1/
      DATA W(64),W(65),W(66),W(67),W(68),W(69),W(70),W(71),W(72),W(73),
     *     W(74),W(75),W(76),W(77),W(78),W(79)/ .449380960292090D-3,
     *  .104581267934035D-2,.164250301866903D-2,.223828843096262D-2,
     *  .283275147145799D-2,.342552604091022D-2,.401625498373864D-2,
     *  .460458425670296D-2,.519016183267633D-2,.577263754286570D-2,
     *  .635166316170719D-2,.692689256689881D-2,.749798192563473D-2,
     *  .806458989048606D-2,.862637779861675D-2,.918300987166087D-2/
      DATA W(80),W(81),W(82),W(83),W(84),W(85),W(86),W(87),W(88),
     *   W(89),W(90),W(91),W(92),W(93),W(94),W(95)/.973415341500681D-2,
     *  .102794790158322D-1,.108186607395031D-1,.113513763240804D-1,
     *  .118773073727403D-1,.123961395439509D-1,.129075627392673D-1,
     *  .134112712886163D-1,.139069641329520D-1,.143943450041668D-1,
     *  .148731226021473D-1,.153430107688651D-1,.158037286593993D-1,
     *  .162550009097852D-1,.166965578015892D-1,.171281354231114D-1/
      DATA W(96),W(97),W(98),W(99),W(100),W(101),W(102),W(103),W(104),
     *     W(105),W(106),W(107),W(108),W(109),W(110),W(111),W(112)/
     *  .175494758271177D-1,.179603271850087D-1,.183604439373313D-1,
     *  .187495869405447D-1,.191275236099509D-1,.194940280587066D-1,
     *  .198488812328309D-1,.201918710421300D-1,.205227924869601D-1,
     *  .208414477807511D-1,.211476464682213D-1,.214412055392085D-1,
     *  .217219495380521D-1,.219897106684605D-1,.222443288937998D-1,
     *  .224856520327450D-1,.227135358502365D-1/
      DATA W(113),W(114),W(115),W(116),W(117),W(118),W(119),W(120),
     *     W(121),W(122),W(123),W(124),W(125),W(126),W(127)/
     *  .229278441436868D-1,.231284488243870D-1,.233152299940628D-1,
     *  .234880760165359D-1,.236468835844476D-1,.237915577810034D-1,
     *  .239220121367035D-1,.240381686810241D-1,.241399579890193D-1,
     *  .242273192228152D-1,.243002001679719D-1,.243585572646906D-1,
     *  .244023556338496D-1,.244315690978500D-1,.244461801962625D-1/
      DATA W(128),W(129),W(130),W(131),W(132),W(133),W(134),W(135),
     *     W(136),W(137),W(138),W(139),W(140),W(141),W(142)/
     *  .112789017822272D-3,.262534944296446D-3,.412463254426176D-3,
     *  .562348954031410D-3,.712154163473321D-3,.861853701420089D-3,
     *  .101142439320844D-2,.116084355756772D-2,.131008868190250D-2,
     *  .145913733331073D-2,.160796713074933D-2,.175655573633073D-2,
     *  .190488085349972D-2,.205292022796614D-2,.220065164983991D-2/
      DATA W(143),W(144),W(145),W(146),W(147),W(148),W(149),W(150),
     *     W(151),W(152),W(153),W(154),W(155),W(156),W(157)/
     *  .234805295632731D-2,.249510203470371D-2,.264177682542749D-2,
     *  .278805532532771D-2,.293391559082972D-2,.307933574119934D-2,
     *  .322429396179420D-2,.336876850731555D-2,.351273770505631D-2,
     *  .365617995814250D-2,.379907374876626D-2,.394139764140883D-2,
     *  .408313028605267D-2,.422425042138154D-2,.436473687796806D-2/
      DATA W(158),W(159),W(160),W(161),W(162),W(163),W(164),W(165),
     *     W(166),W(167),W(168),W(169),W(170),W(171),W(172)/
     *  .450456858144790D-2,.464372455568006D-2,.478218392589269D-2,
     *  .491992592181387D-2,.505692988078684D-2,.519317525086928D-2,
     *  .532864159391593D-2,.546330858864431D-2,.559715603368291D-2,
     *  .573016385060144D-2,.586231208692265D-2,.599358091911534D-2,
     *  .612395065556793D-2,.625340173954240D-2,.638191475210788D-2/
      DATA W(173),W(174),W(175),W(176),W(177),W(178),W(179),W(180),
     *     W(181),W(182),W(183),W(184),W(185),W(186),W(187)/
     *  .650947041505366D-2,.663604959378107D-2,.676163330017380D-2,
     *  .688620269544632D-2,.700973909296982D-2,.713222396107539D-2,
     *  .725363892583391D-2,.737396577381235D-2,.749318645480588D-2,
     *  .761128308454566D-2,.772823794738156D-2,.784403349893971D-2,
     *  .795865236875435D-2,.807207736287350D-2,.818429146643827D-2/
      DATA W(188),W(189),W(190),W(191),W(192),W(193),W(194),W(195),
     *     W(196),W(197),W(198),W(199),W(200),W(201),W(202)/
     *  .829527784623523D-2,.840501985322154D-2,.851350102502249D-2,
     *  .862070508840101D-2,.872661596169881D-2,.883121775724875D-2,
     *  .893449478375821D-2,.903643154866287D-2,.913701276045081D-2,
     *  .923622333095630D-2,.933404837762327D-2,.943047322573775D-2,
     *  .952548341062928D-2,.961906467984073D-2,.971120299526628D-2/
      DATA W(203),W(204),W(205),W(206),W(207),W(208),W(209),W(210),
     *     W(211),W(212),W(213),W(214),W(215),W(216),W(217)/
     *  .980188453525733D-2,.989109569669583D-2,.997882309703491D-2,
     *  .100650535763064D-1,.101497741990949D-1,.102329722564782D-1,
     *  .103146352679340D-1,.103947509832117D-1,.104733073841704D-1,
     *  .105502926865815D-1,.106256953418966D-1,.106995040389798D-1,
     *  .107717077058046D-1,.108422955111148D-1,.109112568660490D-1/
      DATA W(218),W(219),W(220),W(221),W(222),W(223),W(224),W(225),
     *     W(226),W(227),W(228),W(229),W(230),W(231),W(232)/
     *  .109785814257296D-1,.110442590908139D-1,.111082800090098D-1,
     *  .111706345765534D-1,.112313134396497D-1,.112903074958755D-1,
     *  .113476078955455D-1,.114032060430392D-1,.114570935980906D-1,
     *  .115092624770395D-1,.115597048540436D-1,.116084131622531D-1,
     *  .116553800949452D-1,.117005986066207D-1,.117440619140606D-1/
      DATA W(233),W(234),W(235),W(236),W(237),W(238),W(239),W(240),
     *     W(241),W(242),W(243),W(244),W(245),W(246),W(247)/
     *  .117857634973434D-1,.118256971008240D-1,.118638567340711D-1,
     *  .119002366727665D-1,.119348314595636D-1,.119676359049059D-1,
     $  .119986450878058D-1,.120278543565826D-1,.120552593295601D-1,
     *  .120808558957245D-1,.121046402153405D-1,.121266087205273D-1,
     *  .121467581157945D-1,.121650853785355D-1,.121815877594818D-1/
      DATA W(248),W(249),W(250),W(251),W(252),W(253),W(254),W(255)/
     *  .121962627831147D-1,.122091082480372D-1,.122201222273040D-1,
     *  .122293030687103D-1,.122366493950402D-1,.122421601042728D-1,
     *  .122458343697479D-1,.122476716402898D-1/
      SCALE=(B-A)/2.0D0
      M=N/2
      NPLACE=M-1
      DO 1 I=1,M
      S=T(NPLACE+I)
      R=W(NPLACE+I)*SCALE
      TT(I)=(A*(1.0D0+S)+(1.0D0-S)*B)/2.0D0
      TT( N+1-I)=(A*(1.0D0-S)+B*(1.0D0+S))/2.0D0
      WW(I)=R
1     WW(N+1-I)=R
      RETURN
      END

      SUBROUTINE LINSYS(MAT,B,N,MD,IER)
C     -----------------
C
C     THIS ROUTINE SOLVES A SYSTEM OF LINEAR EQUATIONS
C             A*X = B
C     THE METHOD USED IS GAUSSIAN ELIMINATION WITH
C     PARTIAL PIVOTING.
C
C     INPUT:
C     THE COEFFICIENT MATRIX A IS STORED IN THE ARRAY MAT.
C     THE RIGHT SIDE CONSTANTS ARE IN THE ARRAY B.
C     THE ORDER OF THE LINEAR SYSTEM IS N.
C     THE VARIABLE MD IS THE NUMBER OF ROWS THAT MAT
C     IS DIMENSIONED AS HAVING IN THE CALLING PROGRAM.
C
C     OUTPUT:
C     THE ARRAY B CONTAINS THE SOLUTION X.
C     MAT CONTAINS THE UPPER TRIANGULAR MATRIX U
C     OBTAINED BY ELIMINATION. THE ROW MULTIPLIERS
C     USED IN THE ELIMINATION ARE STORED IN THE
C     LOWER TRIANGULAR PART OF MAT.
C     IER=0 MEANS THE MATRIX A WAS COMPUTATIONALLY
C     NONSINGULAR, AND THE GAUSSIAN ELIMINATION
C     WAS COMPLETED SATISFACTORILY.
C     IER=1 MEANS THAT THE MATRIX A WAS
C     COMPUTATIONALLY SINGULAR.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER PIVOT
      DOUBLE PRECISION MULT,MAT
      DIMENSION MAT(MD,*),B(*)
C
C     BEGIN ELIMINATION STEPS.
      DO 40 K=1,N-1
C       CHOOSE PIVOT ROW.
        PIVOT=K
        AMAX=ABS(MAT(K,K))
        DO 10 I=K+1,N
          ABSA=ABS(MAT(I,K))
          IF(ABSA .GT. AMAX) THEN
              PIVOT=I
              AMAX=ABSA
          END IF
10        CONTINUE
C
        IF(AMAX .EQ. 0.0D0) THEN
C           COEFFICIENT MATRIX IS SINGULAR.
            IER=1
            RETURN
        END IF
C
        IF(PIVOT .NE. K) THEN
C           SWITCH ROWS K AND PIVOT.
            TEMP=B(K)
            B(K)=B(PIVOT)
            B(PIVOT)=TEMP
            DO 20 J=K,N
              TEMP=MAT(K,J)
              MAT(K,J)=MAT(PIVOT,J)
20            MAT(PIVOT,J)=TEMP
        END IF
C
C       PERFORM STEP #K OF ELIMINATION.
        DO 30 I=K+1,N
          MULT=MAT(I,K)/MAT(K,K)
          MAT(I,K)=MULT
          B(I)=B(I)-MULT*B(K)
          DO 30 J=K+1,N
30          MAT(I,J)=MAT(I,J)-MULT*MAT(K,J)
40      CONTINUE
C
      IF(MAT(N,N) .EQ. 0.0D0) THEN
C         COEFFICIENT MATRIX IS SINGULAR.
          IER=1
          RETURN
      END IF
C
C     SOLVE FOR SOLUTION X USING BACK SUBSTITUTION.
      DO 60 I=N,1,-1
        SUM=0.0D0
        DO 50 J=I+1,N
50        SUM=SUM+MAT(I,J)*B(J)
60      B(I)=(B(I)-SUM)/MAT(I,I)
      IER=0
      RETURN
      END

      DOUBLE PRECISION FUNCTION D1MACH (I) 
C     --------------------------------
C
C***BEGIN PROLOGUE  D1MACH
C***DATE WRITTEN   750101   (YYMMDD)
C***REVISION DATE  860501   (YYMMDD)
C***CATEGORY NO.  R1
C***KEYWORDS  MACHINE CONSTANTS
C***AUTHOR  FOX, P. A., (BELL LABS)
C           HALL, A. D., (BELL LABS)
C           SCHRYER, N. L., (BELL LABS)
C***PURPOSE  RETURN DOUBLE PRECISION MACHINE DEPENDENT CONSTANTS.
C***DESCRIPTION
C
C     D1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS
C     FOR THE LOCAL MACHINE ENVIRONMENT.  IT IS A FUNCTION
C     SUBPROGRAM WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED
C     AS FOLLOWS, FOR EXAMPLE
C
C          D = D1MACH(I)
C
C     WHERE I=1,...,5.  THE (OUTPUT) VALUE OF D ABOVE IS
C     DETERMINED BY THE (INPUT) VALUE OF I.  THE RESULTS FOR
C     VARIOUS VALUES OF I ARE DISCUSSED BELOW.
C
C  DOUBLE-PRECISION MACHINE CONSTANTS
C  D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
C  D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
C  D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
C  D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
C  D1MACH( 5) = LOG10(B)
C***REFERENCES  FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A
C                 PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL
C                 SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188.
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  D1MACH
C
      INTEGER SMALL(4)
      INTEGER LARGE(4)
      INTEGER RIGHT(4)
      INTEGER DIVER(4)
      INTEGER LOG10(4)
C
      DOUBLE PRECISION DMACH(5)
C
      EQUIVALENCE (DMACH(1),SMALL(1))
      EQUIVALENCE (DMACH(2),LARGE(1))
      EQUIVALENCE (DMACH(3),RIGHT(1))
      EQUIVALENCE (DMACH(4),DIVER(1))
      EQUIVALENCE (DMACH(5),LOG10(1))
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C     DATA SMALL(1) / ZC00800000 /
C     DATA SMALL(2) / Z000000000 /
C
C     DATA LARGE(1) / ZDFFFFFFFF /
C     DATA LARGE(2) / ZFFFFFFFFF /
C
C     DATA RIGHT(1) / ZCC5800000 /
C     DATA RIGHT(2) / Z000000000 /
C
C     DATA DIVER(1) / ZCC6800000 /
C     DATA DIVER(2) / Z000000000 /
C
C     DATA LOG10(1) / ZD00E730E7 /
C     DATA LOG10(2) / ZC77800DC0 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
C
C     DATA SMALL(1) / O1771000000000000 /
C     DATA SMALL(2) / O0000000000000000 /
C
C     DATA LARGE(1) / O0777777777777777 /
C     DATA LARGE(2) / O0007777777777777 /
C
C     DATA RIGHT(1) / O1461000000000000 /
C     DATA RIGHT(2) / O0000000000000000 /
C
C     DATA DIVER(1) / O1451000000000000 /
C     DATA DIVER(2) / O0000000000000000 /
C
C     DATA LOG10(1) / O1157163034761674 /
C     DATA LOG10(2) / O0006677466732724 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
C
C     DATA SMALL(1) / O1771000000000000 /
C     DATA SMALL(2) / O7770000000000000 /
C
C     DATA LARGE(1) / O0777777777777777 /
C     DATA LARGE(2) / O7777777777777777 /
C
C     DATA RIGHT(1) / O1461000000000000 /
C     DATA RIGHT(2) / O0000000000000000 /
C
C     DATA DIVER(1) / O1451000000000000 /
C     DATA DIVER(2) / O0000000000000000 /
C
C     DATA LOG10(1) / O1157163034761674 /
C     DATA LOG10(2) / O0006677466732724 /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C     FOR FTN4
C
C     DATA SMALL(1) / 00564000000000000000B /
C     DATA SMALL(2) / 00000000000000000000B /
C
C     DATA LARGE(1) / 37757777777777777777B /
C     DATA LARGE(2) / 37157777777777777777B /
C
C     DATA RIGHT(1) / 15624000000000000000B /
C     DATA RIGHT(2) / 00000000000000000000B /
C
C     DATA DIVER(1) / 15634000000000000000B /
C     DATA DIVER(2) / 00000000000000000000B /
C
C     DATA LOG10(1) / 17164642023241175717B /
C     DATA LOG10(2) / 16367571421742254654B /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C     FOR FTN5
C
C     DATA SMALL(1) / O"00564000000000000000" /
C     DATA SMALL(2) / O"00000000000000000000" /
C
C     DATA LARGE(1) / O"37757777777777777777" /
C     DATA LARGE(2) / O"37157777777777777777" /
C
C     DATA RIGHT(1) / O"15624000000000000000" /
C     DATA RIGHT(2) / O"00000000000000000000" /
C
C     DATA DIVER(1) / O"15634000000000000000" /
C     DATA DIVER(2) / O"00000000000000000000" /
C
C     DATA LOG10(1) / O"17164642023241175717" /
C     DATA LOG10(2) / O"16367571421742254654" /
C
C     MACHINE CONSTANTS FOR THE CRAY 1
C
C     DATA SMALL(1) / 201354000000000000000B /
C     DATA SMALL(2) / 000000000000000000000B /
C
C     DATA LARGE(1) / 577767777777777777777B /
C     DATA LARGE(2) / 000007777777777777774B /
C
C     DATA RIGHT(1) / 376434000000000000000B /
C     DATA RIGHT(2) / 000000000000000000000B /
C
C     DATA DIVER(1) / 376444000000000000000B /
C     DATA DIVER(2) / 000000000000000000000B /
C
C     DATA LOG10(1) / 377774642023241175717B /
C     DATA LOG10(2) / 000007571421742254654B /
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD -
C     STATIC DMACH(5)
C
C     DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/
C     DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/
C     DATA LOG10/40423K,42023K,50237K,74776K/
C
C     MACHINE CONSTANTS FOR THE HARRIS 220
C
C     DATA SMALL(1), SMALL(2) / '20000000, '00000201 /
C     DATA LARGE(1), LARGE(2) / '37777777, '37777577 /
C     DATA RIGHT(1), RIGHT(2) / '20000000, '00000333 /
C     DATA DIVER(1), DIVER(2) / '20000000, '00000334 /
C     DATA LOG10(1), LOG10(2) / '23210115, '10237777 /
C
C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.
C
C     DATA SMALL(1), SMALL(2) / O402400000000, O000000000000 /
C     DATA LARGE(1), LARGE(2) / O376777777777, O777777777777 /
C     DATA RIGHT(1), RIGHT(2) / O604400000000, O000000000000 /
C     DATA DIVER(1), DIVER(2) / O606400000000, O000000000000 /
C     DATA LOG10(1), LOG10(2) / O776464202324, O117571775714 /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     THREE WORD DOUBLE PRECISION OPTION WITH FTN4
C
C     DATA SMALL(1), SMALL(2), SMALL(3) / 40000B,       0,       1 /
C     DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B /
C     DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B,       0,    265B /
C     DATA DIVER(1), DIVER(2), DIVER(3) / 40000B,       0,    276B /
C     DATA LOG10(1), LOG10(2), LOG10(3) / 46420B,  46502B,  77777B /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     FOUR WORD DOUBLE PRECISION OPTION WITH FTN4
C
C     DATA SMALL(1), SMALL(2) /  40000B,       0 /
C     DATA SMALL(3), SMALL(4) /       0,       1 /
C     DATA LARGE(1), LARGE(2) /  77777B, 177777B /
C     DATA LARGE(3), LARGE(4) / 177777B, 177776B /
C     DATA RIGHT(1), RIGHT(2) /  40000B,       0 /
C     DATA RIGHT(3), RIGHT(4) /       0,    225B /
C     DATA DIVER(1), DIVER(2) /  40000B,       0 /
C     DATA DIVER(3), DIVER(4) /       0,    227B /
C     DATA LOG10(1), LOG10(2) /  46420B,  46502B /
C     DATA LOG10(3), LOG10(4) /  76747B, 176377B /
C
C     MACHINE CONSTANTS FOR THE HP 9000
C
C     D1MACH(1) = 2.8480954D-306
C     D1MACH(2) = 1.40444776D+306
C     D1MACH(3) = 2.22044605D-16
C     D1MACH(4) = 4.44089210D-16
C     D1MACH(5) = 3.01029996D-1
C
C     DATA SMALL(1), SMALL(2) / 00040000000B, 00000000000B /
C     DATA LARGE(1), LARGE(2) / 17737777777B, 37777777777B /
C     DATA RIGHT(1), RIGHT(2) / 07454000000B, 00000000000B /
C     DATA DIVER(1), DIVER(2) / 07460000000B, 00000000000B /
C     DATA LOG10(1), LOG10(2) / 07764642023B, 12047674777B /
C
C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND
C     THE PERKIN ELMER (INTERDATA) 7/32.
C
C     DATA SMALL(1), SMALL(2) / Z00100000, Z00000000 /
C     DATA LARGE(1), LARGE(2) / Z7FFFFFFF, ZFFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / Z33100000, Z00000000 /
C     DATA DIVER(1), DIVER(2) / Z34100000, Z00000000 /
C     DATA LOG10(1), LOG10(2) / Z41134413, Z509F79FF /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
C
C     DATA SMALL(1), SMALL(2) / "033400000000, "000000000000 /
C     DATA LARGE(1), LARGE(2) / "377777777777, "344777777777 /
C     DATA RIGHT(1), RIGHT(2) / "113400000000, "000000000000 /
C     DATA DIVER(1), DIVER(2) / "114400000000, "000000000000 /
C     DATA LOG10(1), LOG10(2) / "177464202324, "144117571776 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
C
C     DATA SMALL(1), SMALL(2) / "000400000000, "000000000000 /
C     DATA LARGE(1), LARGE(2) / "377777777777, "377777777777 /
C     DATA RIGHT(1), RIGHT(2) / "103400000000, "000000000000 /
C     DATA DIVER(1), DIVER(2) / "104400000000, "000000000000 /
C     DATA LOG10(1), LOG10(2) / "177464202324, "476747767461 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1), SMALL(2) /    8388608,           0 /
C     DATA LARGE(1), LARGE(2) / 2147483647,          -1 /
C     DATA RIGHT(1), RIGHT(2) /  612368384,           0 /
C     DATA DIVER(1), DIVER(2) /  620756992,           0 /
C     DATA LOG10(1), LOG10(2) / 1067065498, -2063872008 /
C
C     DATA SMALL(1), SMALL(2) / O00040000000, O00000000000 /
C     DATA LARGE(1), LARGE(2) / O17777777777, O37777777777 /
C     DATA RIGHT(1), RIGHT(2) / O04440000000, O00000000000 /
C     DATA DIVER(1), DIVER(2) / O04500000000, O00000000000 /
C     DATA LOG10(1), LOG10(2) / O07746420232, O20476747770 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1), SMALL(2) /    128,      0 /
C     DATA SMALL(3), SMALL(4) /      0,      0 /
C
C     DATA LARGE(1), LARGE(2) /  32767,     -1 /
C     DATA LARGE(3), LARGE(4) /     -1,     -1 /
C
C     DATA RIGHT(1), RIGHT(2) /   9344,      0 /
C     DATA RIGHT(3), RIGHT(4) /      0,      0 /
C
C     DATA DIVER(1), DIVER(2) /   9472,      0 /
C     DATA DIVER(3), DIVER(4) /      0,      0 /
C
C     DATA LOG10(1), LOG10(2) /  16282,   8346 /
C     DATA LOG10(3), LOG10(4) / -31493, -12296 /
C
C     DATA SMALL(1), SMALL(2) / O000200, O000000 /
C     DATA SMALL(3), SMALL(4) / O000000, O000000 /
C
C     DATA LARGE(1), LARGE(2) / O077777, O177777 /
C     DATA LARGE(3), LARGE(4) / O177777, O177777 /
C
C     DATA RIGHT(1), RIGHT(2) / O022200, O000000 /
C     DATA RIGHT(3), RIGHT(4) / O000000, O000000 /
C
C     DATA DIVER(1), DIVER(2) / O022400, O000000 /
C     DATA DIVER(3), DIVER(4) / O000000, O000000 /
C
C     DATA LOG10(1), LOG10(2) / O037632, O020232 /
C     DATA LOG10(3), LOG10(4) / O102373, O147770 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER
C
C     DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 /
C     DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 /
C     DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 /
C     DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 /
C     DATA LOG10(1), LOG10(2) / O177746420232, O411757177572 /
C
C     MACHINE CONSTANTS FOR VAX 11/780
C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
C     ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***
C     *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C     DATA SMALL(1), SMALL(2) /        128,           0 /
C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
C     DATA RIGHT(1), RIGHT(2) /       9344,           0 /
C     DATA DIVER(1), DIVER(2) /       9472,           0 /
C     DATA LOG10(1), LOG10(2) /  546979738,  -805796613 /
C
C     DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 /
C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 /
C     DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 /
C     DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB /
C
C     MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING)
C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
C     ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***
C     *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C     DATA SMALL(1), SMALL(2) /         16,           0 /
C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
C     DATA RIGHT(1), RIGHT(2) /      15552,           0 /
C     DATA DIVER(1), DIVER(2) /      15568,           0 /
C     DATA LOG10(1), LOG10(2) /  1142112243, 2046775455 /
C
C     DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 /
C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 /
C     DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 /
C     DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F /
C
C     MACHINE CONSTANTS FOR THE ELXSI 6400
C     (ASSUMING REAL*8 IS THE DEFAULT DOUBLE PRECISION)
C
C     DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X /
C     DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF'X /
C     DATA RIGHT(1), RIGHT(2) / '3CB00000'X,'00000000'X /
C     DATA DIVER(1), DIVER(2) / '3CC00000'X,'00000000'X /
C     DATA LOG10(1), DIVER(2) / '3FD34413'X,'509F79FF'X /
C
C     MACHINE CONSTANTS FOR THE IBM PC - MICROSOFT FORTRAN
C
C     DATA SMALL(1), SMALL(2) / #00000000, #00100000 /
C     DATA LARGE(1), LARGE(2) / #FFFFFFFF, #7FEFFFFF /
C     DATA RIGHT(1), RIGHT(2) / #00000000, #3CA00000 /
C     DATA DIVER(1), DIVER(2) / #00000000, #3CB00000 /
C     DATA LOG10(1), LOG10(2) / #509F79FF, #3FD34413 /
C
C     MACHINE CONSTANTS FOR THE IBM PC - PROFESSIONAL FORTRAN
C     AND LAHEY FORTRAN
C
C     DATA SMALL(1), SMALL(2) / Z'00000000', Z'00100000' /
C     DATA LARGE(1), LARGE(2) / Z'FFFFFFFF', Z'7FEFFFFF' /
C     DATA RIGHT(1), RIGHT(2) / Z'00000000', Z'3CA00000' /
C     DATA DIVER(1), DIVER(2) / Z'00000000', Z'3CB00000' /
C     DATA LOG10(1), LOG10(2) / Z'509F79FF', Z'3FD34413' /
C
C     MACHINE CONSTANTS FOR THE PRIME 750/850.   {9/23/86 TSENG}
C
C     DATA SMALL(1), SMALL(2) / :10000000000, :00000100000 /
C     DATA LARGE(1), LARGE(2) / :17777777777, :37777677602 /
C     DATA RIGHT(1), RIGHT(2) / :10000000000, :00000000122 /
C     DATA DIVER(1), DIVER(2) / :10000000000, :00000000123 /
C     DATA LOG10(1), LOG10(2) / :11504046502, :17572000177 /
C
C     MACHINE CONSTANTS FOR THE APOLLO.   {9/23/86 TSENG}
C
      DATA SMALL(1), SMALL(2) / 8#00004000000, 8#00000000000 /
      DATA LARGE(1), LARGE(2) / 8#17773777777, 8#37777777777 /
      DATA RIGHT(1), RIGHT(2) / 8#07450000000, 8#00000000000 /
      DATA DIVER(1), DIVER(2) / 8#07454000000, 8#00000000000 /
      DATA LOG10(1), LOG10(2) / 8#07764642023, 8#12047674777 /
C
C
C***FIRST EXECUTABLE STATEMENT  D1MACH
C     IF (I .LT. 1  .OR.  I .GT. 5)
C    1   CALL XERROR( 'D1MACH -- I OUT OF BOUNDS',25,1,2)
C
      D1MACH = DMACH(I)
      RETURN
C
      END

