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