C**********************************************************************
C
C     Copyright (C) 1992  Roland W. Freund and Noel M. Nachtigal
C     All rights reserved.
C
C     This code is part of a copyrighted package.  For details, see the
C     file "cpyrit.doc" in the top-level directory.
C
C     *****************************************************************
C     ANY USE OF  THIS CODE CONSTITUTES ACCEPTANCE OF  THE TERMS OF THE
C                             COPYRIGHT NOTICE
C     *****************************************************************
C
C**********************************************************************
C
C     This file contains the routine for the QMR-from-BCG algorithm for
C     symmetric matrices.
C
C**********************************************************************
C
      SUBROUTINE DSQBG (NDIM,NLEN,NLIM,VECS,TOL,INFO)
C
C     Purpose:
C     This subroutine  uses the QMR-from-BCG algorithm  to solve linear
C     systems.  It runs the  algorithm to convergence  or until a user-
C     specified  iteration  limit  is reached.  It is  set up  to solve
C     symmetric systems starting with identical starting vectors.
C
C     The  code is  set up  to solve  the system  A x = b  with initial
C     guess x_0 = 0.  Here  A x = b  denotes the preconditioned system,
C     and  it  is  connected with the  original system as follows.  Let
C     B y = c be the original unpreconditioned system to be solved, and
C     let y_0 be an arbitrary initial guess for its solution.  Then:
C          A x = b, where  A = M_1^{-1} B M_2^{-1},
C                          x = M_2 (y - y_0), b = M_1^{-1} (c - B y_0).
C     Here M = M_1 M_2 is the preconditioner.
C
C     To recover the final iterate y_n for  the original system B y = c
C     from the final iterate x_n for the preconditioned system A x = b,
C     set
C               y_n = y_0 + M_2^{-1} x_n.
C
C     The algorithm was described  in the RIACS Technical Report 91.26,
C     `A  Quasi-Minimal  Residual Squared  Algorithm for  Non-Hermitian
C     Linear Systems`, by Freund and Szeto, December 1991.
C
C     Parameters:
C     For a description of  the parameters, see the file "dsqbg.doc" in
C     the current directory.
C
C     External routines used:
C     double precision dlamch(ch)
C        LAPACK routine, computes machine-related constants.
C     double precision dnrm2(n,x,incx)
C        BLAS-1 routine, computes the 2-norm of x.
C     subroutine daxpby(n,z,a,x,b,y)
C        Library routine, computes z = a * x + b * y.
C     double precision ddot(n,x,incx,y,incy)
C        BLAS-1 routine, computes y^H * x.
C     subroutine drandn(n,x,seed)
C        Library routine, fills x with random numbers.
C
C     Noel M. Nachtigal
C     March 9, 1992
C
C**********************************************************************
C
      INTRINSIC DBLE, DSQRT, MAX0
      EXTERNAL DLAMCH, DNRM2, DAXPBY, DDOT, DRANDN
      DOUBLE PRECISION DDOT
      DOUBLE PRECISION DLAMCH, DNRM2
C
      INTEGER INFO(4), NLEN, NDIM, NLIM
      DOUBLE PRECISION VECS(NDIM,6)
      DOUBLE PRECISION TOL
C
C     Miscellaneous parameters.
C
      DOUBLE PRECISION DHUN, DONE, DTEN, DZERO
      PARAMETER (DHUN = 1.0D2,DONE = 1.0D0,DTEN = 1.0D1,DZERO = 0.0D0)
C
C     Local variables, permanent.
C
      INTEGER IERR, N, RETLBL, TF, TRES, VF
      SAVE    IERR, N, RETLBL, TF, TRES, VF
      DOUBLE PRECISION ALPHA, RHO
      SAVE           ALPHA, RHO
      DOUBLE PRECISION QPRD, R0, RESN, TAU, UNRM, VAR
      SAVE             QPRD, R0, RESN, TAU, UNRM, VAR
C
C     Local variables, transient.
C
      INTEGER I, REVCOM
      DOUBLE PRECISION BETA
      DOUBLE PRECISION COS, DTMP
C
C     Initialize some of the permanent variables.
C
      DATA RETLBL /0/
C
C     Check the reverse communication flag to see where to branch.
C        REVCOM   RETLBL      Comment
C           0        0    first call, go to label 10
C           1       30    returning from AXB, go to label 30
C           1       40    returning from AXB, go to label 40
C
      REVCOM  = INFO(2)
      INFO(2) = 0
      IF (REVCOM.EQ.0) THEN
         N = 0
         IF (RETLBL.EQ.0) GO TO 10
      ELSE IF (REVCOM.EQ.1) THEN
         IF (RETLBL.EQ.30) THEN
            GO TO 30
         ELSE IF (RETLBL.EQ.40) THEN
            GO TO 40
         END IF
      END IF
      IERR = 1
      GO TO 70
C
C     Check whether the inputs are valid.
C
 10   IERR = 0
      IF (NDIM.LT.1)    IERR = 2
      IF (NLEN.LT.1)    IERR = 2
      IF (NLIM.LT.1)    IERR = 2
      IF (NLEN.GT.NDIM) IERR = 2
      IF (IERR.NE.0) GO TO 70
C
C     Extract from INFO the output units TF and VF, the true residual
C     flag TRES.
C
      VF   = MAX0(INFO(1),0)
      I    = VF / 100000
      VF   = VF - I * 100000
      TRES = VF / 10000
      VF   = VF - TRES * 10000
      TF   = VF / 100
      VF   = VF - TF * 100
C
C     Check the covergence tolerance.
C
      IF (TOL.LE.DZERO) TOL = DSQRT(DLAMCH('E'))
C
C     Start the trace messages and convergence history.
C
      IF (VF.NE.0) WRITE (VF,'(I8,2E11.4)') 0, DONE, DONE
      IF (TF.NE.0) WRITE (TF,'(I8,2E11.4)') 0, DONE, DONE
C
C     Set x_0 = 0 and compute the norm of the initial residual.
C
      CALL DAXPBY (NLEN,VECS(1,3),DONE,VECS(1,2),DZERO,VECS(1,3))
      CALL DAXPBY (NLEN,VECS(1,1),DZERO,VECS(1,1),DZERO,VECS(1,1))
      R0 = DNRM2(NLEN,VECS(1,3),1)
      IF ((TOL.GE.DONE).OR.(R0.EQ.DZERO)) GO TO 70
C
C     Initialize the variables.
C
      N    = 1
      QPRD = R0
      RESN = DONE
      VAR  = DZERO
      TAU  = R0 * R0
      RHO  = DDOT(NLEN,VECS(1,3),1,VECS(1,3),1)
      CALL DAXPBY (NLEN,VECS(1,4),DONE,VECS(1,3),DZERO,VECS(1,4))
C
C     This is one step of the QMR-BCG algorithm.
C
C     Have the caller carry out AXB, then return here.
C        CALL AXB (VECS(1,4),VECS(1,6))
C
 20   INFO(2) = 1
      INFO(3) = 4
      INFO(4) = 6
      RETLBL  = 30
      RETURN
C
C     Compute \sigma and check for breakdown.
C
 30   DTMP = DDOT(NLEN,VECS(1,4),1,VECS(1,6),1)
      IF (DTMP.EQ.DZERO) THEN
         IERR = 8
         GO TO 70
      END IF
C
C     Compute \alpha and the updated BCG residual.
C
      ALPHA = RHO / DTMP
      CALL DAXPBY (NLEN,VECS(1,3),DONE,VECS(1,3),-ALPHA,VECS(1,6))
C
C     Compute the updated BCG residual norm.
C
      DTMP = DNRM2(NLEN,VECS(1,3),1)
      DTMP = DTMP * DTMP
C
C     Compute the QMR quantities.
C
      BETA  = VAR
      VAR   = DTMP / TAU
      COS   = DONE / (DONE + VAR )
      TAU   = DTMP * COS
      QPRD  = QPRD * DSQRT(DONE - COS)
      UNRM  = QPRD * DSQRT(DBLE(N+1))
C
C     Compute \hat{p} and the QMR iterate x_n.
C
      DTMP = ALPHA * COS
      BETA = BETA  * COS
      CALL DAXPBY (NLEN,VECS(1,5),BETA,VECS(1,5),DTMP,VECS(1,4))
      CALL DAXPBY (NLEN,VECS(1,1),DONE,VECS(1,1),DONE,VECS(1,5))
C
C     Compute the updated residual norm.
C     If the updated residual norm is  within one order of magnitude of
C     the target convergence norm, compute the true residual norm.
C
      IF ((TRES.EQ.0).AND.(UNRM/TOL.GT.DTEN).AND.(N.LT.NLIM)) GO TO 50
C
C     Have the caller carry out AXB, then return here.
C        CALL AXB (VECS(1,1),VECS(1,6))
C
      INFO(2) = 1
      INFO(3) = 1
      INFO(4) = 6
      RETLBL  = 40
      RETURN
 40   CALL DAXPBY (NLEN,VECS(1,6),DONE,VECS(1,2),-DONE,VECS(1,6))
      RESN = DNRM2(NLEN,VECS(1,6),1) / R0
C
C     Output the convergence history.
C
 50   IF (VF.NE.0) WRITE (VF,'(I8,2E11.4)') N, UNRM, RESN
      IF (TF.NE.0) WRITE (TF,'(I8,2E11.4)') N, UNRM, RESN
C
C     Check for convergence or termination.  Stop if:
C         1. algorithm converged;
C         2. the updated residual norm is smaller than the computed
C     residual norm by a factor of at least 100;
C         3. algorithm exceeded the iterations limit.
C
      IF (RESN.LE.TOL) THEN
         IERR = 0
         GO TO 70
      ELSE IF (UNRM.LT.RESN/DHUN) THEN
         IERR = 4
         GO TO 70
      ELSE IF (N.GE.NLIM) THEN
         IERR = 4
         GO TO 70
      END IF
C
C     Check for breakdown due to \rho = 0.
C
      IF (RHO.EQ.DZERO) THEN
         IERR = 8
         GO TO 70
      END IF
C
C     Compute \beta and \rho.
C
      DTMP = DDOT(NLEN,VECS(1,3),1,VECS(1,3),1)
      BETA = DTMP / RHO
      RHO  = DTMP
C
C     Compute q, \tilde{q}, and A q.
C
      CALL DAXPBY (NLEN,VECS(1,4),BETA,VECS(1,4),DONE,VECS(1,3))
C
C     Update the running counter.
C
      N = N + 1
      GO TO 20
C
C     Done.
C
 70   NLIM    = N
      RETLBL  = 0
      INFO(1) = IERR
C
      RETURN
      END
C
C**********************************************************************