share/odepack/fortran/dorthog.f

   1 *DECK DORTHOG
   2       SUBROUTINE DORTHOG (VNEW, V, HES, N, LL, LDHES, KMP, SNORMW)
   3       INTEGER N, LL, LDHES, KMP
   4       DOUBLE PRECISION VNEW, V, HES, SNORMW
   5       DIMENSION VNEW(*), V(N,*), HES(LDHES,*)
   6 C-----------------------------------------------------------------------
   7 C This routine orthogonalizes the vector VNEW against the previous
   8 C KMP vectors in the V array.  It uses a modified Gram-Schmidt
   9 C orthogonalization procedure with conditional reorthogonalization.
  10 C This is the version of 28 may 1986.
  11 C-----------------------------------------------------------------------
  12 C
  13 C      On entry
  14 C
  15 C         VNEW = the vector of length N containing a scaled product
  16 C                of the Jacobian and the vector V(*,LL).
  17 C
  18 C         V    = the N x l array containing the previous LL
  19 C                orthogonal vectors v(*,1) to v(*,LL).
  20 C
  21 C         HES  = an LL x LL upper Hessenberg matrix containing,
  22 C                in HES(i,k), k.lt.LL, scaled inner products of
  23 C                A*V(*,k) and V(*,i).
  24 C
  25 C        LDHES = the leading dimension of the HES array.
  26 C
  27 C         N    = the order of the matrix A, and the length of VNEW.
  28 C
  29 C         LL   = the current order of the matrix HES.
  30 C
  31 C          KMP = the number of previous vectors the new vector VNEW
  32 C                must be made orthogonal to (KMP .le. MAXL).
  33 C
  34 C
  35 C      On return
  36 C
  37 C         VNEW = the new vector orthogonal to V(*,i0) to V(*,LL),
  38 C                where i0 = MAX(1, LL-KMP+1).
  39 C
  40 C         HES  = upper Hessenberg matrix with column LL filled in with
  41 C                scaled inner products of A*V(*,LL) and V(*,i).
  42 C
  43 C       SNORMW = L-2 norm of VNEW.
  44 C
  45 C-----------------------------------------------------------------------
  46       INTEGER I, I0
  47       DOUBLE PRECISION ARG, DDOT, DNRM2, SUMDSQ, TEM, VNRM
  48 C
  49 C Get norm of unaltered VNEW for later use. ----------------------------
  50       VNRM = DNRM2 (N, VNEW, 1)
  51 C-----------------------------------------------------------------------
  52 C Do modified Gram-Schmidt on VNEW = A*v(LL).
  53 C Scaled inner products give new column of HES.
  54 C Projections of earlier vectors are subtracted from VNEW.
  55 C-----------------------------------------------------------------------
  56       I0 = MAX(1,LL-KMP+1)
  57       DO 10 I = I0,LL
  58         HES(I,LL) = DDOT (N, V(1,I), 1, VNEW, 1)
  59         TEM = -HES(I,LL)
  60         CALL DAXPY (N, TEM, V(1,I), 1, VNEW, 1)
  61  10     CONTINUE
  62 C-----------------------------------------------------------------------
  63 C Compute SNORMW = norm of VNEW.
  64 C If VNEW is small compared to its input value (in norm), then
  65 C reorthogonalize VNEW to V(*,1) through V(*,LL).
  66 C Correct if relative correction exceeds 1000*(unit roundoff).
  67 C finally, correct SNORMW using the dot products involved.
  68 C-----------------------------------------------------------------------
  69       SNORMW = DNRM2 (N, VNEW, 1)
  70       IF (VNRM + 0.001D0*SNORMW .NE. VNRM) RETURN
  71       SUMDSQ = 0.0D0
  72       DO 30 I = I0,LL
  73         TEM = -DDOT (N, V(1,I), 1, VNEW, 1)
  74         IF (HES(I,LL) + 0.001D0*TEM .EQ. HES(I,LL)) GO TO 30
  75         HES(I,LL) = HES(I,LL) - TEM
  76         CALL DAXPY (N, TEM, V(1,I), 1, VNEW, 1)
  77         SUMDSQ = SUMDSQ + TEM**2
  78  30     CONTINUE
  79       IF (SUMDSQ .EQ. 0.0D0) RETURN
  80       ARG = MAX(0.0D0,SNORMW**2 - SUMDSQ)
  81       SNORMW = SQRT(ARG)
  82 C
  83       RETURN
  84 C----------------------- End of Subroutine DORTHOG ---------------------
  85       END