SUBROUTINE ADIP ( A, B, C, D, IOPT, N, ITNS, DMU, OMEH,
     &  OMEV, IER )

c*********************************************************************72
c
C  GIVEN A MATRIX EQUATION GZ=S, WHERE G IS A REAL POSITIVE
C  DEFINITE MATRIX, S IS A KNOWN, AND Z THE UNKNOWN, VECTOR,
C  LET H AND V BE SYMMETRIC COMMUTING MATRICES SUCH THAT
C  G=H+V.  BEGINNING WITH AN INITIAL APPROXIMATION Z(0), LET
C  Z(K+1/2)=(H+OMEH(K)*I)**(-1)*(S-(V-OMEH(K)*I)*Z(K),
C  Z(K+1)  =(V+OMEGV(K)*I)**(-1)*(S-(H-OMEV(K)*I)*Z(K+1/2),
C  WHERE I IS THE IDENTITY MATRIX.  FINALLY, LET ONE OF ITNS
C  AND DMU BE GIVEN.  THEN THIS SUBROUTINE COMPUTES THE PAR-
C  AMETERS OMEH(K), OMEV(K) THAT MINIMIZE THE VALUE OF 
C  DMU AND ITNS WHICH IS NOT GIVEN WHILE SATISFYING THE
C  INEQUALITY /Z-Z(ITNS)/.LE.DMU*/Z-Z(0)/, WHERE // DENOTES
C  THE EUCLIDEAN NORM.
C  THE SUBROUTINE ARGUMENTS HAVE THE FOLOWING MEANING.
C  A AND D ARE LOWER AND UPPER BOUNDS, RESPECTIVELY, ON THE
C  EIGENVALUES OF H.  C AND D ARE LOWER AND UPPER BOUNDS,
C  RESPECTIVELY, ON THE EIGENVALUES OF V.  THE VALUES OF A, B,
C  C AND D MUST SATISFY THE INEQUALITIES 0.LT.A.LE.B AND
C  0.LT.C.LE.D.
C  IOPT DENOTES THE INPUT OPTION.  IF IOPT=1, THEN THE VALUE OF
C  ITNS MUST BE SPECIFIED ON ENTRY AND DMU WILL BE COMPUTED.
C  IF IOPT = 2 THEN THE VALULE OF DMU MUST BE SPECIFIED ON ENTRY
C  AND ITNS WILL BE COMPUTED.
C  IF IOPT=1 THEN THE INEQULAITY 1.LE.ITNS.LE.N MUST BE
C  SATISFIED, WHILE IF IOPT=2 THEN THE INEQUALITIES N.GE.1
C  AND DMU.GT.0 MUST BE SATISFIED.
C  N IS THE DIMENSION OF THE ARRAYS OMEV AND OMEH.
C  ITNS IS THE NUMBER OF ITERATIONS TO BE PERFORMED.
C  DMU IS A BOUND ON THE SPECTRAL NORM OF THE ITERATION
C  MATRIX TO THE ITNS POWER.  IF IOPT=2 A VALUE FOR DMU MUST
C  BE SPECIFIED ON ENTRY, AND THIS VALUE MAY BE CHANGED BY
C  ADIP (SEE IER, BELOW).
C  THE VALUES OF THE REQUIRED PARAMETERS ARE CONTAINED IN THE
C  LOCATIONS OMEV(K), OMEH(K), K=1,...,ITNS ON EXIT FROM
C  ADIP.
C  IER IS A VARIABLE WHOSE VALUE ON EXIT FROM ADIP HAS THE
C  FOLLOWING MEANING
C  IER=0 SIGNIFIES COMPUTATION OF THE PARAMETERS HAS BEEN
C  PERFORMED WITH NO CHANGED OF THE VALUES SPECIFIED ON ENTRY.
C  IER=1 SIGNIFIES THAT SOME INPUT VALUE VIOLATES THE
C  CONSTRAINTS GIVEN ABOVE, AND HENCE THE PARAMETERS HAVE NOT
C  BEEN COMPUTED.
C  IER=2 (POSSIBLE ONLY IF IOPT=2) SIGNIFIES THAT FOR THE
C  GIVEN VALUE OF DMU, THE COMPUTED VALUE OF ITNS WOULD BE
C  GREATER THAN N, SO THAT ITNS HAS BEEN SET EQUAL TO N AND
C  DMU HAS BEEN RECOMPUTED AS FOR IOPT=1.
C  TEST INPUT VALUES FOR RANGE.
c
      DOUBLE PRECISION A, ALFA, B, BETA, BMD, BPD, C, CMA, CPA,
     &  D, DEL, DEXP, DKPR, DLOG, DM, DMU, DRJ, DSQRT, OJ,
     &  OMEH(N), OMEV(N), PISQ, TEMP, TEMPA, TEMPB, TEMPC
      DATA PISQ / 9.869604401089359D0 /

      IER = 1
      IF ( .NOT. ( A.GT.0.D0 .AND. A.LE.B .AND. C.GT.0.D0 .AND.
     &  C .LE. D ) ) GO TO 90
      IF ( .NOT. ( IOPT .EQ. 1 .OR. IOPT .EQ. 2 ) ) GO TO 90
      IF ( IOPT .EQ. 2 ) GO TO 10
      IF ( .NOT. ( ITNS .GE. 1 .AND. ITNS .LE. N ) ) GO TO 90
      GO TO 20
10    IF ( .NOT. (N.GE.1 .AND. ITNS.LE.N ) ) GO TO 90
C  STATE 1 - PRELIMINARY COMPUTATIONS COMMON TO BOTH OPTIONS.
20    IER = 0
      BPD = B + D
      BMD = B - D
      CPA = C + A
      CMA = C - A
      DM = 2.D0 * ( ( D - C ) * ( B - A ) ) / ( BPD * CPA )
      DKPR = 1.D0 / ( 1.D0 + DM + DSQRT ( DM * ( DM + 2.D0 ) ) )
      DEL = 0.D0
      IF ( BMD .EQ. 0.D0 .AND. CMA .EQ. 0.D0 ) GO TO 30
      TEMP = BPD * DKPR
      DEL = 2.D0 * ( TEMP - CPA ) / ( CPA * BMD + TEMP * CMA )
30    ALFA = DKPR * ( CMA + 2.D0 * DEL * A * C ) / CPA
      BETA = ( 2.D0 + DEL * BMD ) / BPD
      TEMP = DKPR / 4.D0
C  END OF STAGE 1 - COMPUTE ITNS FOR OPTION 2.
      IF ( IOPT .EQ. 1 ) GO TO 40
      ITNS = ( DLOG ( DMU / 4.D0 ) * DLOG ( TEMP ) ) / PISQ + 1.D0
      IF ( ITNS .LE. N ) GO TO 40
      ITNS = N
      IER = 2
C  STAGE 2 - COMPUTATION OF THE OPTIMAL PARAMETERS.
40    TEMPA = 2 * ITNS
      TEMPB = TEMP * TEMP
      DO 50 J = 1, ITNS
        DRJ = 2 * J - 1
        DRJ = DRJ / TEMPA
        TEMPC = TEMP**DRJ
        OJ = 2.D0 + ( TEMPC+TEMPB/TEMPC ) / ( 1.D0 + TEMPC * TEMPC )
        TEMPC = DEL + OJ
        OMEV(J) = ( OJ - ALFA ) / ( BETA - TEMPC )
        OMEH(J) = ( OJ + ALFA ) / ( BETA + TEMPC )
50    CONTINUE
      IF ( IOPT .EQ. 2 .AND. IER .EQ. 0 ) GO TO 90
C  END OF STAGE 2 - COMPUTE DMU FOR OPTION 1.
      TEMPA = ITNS
      TEMP = PISQ * TEMPA / DLOG ( TEMPB * ( 1.D0 + 8.D0 * TEMPB ) )
C  CHOOSE PROPER FORMULA TO AVOID UNDERFLOW OR OVERFLOW.
      IF ( TEMP .LE. -90.D0 .OR. TEMP .GE. 30.D0 ) GO TO 60
      IF ( TEMP .LE. -10.D0 ) GO TO 70
      IF ( TEMP .LT. 10.D0 ) GO TO 80
      DMU = DEXP ( -6.D0 * TEMP )
      GO TO 90
60    DMU = 0.D0
      GO TO 90
70    DMU = 4.D0 * DEXP ( 2.0 * TEMP )
      GO TO 90
80    TEMP = DEXP ( TEMP )
      DMU = ( ( 2.D0 * TEMP ) / ( 1.D0 + 2.0D0 * TEMP**4 ) )**2
90    RETURN
      END