C*********************************************************************
C      Subroutine GLS by www.boundary-element-method.com             *              
C*********************************************************************
C
C  GLS version 2 (October 2014) by Stephen Kirkup
C  replaces GLS verstion 1 (1998)
C  [ GLS version 1 is a avalable on
C  www.boundary-element-method.com/fortran/GLS.FOR ]
C
C  School of Computing Engineering and Physical Sciences
C  University of Central Lancashire - www.uclan.ac.uk 
C  smkirkup@uclan.ac.uk
C  https://www.researchgate.net/profile/Stephen_Kirkup
C
C  This open source code can be found at
C   www.boundary-element-method.com/fortran/GLS2.FOR 
C  The user guide for this code can be found on 
C   www.bounday-element-method.com/fortran/GLS_FOR.pdf
C  Their are two subroutines that GLS depends on and the 
C   source codes can be downloaded from
C   www.numerical-methods.com/fortran/LUFAC.FOR
C   www.numerical-methods.com/fortran/LUFBSUB.FOR
C
C Licence: This is 'open source'; the software may be used and 
C  applied within other systems as long as its provenance is 
C  appropriately acknowledged. See the GNU Licence 
C  http://www.gnu.org/licenses/lgpl.txt
C
C**********************************************************************
C
C Subroutine GLS returns the solution x,y to a problem of the
C  form  
C                             A x = B y
C
C where A and B are n by n real matrices under the condition(s)
C                         
C             {\alpha}_i x_i + {\beta}_i y_i = f_i   for i=1..n.
C
C Clearly only one of {\alpha}_i or {\beta}_i can be zero for#
C  each i.
C
C The method employed involves forming a linear system of the 
C form Cz=d where the n by n matrix C and the vector d can be 
C determined from A,B and the {\alpha}_i and {\beta}_i. A 
C standard LU factorisation solution method is then employed 
C to return z. From z the actual solutions x,y can be 
C determined.
C
C *********************************************************************


      SUBROUTINE GLS(MAXN,N,A,B,C,ALPHA,BETA,F,X,Y,LFAIL,
     * PERM, XORY, WKSPC1, WKSPC2, WKSPC)

C Input parameters
C ----------------
C The limit on the dimension of the matrices A and B 
      INTEGER    MAXN
C The dimension of the matrices and vectors
      INTEGER    N
C The matrix A
      REAL*8     A(MAXN,MAXN)
C The matrix B
      REAL*8     B(MAXN,MAXN)
C The vector c
      REAL*8     C(MAXN)
C The {\alpha}_i
      REAL*8     ALPHA(MAXN)
C The {\beta}_i
      REAL*8     BETA(MAXN)
C The f_i
      REAL*8     F(MAXN)

C Output parameters
C -----------------
C The solution vector x
      REAL*8     X(MAXN)
C The solution vector y
      REAL*8     Y(MAXN)
C The following parameters are 'output', but only needed if new 
C  'boundary conditions' are to be applied to the system
C The matrix A is overwritten by the LU factorisation and B
C  is altered by row swaps with A
C The permutation matrix (resulting from the LU factorisation)
      INTEGER    PERM(MAXN)
C The record of column exchanges
      LOGICAL    XORY(MAXN)

C Work space
C ----------
      REAL*8     WKSPC1(MAXN)
      REAL*8     WKSPC2(MAXN)
      REAL*8     WKSPC(MAXN)
      


C Local variables
C ---------------
      LOGICAL    LFAIL,LERROR
      REAL*8     ANORM,BNORM
      REAL*8     EPS
      LOGICAL    A_DIAGONAL
      REAL*8     TEMP

C Initialisation
      EPS=1.0D-10
      LFAIL=.FALSE.

C Validation
      IF (MAXN.LT.1.OR.N.LT.1.OR.N.GT.MAXN) THEN
        WRITE(*,*) 'ERROR(GLS) - Check N,MAXN parameters'
        LFAIL=.TRUE.
      END IF

      LERROR=.FALSE.
      DO 10 I=1,N
        IF (ABS(ALPHA(I)).LT.EPS.AND.ABS(BETA(I)).LT.EPS) THEN
          LERROR=.TRUE.
        END IF
10    CONTINUE
      IF (LERROR) THEN
        WRITE(*,*) 'ERROR(GLS) - ALPHA(i) and BETA(i) must not both'
        WRITE(*,*) ' be zero for any value of i in 1..N'
        LFAIL=.TRUE.
      END IF
  
      IF (LFAIL) GOTO 999


C Compute the column (1-)norms of the matrices A and B
C  and the 1-norms of the matrices A and B
C Compute the column (1-)norms of the matrices A and B
C  and the 1-norms of the matrices A and B
      ANORM = 0.0D0
      BNORM = 0.0D0
      DO 100 J=1,N
        WKSPC1(J)=0.0D0
        WKSPC2(J)=0.0D0
        DO 110 i=1,N
          WKSPC1(J) = WKSPC1(J) + ABS(A(I,J))
          WKSPC2(J) = WKSPC2(J) + ABS(B(I,J))
110     CONTINUE
      ANORM=ANORM+WKSPC1(J)
      BNORM=BNORM+WKSPC2(J)
100   CONTINUE     
          


C Check that each alpha_i and beta_i are not both zero. If they are then the
C  subroutine is abandoned
      DO 200 I=1,N
        IF (ABS(ALPHA(I)).LT.(EPS*ANORM)) THEN
          IF ((ABS(BETA(I)).LT.EPS*BNORM)) THEN
            WRITE(*,*) 'Both ALPHA(i) and BETA(i) are zero. '
            WRITE(*,*) 'GLS is abandoned'
            LFAIL=.TRUE.
          END IF
        END IF
200    CONTINUE

       IF (LFAIL) GOTO 999
              
C Set XORY. This states whether the x_i or y_i is initially found
       DO 300 I=1,N
         XORY(I)=.FALSE.
         IF (BETA(I)**2*WKSPC1(I).GT.ALPHA(I)**2*WKSPC2(I)) THEN 
           XORY(I)=.TRUE.
         END IF
300    CONTINUE
  
C The matrices A and B are then over-written
C  First the columns of the matrices are exchanged when xory_j if 'false'
      DO 400 J=1,N
        IF (.NOT.XORY(J)) THEN
          DO 410 I=1,N
            TEMP=A(I,J)
            A(I,J)=-B(I,J)
            B(I,J)=-TEMP
410       CONTINUE
        END IF
400   CONTINUE


C  The 'B' matrix is divided by alpha_j or beta_j, depending on xory_j
      DO 500 J=1,N
        IF (XORY(J)) THEN
          DO 510 I=1,N
            B(I,J)=B(I,J)/BETA(J)
510       CONTINUE
        ELSE
          DO 520 I=1,N
            B(I,J)=B(I,J)/ALPHA(J)
520       CONTINUE
        END IF
500   CONTINUE

C The terms are subtracted from 'A'
      DO 540 J=1,N
        IF (XORY(J)) THEN
          DO 550 I=1,N
            A(I,J)=A(I,J)+ALPHA(J)*B(I,J)
550       CONTINUE
        ELSE
          DO 560 I=1,N
            A(I,J)=A(I,J)+BETA(J)*B(I,J)
560       CONTINUE
        END IF
540    CONTINUE
            
C Prepare the right hand side for the carrying out forward and backward substitution

       DO 600 I=1,N
         DO 610 J=1,N
           C(I)=C(I)+B(I,J)*F(J)
610      CONTINUE
600    CONTINUE

C Carry out an LU factorisation and solve the system
C  Check for special case that 'a' is diagonal
C   Check to see if 'a' is diagonal
       A_DIAGONAL=.TRUE.
       DO 700 I=1,N
         DO 710 J=1,N
           IF(I.NE.J) THEN
             IF (ABS(A(I,J)).GT.EPS*ANORM) THEN 
               A_DIAGONAL=.FALSE.
             END IF
           END IF
710      CONTINUE
700    CONTINUE

C   Solve the system
       IF (A_DIAGONAL) THEN
C   Set perm to its default
         DO 800 I=1,N
           PERM(I)=I
800      CONTINUE
C   Solve the diagonal system
         DO 820 I=1,N
           C(I)=C(I)/A(I,I)
820      CONTINUE       
       ELSE


 
         
C   Carry out an LU factorisation of the matrix a using the LUfac subroutine.
C    A is overwritten by its LU factors
        CALL LUFAC(MAXN,N,A,PERM,LFAIL,WKSPC)
C   Carry out forward and back substitution
C    c is overwritten by the solution (each term dependent on xory)
        CALL LUFBSUB(MAXN,N,A,PERM,C,WKSPC)
      END IF

C Complete the solution
      DO 900 I=1,N
        IF (XORY(I)) THEN
          X(I)=C(I)
          Y(I)=(F(I)-ALPHA(I)*X(I))/BETA(I)
        ELSE
          Y(I)=C(I)
          X(I)=(F(I)-BETA(I)*Y(I))/ALPHA(I)
        END IF
900   CONTINUE

999   CONTINUE

      END

