GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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
subroutine dorth(VNEW, V, HES, N, LL, LDHES, KMP, SNORMW)
Definition: dorth.f:5
std::string dimension(void) const
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:233
octave_value sqrt(void) const
Definition: ov.h:1200