3911
|
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 |