C***************************************************************
C *   Test program for subroutine H3LC by Stephen Kirkup      
C***************************************************************
C
C  Copyright 1998- Stephen Kirkup
C  School of Computing Engineering and Physical Sciences
C  University of Central Lancashire - www.uclan.ac.uk 
C  United Kingdom
C  smkirkup@uclan.ac.uk
C  http://www.researchgate.net/profile/Stephen_Kirkup
C
C  This open source code can be found at
C   www.boundary-element-method.com/fortran/H3LC_T.FOR
C
C  Issued under the GNU General Public License 2007, see gpl.txt
C
C  Part of the the author's open source BEM packages. 
C  All codes and manuals can be downloaded from 
C  www.boundary-element-method.com
C
C***************************************************************
C
C
C This program is a test for subroutine H3LC. The program computes the
C  solution to the Helmholtz problem  exterior to a cube via the integral
C  equation of Burton & Miller. In order to use H3LC, the cube is
C  represented by 24 uniform triangles. The boundary functions are
C  approximated by a constant on each element.
C Simple changes in the program allow the set up of the following data:
C  (1) the choice of quadrature rule,
C  (2) the geometry of the surface (of the approximating polyhedron),
C  (3) the set of wavenumbers,
C  (4) the parameter(s) in the Burton & Miller integral equation,
C  (5) the position of one source point which determines the boundary
C       condition and the exact solution,
C  (6) the points in the exterior where the solution is sought.
C During execution the program gives the solution at the collocation
C  points (the points at the centre of each element) and the solution at
C  the selected exterior point. The program also give the exact
C  solution at the same points so that computed and exact solutions
C  may be compared.

C The PARAMETER statement
C =======================
C There are five components in the PARAMETER statement.
C integer LIMK  : The limit on the number of wavenumbers.
C integer LIMN  : The limit on the number of elements.
C integer LIMNV : The limit on the number of vertices.
C integer LIMNQ : The limit on the number of quadrature points allowed
C                 in the quadrature rules.
C integer LIMNPE: The limit on the number of exterior points.

C External modules related to the test program
C ============================================
C Subroutine CLINSL: This computes the solution to a linear system of
C  equations.
C real function FNSQRT: Returns the square root.
C complex function FNEXP: Returns the exponential.

C External modules related to the package
C =======================================
C Subroutine H3LC: This computes the discrete Helmholtz integral
C  operators.

      PROGRAM H3LCT
C Declaration of variables
C  Limits on the size of the arrays
C   Limit on the number of wavenumbers
      INTEGER    LIMK
C   Limit on the number of elements
      INTEGER    LIMN
C   Limit on the number of vertices
      INTEGER    LIMNV
C   Limit on the number of quadrature points in any quadrature rule
      INTEGER    LIMNQ
C   Limit on the number of exterior points where a solution is sought
      INTEGER    LIMNPE
C  Set PARAMETERs
      PARAMETER (LIMK=4,LIMN=24,LIMNV=20,LIMNQ=30,LIMNPE=4)
C  Constants
      REAL*8     ZERO,ONE,TWO,FOUR,EIGHT,THIRD,THREE,PI,ROOT2
      COMPLEX*16 CIMAG
C  Number of elements
      INTEGER*4  N
C  Wavenumber information
      INTEGER    NOK
      COMPLEX*16 K,KVAL(LIMK)
C  Geometrical information
      REAL*8     VERTPT(LIMNV,3),ELNORM(LIMN,3)
      INTEGER    ELVERT(LIMN,3),NOVERT,NPEXT
      REAL*8     P(3),NORMP(3),Q(3),QA(3),QB(3),QC(3),NORMQ(3)
      REAL*8     PEXT(LIMNPE,3),PP(3),RR(3),R,DRBDNQ,VECP(3)
      LOGICAL    LPONEL
C  Error flag used for checking the set-up information
      LOGICAL    LERROR
C  Quadrature rules
C   General quadrature rule
C    Number of points
      INTEGER    NQ
C    Points (x and y coordinates) and weights
      REAL*8     XQ(7),YQ(7),WQ(7)
C   Particular quadrature rule
C    Number of points
      INTEGER    NQT
C    Points (x and y coordinates) and weights
      REAL*8     XQT(LIMNQ),YQT(LIMNQ),WQT(LIMNQ)
C  Validation and control variables in H3LC
      LOGICAL    LVALID,LFAIL
      REAL*8     EK,EGEOM,EQRULE
      LOGICAL    LLK,LMK,LMKT,LNK
C  Storage of discrete integral operators on H3LC
      COMPLEX*16 DISLK,DISMK,DISMKT,DISNK
C  Matrices corresponding to the Helmholtz integral operators
      COMPLEX*16 LK(LIMN,LIMN),MK(LIMN,LIMN)
      COMPLEX*16 MKT(LIMN,LIMN),NK(LIMN,LIMN)
C  Parameters in the Burton and Miller integral equation for each
C   wavenumber
      COMPLEX*16 ALPHA(LIMK),BETA(LIMK)
C  Composite matrices used in the Burton and Miller equation
      COMPLEX*16 AK(LIMN,LIMN),BK(LIMN,LIMN)
C  Exact Neumann boundary condition, exact solution and computed 
C   solution
      COMPLEX*16 VEL(LIMN),PHI(LIMN),APPPHI(LIMN)
C  Solution at the exterior point
      COMPLEX*16 PHIPT
C  Temporary variable used in computing the solutions
      COMPLEX*16 BVEL(LIMN),SUM
C  Counters and indices
      INTEGER*4 I,IQ,J,I1,I2,I3
C  Variables used in determining the test problem
      COMPLEX*16 E
      COMPLEX*16 KR,IKR
C  Counters
      INTEGER    IPE,IK

C  Set up constants
      ZERO=0.0D0
      THIRD=1.0D0/3.0D0
      ONE=1.0D0
      TWO=2.0D0
      THREE=3.0D0
      FOUR=4.0D0
      EIGHT=8.0D0
      CIMAG=CMPLX(ZERO,ONE)
      PI=FOUR*ATAN(ONE)
      ROOT2=SQRT(TWO)

C  Set up control and tolerances and open file for storing the possible
C   error messages
      LVALID=.TRUE.
      EK=1D-08
      EGEOM=1D-08
      EQRULE=1D-08
      OPEN(UNIT=10,FILE='H3LC.ERR',STATUS='UNKNOWN')

C  Set up standard Gaussian Quadrature rule. Stores the values
C   for the 7 point Gaussian quadrature over the standard triangle.
C   Taken from "Some criteria for numerically integrated matrices and
C   quadrature formulas for triangles" by M.E.Laursen and M.Gellert in
C   International Journal for numerical methods in engineering Vol 12 
C   67-76 1978.

      NQ=7
      DATA WQ/ 0.225000000000000D0,
     *         0.125939180544827D0,
     *         0.125939180544827D0,
     *         0.125939180544827D0,
     *         0.132394152788506D0,
     *         0.132394152788506D0,
     *         0.132394152788506D0/

      DATA XQ/ 0.333333333333333D0,
     *         0.797426985353087D0,
     *         0.101286507323456D0,
     *         0.101286507323456D0,
     *         0.470142064105115D0,
     *         0.470142064105115D0,
     *         0.059715871789770D0/

      DATA YQ/ 0.333333333333333D0,
     *         0.101286507323456D0,
     *         0.797426985353087D0,
     *         0.101286507323456D0,
     *         0.470142064105115D0,
     *         0.059715871789770D0,
     *         0.470142064105115D0/



C  Set up the geometry of the structure  by completing the data
C   structures VERTPT and ELVERT.

C   Set up VERTPT, the coordinates of the vertices of the elements
      NOVERT=14
      DATA ((VERTPT(I,J),J=1,3),I=1,14)
     * /-1.0D0,-1.0D0, 1.0D0,
     *  -1.0D0, 1.0D0, 1.0D0,
     *  -1.0D0, 1.0D0,-1.0D0,
     *  -1.0D0,-1.0D0,-1.0D0,
     *   1.0D0,-1.0D0, 1.0D0,
     *   1.0D0, 1.0D0, 1.0D0,
     *   1.0D0, 1.0D0,-1.0D0,
     *   1.0D0,-1.0D0,-1.0D0,
     *   0.0D0, 0.0D0, 1.0D0,
     *   0.0D0, 1.0D0, 0.0D0,
     *  -1.0D0, 0.0D0, 0.0D0,
     *   0.0D0, 0.0D0,-1.0D0,
     *   1.0D0, 0.0D0, 0.0D0,
     *   0.0D0,-1.0D0, 0.0D0/

C   Set the number of elements (N)
      N=24

C   Set up ELEMVERT, the list of elements and the indices of the 
C    vertices that describe the elements. Since the integral equation 
C    that is applied assumes the normal to be outward, the vertices are
C    enumerated anticlockwise about the element when viewed from a 
C    point near to the element but exterior to the surface. This is
C    so that the normals to the elements are implicitly defined as
C    being outward within the subroutine H3LC.
      DATA ((ELVERT(I,J),J=1,3),I=1,24) /
     *   9, 1, 5,    9, 5, 6,    9, 6, 2,    9, 2, 1,
     *  10, 2, 6,   10, 6, 7,   10, 7, 3,   10, 3, 2,
     *  11, 2, 3,   11, 1, 2,   11, 4, 1,   11, 3, 4,
     *  12, 3, 7,   12, 4, 3,   12, 8, 4,   12, 7, 8,
     *  13, 7, 6,   13, 8, 7,   13, 5, 8,   13, 6, 5,
     *  14, 4, 8,   14, 1, 4,   14, 5, 1,   14, 8, 5 /


C   The outward normals to the elements are stated explicitly. Equally,
C    they could have been computed from the description of each element
      DATA ((ELNORM(I,J),J=1,3),I=1,12) /
     *   0.0D0, 0.0D0, 1.0D0,
     *   0.0D0, 0.0D0, 1.0D0,
     *   0.0D0, 0.0D0, 1.0D0,
     *   0.0D0, 0.0D0, 1.0D0,
     *   0.0D0, 1.0D0, 0.0D0,
     *   0.0D0, 1.0D0, 0.0D0,
     *   0.0D0, 1.0D0, 0.0D0,
     *   0.0D0, 1.0D0, 0.0D0,
     *  -1.0D0, 0.0D0, 0.0D0,
     *  -1.0D0, 0.0D0, 0.0D0,
     *  -1.0D0, 0.0D0, 0.0D0,
     *  -1.0D0, 0.0D0, 0.0D0/
      DATA ((ELNORM(I,J),J=1,3),I=13,24) /
     *   0.0D0, 0.0D0,-1.0D0,
     *   0.0D0, 0.0D0,-1.0D0,
     *   0.0D0, 0.0D0,-1.0D0,
     *   0.0D0, 0.0D0,-1.0D0,
     *   1.0D0, 0.0D0, 0.0D0,
     *   1.0D0, 0.0D0, 0.0D0,
     *   1.0D0, 0.0D0, 0.0D0,
     *   1.0D0, 0.0D0, 0.0D0,
     *   0.0D0,-1.0D0, 0.0D0,
     *   0.0D0,-1.0D0, 0.0D0,
     *   0.0D0,-1.0D0, 0.0D0,
     *   0.0D0,-1.0D0, 0.0D0/

C Set the number and values of wavenumber required (NOK and KVAL)
      NOK=4
      DATA KVAL/(0.0D0,0.0D0),(1.0D0,0.0D0),
     *          (0.0D0,1.0D0),(1.0D0,1.0D0)/

C Set the parameters of the Burton and Miller Equation corresponding to 
C  the wavenumbers (ALPHA and BETA)
      DATA ALPHA/(1.0D0,0.0D0),(1.0D0,0.0D0),
     *           (1.0D0,0.0D0),(1.0D0,0.0D0)/
      DATA  BETA/(0.0D0,1.0D0),(0.0D0,1.0D0),
     *           (0.0D0,1.0D0),(0.0D0,1.0D0)/

C Set interior point that determine the Neumann boundary condition and
C  exact solution
      PP(1)=0.0D0
      PP(2)=0.0D0
      PP(3)=0.0D0

C Set exterior points, where an approximation the the solution phi is 
C  sought
      NPEXT=1
      PEXT(1,1)=0.0D0
      PEXT(1,2)=0.0D0
      PEXT(1,3)=2.0D0

C Check that the set-up data is satisfactory
      LERROR=.FALSE.
      IF (NQ.GT.LIMNQ) THEN
        WRITE(*,*) 'ERROR - NQ must be less than or equal to LIMNQ'
        LERROR=.TRUE.
      END IF
      IF (N.GT.LIMN) THEN
        WRITE(*,*) 'ERROR - N must be less than or equal to LIMN'
        LERROR=.TRUE.
      END IF
      IF (NOK.GT.LIMK) THEN
        WRITE(*,*) 'ERROR - NOK must be less than or equal to LIMK'
        LERROR=.TRUE.
      END IF
      IF (ABS(PP(1)).GE.ONE.OR.ABS(PP(2)).GE.ONE.
     * OR.ABS(PP(3)).GE.ONE) THEN
        WRITE(*,*) 'ERROR - PP must be interior to the cube'
        LERROR=.TRUE.
      END IF
      IF (NPEXT.GT.LIMNPE) THEN
        WRITE(*,*) 'ERROR - NPEXT must be less than or equal to LIMNPE'
        LERROR=.TRUE.
      END IF
      DO 50 IPE=1,NPEXT
        IF (ABS(PEXT(IPE,1)).LE.ONE.AND.ABS(PEXT(IPE,2)).LE.ONE.
     *   AND.ABS(PEXT(IPE,3)).LE.ONE) THEN
          WRITE(*,*) 'ERROR - PEXT must be exterior to the circle'
          LERROR=.TRUE.
        END IF
50    CONTINUE
      IF (LERROR) STOP

C Output description of the test problem
      WRITE(*,*) 'TEST FOR SUBROUTINE H3LC'
      WRITE(*,*) '========================'
      WRITE(*,*) 'Test problem is that of a cube with sides'
      WRITE(*,*) 'of length 2.0 and centred at the origin with'
      WRITE(*,*) 'a Neumann boundary condition. Solution is via'
      WRITE(*,*) 'the integral equation of Burton and Miller'
      WRITE(*,*) 'Number of elements is',N

C Loop(IK) through each wavenumber
      DO 250 IK=1,NOK

C  Set wavenumber K
        K=KVAL(IK)

C  Output the value of the wavenumber
        WRITE(*,*) 'Wavenumber = ',K

C  Loop(I) through each collocation point
      DO 100 I=1,N

C   Set P(1..3) (the collocation point) and NORMP(1..3) the normal at 
C    P(1..3)
        I1=ELVERT(I,1)
        I2=ELVERT(I,2)
        I3=ELVERT(I,3)
        P(1)=(VERTPT(I1,1)+VERTPT(I2,1)+VERTPT(I3,1))/THREE
        P(2)=(VERTPT(I1,2)+VERTPT(I2,2)+VERTPT(I3,2))/THREE
        P(3)=(VERTPT(I1,3)+VERTPT(I2,3)+VERTPT(I3,3))/THREE
        NORMP(1)=ELNORM(I,1)
        NORMP(2)=ELNORM(I,2)
        NORMP(3)=ELNORM(I,3)

C   Loop(J) through each element
        DO 110 J=1,N

C    Set QA(1..3), QB(1..3) and QC(1..3), the vertices of the element.
C     The vertices are enumerated in anti-clockwise direction (when 
C     viewed from a point near to the element and exterior to the 
C     surface). Hence the normal will be implicitly defined as outward
C     from the boundary in subroutine H3LC.
          QA(1)=VERTPT(ELVERT(J,1),1)
          QA(2)=VERTPT(ELVERT(J,1),2)
          QA(3)=VERTPT(ELVERT(J,1),3)
          QB(1)=VERTPT(ELVERT(J,2),1)
          QB(2)=VERTPT(ELVERT(J,2),2)
          QB(3)=VERTPT(ELVERT(J,2),3)
          QC(1)=VERTPT(ELVERT(J,3),1)
          QC(2)=VERTPT(ELVERT(J,3),2)
          QC(3)=VERTPT(ELVERT(J,3),3)

C    Set LPONEL
          IF (I.EQ.J) THEN
            LPONEL=.TRUE.
          ELSE
            LPONEL=.FALSE.
          END IF

C   Set up quadrature rule data. If LPONEL is false then use the standard
C    Gaussian quadrature rule above. If LPONEL is true then then a
C    quadrature rule with 3 times as many points is used, this is made
C    up from three standard quadrature rules with the quadrature points
C    translated to the three triangles that each have the cetroid and two
C    of the original vertices as its vertices.
          IF (LPONEL) THEN
            NQT=21
            DO 120 IQ=1,7
              XQT(IQ)=XQ(IQ)*THIRD+YQ(IQ)
              YQT(IQ)=XQ(IQ)*THIRD
              WQT(IQ)=WQ(IQ)/THREE
              XQT(IQ+7)=XQ(IQ)*THIRD
              YQT(IQ+7)=XQ(IQ)*THIRD+YQ(IQ)
              WQT(IQ+7)=WQ(IQ)/THREE
              XQT(IQ+14)=THIRD*(ONE+TWO*XQ(IQ)-YQ(IQ))
              YQT(IQ+14)=THIRD*(ONE-XQ(IQ)+TWO*YQ(IQ))
              WQT(IQ+14)=WQ(IQ)/THREE
120         CONTINUE
          ELSE
            NQT=7
            DO 130 IQ=1,7
              XQT(IQ)=XQ(IQ)
              YQT(IQ)=YQ(IQ)
              WQT(IQ)=WQ(IQ)
130         CONTINUE
          END IF

C   For primary stage of method all integral operators are required
C    in general (for the Burton and Miller equation)
          LLK=.TRUE.
          LMK=.TRUE.
          LMKT=.TRUE.
          LNK=.TRUE.

C   Call of H3LC subroutine to compute [Lk], [Mk], [Mkt], [Nk]
          CALL H3LC(K,P,NORMP,QA,QB,QC,LPONEL,
     *     LIMNQ,NQT,XQT,YQT,WQT,
     *     LVALID,EK,EGEOM,EQRULE,LFAIL,
     *     LLK,LMK,LMKT,LNK,DISLK,DISMK,DISMKT,DISNK)

C   Storage of results in matrices
          LK(I,J)=DISLK
          MK(I,J)=DISMK
          MKT(I,J)=DISMKT
          NK(I,J)=DISNK

C   End loop(J) through elements
110     CONTINUE

C  End loop(I) through collocation points
100   CONTINUE

C  Output the position of the source point that prescribes the boundary
C   condition and exact solution
      WRITE(*,*) 'Neumann boundary condition is that produced by'
      WRITE(*,*) 'a point source at ',PP

C  Set up test Neumann boundary condition (VEL) and the exact solution
C   on the boundary (PHI)
      DO 140 I=1,N
        I1=ELVERT(I,1)
        I2=ELVERT(I,2)
        I3=ELVERT(I,3)
        Q(1)=(VERTPT(I1,1)+VERTPT(I2,1)+VERTPT(I3,1))/THREE
        Q(2)=(VERTPT(I1,2)+VERTPT(I2,2)+VERTPT(I3,2))/THREE
        Q(3)=(VERTPT(I1,3)+VERTPT(I2,3)+VERTPT(I3,3))/THREE
        NORMQ(1)=ELNORM(I,1)
        NORMQ(2)=ELNORM(I,2)
        NORMQ(3)=ELNORM(I,3)
        RR(1)=Q(1)-PP(1)
        RR(2)=Q(2)-PP(2)
        RR(3)=Q(3)-PP(3)
        R=SQRT(RR(1)*RR(1)+RR(2)*RR(2)+RR(3)*RR(3))
        DRBDNQ=(RR(1)*NORMQ(1)+RR(2)*NORMQ(2)+RR(3)*NORMQ(3))/R
        KR=K*R
        IKR=CIMAG*KR
        E=EXP(IKR)
        PHI(I)=E/R
        VEL(I)=E*DRBDNQ*(IKR-ONE)/R/R
140   CONTINUE

C  Construct composite matrices AK and BK of the Burton and Miller 
C   integral equation
      DO 150 I=1,N
        DO 160 J=1,N
          AK(I,J)=ALPHA(IK)*MK(I,J)+BETA(IK)*NK(I,J)
          BK(I,J)=ALPHA(IK)*LK(I,J)+BETA(IK)*MKT(I,J)
160     CONTINUE
        AK(I,I)=AK(I,I)-ALPHA(IK)/TWO
        BK(I,I)=BK(I,I)+BETA(IK)/TWO
150   CONTINUE

C  Solve linear system AK APPPHI = BK VEL to give APPPHI, the 
C   approximation to PHI.
      DO 170 I=1,N
        BVEL(I)=ZERO
        DO 180 J=1,N
          BVEL(I)=BVEL(I)+BK(I,J)*VEL(J)
180     CONTINUE
170   CONTINUE

C   Call the subroutine for solving the linear system of equations
      CALL CLINSL(LIMN,N,AK,APPPHI,BVEL)

C  Output of solution on boundary
999   FORMAT(2F12.6,'   ',2F12.6)
      WRITE(*,*) '      EXACT PHI ON S            APPROX PHI ON S'
      DO 190 I=1,N
        WRITE(*,999) DBLE(PHI(I)),DIMAG(PHI(I)),
     *   DBLE(APPPHI(I)),DIMAG(APPPHI(I))
190   CONTINUE


C  Loop (IPE) through the points in the exterior
      DO 300 IPE=1,NPEXT
        P(1)=PEXT(IPE,1)
        P(2)=PEXT(IPE,2)
        P(3)=PEXT(IPE,3)

C   Output the coordinates of the point in the exterior
        WRITE(*,*) 'Point in the exterior is',P

C   Compute the exact solution PHIPT at the exterior point
        RR(1)=P(1)-PP(1)
        RR(2)=P(2)-PP(2)
        RR(3)=P(3)-PP(3)
        R=SQRT(RR(1)*RR(1)+RR(2)*RR(2)+RR(3)*RR(3))
        KR=K*R
        IKR=CIMAG*KR
        E=EXP(IKR)
        PHIPT=E/R

C   Initialise SUM, this ultimately stores the value of the summation
C    that approximates to phi at the exterior point
        SUM=ZERO

C   Loop(J) through the elements
        DO 200 J=1,N

C    Set QA(1..3), QB(1..3) and QC(1..3) the vertices of the element.
          QA(1)=VERTPT(ELVERT(J,1),1)
          QA(2)=VERTPT(ELVERT(J,1),2)
          QA(3)=VERTPT(ELVERT(J,1),3)
          QB(1)=VERTPT(ELVERT(J,2),1)
          QB(2)=VERTPT(ELVERT(J,2),2)
          QB(3)=VERTPT(ELVERT(J,2),3)
          QC(1)=VERTPT(ELVERT(J,3),1)
          QC(2)=VERTPT(ELVERT(J,3),2)
          QC(3)=VERTPT(ELVERT(J,3),3)
       
C    Set up quadrature rule
          NQT=7
          DO 210 IQ=1,7
            WQT(IQ)=WQ(IQ)
            XQT(IQ)=XQ(IQ)
            YQT(IQ)=YQ(IQ)
210       CONTINUE

C    For secondary stage of method Lk and Mk are required.
          LLK=.TRUE.
          LMK=.TRUE.
          LMKT=.FALSE.
          LNK=.FALSE.

C    Since P is in the exterior then LPONEL is always false
          LPONEL=.FALSE.

C    Value of VECP is not relevant

C    Call of H3LC routine to compute [Lk], [Mk]
          CALL H3LC(K,P,VECP,QA,QB,QC,LPONEL,
     *     LIMNQ,NQT,XQT,YQT,WQT,
     *     LVALID,EK,EGEOM,EQRULE,LFAIL,
     *     LLK,LMK,LMKT,LNK,DISLK,DISMK,DISMKT,DISNK)

C    Accumulate SUM
          SUM=SUM+DISMK*APPPHI(J)-DISLK*VEL(J)

C   End loop(J) through elements
200     CONTINUE

C   Ouput the solution at the exterior point
        WRITE(*,*) 'EXACT  PHI AT POINT = ',PHIPT
        WRITE(*,*) 'APPROX PHI AT POINT = ',SUM

C  End loop(IPE) through the exterior points
300     CONTINUE

C End loop(IK) through the wavenumbers
250   CONTINUE

      CLOSE(10)

      END



C Subroutine CLINSL
C =================
C This subroutine calculates the solution to the linear system
C of complex equations Ax=b.
C Input parameters
C integer MAXN : The maximum dimension of the matrix. Its value must not
C  be changed between calls of CLINSL.
C integer N    : The dimension of the matrix.
C complex*16 A(MAXN,MAXN) : The matrix `A'.
C complex*16 B(MAXN) : The vector `b'.

C Output parameters
C complex*16 X(MAXN) : The vector `x'.

        SUBROUTINE CLINSL(MAXN,N,A,X,B)
        INTEGER    MAXN,N
        COMPLEX*16 A(MAXN,MAXN),X(MAXN),B(MAXN)
        COMPLEX*16 TEMP,CSUM,RATIO
        REAL*8     COLMAX,AA
        INTEGER    I,J,II,IIMAX
        IF (MAXN.LT.N) THEN
         WRITE(6,*) 'ERROR(CLINSL) - Must have MAXN>=N'
         STOP        
        ENDIF
        DO 10 I=1,N-1
          COLMAX=0
          DO 20 II=I,N
            AA=ABS(A(II,I))
            IF (AA.GT.COLMAX) THEN
              COLMAX=AA
              IIMAX=II
            ENDIF
20        CONTINUE
          DO 30 J=I,N
            TEMP=A(I,J)
            A(I,J)=A(IIMAX,J)
            A(IIMAX,J)=TEMP
30        CONTINUE
          TEMP=B(I)
          B(I)=B(IIMAX)
          B(IIMAX)=TEMP
          DO 40 II=I+1,N
            RATIO=A(II,I)/A(I,I)
            DO 50 J=I+1,N
              A(II,J)=A(II,J)-RATIO*A(I,J)
50          CONTINUE
          B(II)=B(II)-RATIO*B(I)
40        CONTINUE
10      CONTINUE
        X(N)=B(N)/A(N,N)
        DO 60 I=N-1,1,-1
          CSUM=B(I)
          DO 70 J=I+1,N
            CSUM=CSUM-A(I,J)*X(J)
70        CONTINUE
          X(I)=CSUM/A(I,I)
60      CONTINUE
        END


      REAL*8 FUNCTION FNSQRT(X)
      REAL*8 X
      FNSQRT=SQRT(X)
      END

      COMPLEX*16 FUNCTION FNEXP(Z)
      COMPLEX*16 Z
      FNEXP=EXP(Z)
      END