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