GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dorth.f
Go to the documentation of this file.
1 C Work performed under the auspices of the U.S. Department of Energy
2 C by Lawrence Livermore National Laboratory under contract number
3 C W-7405-Eng-48.
4 C
5  SUBROUTINE dorth (VNEW, V, HES, N, LL, LDHES, KMP, SNORMW)
6 C
7 C***BEGIN PROLOGUE DORTH
8 C***DATE WRITTEN 890101 (YYMMDD)
9 C***REVISION DATE 900926 (YYMMDD)
10 C
11 C
12 C-----------------------------------------------------------------------
13 C***DESCRIPTION
14 C
15 C This routine orthogonalizes the vector VNEW against the previous
16 C KMP vectors in the V array. It uses a modified Gram-Schmidt
17 C orthogonalization procedure with conditional reorthogonalization.
18 C
19 C On entry
20 C
21 C VNEW = The vector of length N containing a scaled product
22 C OF The Jacobian and the vector V(*,LL).
23 C
24 C V = The N x LL array containing the previous LL
25 C orthogonal vectors V(*,1) to V(*,LL).
26 C
27 C HES = An LL x LL upper Hessenberg matrix containing,
28 C in HES(I,K), K.LT.LL, scaled inner products of
29 C A*V(*,K) and V(*,I).
30 C
31 C LDHES = The leading dimension of the HES array.
32 C
33 C N = The order of the matrix A, and the length of VNEW.
34 C
35 C LL = The current order of the matrix HES.
36 C
37 C KMP = The number of previous vectors the new vector VNEW
38 C must be made orthogonal to (KMP .LE. MAXL).
39 C
40 C
41 C On return
42 C
43 C VNEW = The new vector orthogonal to V(*,I0),
44 C where I0 = MAX(1, LL-KMP+1).
45 C
46 C HES = Upper Hessenberg matrix with column LL filled in with
47 C scaled inner products of A*V(*,LL) and V(*,I).
48 C
49 C SNORMW = L-2 norm of VNEW.
50 C
51 C-----------------------------------------------------------------------
52 C***ROUTINES CALLED
53 C DDOT, DNRM2, DAXPY
54 C
55 C***END PROLOGUE DORTH
56 C
57  INTEGER N, LL, LDHES, KMP
58  DOUBLE PRECISION VNEW, V, HES, SNORMW
59  dimension vnew(*), v(n,*), hes(ldhes,*)
60  INTEGER I, I0
61  DOUBLE PRECISION ARG, DDOT, DNRM2, SUMDSQ, TEM, VNRM
62 C
63 C-----------------------------------------------------------------------
64 C Get norm of unaltered VNEW for later use.
65 C-----------------------------------------------------------------------
66  vnrm = dnrm2(n, vnew, 1)
67 C-----------------------------------------------------------------------
68 C Do Modified Gram-Schmidt on VNEW = A*V(LL).
69 C Scaled inner products give new column of HES.
70 C Projections of earlier vectors are subtracted from VNEW.
71 C-----------------------------------------------------------------------
72  i0 = max0(1,ll-kmp+1)
73  DO 10 i = i0,ll
74  hes(i,ll) = ddot(n, v(1,i), 1, vnew, 1)
75  tem = -hes(i,ll)
76  CALL daxpy (n, tem, v(1,i), 1, vnew, 1)
77  10 CONTINUE
78 C-----------------------------------------------------------------------
79 C Compute SNORMW = norm of VNEW.
80 C If VNEW is small compared to its input value (in norm), then
81 C Reorthogonalize VNEW to V(*,1) through V(*,LL).
82 C Correct if relative correction exceeds 1000*(unit roundoff).
83 C Finally, correct SNORMW using the dot products involved.
84 C-----------------------------------------------------------------------
85  snormw = dnrm2(n, vnew, 1)
86  IF (vnrm + 0.001d0*snormw .NE. vnrm) RETURN
87  sumdsq = 0.0d0
88  DO 30 i = i0,ll
89  tem = -ddot(n, v(1,i), 1, vnew, 1)
90  IF (hes(i,ll) + 0.001d0*tem .EQ. hes(i,ll)) GO TO 30
91  hes(i,ll) = hes(i,ll) - tem
92  CALL daxpy (n, tem, v(1,i), 1, vnew, 1)
93  sumdsq = sumdsq + tem**2
94  30 CONTINUE
95  IF (sumdsq .EQ. 0.0d0) RETURN
96  arg = max(0.0d0,snormw**2 - sumdsq)
97  snormw = sqrt(arg)
98  RETURN
99 C
100 C------END OF SUBROUTINE DORTH------------------------------------------
101  END
max
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:230
dorth
subroutine dorth(VNEW, V, HES, N, LL, LDHES, KMP, SNORMW)
Definition: dorth.f:6