c----------- cgcommon package -------------------------------------------
c these routines are required for use by the conjugate gradient machinery
c history
c 03.01.28 (BTD) moved routines CAXPY and CSWAP to blas.f package of
c                basic linear algebra subroutines
c end history
c------------------------------------------------------------------------

      REAL FUNCTION SCSETRHSSTOP(B,WRK,EPSILON,IPAR,PRECONL,PSCNRM)
      IMPLICIT NONE

*     .. Scalar Arguments ..
      REAL EPSILON
*     ..
*     .. Array Arguments ..
      COMPLEX B(*),WRK(*)
      INTEGER IPAR(*)
*     ..
*     .. Function Arguments ..
      REAL PSCNRM
      EXTERNAL PSCNRM
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL PRECONL
*     ..
*     .. Local Scalars ..
      INTEGER LOCLEN,STOPTYPE
*     ..
      LOCLEN = IPAR(4)
      STOPTYPE = IPAR(9)
      IF ((STOPTYPE.EQ.1) .OR. (STOPTYPE.EQ.4) .OR.
     +    (STOPTYPE.EQ.7)) THEN
*  ||r||<epsilon or ||Q1r||<epsilon ||x(k)-x(k-1)||<epsilon
          SCSETRHSSTOP = EPSILON

      ELSE IF ((STOPTYPE.EQ.2) .OR. (STOPTYPE.EQ.3) .OR.
     +         (STOPTYPE.EQ.5)) THEN
*  ||r||<epsilon||b|| or sqrt(r(Q1r))<epsilon||b|| or
*  ||Q1r||<epsilon||b||
          SCSETRHSSTOP = EPSILON*PSCNRM(LOCLEN,B)

      ELSE IF (STOPTYPE.EQ.6) THEN
*  ||Q1r||<epsilon||Q1b||
          CALL PRECONL(B,WRK,IPAR)
          SCSETRHSSTOP = EPSILON*PSCNRM(LOCLEN,WRK)
      END IF

      RETURN

      END
      SUBROUTINE STOPCRIT(B,R,RTRUE,X,XOLD,WRK,RHSSTOP,CNVRTX,EXITNORM,
     +                    STATUS,IPAR,MATVEC,TMATVEC,PRECONR,PCSUM,
     +                    PSCNRM)
      IMPLICIT NONE

*     .. Scalar Arguments ..
      REAL EXITNORM,RHSSTOP
      INTEGER CNVRTX,STATUS
*     ..
*     .. Array Arguments ..
      COMPLEX B(*),R(*),RTRUE(*),WRK(*),X(*),XOLD(*)
      INTEGER IPAR(*)
*     ..
*     .. Function Arguments ..
      REAL PSCNRM
      EXTERNAL PSCNRM
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PCSUM,PRECONR,TMATVEC
*     ..
*     .. Local Scalars ..
      INTEGER LOCLEN,PRECONTYPE,STOPTYPE
*     ..
*     .. Local Arrays ..
      COMPLEX DOTS(1)
*     ..
*     .. External Functions ..
C 98.04.27 BTD: CDOTC is defined as a complex function!
C 98.12.03 (BTD) this error also pointed out by Rodolphe Vaillon
C      REAL CDOTC
      COMPLEX CDOTC
C------------------------------------------------
      EXTERNAL CDOTC
*     ..
*     .. External Subroutines ..
      EXTERNAL CAXPY,CCOPY
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS,SQRT
*     ..
*     .. Parameters ..
      COMPLEX CONE
      PARAMETER (CONE= (1.0E0,0.0E0))
*     ..
      LOCLEN = IPAR(4)
      PRECONTYPE = IPAR(8)
      STOPTYPE = IPAR(9)

      IF ((STOPTYPE.EQ.1) .OR. (STOPTYPE.EQ.2) .OR.
     +    (STOPTYPE.EQ.3)) THEN

*  Compute true residual if needed
          CALL CCOPY(LOCLEN,B,1,RTRUE,1)

          IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
              CALL PRECONR(X,WRK,IPAR)
              IF (CNVRTX.EQ.1) THEN
                  CALL TMATVEC(WRK,XOLD,IPAR)
                  CALL MATVEC(XOLD,WRK,IPAR)
                  CALL CAXPY(LOCLEN,-CONE,WRK,1,RTRUE,1)
              ELSE
                  CALL MATVEC(WRK,XOLD,IPAR)
                  CALL CAXPY(LOCLEN,-CONE,XOLD,1,RTRUE,1)
              END IF
          ELSE IF (CNVRTX.EQ.1) THEN
              CALL TMATVEC(X,XOLD,IPAR)
              CALL MATVEC(XOLD,WRK,IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK,1,RTRUE,1)
          ELSE
              CALL MATVEC(X,WRK,IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK,1,RTRUE,1)
          END IF
      END IF

      IF ((STOPTYPE.EQ.1) .OR. (STOPTYPE.EQ.2)) THEN

*  ||r|| < epsilon or ||r|| < epsilon ||b||

          EXITNORM = PSCNRM(LOCLEN,RTRUE)
          IF (EXITNORM.LT.RHSSTOP) THEN
              STATUS = 0
          ELSE
              STATUS = -99
          END IF

      ELSE IF (STOPTYPE.EQ.3) THEN

*  sqrt(r(Q1r))<epsilon||b||
C
C BTD 98.04.27 CDOTC is defined as a complex function!
C          DOTS(1) = CDOTC(LOCLEN,RTRUE,1,R,1)
C
          DOTS(1) = REAL(CDOTC(LOCLEN,RTRUE,1,R,1))
          CALL PCSUM(1,DOTS(1))

          EXITNORM = ABS(SQRT(DOTS(1)))
          IF (EXITNORM.LT.RHSSTOP) THEN
              STATUS = 0
          ELSE
              STATUS = -99
          END IF

      ELSE IF ((STOPTYPE.EQ.4) .OR. (STOPTYPE.EQ.5) .OR.
     +         (STOPTYPE.EQ.6)) THEN

*  ||Q1r|| < epsilon or ||Q1r|| < epsilon||b|| or ||Q1r|| < epsilon||Q1b||

          EXITNORM = PSCNRM(LOCLEN,R)
          IF (EXITNORM.LT.RHSSTOP) THEN
              STATUS = 0
          ELSE
              STATUS = -99
          END IF

      ELSE IF (STOPTYPE.EQ.7) THEN

*  ||x-x0||<epsilon
          CALL CCOPY(LOCLEN,X,1,WRK,1)
          CALL CAXPY(LOCLEN,-CONE,XOLD,1,WRK,1)
          EXITNORM = PSCNRM(LOCLEN,WRK)
          IF (EXITNORM.LT.RHSSTOP) THEN
              STATUS = 0
          ELSE
              STATUS = -99
          END IF
      END IF
      
      RETURN

      END

      SUBROUTINE PROGRESS(LOCLEN,ITNO,NORMRES,X,RES,TRUERES)
C part of PIM
C history
C 95.08.14 (BTD): modified to include COMMON/NORMERR/ERRSCAL to allow 
C                 desired normalization of error.
C 96.11.21 (BTD): added IDVOUT2 to COMMON/NORMERR/ to allow passing
C                 of device number for standard output (so that it
C                 doesn't need to be hardwired here).
C 98.10.08 (BTD): modified to keep track of minimum error, and to
C                 note when new error minimum is attained.  This is
C                 to observe convergence behavior of PBCGST.
C 98.10.26 (BTD): modified to PROGRESS to also print information on
C                 rate of convergence (variable RATE).
C 03.04.13 (BTD): added ITERMX and ITERN to COMMON/NORMERR/
C                 to communicate number of iterations and max number
C                 of iterations allowed.
C end history
C
C Note: when used with STOPTYPE=5, NORMRES=|Ax-b|.
C If quantity ERRSCAL is set to |b| elsewhere, then NORMRES/ERRSCAL
C is a measure of the fractional error in the solution vector x.

      IMPLICIT NONE
C Arguments:
      REAL NORMRES
      INTEGER ITNO,LOCLEN
      COMPLEX RES(*),TRUERES(*),X(*)
C
C Common:
      INTEGER IDVOUT2,ITERMX,ITERN
      REAL ERRSCAL
C-----------------------------------------------------------------------
      COMMON/NORMERR/ERRSCAL,IDVOUT2,ITERMX,ITERN
C-----------------------------------------------------------------------
C
C Local:
      INTEGER ITNOL
      REAL ERRMIN,NORMERR,RATE
      SAVE ERRMIN,ITNOL
C
      ITERN=ITNO
      NORMERR=NORMRES/ERRSCAL
      IF(ITNO.LE.2)THEN
         ERRMIN=NORMERR
         ITNOL=ITNO
         WRITE(IDVOUT2,FMT=9000)ITNO,NORMERR
      ELSE
         IF(NORMERR.LT.ERRMIN)THEN
            RATE=ALOG(ERRMIN/NORMERR)/(ITNO-ITNOL)
            ERRMIN=NORMERR
            ITNOL=ITNO
            WRITE(IDVOUT2,FMT=9001)ITNO,NORMERR,ERRMIN,RATE
         ELSE
            WRITE(IDVOUT2,FMT=9000)ITNO,NORMERR
         ENDIF
      ENDIF
c*** btd 98.12.07 experiment
      RETURN
 9000 FORMAT(1X,'iter= ',I4,' frac.err= ',0PF11.7)
 9001 FORMAT(1X,'iter= ',I4,' frac.err= ',0PF11.7,' min.err=',0PF11.7,
     &      ' rate=',0PF8.6)
      END

      SUBROUTINE CVPROD(N,CX,INCX,CY,INCY)
C part of PIM
      IMPLICIT NONE
*
*     Modified from saxpy level 1 BLAS
*     element-wise vector multiplication, y<-x*y
*     Rudnei Dias da Cunha, 16/6/93
*
*     constant times a vector plus a vector.
*     uses unrolled loops for increments equal to one.
*     jack dongarra, linpack, 3/11/78.
*
*
*     .. Scalar Arguments ..
      INTEGER INCX,INCY,N
*     ..
*     .. Array Arguments ..
      COMPLEX CX(*),CY(*)
*     ..
*     .. Local Scalars ..
      INTEGER I,IX,IY,M,MP1
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC MOD
*     ..
      IF (N.LE.0) RETURN
      IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20
*
*        code for unequal increments or equal increments
*          not equal to 1
*
      IX = 1
      IY = 1
      IF (INCX.LT.0) IX = (-N+1)*INCX + 1
      IF (INCY.LT.0) IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
          CY(IY) = CY(IY)*CX(IX)
          IX = IX + INCX
          IY = IY + INCY
   10 CONTINUE
      RETURN
*
*        code for both increments equal to 1
*
*
*        clean-up loop
*
   20 M = MOD(N,4)
      IF (M.EQ.0) GO TO 40
      DO 30 I = 1,M
          CY(I) = CY(I)*CX(I)
   30 CONTINUE
      IF (N.LT.4) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,4
          CY(I) = CY(I)*CX(I)
          CY(I+1) = CY(I+1)*CX(I+1)
          CY(I+2) = CY(I+2)*CX(I+2)
          CY(I+3) = CY(I+3)*CX(I+3)
   50 CONTINUE
      RETURN

      END
      SUBROUTINE PRINTV(N,U)
C part of PIM
      IMPLICIT NONE
*     .. Scalar Arguments ..
      INTEGER N
*     ..
*     .. Array Arguments ..
      COMPLEX U(*)
*     ..
*     .. Local Scalars ..
      INTEGER I
*     ..
      WRITE (6,FMT=9000) (U(I),I=1,N)
      RETURN
 9000 FORMAT (8 (E14.8,1X))
      END

      SUBROUTINE CINIT(N,ALPHA,CX,INCX)
C part of PIM
      IMPLICIT NONE

*     Initialises a vector x with a scalar alpha.
*     Modified from scopy, BLAS Level 1.
*     Rudnei Dias da Cunha, 14/6/93.
*

*     copies a vector, x, to a vector, y.
*     uses unrolled loops for increments equal to one.
*     jack dongarra, linpack, 3/11/78.
*
*
*     .. Scalar Arguments ..
      COMPLEX ALPHA
      INTEGER INCX,N
*     ..
*     .. Array Arguments ..
      COMPLEX CX(*)
*     ..
*     .. Local Scalars ..
      INTEGER I,IX,M,MP1
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC MOD
*     ..
      IF (N.LE.0) RETURN
      IF (INCX.EQ.1) GO TO 20
*
*        code for unequal increments or equal increments
*          not equal to 1
*
      IX = 1
      IF (INCX.LT.0) IX = (-N+1)*INCX + 1
      DO 10 I = 1,N
          CX(IX) = ALPHA
          IX = IX + INCX
   10 CONTINUE
      RETURN
*
*        code for both increments equal to 1
*
*
*        clean-up loop
*
   20 M = MOD(N,7)
      IF (M.EQ.0) GO TO 40
      DO 30 I = 1,M
          CX(I) = ALPHA
   30 CONTINUE
      IF (N.LT.7) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,7
          CX(I) = ALPHA
          CX(I+1) = ALPHA
          CX(I+2) = ALPHA
          CX(I+3) = ALPHA
          CX(I+4) = ALPHA
          CX(I+5) = ALPHA
          CX(I+6) = ALPHA
   50 CONTINUE
      RETURN

      END

      COMPLEX FUNCTION CSIGN(X)
C part of PIM
      IMPLICIT NONE
*     .. Scalar Arguments ..
      COMPLEX X
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS
*     ..
      CSIGN = X/ABS(X)
      RETURN

      END

      SUBROUTINE DECODE(RHO,C,S)
C part of PIM
      IMPLICIT NONE
*     .. Scalar Arguments ..
      COMPLEX C,RHO,S
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS,SQRT
*     ..
*     .. Parameters ..
      REAL ONE
      PARAMETER (ONE=1.0)
      COMPLEX CZERO
      PARAMETER (CZERO= (0.0,0.0))
      COMPLEX CONE
      PARAMETER (CONE= (1.0,0.0))
*     ..
      IF (RHO.EQ.CONE) THEN
          C = CZERO
          S = CONE

      ELSE IF (ABS(RHO).LT.ONE) THEN
          S = 2.0*RHO
          C = SQRT(CONE-S**2)

      ELSE
          C = 2.0/RHO
          S = SQRT(CONE-C**2)
      END IF

      RETURN

      END

      SUBROUTINE ENCODE(RHO,C,S)
C part of PIM
      IMPLICIT NONE
*     .. Scalar Arguments ..
      COMPLEX C,RHO,S
*     ..
*     .. External Functions ..
      COMPLEX CSIGN
      EXTERNAL CSIGN
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS
*     ..
*     .. Parameters ..
      COMPLEX CZERO
      PARAMETER (CZERO= (0.0,0.0))
      COMPLEX CONE
      PARAMETER (CONE= (1.0,0.0))
*     ..
      IF (C.EQ.CZERO) THEN
          RHO = CONE

      ELSE IF (ABS(S).LT.ABS(C)) THEN
          RHO = CSIGN(C)*S/2.0

      ELSE
          RHO = 2.0*CSIGN(S)/C
      END IF

      RETURN

      END

      SUBROUTINE GIVENS(A,B,C,S)
C part of PIM
      IMPLICIT NONE
*     .. Scalar Arguments ..
      COMPLEX A,B,C,S
*     ..
*     .. Local Scalars ..
      COMPLEX TAU
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS,SQRT
*     ..
*     .. Parameters ..
      COMPLEX CZERO
      PARAMETER (CZERO= (0.0,0.0))
      COMPLEX CONE
      PARAMETER (CONE= (1.0,0.0))
*     ..
      IF (B.EQ.CZERO) THEN
          C = CONE
          S = CZERO

      ELSE IF (ABS(B).GT.ABS(A)) THEN
          TAU = -A/B
          S = CONE/SQRT(CONE+TAU**2)
          C = S*TAU

      ELSE
          TAU = -B/A
          C = CONE/SQRT(CONE+TAU**2)
          S = C*TAU
      END IF

      RETURN

      END
      REAL FUNCTION PSCNRM2(LOCLEN,U)
C part of PIM package
*     .. Scalar Arguments ..
      INTEGER LOCLEN
*     ..
*     .. Array Arguments ..
      COMPLEX U(*)
*     ..
*     .. External Functions ..
      REAL SCNRM2
      EXTERNAL SCNRM2
*     ..
      PSCNRM2 = SCNRM2(LOCLEN,U,1)
      RETURN

      END

      SUBROUTINE PCSUM(ISIZE,X)
*     .. Scalar Arguments ..
      INTEGER ISIZE
*     ..
*     .. Array Arguments ..
      COMPLEX X(*)
*     ..
      RETURN

      END

      subroutine  ccopy(n,cx,incx,cy,incy)
c
c     copies a vector, x, to a vector, y.
c     jack dongarra, linpack, 4/11/78.
c
      complex cx(1),cy(1)
      integer i,incx,incy,ix,iy,n
c
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        cy(iy) = cx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
   20 do 30 i = 1,n
        cy(i) = cx(i)
   30 continue
      return
      end
      complex function cdotc(n,cx,incx,cy,incy)
c
c     forms the dot product of a vector.
c     jack dongarra, 3/11/78.
c
      complex cx(1),cy(1),ctemp
      integer i,ix,iy,n,incx,incy
      ctemp = (0.0e0,0.0e0)
      cdotc = (0.0e0,0.0e0)
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        ctemp = ctemp + conjg(cx(ix))*cy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      cdotc = ctemp
      return
c
c        code for both increments equal to 1
c
   20 do 30 i = 1,n
        ctemp = ctemp + conjg(cx(i))*cy(i)
   30 continue
      cdotc = ctemp
      return
      end
      complex function cdotu(n,cx,incx,cy,incy)
c
c     forms the dot product of a vector.
c     jack dongarra, 3/11/78.
c
      complex cx(1),cy(1),ctemp
      integer i,ix,iy,n,incx,incy
      ctemp = (0.0e0,0.0e0)
      cdotu = (0.0e0,0.0e0)
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        ctemp = ctemp + cx(ix)*cy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      cdotu = ctemp
      return
c
c        code for both increments equal to 1
c
   20 do 30 i = 1,n
        ctemp = ctemp + cx(i)*cy(i)
   30 continue
      cdotu = ctemp
      return
      end
      subroutine crotg(ca,cb,c,s)
      complex ca,cb,s
      real c
      real norm,scale
      complex alpha
      if (cabs(ca) .ne. 0.0e0) go to 10
         c = 0.0e0
         s = (1.0e0,0.0e0)
         ca = cb
         go to 20
   10 continue
         scale = cabs(ca) + cabs(cb)
         norm = scale*sqrt((cabs(ca/cmplx(scale,0.0e0)))**2 +
     *                      (cabs(cb/cmplx(scale,0.0e0)))**2)
         alpha = ca /cabs(ca)
         c = cabs(ca) / norm
         s = alpha * conjg(cb) / norm
         ca = alpha * norm
   20 continue
      return
      end
      subroutine  csrot (n,cx,incx,cy,incy,c,s)
c
c     applies a plane rotation, where the cos and sin (c and s) are
c     real and the vectors cx and cy are complex.
c     jack dongarra, linpack, 3/11/78.
c
      complex cx(1),cy(1),ctemp
      real c,s
      integer i,incx,incy,ix,iy,n
c
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c       code for unequal increments or equal increments not equal
c         to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        ctemp = c*cx(ix) + s*cy(iy)
        cy(iy) = c*cy(iy) - s*cx(ix)
        cx(ix) = ctemp
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c       code for both increments equal to 1
c
   20 do 30 i = 1,n
        ctemp = c*cx(i) + s*cy(i)
        cy(i) = c*cy(i) - s*cx(i)
        cx(i) = ctemp
   30 continue
      return
      end
      REAL FUNCTION SCNRM2(N,CX,INCX)
C Arguments:
      INTEGER INCX,N
      COMPLEX CX(1)
C Local variables:
      INTEGER I,J
      DOUBLE PRECISION SUM
C***********************************************************************
C Returns SCNRM2=unitary norm of complex n-vector stored in CX with
C storage increment INCX.
C Written to replace SCNRM2 written by C.L.Lawson, as g77 compiler 
C produces bad code when optimizing Lawson's routine.
C In any event, Lawson's routine appears to be unnecessarily complicated
C for needs of DDSCAT.
C
C History:
C 98.10.07 (BTD) Created by B.T. Draine, Princeton Univ. Observatory,
C                for use by DDSCAT
C 00.06.13 (BTD) Make data conversion explicit.
C end history
C***********************************************************************
      SUM=0.
      DO I=1,N
         J=1+(I-1)*INCX
         SUM=SUM+DBLE(CX(J)*CONJG(CX(J)))
      ENDDO
      SCNRM2=REAL(SQRT(SUM))
      RETURN
      END
C-----------------------------------------------------------------------
C Following code due to C.L. Lawson has been replaced by above module
C because it generated bad code when compiled with g77.
C For use by DDSCAT the overflow/underflow avoidance strategies
C used by this routine do not appear to be necessary, so they can
C be omitted in interests of speed and clean code.
C 98.10.07 BTD
C
C      real function scnrm2( n, cx, incx)
C      logical imag, scale
C      integer i, incx, ix, n, next
C      real cutlo, cuthi, hitest, sum, xmax, absx, zero, one
C      complex      cx(1)
C      real real,aimag
C      complex zdumr,zdumi
C      real(zdumr) = zdumr
C      aimag(zdumi) = (0.0e0,-1.0e0)*zdumi
C      data         zero, one /0.0e0, 1.0e0/
Cc
Cc     unitary norm of the complex n-vector stored in cx() with storage
Cc     increment incx .
Cc     if    n .le. 0 return with result = 0.
Cc     if n .ge. 1 then incx must be .ge. 1
Cc
Cc           c.l.lawson , 1978 jan 08
Cc     modified to correct problem with negative increment, 8/21/90.
Cc
Cc     four phase method     using two built-in constants that are
Cc     hopefully applicable to all machines.
Cc         cutlo = maximum of  sqrt(u/eps)  over all known machines.
Cc         cuthi = minimum of  sqrt(v)      over all known machines.
Cc     where
Cc         eps = smallest no. such that eps + 1. .gt. 1.
Cc         u   = smallest positive no.   (underflow limit)
Cc         v   = largest  no.            (overflow  limit)
Cc
Cc     brief outline of algorithm..
Cc
Cc     phase 1    scans zero components.
Cc     move to phase 2 when a component is nonzero and .le. cutlo
Cc     move to phase 3 when a component is .gt. cutlo
Cc     move to phase 4 when a component is .ge. cuthi/m
Cc     where m = n for x() real and m = 2*n for complex.
Cc
Cc     values for cutlo and cuthi..
Cc     from the environmental parameters listed in the imsl converter
Cc     document the limiting values are as follows..
Cc     cutlo, s.p.   u/eps = 2**(-102) for  honeywell.  close seconds are
Cc                   univac and dec at 2**(-103)
Cc                   thus cutlo = 2**(-51) = 4.44089e-16
Cc     cuthi, s.p.   v = 2**127 for univac, honeywell, and dec.
Cc                   thus cuthi = 2**(63.5) = 1.30438e19
Cc     cutlo, d.p.   u/eps = 2**(-67) for honeywell and dec.
Cc                   thus cutlo = 2**(-33.5) = 8.23181d-11
Cc     cuthi, d.p.   same as s.p.  cuthi = 1.30438d19
Cc     data cutlo, cuthi / 8.232d-11,  1.304d19 /
Cc     data cutlo, cuthi / 4.441e-16,  1.304e19 /
C      data cutlo, cuthi / 8.232d-11,  1.304d19 /
Cc
C      if(n .gt. 0) go to 10
C         scnrm2  = zero
C         go to 300
Cc
C   10 assign 30 to next
C      sum = zero
C      i = 1
C      if( incx .lt. 0 )i = (-n+1)*incx + 1
Cc                                                 begin main loop
C      do 220 ix = 1,n
C         absx = abs(real(cx(i)))
C         imag = .false.
C         go to next,(30, 50, 70, 90, 110)
C   30 if( absx .gt. cutlo) go to 85
C      assign 50 to next
C      scale = .false.
Cc
Cc                        phase 1.  sum is zero
Cc
C   50 if( absx .eq. zero) go to 200
C      if( absx .gt. cutlo) go to 85
Cc
Cc                                prepare for phase 2.
C      assign 70 to next
C      go to 105
Cc
Cc                                prepare for phase 4.
Cc
C  100 assign 110 to next
C      sum = (sum / absx) / absx
C  105 scale = .true.
C      xmax = absx
C      go to 115
Cc
Cc                   phase 2.  sum is small.
Cc                             scale to avoid destructive underflow.
Cc
C   70 if( absx .gt. cutlo ) go to 75
Cc
Cc                     common code for phases 2 and 4.
Cc                     in phase 4 sum is large.  scale to avoid overflow.
Cc
C  110 if( absx .le. xmax ) go to 115
C         sum = one + sum * (xmax / absx)**2
C         xmax = absx
C         go to 200
Cc
C  115 sum = sum + (absx/xmax)**2
C      go to 200
Cc
Cc
Cc                  prepare for phase 3.
Cc
C   75 sum = (sum * xmax) * xmax
Cc
C   85 assign 90 to next
C      scale = .false.
Cc
Cc     for real or d.p. set hitest = cuthi/n
Cc     for complex      set hitest = cuthi/(2*n)
Cc
C      hitest = cuthi/float( 2*n )
Cc
Cc                   phase 3.  sum is mid-range.  no scaling.
Cc
C   90 if(absx .ge. hitest) go to 100
C         sum = sum + absx**2
C  200 continue
Cc                  control selection of real and imaginary parts.
Cc
C      if(imag) go to 210
C         absx = abs(aimag(cx(i)))
C         imag = .true.
C      go to next,(  50, 70, 90, 110 )
Cc
C  210 continue
C      i = i + incx
C  220 continue
Cc
Cc              end of main loop.
Cc              compute square root and adjust for scaling.
Cc
C      scnrm2 = sqrt(sum)
C      if(scale) scnrm2 = scnrm2 * xmax
C  300 continue
C      return
C      end
C-----------------------------------------------------------------------
      subroutine srotg(da,db,c,s)
c
c     construct givens plane rotation.
c     jack dongarra, linpack, 3/11/78.
c                    modified 9/27/86.
c
      real da,db,c,s,roe,scale,r,z
c
      roe = db
      if( abs(da) .gt. abs(db) ) roe = da
      scale = abs(da) + abs(db)
      if( scale .ne. 0.0e0 ) go to 10
         c = 1.0e0
         s = 0.0e0
         r = 0.0e0
         go to 20
   10 r = scale*sqrt((da/scale)**2 + (db/scale)**2)
      r = sign(1.0e0,roe)*r
      c = da/r
      s = db/r
   20 z = s
      if( abs(c) .gt. 0.0e0 .and. abs(c) .le. s ) z = 1.0e0/c
      da = r
      db = z
      return
      end

      SUBROUTINE DMACHCONS(WHAT,RESULT)

* These values are for IEEE-754 arithmetic

*     .. Parameters ..
      DOUBLE PRECISION MACHEPS
      PARAMETER (MACHEPS=2.2204460492503D-16)
      DOUBLE PRECISION UNDERFLOW
      PARAMETER (UNDERFLOW=0.22250739D-307)
      DOUBLE PRECISION OVERFLOW
      PARAMETER (OVERFLOW=1.7976313D+308)
*     ..
*     .. Scalar Arguments ..
      DOUBLE PRECISION RESULT
      CHARACTER WHAT
*     ..
      IF ((WHAT.EQ.'M') .OR. (WHAT.EQ.'m')) THEN
          RESULT = MACHEPS
      ELSE IF ((WHAT.EQ.'U') .OR. (WHAT.EQ.'u')) THEN
          RESULT = UNDERFLOW
      ELSE IF ((WHAT.EQ.'O') .OR. (WHAT.EQ.'o')) THEN
          RESULT = OVERFLOW
      END IF

      RETURN

      END
      SUBROUTINE PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,
     +                      NPROCS,PROCID,PRECONTYPE,STOPTYPE,MAXIT,
     +                      ITNO,STATUS,STEPERR,EPSILON,EXITNORM)

c           PIM -- The Parallel Iterative Methods package
c           ---------------------------------------------
c
c                      Rudnei Dias da Cunha
c               Centro de Processamento de Dados,
c         Universidade Federal do Rio Grande do Sul, Brasil
c                              and
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c                          Tim Hopkins
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c ----------------------------------------------------------------------
c
c  Description of parameter arrays
c   IPAR (INPUT)  : integer
c     ipar( 1): lda    (Leading dimension of a)
c           2 : n      (Number of rows/columns of a)
c           3 : blksz  (Size of block of data; used when data is
c                       partitioned using cyclic mode)
c           4 : loclen (Number of elements stored locally;
c                       *PARALLEL: Equal to at least m/nprocs or
c                                  n/procs depending if row or
c                                  column partitioning is used or,
c                                  in the case of cyclic partitioning,
c                                  it is a multiple of either
c                                  m/(blksz*nprocs) or n/(blksz*nprocs).
c                       *SEQUENTIAL: equal to n)
c           5 : basisdim (Dimension of orthogonal basis, used in
c                       GMRES)
c           6 : nprocs (Number of processors)
c           7 : procid (Processor identification)
c           8 : precontype (Type of preconditioning; one of
c                           0 : No preconditioning,
c                           1 : Left preconditioning,
c                           2 : Right preconditioning,
c                           3 : Symmetric preconditioning)
c           9 : stoptype (Type of stopping criteria used)
c          10 : maxit  (Maximum number of iterations allowed)
c
c   IPAR (OUTPUT) : integer
c     ipar(11): itno   (Number of iterations executed)
c          12 : status (On exit of iterative method, one of
c                        0: converged to solution
c                       -1: no convergence has been achieved
c                       -2: "soft"-breakdown, solution may have
c                           been found
c                       -3: "hard"-breakdown, no solution)
c                       -4: conflict in preconditioner and stopping
c                           criterion selected
c                       -5: error in stopping criterion 3, r^{T}z<0)
c          13 : steperr (If status is either -2 or -3, it gives
c                         the step number of the respective algorithm
c                         where a breakdown has occurred. Refer to the
c                         User's Guide for further information)
c
c   RPAR/DPAR (INPUT)  : real/double precision
c     dpar( 1): epsilon (The value of epsilon for use in the
c                        stopping criterion)
c
c   RPAR/DPAR (OUTPUT) : real/double precision
c     dpar( 2): exitnorm (The value of a norm of either the residual
c                         vector or the difference between two 
c                         successive solution estimates according to
c                         the value of stoptype)
c           3,
c           4 : smallest and largest eigenvalues of Q1AQ2 (in the
c               symmetric case) OR smallest and largest real parts
c               (in the nonsymmetric case)
c           5,
c           6 : smallest and largest imaginary parts (only in the
c               nonsymmetric case)
c

C     .. Parameters ..
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=6)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION EPSILON,EXITNORM
      INTEGER BASISDIM,BLKSZ,ITNO,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,
     +        PROCID,STATUS,STEPERR,STOPTYPE
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION DPAR(DPARSIZ)
      INTEGER IPAR(IPARSIZ)
C     ..
      LDA = IPAR(1)
      N = IPAR(2)
      BLKSZ = IPAR(3)
      LOCLEN = IPAR(4)
      BASISDIM = IPAR(5)
      NPROCS = IPAR(6)
      PROCID = IPAR(7)
      PRECONTYPE = IPAR(8)
      STOPTYPE = IPAR(9)
      MAXIT = IPAR(10)
      ITNO = IPAR(11)
      STATUS = IPAR(12)
      STEPERR = IPAR(13)

      EPSILON = DPAR(1)
      EXITNORM = DPAR(2)

      RETURN

      END
      SUBROUTINE PIMDPRTPAR(IPAR,DPAR)

c           PIM -- The Parallel Iterative Methods package
c           ---------------------------------------------
c
c                      Rudnei Dias da Cunha
c               Centro de Processamento de Dados,
c         Universidade Federal do Rio Grande do Sul, Brasil
c                              and
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c                          Tim Hopkins
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c ----------------------------------------------------------------------
c
c  Description of parameter arrays
c   IPAR (INPUT)  : integer
c     ipar( 1): lda    (Leading dimension of a)
c           2 : n      (Number of rows/columns of a)
c           3 : blksz  (Size of block of data; used when data is
c                       partitioned using cyclic mode)
c           4 : loclen (Number of elements stored locally;
c                       *PARALLEL: Equal to at least m/nprocs or
c                                  n/procs depending if row or
c                                  column partitioning is used or,
c                                  in the case of cyclic partitioning,
c                                  it is a multiple of either
c                                  m/(blksz*nprocs) or n/(blksz*nprocs).
c                       *SEQUENTIAL: equal to n)
c           5 : basisdim (Dimension of orthogonal basis, used in
c                       GMRES)
c           6 : nprocs (Number of processors)
c           7 : procid (Processor identification)
c           8 : precontype (Type of preconditioning; one of
c                           0 : No preconditioning,
c                           1 : Left preconditioning,
c                           2 : Right preconditioning,
c                           3 : Symmetric preconditioning)
c           9 : stoptype (Type of stopping criteria used)
c          10 : maxit  (Maximum number of iterations allowed)
c
c   IPAR (OUTPUT) : integer
c     ipar(11): itno   (Number of iterations executed)
c          12 : status (On exit of iterative method, one of
c                        0: converged to solution
c                       -1: no convergence has been achieved
c                       -2: "soft"-breakdown, solution may have
c                           been found
c                       -3: "hard"-breakdown, no solution)
c                       -4: conflict in preconditioner and stopping
c                           criterion selected
c                       -5: error in stopping criterion 3, r^{T}z<0)
c          13 : steperr (If status is either -2 or -3, it gives
c                         the step number of the respective algorithm
c                         where a breakdown has occurred. Refer to the
c                         User's Guide for further information)
c
c   RPAR/DPAR (INPUT)  : real/double precision
c     dpar( 1): epsilon (The value of epsilon for use in the
c                        stopping criterion)
c
c   RPAR/DPAR (OUTPUT) : real/double precision
c     dpar( 2): exitnorm (The value of a norm of either the residual
c                         vector or the difference between two 
c                         successive solution estimates according to 
c                         the value of stoptype)
c           3,
c           4 : smallest and largest eigenvalues of Q1AQ2 (in the
c               symmetric case) OR smallest and largest real parts
c               (in the nonsymmetric case)
c           5,
c           6 : smallest and largest imaginary parts (only in the
c               nonsymmetric case)
c

C     .. Parameters ..
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=6)
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION DPAR(DPARSIZ)
      INTEGER IPAR(IPARSIZ)
C     ..
C     .. Local Scalars ..
      INTEGER I
C     ..
      WRITE (6,FMT=10) (IPAR(I),I=1,IPARSIZ)

   10 FORMAT ('lda=',i6,/,'n=',i6,/,'blksz=',i6,/,'loclen=',i4,/,
     +       'basisdim=',i4,/,'nprocs=',i4,/,'procid=',i4,/,
     +       'precontype=',i4,/,'stoptype=',i4,/,'maxit=',i4,/,'itno=',
     +       i4,/,'status=',i4,/,'steperr=',i4)

      WRITE (6,FMT=20) (DPAR(I),I=1,DPARSIZ)

   20 FORMAT ('epsilon=',d20.12,/,'exitnorm=',d20.12,/,
     +       'eigenvalues region:',/,4(d20.12,1x))

      RETURN

      END
      SUBROUTINE PIMDSETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,
     +                      NPROCS,PROCID,PRECONTYPE,STOPTYPE,MAXIT,
     +                      EPSILON)

c           PIM -- The Parallel Iterative Methods package
c           ---------------------------------------------
c
c                      Rudnei Dias da Cunha
c               Centro de Processamento de Dados,
c         Universidade Federal do Rio Grande do Sul, Brasil
c                              and
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c                          Tim Hopkins
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c ----------------------------------------------------------------------
c
c  Description of parameter arrays
c   IPAR (INPUT)  : integer
c     ipar( 1): lda    (Leading dimension of a)
c           2 : n      (Number of rows/columns of a)
c           3 : blksz  (Size of block of data; used when data is
c                       partitioned using cyclic mode)
c           4 : loclen (Number of elements stored locally;
c                       *PARALLEL: Equal to at least m/nprocs or
c                                  n/procs depending if row or
c                                  column partitioning is used or,
c                                  in the case of cyclic partitioning,
c                                  it is a multiple of either
c                                  m/(blksz*nprocs) or n/(blksz*nprocs).
c                       *SEQUENTIAL: equal to n)
c           5 : basisdim (Dimension of orthogonal basis, used in
c                       GMRES)
c           6 : nprocs (Number of processors)
c           7 : procid (Processor identification)
c           8 : precontype (Type of preconditioning; one of
c                           0 : No preconditioning,
c                           1 : Left preconditioning,
c                           2 : Right preconditioning,
c                           3 : Symmetric preconditioning)
c           9 : stoptype (Type of stopping criteria used)
c          10 : maxit  (Maximum number of iterations allowed)
c
c   IPAR (OUTPUT) : integer
c     ipar(11): itno   (Number of iterations executed)
c          12 : status (On exit of iterative method, one of
c                        0: converged to solution
c                       -1: no convergence has been achieved
c                       -2: "soft"-breakdown, solution may have
c                           been found
c                       -3: "hard"-breakdown, no solution)
c                       -4: conflict in preconditioner and stopping
c                           criterion selected
c                       -5: error in stopping criterion 3, r^{T}z<0)
c          13 : steperr (If status is either -2 or -3, it gives
c                         the step number of the respective algorithm
c                         where a breakdown has occurred. Refer to the
c                         User's Guide for further information)
c
c   RPAR/DPAR (INPUT)  : real/double precision
c     dpar( 1): epsilon (The value of epsilon for use in the
c                        stopping criterion)
c
c   RPAR/DPAR (OUTPUT) : real/double precision
c     dpar( 2): exitnorm (The value of a norm of either the residual
c                         vector or the difference between two
c                         successive solution estimates according to
c                         the value of stoptype)
c           3,
c           4 : smallest and largest eigenvalues of Q1AQ2 (in the
c               symmetric case) OR smallest and largest real parts
c               (in the nonsymmetric case)
c           5,
c           6 : smallest and largest imaginary parts (only in the
c               nonsymmetric case)
c

C     .. Parameters ..
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=6)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION EPSILON
      INTEGER BASISDIM,BLKSZ,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,
     +        PROCID,STOPTYPE
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION DPAR(DPARSIZ)
      INTEGER IPAR(IPARSIZ)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION EXITNORM
      INTEGER ITNO,STATUS,STEPERR
C     ..
      IPAR(1) = LDA
      IPAR(2) = N
      IPAR(3) = BLKSZ
      IPAR(4) = LOCLEN
      IPAR(5) = BASISDIM
      IPAR(6) = NPROCS
      IPAR(7) = PROCID
      IPAR(8) = PRECONTYPE
      IPAR(9) = STOPTYPE
      IPAR(10) = MAXIT
      IPAR(11) = -1
      IPAR(12) = -1
      IPAR(13) = -1

      DPAR(1) = EPSILON
      DPAR(2) = -ONE

      RETURN

      END
      SUBROUTINE PIMSGETPAR(IPAR,SPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,
     +                      NPROCS,PROCID,PRECONTYPE,STOPTYPE,MAXIT,
     +                      ITNO,STATUS,STEPERR,EPSILON,EXITNORM)

c           PIM -- The Parallel Iterative Methods package
c           ---------------------------------------------
c
c                      Rudnei Dias da Cunha
c               Centro de Processamento de Dados,
c         Universidade Federal do Rio Grande do Sul, Brasil
c                              and
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c                          Tim Hopkins
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c ----------------------------------------------------------------------
c
c  Description of parameter arrays
c   IPAR (INPUT)  : integer
c     ipar( 1): lda    (Leading dimension of a)
c           2 : n      (Number of rows/columns of a)
c           3 : blksz  (Size of block of data; used when data is
c                       partitioned using cyclic mode)
c           4 : loclen (Number of elements stored locally;
c                       *PARALLEL: Equal to at least m/nprocs or
c                                  n/procs depending if row or
c                                  column partitioning is used or,
c                                  in the case of cyclic partitioning,
c                                  it is a multiple of either
c                                  m/(blksz*nprocs) or n/(blksz*nprocs).
c                       *SEQUENTIAL: equal to n)
c           5 : basisdim (Dimension of orthogonal basis, used in
c                       GMRES)
c           6 : nprocs (Number of processors)
c           7 : procid (Processor identification)
c           8 : precontype (Type of preconditioning; one of
c                           0 : No preconditioning,
c                           1 : Left preconditioning,
c                           2 : Right preconditioning,
c                           3 : Symmetric preconditioning)
c           9 : stoptype (Type of stopping criteria used)
c          10 : maxit  (Maximum number of iterations allowed)
c
c   IPAR (OUTPUT) : integer
c     ipar(11): itno   (Number of iterations executed)
c          12 : status (On exit of iterative method, one of
c                        0: converged to solution
c                       -1: no convergence has been achieved
c                       -2: "soft"-breakdown, solution may have
c                           been found
c                       -3: "hard"-breakdown, no solution)
c                       -4: conflict in preconditioner and stopping
c                           criterion selected
c                       -5: error in stopping criterion 3, r^{T}z<0)
c          13 : steperr (If status is either -2 or -3, it gives
c                         the step number of the respective algorithm
c                         where a breakdown has occurred. Refer to the
c                         User's Guide for further information)
c
c   RPAR/DPAR (INPUT)  : real/real
c     spar( 1): epsilon (The value of epsilon for use in the
c                        stopping criterion)
c
c   RPAR/DPAR (OUTPUT) : real/real
c     spar( 2): exitnorm (The value of a norm of either the residual
c                         vector or the difference between two 
c                         successive solution estimates according to 
c                         the value of stoptype)
c           3,
c           4 : smallest and largest eigenvalues of Q1AQ2 (in the
c               symmetric case) OR smallest and largest real parts
c               (in the nonsymmetric case)
c           5,
c           6 : smallest and largest imaginary parts (only in the
c               nonsymmetric case)
c

C     .. Parameters ..
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER SPARSIZ
      PARAMETER (SPARSIZ=6)
C     ..
C     .. Scalar Arguments ..
      REAL EPSILON,EXITNORM
      INTEGER BASISDIM,BLKSZ,ITNO,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,
     +        PROCID,STATUS,STEPERR,STOPTYPE
C     ..
C     .. Array Arguments ..
      REAL SPAR(SPARSIZ)
      INTEGER IPAR(IPARSIZ)
C     ..
      LDA = IPAR(1)
      N = IPAR(2)
      BLKSZ = IPAR(3)
      LOCLEN = IPAR(4)
      BASISDIM = IPAR(5)
      NPROCS = IPAR(6)
      PROCID = IPAR(7)
      PRECONTYPE = IPAR(8)
      STOPTYPE = IPAR(9)
      MAXIT = IPAR(10)
      ITNO = IPAR(11)
      STATUS = IPAR(12)
      STEPERR = IPAR(13)

      EPSILON = SPAR(1)
      EXITNORM = SPAR(2)

      RETURN

      END
      SUBROUTINE PIMSPRTPAR(IPAR,SPAR)

c           PIM -- The Parallel Iterative Methods package
c           ---------------------------------------------
c
c                      Rudnei Dias da Cunha
c               Centro de Processamento de Dados,
c         Universidade Federal do Rio Grande do Sul, Brasil
c                              and
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c                          Tim Hopkins
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c ----------------------------------------------------------------------
c
c  Description of parameter arrays
c   IPAR (INPUT)  : integer
c     ipar( 1): lda    (Leading dimension of a)
c           2 : n      (Number of rows/columns of a)
c           3 : blksz  (Size of block of data; used when data is
c                       partitioned using cyclic mode)
c           4 : loclen (Number of elements stored locally;
c                       *PARALLEL: Equal to at least m/nprocs or
c                                  n/procs depending if row or
c                                  column partitioning is used or,
c                                  in the case of cyclic partitioning,
c                                  it is a multiple of either
c                                  m/(blksz*nprocs) or n/(blksz*nprocs).
c                       *SEQUENTIAL: equal to n)
c           5 : basisdim (Dimension of orthogonal basis, used in
c                       GMRES)
c           6 : nprocs (Number of processors)
c           7 : procid (Processor identification)
c           8 : precontype (Type of preconditioning; one of
c                           0 : No preconditioning,
c                           1 : Left preconditioning,
c                           2 : Right preconditioning,
c                           3 : Symmetric preconditioning)
c           9 : stoptype (Type of stopping criteria used)
c          10 : maxit  (Maximum number of iterations allowed)
c
c   IPAR (OUTPUT) : integer
c     ipar(11): itno   (Number of iterations executed)
c          12 : status (On exit of iterative method, one of
c                        0: converged to solution
c                       -1: no convergence has been achieved
c                       -2: "soft"-breakdown, solution may have
c                           been found
c                       -3: "hard"-breakdown, no solution)
c                       -4: conflict in preconditioner and stopping
c                           criterion selected
c                       -5: error in stopping criterion 3, r^{T}z<0)
c          13 : steperr (If status is either -2 or -3, it gives
c                         the step number of the respective algorithm
c                         where a breakdown has occurred. Refer to the
c                         User's Guide for further information)
c
c   RPAR/DPAR (INPUT)  : real/real
c     spar( 1): epsilon (The value of epsilon for use in the
c                        stopping criterion)
c
c   RPAR/DPAR (OUTPUT) : real/real
c     spar( 2): exitnorm (The value of a norm of either the residual
c                         vector or the difference between two 
c                         successive solution estimates according to 
c                         the value of stoptype)
c           3,
c           4 : smallest and largest eigenvalues of Q1AQ2 (in the
c               symmetric case) OR smallest and largest real parts
c               (in the nonsymmetric case)
c           5,
c           6 : smallest and largest imaginary parts (only in the
c               nonsymmetric case)
c

C     .. Parameters ..
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER SPARSIZ
      PARAMETER (SPARSIZ=6)
C     ..
C     .. Array Arguments ..
      REAL SPAR(SPARSIZ)
      INTEGER IPAR(IPARSIZ)
C     ..
C     .. Local Scalars ..
      INTEGER I
C     ..
      WRITE (6,FMT=10) (IPAR(I),I=1,IPARSIZ)

   10 FORMAT ('lda=',i6,/,'n=',i6,/,'blksz=',i6,/,'loclen=',i4,/,
     +       'basisdim=',i4,/,'nprocs=',i4,/,'procid=',i4,/,
     +       'precontype=',i4,/,'stoptype=',i4,/,'maxit=',i4,/,'itno=',
     +       i4,/,'status=',i4,/,'steperr=',i4)

      WRITE (6,FMT=20) (SPAR(I),I=1,SPARSIZ)

   20 FORMAT ('epsilon=',e20.12,/,'exitnorm=',e20.12,
     +       'eigenvalues region:',/,4(e20.12,1x))


      RETURN

      END
      SUBROUTINE PIMSSETPAR(IPAR,SPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,
     +                      NPROCS,PROCID,PRECONTYPE,STOPTYPE,MAXIT,
     +                      EPSILON)

c           PIM -- The Parallel Iterative Methods package
c           ---------------------------------------------
c
c                      Rudnei Dias da Cunha
c               Centro de Processamento de Dados,
c         Universidade Federal do Rio Grande do Sul, Brasil
c                              and
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c                          Tim Hopkins
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c ----------------------------------------------------------------------
c
c  Description of parameter arrays
c   IPAR (INPUT)  : integer
c     ipar( 1): lda    (Leading dimension of a)
c           2 : n      (Number of rows/columns of a)
c           3 : blksz  (Size of block of data; used when data is
c                       partitioned using cyclic mode)
c           4 : loclen (Number of elements stored locally;
c                       *PARALLEL: Equal to at least m/nprocs or
c                                  n/procs depending if row or
c                                  column partitioning is used or,
c                                  in the case of cyclic partitioning,
c                                  it is a multiple of either
c                                  m/(blksz*nprocs) or n/(blksz*nprocs).
c                       *SEQUENTIAL: equal to n)
c           5 : basisdim (Dimension of orthogonal basis, used in
c                       GMRES)
c           6 : nprocs (Number of processors)
c           7 : procid (Processor identification)
c           8 : precontype (Type of preconditioning; one of
c                           0 : No preconditioning,
c                           1 : Left preconditioning,
c                           2 : Right preconditioning,
c                           3 : Symmetric preconditioning)
c           9 : stoptype (Type of stopping criteria used)
c          10 : maxit  (Maximum number of iterations allowed)
c
c   IPAR (OUTPUT) : integer
c     ipar(11): itno   (Number of iterations executed)
c          12 : status (On exit of iterative method, one of
c                        0: converged to solution
c                       -1: no convergence has been achieved
c                       -2: "soft"-breakdown, solution may have
c                           been found
c                       -3: "hard"-breakdown, no solution)
c                       -4: conflict in preconditioner and stopping
c                           criterion selected
c                       -5: error in stopping criterion 3, r^{T}z<0)
c          13 : steperr (If status is either -2 or -3, it gives
c                         the step number of the respective algorithm
c                         where a breakdown has occurred. Refer to the
c                         User's Guide for further information)
c
c   RPAR/DPAR (INPUT)  : real/real
c     spar( 1): epsilon (The value of epsilon for use in the
c                        stopping criterion)
c
c   RPAR/DPAR (OUTPUT) : real/real
c     spar( 2): exitnorm (The value of a norm of either the residual
c                         vector or the difference between two
c                         successive solution estimates according to
c                         the value of stoptype)
c           3,
c           4 : smallest and largest eigenvalues of Q1AQ2 (in the
c               symmetric case) OR smallest and largest real parts
c               (in the nonsymmetric case)
c           5,
c           6 : smallest and largest imaginary parts (only in the
c               nonsymmetric case)
c

C     .. Parameters ..
      REAL ONE
      PARAMETER (ONE=1.0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER SPARSIZ
      PARAMETER (SPARSIZ=6)
C     ..
C     .. Scalar Arguments ..
      REAL EPSILON
      INTEGER BASISDIM,BLKSZ,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,
     +        PROCID,STOPTYPE
C     ..
C     .. Array Arguments ..
      REAL SPAR(SPARSIZ)
      INTEGER IPAR(IPARSIZ)
C     ..
C     .. Local Scalars ..
      REAL EXITNORM
      INTEGER ITNO,STATUS,STEPERR
C     ..
      IPAR(1) = LDA
      IPAR(2) = N
      IPAR(3) = BLKSZ
      IPAR(4) = LOCLEN
      IPAR(5) = BASISDIM
      IPAR(6) = NPROCS
      IPAR(7) = PROCID
      IPAR(8) = PRECONTYPE
      IPAR(9) = STOPTYPE
      IPAR(10) = MAXIT
      IPAR(11) = -1
      IPAR(12) = -1
      IPAR(13) = -1

      SPAR(1) = EPSILON
      SPAR(2) = -ONE

      RETURN

      END
      SUBROUTINE SMACHCONS(WHAT,RESULT)

* These values are for IEEE-754 arithmetic

*     .. Parameters ..
      REAL MACHEPS
      PARAMETER (MACHEPS=1.19209E-07)
      REAL UNDERFLOW
      PARAMETER (UNDERFLOW=0.11754945E-37)
      REAL OVERFLOW
C*****************
C 98.10.27 BTD
C following line generated complaint from Solaris compiler, with
C OVERFLOW being instead set equal to "IEEE's +Inf" instead
C Therefore try reducing OVERFLOW to 2**127
C      PARAMETER (OVERFLOW=3.40282347E38)
      PARAMETER (OVERFLOW=1.7014118E38)
C****************
*     ..
*     .. Scalar Arguments ..
      REAL RESULT
      CHARACTER WHAT
*     ..
      IF ((WHAT.EQ.'M') .OR. (WHAT.EQ.'m')) THEN
          RESULT = MACHEPS
      ELSE IF ((WHAT.EQ.'U') .OR. (WHAT.EQ.'u')) THEN
          RESULT = UNDERFLOW
      ELSE IF ((WHAT.EQ.'O') .OR. (WHAT.EQ.'o')) THEN
          RESULT = OVERFLOW
      END IF

      RETURN

      END