C **********************************************************************
C Subroutine CREGLS.FOR by www.boundary-element-method.com *
C **********************************************************************
C
C CREGLS version 1 (June 2015) by Stephen Kirkup
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 Subroutine REGLS returns the solution x,y to a general linear
C system. The code is normally used in conjunction with GLS and
C it computed a further solution, but with information taken from
C a prior run of GLS.
C
C The A,B,alpha and beta have been pre-determined and presented to
C the subroutine GLS.FOR. GLS has carried our row and column
C rearrangements and returned matrices A_GLS, B_GLS, XORY and PERM.
C
C REGLS is then called with 'new' values of c and f (C_NEW and F_NEW)
C and the new solutions are returned in X and Y.
C
C Source of the code:
C http://www.boundary-element-method.com/fortran/CREGLS.FOR
C Source of the user-guide:
C http://www.boundary-element-method.com/fortran/CREGLS_FOR.pdf
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
SUBROUTINE CREGLS(MAXN,N, A_GLS, B_GLS, XORY, PERM, ALPHA, BETA,
* C_NEW, F_NEW, X, Y, 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
INTEGER N
C The matrix A_GLS
COMPLEX*16 A_GLS(MAXN,MAXN)
C The matrix B_GLS
COMPLEX*16 B_GLS(MAXN,MAXN)
C The record of column exchanges
LOGICAL XORY(MAXN)
C The vector PERM
INTEGER PERM(MAXN)
C The {\alpha}_i
COMPLEX*16 ALPHA(MAXN)
C The {\beta}_i
COMPLEX*16 BETA(MAXN)
C The vector C_NEW
COMPLEX*16 C_NEW(MAXN)
C The vector F_NEW
COMPLEX*16 F_NEW(MAXN)
C Output parametrs
C The vector X
COMPLEX*16 X(MAXN)
C The vector Y
COMPLEX*16 Y(MAXN)
C Workspace
COMPLEX*16 WKSPC(MAXN)
C local variables
INTEGER I,J
REAL*8 ANORM
LOGICAL DIAGONAL
REAL*8 TINY
TINY=1D-06
C Carry out forward and back substitution
C Prepare the right hand side for the carrying out forward and backward substitution
DO 200 I=1,N
DO 210 J=1,N
C_NEW(I)=C_NEW(I)+B_GLS(I,J)*F_NEW(J)
210 CONTINUE
200 CONTINUE
C Check for special case that 'A' is diagonal
C First find the norm of 'A'
ANORM=0.0D0
DO 300 I=1,N
DO 310 J=1,N
ANORM=ANORM+ABS(A_GLS(I,J))
310 CONTINUE
300 CONTINUE
C Check to see if 'A' is diagonal
DIAGONAL=.TRUE.
DO 400 I=1,N
DO 410 J=1,N
IF (I.NE.J) THEN
IF (ABS(A_GLS(I,J)).GT.TINY*ANORM) THEN
DIAGONAL=.FALSE.
END IF
ENDIF
410 CONTINUE
400 CONTINUE
C Solve the system
IF (DIAGONAL) THEN
C Solve the diagonal system
DO 500 I=1,N
C_NEW(I)=C_NEW(I)/A_GLS(I,I)
500 CONTINUE
ELSE
C Carry out forward and back substitution
C C is overwritten by the solution (each term dependent on XORY)
CALL CLUFBSUB(MAXN, N, A_GLS, PERM, C_NEW, WKSPC)
END IF
C Complete the solution
DO 600 I=1,N
IF (XORY(I)) THEN
X(I)=C_NEW(I)
Y(I)=(F_NEW(I)-ALPHA(I)*X(I))/BETA(I)
ELSE
Y(I)=C_NEW(I)
X(I)=(F_NEW(I)-BETA(I)*Y(I))/ALPHA(I)
END IF
600 CONTINUE
END