SUBROUTINE GMRES(N, KDMAX, ITMAX, RHS, X, KVAL,
& PRECON, IFLAG, ROWPOSP, COLPOSP)
C
C THIS ROUTINE IS AN EXTENSION OF THE STANDARD RESTARTED GMRES METHOD OF
C Y. SAAD AND M. SCHULTZ, "GMRES: A GENERALIZED MINIMAL RESIDUAL METHOD
C FOR SOLVING NONSYMMETRIC LINEAR SYSTEMS", SIAM J. SCI. STAT. COMPUT., 7
C (1986), PP. 856-869, THAT ALLOWS THE MAXIMUM KRYLOV SUBSPACE DIMENSION
C TO INCREASE (UP TO A MAXIMUM SPECIFIED PARAMETER VALUE) IF STAGNATION
C OCCURS. THE ARNOLDI BASIS VECTORS ARE GENERATED USING ORTHOGONALIZATION
C WITH HOUSEHOLDER TRANSFORMATIONS.
C ON RESTARTS, RESIDUAL VECTORS ARE COMPUTED BY DIRECT EVALUATION.
C CONDITIONING OF THE GMRES LEAST-SQUARES PROBLEM IS MONITORED USING
C LAPACK ROUTINE DLAIC1 (SEE P. N. BROWN AND H. F. WALKER, "GMRES ON
C (NEARLY) SINGULAR SYSTEMS", UTAH STATE UNIV. C MATH. STAT. DEPT. RES.
C REPORT 2/94/73, FEBRUARY, 1994).
C
C
C INPUT VARIABLES:
C
C N IS THE DIMENSION OF THE LINEAR SYSTEM AA*X = RHS. THE SPARSE
C MATRIX DATA STRUCTURE FOR AA IS STORED IN THE MODULE HOMOTOPY.
C
C KDMAX IS THE MAXIMUM KRYLOV SUBSPACE DIMENSION ALLOWED BEFORE
C STAGNATION DETECTION BEGINS.
C
C ITMAX IS THE MAXIMUM OVERALL NUMBER OF GMRES ITERATIONS ALLOWED.
C
C RHS IS THE RIGHT HAND SIDE OF THE LINEAR SYSTEM.
C
C X IS THE INITIAL APPROXIMATE SOLUTION OF THE LINEAR SYSTEM.
C
C KVAL IS AN INDEX INTO X REQUIRED FOR MATRIX-VECTOR MULTIPLICATION.
C
C PRECON IS A ONE-DIMENSIONAL REAL ARRAY CONTAINING A PRECONDITIONING
C MATRIX Q FOR AA.
C
C IFLAG INDICATES WHAT DATA IS USED IN MATRIX-VECTOR MULTIPLICATION.
C
C ROWPOSP IS AN OPTIONAL INTEGER ARRAY USED IN THE DATA STRUCTURE
C DESCRIBING THE SPARSE PRECONDITIONING MATRIX Q.
C
C COLPOSP IS AN OPTIONAL INTEGER ARRAY USED IN THE DATA STRUCTURE
C DESCRIBING THE SPARSE PRECONDITIONING MATRIX Q.
C
C OUTPUT VARIABLES:
C
C KDMAX IS THE INCREASED MAXIMUM KRYLOV SUBSPACE DIMENSION ON
C TERMINATION.
C
C IFLAG IS UNCHANGED ON NORMAL TERMINATION (DESIRED TOLERANCE MET);
C = 4, IF GMRES ITERATION FAILS TO CONVERGE BECAUSE
C - DESIRED TOLERANCE NOT MET AFTER ITMAX ITERATIONS;
C - DANGEROUS NEAR-SINGULARITY DETECTED;
C - LIMITS OF NUMERICAL ACCURACY REACHED;
C - STAGNATION OCCURS.
C
C ITMAX IS THE TOTAL NUMBER OF GMRES ITERATIONS PERFORMED.
C
C X IS THE FINAL APPROXIMATE SOLUTION.
C
C
C CALLS DNRM2, ILUSOLVDS, MULTDS, MULT2DS, SOLVDS, AND INTERNAL
C SUBROUTINES MULM1 AND MULM2 AND INTERNAL FUNCTIONS APPL_HOUSE, HOUSE.
C
USE HOMOTOPY, AA=>QRSPARSE
USE REAL_PRECISION
INTEGER, INTENT(IN):: KVAL, N
INTEGER, INTENT(IN OUT):: IFLAG, ITMAX, KDMAX
REAL (KIND=R8), INTENT(IN):: PRECON(:), RHS(N)
REAL (KIND=R8), INTENT(IN OUT):: X(N)
INTEGER, INTENT(IN), OPTIONAL:: COLPOSP(:), ROWPOSP(:)
C
C IPRINT IS A PARAMETER FOR THE FORTRAN UNIT NUMBER, WHICH IF POSITIVE
C ACTIVATES TRACING OF THE GMRES ITERATIONS.
C M IS A PARAMETER USED IN THE ADAPTIVE STRATEGY FOR INCREMENTING THE
C KRYLOV SUBSPACE DIMENSION (VALUE KDMAX).
C STRONG_VERSION IS A LOGICAL PARAMETER, WHICH IF TRUE CAUSES GMRES TO
C STOP IN TWO CASES: WHEN THE ESTIMATED CONDITION NUMBER IS GREATER
C THAN 1.0/(50*EPSILON);
C WHEN THE TRUE RESIDUAL NORM IS LARGE AND IT INCREASED DURING THE LAST
C RESTART CYCLE.
C
INTEGER, PARAMETER:: IPRINT=0, M=4
LOGICAL, PARAMETER:: STRONG_VERSION=.TRUE.
C
C LOCAL VARIABLES.
C
INTEGER:: I, IFLAGC, IFLAGI, IQUIT, ITNO, KD,
& KDP1, KDLIMIT, LENAA
REAL (KIND=R8), ALLOCATABLE:: C(:), R(:,:), S(:), SVBIG(:),
& SVSML(:), V(:,:), W(:)
REAL (KIND=R8):: VTEMP(N), VTEMP2(N)
REAL (KIND=R8):: BIG, BIGCND, CC, CNDMAX,
& RSN, RSNOLD, SESTPR, SMALL, SS, TEMP, TOL
INTERFACE
FUNCTION DNRM2(N,X,STRIDE)
USE REAL_PRECISION
INTEGER:: N,STRIDE
REAL (KIND=R8):: DNRM2,X(N)
END FUNCTION DNRM2
END INTERFACE
C
C DEFINE MAXIMUM ALLOWED KRYLOV SUBSPACE DIMENSION.
C
KDLIMIT = MIN(MAX(200, INT(SQRT(REAL(N)))), (N+1)/2)
KDMAX = MIN(KDMAX,KDLIMIT)
C
C ALLOCATE LOCAL ARRAYS.
C
IF(.NOT. ALLOCATED(C)) ALLOCATE(C(KDLIMIT))
IF(.NOT. ALLOCATED(R)) ALLOCATE(R(KDLIMIT,KDLIMIT))
IF(.NOT. ALLOCATED(S)) ALLOCATE(S(KDLIMIT))
IF(.NOT. ALLOCATED(SVBIG)) ALLOCATE(SVBIG(KDLIMIT))
IF(.NOT. ALLOCATED(SVSML)) ALLOCATE(SVSML(KDLIMIT))
IF(.NOT. ALLOCATED(V)) ALLOCATE(V(N,KDLIMIT+1))
IF(.NOT. ALLOCATED(W)) ALLOCATE(W(KDLIMIT+1))
C
C PERFORM INITIALIZATIONS.
C
ITNO = 0
IFLAGC = 0
IFLAGI = 1
LENAA = ROWPOS(N)-1
TOL = MAX(100.0, 1.01*REAL(LENAA)/REAL(N))*EPSILON(1.0_R8)
VTEMP = X
IF (PRESENT(ROWPOSP) .AND. PRESENT( COLPOSP)) THEN
CALL MULM2(V(:,1),VTEMP) ! MODE=2
ELSE
CALL MULM1(V(:,1),VTEMP) ! MODE=1
END IF
V(:,1) = RHS - V(:,1)
RSN = DNRM2(N, V(:,1), 1)
TEMP = DNRM2(N, RHS, 1)
TOL = MAX(RSN,TEMP)*TOL
IF (RSN .LE. TOL) THEN
ITMAX = ITNO
CALL CLEANUP
RETURN
ELSE
IF (ITMAX .EQ. 0) THEN
ITMAX = ITNO
IFLAG = 4
CALL CLEANUP
RETURN
ENDIF
ENDIF
CNDMAX = 1.0/(50.0*EPSILON(1.0_R8))
IQUIT = 0
BIGCND = 0.0
C
C FOR PRINTING:
C
IF (IPRINT .GT. 0) THEN
WRITE(IPRINT, 19) ITNO, RSN
19 FORMAT(//,9X,'GMRES ITERATIONS:',//,T4,'ITERATION',TR4,
& 'RESIDUAL NORM',TR5,'CONDITION NUMBER',/,T4,I6,TR7,ES16.8)
ENDIF
C
C BEGIN THE OUTER LOOP.
C
OUTER: DO
C
C FOR PRINTING:
C
IF (IPRINT .GT. 0 ) THEN
IF (IFLAGI .NE. 0) THEN
WRITE (IPRINT, 21) KDMAX
21 FORMAT(/' RESTART WITH KRYLOV SUBSPACE DIMENSION =',I2)
ELSE
WRITE(IPRINT, 22)
22 FORMAT(/' RESTART')
ENDIF
ENDIF
C
KD = 0
RSNOLD = RSN
IQUIT = 0
IFLAGI = 0
C
C FIND HOUSEHOLDER VECTOR V(:,1).
C
W(1)=V(1,1)
CALL HREFG(N,V(2:N,1),V(1,1),W(1))
C
C BEGIN THE INNER LOOP.
C
200 INNER: DO
KD = KD + 1
KDP1 = KD + 1
ITNO = ITNO + 1
IQUIT = 0
C
C FIND THE KD-TH ORTHOGONAL VECTOR.
C
VTEMP = 0.0
VTEMP(KD) = 1.0
DO I=KD,1,-1
TEMP=V(I,I)
V(I,I)=1.0
CALL HREFX(N-I+1,V(I:N,I),TEMP,VTEMP(I:N))
V(I,I)=TEMP
END DO
C
C EVALUATE AA*(KD-TH KRYLOV SUBSPACE BASIS VECTOR).
C
IF (PRESENT(ROWPOSP) .AND. PRESENT(COLPOSP)) THEN ! MODE = 2
CALL ILUSOLVDS(N, PRECON(1:ROWPOSP(N+1)-1),
& ROWPOSP(N+1)-1, ROWPOSP(1:N+1),
& COLPOSP(1:ROWPOSP(N+1)-1), VTEMP(1:N))
CALL MULM2(VTEMP2, VTEMP(1:N))
ELSE ! MODE = 1
CALL SOLVDS(N, PRECON, ROWPOS(N+1)-1, ROWPOS(1:N+1),
& VTEMP(1:N))
CALL MULM1(VTEMP2, VTEMP(1:N))
ENDIF
C
C FIND HOUSEHOLDER VECTOR V(:,KDP1).
C
VTEMP(1:N)=VTEMP2(1:N)
DO I=1,KD
TEMP=V(I,I)
V(I,I)=1.0
CALL HREFX(N-I+1,V(I:N,I),TEMP,VTEMP(I:N))
V(I,I)=TEMP
END DO
IF ( MAXVAL(ABS(VTEMP(KDP1:N))) .NE. 0.0) THEN
V(KDP1+1:N, KDP1) = VTEMP(KDP1+1:N)
CALL HREFG(N-KDP1+1,V(KDP1+1:N, KDP1),V(KDP1,KDP1),
& VTEMP(KDP1))
ENDIF
C
C IF KD .GT. 1, APPLY THE PREVIOUS GIVENS ROTATIONS.
C
DO I = 1, KD-1
TEMP = VTEMP(I)
VTEMP(I) = C(I)*TEMP + S(I)*VTEMP(I+1)
VTEMP(I+1) = -S(I)*TEMP + C(I)*VTEMP(I+1)
END DO
C
C DETERMINE AND APPLY THE NEXT ROTATION.
C
IF (VTEMP(KDP1) .NE. 0.0) THEN
TEMP = VTEMP(KD)
R(KD,KD) = SQRT(VTEMP(KD)**2 + VTEMP(KDP1)**2)
C(KD) = TEMP/R(KD,KD)
S(KD) = VTEMP(KDP1)/R(KD,KD)
VTEMP(KD) = C(KD)*TEMP + S(KD)*VTEMP(KDP1)
END IF
R(1:KD, KD) = VTEMP(1:KD)
C
C COMPUTE AND TEST INCREMENTAL CONDITION NUMBER.
C
IF (KD .EQ. 1) THEN
BIG = R(KD,KD)
SMALL = BIG
SVBIG(1) = 1.0
SVSML(1) = 1.0
ELSE
I = 1
CALL DLAIC1(I, KD-1, SVBIG, BIG, R(1,KD), R(KD,KD),
& SESTPR, SS, CC)
BIG = SESTPR
SVBIG(1:KD-1) = SS*SVBIG(1:KD-1)
SVBIG(KD) = CC
I = 2
CALL DLAIC1(I, KD-1, SVSML, SMALL, R(1,KD), R(KD,KD),
& SESTPR, SS, CC)
SMALL = SESTPR
SVSML(1:KD-1) = SS*SVSML(1:KD-1)
SVSML(KD) = CC
ENDIF
IF (STRONG_VERSION) THEN
IF (BIG .GE. SMALL*CNDMAX) THEN
IF (IPRINT .GT. 0) THEN
WRITE (IPRINT, 230) CNDMAX
230 FORMAT(/,4X, 'IN GMRES CONDITION NUMBER .GE.',
& ES16.8)
ENDIF
ITNO = ITNO - 1
IFLAGC = 4
EXIT OUTER
ENDIF
ENDIF
TEMP = BIG/SMALL
IF (TEMP .GT. BIGCND) BIGCND = TEMP
C
C UPDATE W AND THE RESIDUAL NORM.
C
IF (VTEMP(KDP1) .NE. 0.0) THEN
W(KDP1)= -S(KD)*W(KD)
W(KD) = C(KD)*W(KD)
ELSE
W(KDP1) = 0.0
ENDIF
RSN = ABS(W(KDP1))
C
C FOR PRINTING:
C
IF (IPRINT .GT. 0 ) THEN
WRITE (IPRINT, 240) ITNO, RSN, BIG/SMALL
240 FORMAT(T4,I6, TR7,2ES16.8)
ENDIF
C
C TEST FOR TERMINATION OF THE INNER LOOP.
C
IF (RSN .LE. TOL .OR. KD .EQ. KDMAX .OR. ITNO .EQ. ITMAX)
& EXIT INNER
END DO INNER
C
C IF TERMINATING THE INNER LOOP, TEST FOR TERMINATION OF THE OUTER LOOP,
C COMPUTE THE CORRECTION, AND UPDATE THE APPROXIMATE SOLUTION.
C
C TEST FOR TERMINATION OF THE OUTER LOOP.
C
IF (RSN .LE. TOL .OR. ITNO .EQ. ITMAX ) THEN
IQUIT = 1
C
C TEST FOR STAGNATION.
C
ELSE
TEMP = KD*LOG(TOL/RSN)/
& LOG(RSN/((1.0 + 10.0*EPSILON(1.0_R8))*RSNOLD))
IF (TEMP .GE. 40.0*(ITMAX - ITNO)) THEN
IF (KDMAX .LE. KDLIMIT-M) THEN
IFLAGI = 3
KDMAX = KDMAX + M ! INCREASE KDMAX BY M
GO TO 200
ELSE IF (KDMAX .NE. KDLIMIT) THEN
IFLAGI = 3
KDMAX = KDLIMIT
GO TO 200
END IF
ENDIF
ENDIF
C
C COMPUTE THE CORRECTION IN V(:,1).
C
DO I = KD, 1, -1
W(I) = W(I)/R(I,I)
IF (I .GT. 1) W(1:I-1) = W(1:I-1)-W(I)*R(1:I-1,I)
END DO
C
C COMPUTE KD ORTHOGONAL VECTORS FROM HOUSEHOLDER VECTORS.
C
VTEMP = 0.0
VTEMP(1:KD)=W(1:KD)
DO I=KD,1,-1
TEMP=V(I,I)
V(I,I)=1.0
CALL HREFX(N-I+1,V(I:N,I), TEMP, VTEMP(I:N))
V(I,I)=TEMP
END DO
C
C ADD THE CORRECTION TO THE APPROXIMATE SOLUTION.
C
IF (PRESENT(ROWPOSP) .AND. PRESENT(COLPOSP)) THEN ! MODE = 2
CALL ILUSOLVDS(N, PRECON(1:ROWPOSP(N+1)-1),
& ROWPOSP(N+1)-1, ROWPOSP(1:N+1),
& COLPOSP(1:ROWPOSP(N+1)-1), VTEMP(1:N))
ELSE ! MODE = 1
CALL SOLVDS(N, PRECON, ROWPOS(N+1)-1, ROWPOS(1:N+1),
& VTEMP(1:N))
ENDIF
X = X+VTEMP
C
C UPDATE THE RESIDUAL VECTOR BY DIRECT EVALUATION IN (V:,1) AND RECOMPUTE
C THE RESIDUAL NORM.
C
IF (IQUIT .EQ. 0 ) THEN
VTEMP=X
IF (PRESENT(ROWPOSP) .AND. PRESENT( COLPOSP)) THEN
CALL MULM2(V(:,1),VTEMP) ! MODE=2
ELSE
CALL MULM1(V(:,1),VTEMP) ! MODE=1
END IF
V(:,1) = RHS - V(:,1)
RSN = DNRM2(N, V(:,1), 1)
ENDIF
C
C TERMINATE, OR GO TO THE TOP OF THE OUTER LOOP.
C
IF (RSN .LE. TOL) THEN
IFLAGC = 0
EXIT OUTER
ENDIF
IF (ITNO .EQ. ITMAX) THEN
IF (IPRINT .GT. 0) THEN
WRITE (IPRINT, 250) ITMAX
250 FORMAT(/,4X, 'MAXIMUM NUMBER OF ITERATIONS REACHED:',I7)
ENDIF
IFLAGC = 4
EXIT OUTER
END IF
IF (IQUIT .EQ. 0 ) THEN
IF (RSN .GT. RSNOLD ) THEN
IF (IPRINT .GT. 0) THEN
WRITE (IPRINT, 260) RSN
260 FORMAT(/,4X, '(TRUE) RESIDUAL NORM INCREASED TO',ES16.8)
ENDIF
IF (RSN .LE. TOL**(2.0/3.0)) THEN
IFLAGC = 0
ELSE IF (STRONG_VERSION) THEN
IF (IPRINT .GT. 0) THEN
WRITE (IPRINT, 270) TOL**(2.0/3.0)
270 FORMAT(/,4X,'(TRUE) RESIDUAL NORM IS LARGER THAN',ES16.8)
ENDIF
IFLAGC = 4
END IF
EXIT OUTER
END IF
C
C TEST FOR STAGNATION USING TRUE RESIDUAL NORM.
C
TEMP = KD*LOG(TOL/RSN)/
& LOG(RSN/((1.0 + 10.0*EPSILON(1.0_R8))*RSNOLD))
IF (TEMP .GE. 70.0*(ITMAX - ITNO)) THEN
IFLAGI = 3
IF ( KDMAX .LE. KDLIMIT-M ) THEN
KDMAX = KDMAX + M ! INCREASE KDMAX
ELSE IF (KDMAX .NE. KDLIMIT) THEN
KDMAX = KDLIMIT
END IF
ELSE IF (TEMP .GE. 1000.0*(ITMAX - ITNO)) THEN
IFLAGC = 4
EXIT OUTER
END IF
END IF
END DO OUTER
C
C RETURN.
C
ITMAX = ITNO
IF (IFLAGC .NE. 0) IFLAG = IFLAGC
CALL CLEANUP
RETURN
C
C INTERNAL SUBROUTINES.
C
CONTAINS
SUBROUTINE MULM1(Y,X)
C
C MATRIX-VECTOR MULTIPLY FOR MODE = 1.
C
REAL (KIND=R8), INTENT(IN):: X(N)
REAL (KIND=R8), INTENT(OUT):: Y(N)
Y = 0.0
CALL MULTDS(Y(1:N-1), AA, X(1:N-1), ROWPOS(1:N),
& N-1, ROWPOS(N)-1)
C
C RESULT MODIFIED ACCORDING TO KVAL.
C
Y(KVAL) = Y(KVAL)+X(N)
IF (KVAL .LT. N)
& Y(N) = X(KVAL) + X(N)*(1.0/ABS(AA(ROWPOS(KVAL)))+1.0)
RETURN
END SUBROUTINE MULM1
C
SUBROUTINE MULM2(Y,X)
C
C MATRIX-VECTOR MULTIPLY FOR MODE = 2.
C
REAL (KIND=R8), INTENT(IN):: X(N)
REAL (KIND=R8), INTENT(OUT):: Y(N)
INTERFACE
SUBROUTINE MULT2DS(Y, B, X, ROWPOS, COLPOS, N, LENB)
USE REAL_PRECISION
INTEGER, INTENT (IN):: LENB, N, ROWPOS(N+1),
& COLPOS(LENB)
REAL (KIND=R8), INTENT(IN):: X(:), B(LENB)
REAL (KIND=R8), INTENT (OUT):: Y(N)
END SUBROUTINE MULT2DS
END INTERFACE
Y = 0.0
IF (IFLAG .NE. -2) THEN
CALL MULT2DS(Y(1:N-1), AA, X(1:N-1), ROWPOS(1:N),
& COLPOS(1:LENAA), N-1, LENAA)
Y(1:N-1) = Y(1:N-1)-X(N)*PP(1:N-1)
ELSE
CALL MULT2DS(Y(1:N-1), AA, X(1:N), ROWPOS(1:N),
& COLPOS(1:LENAA), N-1, LENAA)
END IF
C
C RESULT MODIFIED ACCORDING TO KVAL.
C
Y(N) = X(KVAL)
RETURN
END SUBROUTINE MULM2
SUBROUTINE CLEANUP ! DEALLOCATE TEMPORARY LOCAL STORAGE.
IF(ALLOCATED(C)) DEALLOCATE(C)
IF(ALLOCATED(R)) DEALLOCATE(R)
IF(ALLOCATED(S)) DEALLOCATE(S)
IF(ALLOCATED(SVBIG)) DEALLOCATE(SVBIG)
IF(ALLOCATED(SVSML)) DEALLOCATE(SVSML)
IF(ALLOCATED(V)) DEALLOCATE(V)
IF(ALLOCATED(W)) DEALLOCATE(W)
END SUBROUTINE CLEANUP
C
SUBROUTINE HREFG(N, X, TAU, ALPHA)
C
C HOUSEHOLDER VECTOR CONSTRUCTION. HREFG IS THE LAPACK AUXILIARY
C ROUTINE DLARFG WITH TRIVIAL MODIFICATIONS FOR COMPATIBILITY WITH
C HOMPACK90.
C
INTEGER, INTENT(IN):: N
REAL (KIND=R8), INTENT(OUT):: ALPHA, TAU
REAL (KIND=R8), INTENT(IN OUT):: X(N-1)
!
! Purpose
! =======
!
! HREFG generates a real elementary reflector H of order N, such
! that
!
! H * ( ALPHA ) = ( BETA ), H' * H = I.
! ( X ) ( 0 )
!
! where ALPHA and BETA are scalars, and X is an (N-1)-element real
! vector. H is represented in the form
!
! H = I - TAU * ( 1 ) * ( 1 V' ) ,
! ( V )
!
! where TAU is a real scalar and V is a real (N-1)-element
! vector.
!
! If the elements of X are all zero, then TAU = 0 and H is taken to be
! the unit matrix.
!
! Otherwise 1 <= TAU <= 2.
!
! Arguments
! =========
!
! N (input)
! The order of the elementary reflector.
!
!
! X (input/output)
! On entry, the vector X.
! On exit, it is overwritten with the vector V.
!
! TAU (output)
! The value TAU.
!
! ALPHA (INPUT/OUTPUT)
! On entry, the value ALPHA.
! On exit, it is overwritten with the value BETA.
!
! =====================================================================
!
! .. Parameters ..
REAL (KIND=R8), PARAMETER:: ONE = 1.0_R8, ZERO = 0.0_R8
! .. Local Scalars ..
INTEGER:: J, KNT
REAL (KIND=R8):: BETA, RSAFMN, SAFMIN, XNORM
!
! .. External Functions ..
INTERFACE
FUNCTION DLAMCH(CMACH)
USE REAL_PRECISION
CHARACTER (LEN=1):: CMACH
REAL (KIND=R8):: DLAMCH
END FUNCTION DLAMCH
FUNCTION DLAPY2(X,Y)
USE REAL_PRECISION
REAL (KIND=R8):: DLAPY2,X,Y
END FUNCTION DLAPY2
FUNCTION DNRM2(N,X,STRIDE)
USE REAL_PRECISION
INTEGER:: N,STRIDE
REAL (KIND=R8):: DNRM2,X(N)
END FUNCTION DNRM2
END INTERFACE
!
! .. Executable Statements ..
!
IF( N.LE.1 ) THEN
TAU = ZERO
RETURN
END IF
XNORM = DNRM2( N-1, X, 1)
IF( XNORM.EQ.ZERO ) THEN
!
! H = I
!
TAU = ZERO
ELSE
!
! General case.
!
BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
SAFMIN = DLAMCH( 'S' )
IF( ABS( BETA ).LT.SAFMIN ) THEN
!
! XNORM, BETA may be inaccurate; scale X and recompute them.
!
RSAFMN = ONE / SAFMIN
KNT = 0
DO WHILE (ABS(BETA) < SAFMIN)
KNT = KNT + 1
X=RSAFMN*X
BETA = BETA*RSAFMN
ALPHA = ALPHA*RSAFMN
END DO
!
! New BETA is at most 1, at least SAFMIN.
!
XNORM = DNRM2( N-1, X, 1 )
BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
TAU = ( BETA-ALPHA ) / BETA
X = X/( ALPHA-BETA )
!
! If ALPHA is subnormal, it may lose relative accuracy.
!
ALPHA = BETA
DO J = 1, KNT
ALPHA = ALPHA*SAFMIN
END DO
ELSE
TAU = ( BETA-ALPHA) / BETA
X = X/( ALPHA-BETA )
ALPHA = BETA
END IF
END IF
RETURN
END SUBROUTINE HREFG
C
SUBROUTINE HREFX(M, V, TAU, C)
C
C HOUSEHOLDER VECTOR APPLICATION. HREFX IS A FORTRAN 90 VERSION OF
C THE LAPACK ROUTINE DLARFX WITH MODIFICATIONS FOR
C COMPATIBILITY WITH HOMPACK 90.
C
INTEGER, INTENT(IN):: M
REAL (KIND=R8), INTENT(IN):: TAU, V(M)
REAL (KIND=R8), INTENT(IN OUT):: C(M)
!
! Purpose
! =======
!
! HREFX applies a real elementary reflector H to a real M-element
! vector C. H is represented in the form
!
! H = I - TAU * V * V'
!
! where TAU is a real scalar and V is a real vector.
!
! If TAU = 0, then H is taken to be the unit matrix.
!
! This version uses inline code if H has order < 11.
!
! Arguments
! =========
!
! M (input)
! The order of the vector C.
!
! V (input)
! The vector V in the representation of H.
!
! TAU (input)
! The value TAU in the representation of H.
!
! C (input/output)
! On entry, the M-element vector C.
! On exit, C is overwritten by H * C.
!
! =====================================================================
!
! .. Parameters ..
REAL (KIND=R8), PARAMETER:: ONE=1.0_R8, ZERO=0.0_R8
!
! .. Local Scalars ..
REAL (KIND=R8):: SUM, T1, T10, T2, T3, T4, T5, T6, T7, T8, T9,
& V1, V10, V2, V3, V4, V5, V6, V7, V8, V9
!
! .. Executable Statements ..
!
IF( TAU.EQ.ZERO ) RETURN
!
! Form H * C, where H has order M.
!
SELECT CASE(M)
!
! Code for general M.
!
CASE (11:)
!
! C := C - TAU * V * V' * C
!
C = C - (TAU * DOT_PRODUCT(V,C)) * V
!
! Special code for 1 x 1 Householder.
!
CASE(1)
T1 = ONE - TAU*V( 1 )*V( 1 )
C(1) = T1*C(1)
!
! Special code for 2 x 2 Householder.
!
CASE(2)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
SUM = V1*C(1) + V2*C(2)
C(1) = C(1) - SUM*T1
C(2) = C(2) - SUM*T2
!
! Special code for 3 x 3 Householder.
!
CASE(3)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
SUM = V1*C(1) + V2*C(2) + V3*C(3)
C(1) = C(1) - SUM*T1
C(2) = C(2) - SUM*T2
C(3) = C(3) - SUM*T3
!
! Special code for 4 x 4 Householder.
!
CASE(4)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
V4 = V( 4 )
T4 = TAU*V4
SUM = V1*C(1) + V2*C(2) + V3*C(3) + V4*C(4)
C( 1 ) = C( 1 ) - SUM*T1
C( 2 ) = C( 2 ) - SUM*T2
C( 3 ) = C( 3 ) - SUM*T3
C( 4 ) = C( 4 ) - SUM*T4
!
! Special code for 5 x 5 Householder.
!
CASE(5)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
V4 = V( 4 )
T4 = TAU*V4
V5 = V( 5 )
T5 = TAU*V5
SUM = V1*C( 1 ) + V2*C( 2 ) + V3*C( 3 ) +
& V4*C( 4 ) + V5*C( 5 )
C( 1 ) = C( 1 ) - SUM*T1
C( 2 ) = C( 2 ) - SUM*T2
C( 3 ) = C( 3 ) - SUM*T3
C( 4 ) = C( 4 ) - SUM*T4
C( 5 ) = C( 5 ) - SUM*T5
!
! Special code for 6 x 6 Householder.
!
CASE(6)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
V4 = V( 4 )
T4 = TAU*V4
V5 = V( 5 )
T5 = TAU*V5
V6 = V( 6 )
T6 = TAU*V6
SUM = V1*C( 1 ) + V2*C( 2 ) + V3*C( 3 ) +
& V4*C( 4 ) + V5*C( 5 ) + V6*C( 6 )
C( 1 ) = C( 1 ) - SUM*T1
C( 2 ) = C( 2 ) - SUM*T2
C( 3 ) = C( 3 ) - SUM*T3
C( 4 ) = C( 4 ) - SUM*T4
C( 5 ) = C( 5 ) - SUM*T5
C( 6 ) = C( 6 ) - SUM*T6
!
! Special code for 7 x 7 Householder.
!
CASE(7)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
V4 = V( 4 )
T4 = TAU*V4
V5 = V( 5 )
T5 = TAU*V5
V6 = V( 6 )
T6 = TAU*V6
V7 = V( 7 )
T7 = TAU*V7
SUM = V1*C( 1 ) + V2*C( 2 ) + V3*C( 3 ) +
& V4*C( 4 ) + V5*C( 5 ) + V6*C( 6 ) +
& V7*C( 7 )
C( 1 ) = C( 1 ) - SUM*T1
C( 2 ) = C( 2 ) - SUM*T2
C( 3 ) = C( 3 ) - SUM*T3
C( 4 ) = C( 4 ) - SUM*T4
C( 5 ) = C( 5 ) - SUM*T5
C( 6 ) = C( 6 ) - SUM*T6
C( 7 ) = C( 7 ) - SUM*T7
!
! Special code for 8 x 8 Householder.
!
CASE(8)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
V4 = V( 4 )
T4 = TAU*V4
V5 = V( 5 )
T5 = TAU*V5
V6 = V( 6 )
T6 = TAU*V6
V7 = V( 7 )
T7 = TAU*V7
V8 = V( 8 )
T8 = TAU*V8
SUM = V1*C( 1 ) + V2*C( 2 ) + V3*C( 3 ) +
& V4*C( 4 ) + V5*C( 5 ) + V6*C( 6 ) +
& V7*C( 7 ) + V8*C( 8 )
C( 1 ) = C( 1 ) - SUM*T1
C( 2 ) = C( 2 ) - SUM*T2
C( 3 ) = C( 3 ) - SUM*T3
C( 4 ) = C( 4 ) - SUM*T4
C( 5 ) = C( 5 ) - SUM*T5
C( 6 ) = C( 6 ) - SUM*T6
C( 7 ) = C( 7 ) - SUM*T7
C( 8 ) = C( 8 ) - SUM*T8
!
! Special code for 9 x 9 Householder.
!
CASE(9)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
V4 = V( 4 )
T4 = TAU*V4
V5 = V( 5 )
T5 = TAU*V5
V6 = V( 6 )
T6 = TAU*V6
V7 = V( 7 )
T7 = TAU*V7
V8 = V( 8 )
T8 = TAU*V8
V9 = V( 9 )
T9 = TAU*V9
SUM = V1*C( 1 ) + V2*C( 2 ) + V3*C( 3 ) +
& V4*C( 4 ) + V5*C( 5 ) + V6*C( 6 ) +
& V7*C( 7 ) + V8*C( 8 ) + V9*C( 9 )
C( 1 ) = C( 1 ) - SUM*T1
C( 2 ) = C( 2 ) - SUM*T2
C( 3 ) = C( 3 ) - SUM*T3
C( 4 ) = C( 4 ) - SUM*T4
C( 5 ) = C( 5 ) - SUM*T5
C( 6 ) = C( 6 ) - SUM*T6
C( 7 ) = C( 7 ) - SUM*T7
C( 8 ) = C( 8 ) - SUM*T8
C( 9 ) = C( 9 ) - SUM*T9
!
! Special code for 10 x 10 Householder.
!
CASE (10)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
V4 = V( 4 )
T4 = TAU*V4
V5 = V( 5 )
T5 = TAU*V5
V6 = V( 6 )
T6 = TAU*V6
V7 = V( 7 )
T7 = TAU*V7
V8 = V( 8 )
T8 = TAU*V8
V9 = V( 9 )
T9 = TAU*V9
V10 = V( 10 )
T10 = TAU*V10
SUM = V1*C( 1 ) + V2*C( 2 ) + V3*C( 3 ) +
& V4*C( 4 ) + V5*C( 5 ) + V6*C( 6 ) +
& V7*C( 7 ) + V8*C( 8 ) + V9*C( 9 ) +
& V10*C( 10 )
C( 1 ) = C( 1 ) - SUM*T1
C( 2 ) = C( 2 ) - SUM*T2
C( 3 ) = C( 3 ) - SUM*T3
C( 4 ) = C( 4 ) - SUM*T4
C( 5 ) = C( 5 ) - SUM*T5
C( 6 ) = C( 6 ) - SUM*T6
C( 7 ) = C( 7 ) - SUM*T7
C( 8 ) = C( 8 ) - SUM*T8
C( 9 ) = C( 9 ) - SUM*T9
C( 10 ) = C( 10 ) - SUM*T10
END SELECT
RETURN
END SUBROUTINE HREFX
END SUBROUTINE GMRES