7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
23
24 implicit none
25
26
27
28 integer, intent(in) :: innLoop
29
30 real(r8), dimension(innLoop,innLoop), intent(inout) :: a, ad_a
31 real(r8), dimension(innLoop), intent(inout) :: tau, ad_tau
32 real(r8), dimension(innLoop), intent(inout) :: y, ad_y
33
34
35
36 integer :: i, j, ii, jj, m, kk, iflag
37
38 real(r8) :: znorm, zbeta, zaii, ztemp, zbetas
39 real(r8) :: adfac, ad_znorm, ad_zbeta, ad_zaii, ad_ztemp
40
41 real(r8), dimension(innLoop,innLoop) :: as
42 real(r8), dimension(innLoop) :: arow
43
44
45
46
47
48 ad_znorm=0.0_r8
49 ad_zbeta=0.0_r8
50 ad_zaii=0.0_r8
51 ad_ztemp=0.0_r8
52 DO i=1,innloop
53 ad_y(i)=0.0_r8
54 arow(i)=0.0_r8
55 END DO
56
57
58
59 DO j=1,innloop
60 DO i=1,innloop
61 as(i,j)=a(i,j)
62 END DO
63 END DO
64
65 DO ii=innloop-1,1,-1
66
67
68 ad_zaii=ad_zaii+ad_a(ii,ii)
69 ad_a(ii,ii)=0.0_r8
70
71
72
73 DO j=1,innloop
74 DO i=1,innloop
75 a(i,j)=as(i,j)
76 END DO
77 END DO
78
79 iflag=1
80 kk=ii
81 m=innloop
82 call reclqbs (iflag, kk, m, a, tau, y)
83
84 IF (tau(ii).ne.0.0_r8) THEN
85 jj=ii-1
86 DO j=1,innloop-jj
87 ztemp=-tau(ii)*a(ii,j+jj)
88 DO i=1,innloop-ii
89
90
91
92 ad_ztemp=ad_ztemp+y(i)*ad_a(ii+i,j+jj)
93 ad_y(i)=ad_y(i)+ztemp*ad_a(ii+i,j+jj)
94 END DO
95
96
97 ad_tau(ii)=ad_tau(ii)-a(ii,j+jj)*ad_ztemp
98 ad_a(ii,j+jj)=ad_a(ii,j+jj)-tau(ii)*ad_ztemp
99 ad_ztemp=0.0_r8
100 END DO
101
102 DO j=1,innloop-jj
103 ztemp=a(ii,j+jj)
104 DO i=1,innloop-ii
105
106
107
108 ad_ztemp=ad_ztemp+a(ii+i,j+jj)*ad_y(i)
109 ad_a(ii+i,j+jj)=ad_a(ii+i,j+jj)+ztemp*ad_y(i)
110 END DO
111
112
113 ad_a(ii,j+jj)=ad_a(ii,j+jj)+ad_ztemp
114 ad_ztemp=0.0_r8
115 END DO
116
117 DO i=1,innloop
118
119 ad_y(i)=0.0_r8
120 END DO
121 END IF
122
123
124 ad_a(ii,ii)=0.0_r8
125
126
127 ad_a(ii,ii)=ad_a(ii,ii)+ad_zaii
128 ad_zaii=0.0_r8
129
130
131 ad_zbeta=ad_zbeta+ad_a(ii,ii)
132 ad_a(ii,ii)=0.0_r8
133
134
135
136 DO j=1,innloop
137 DO i=1,innloop
138 a(i,j)=as(i,j)
139 END DO
140 END DO
141
142 iflag=2
143 kk=ii
144 m=innloop
145 call reclqbs (iflag, kk, m, a, tau, y)
146
147 znorm=0.0_r8
148 DO j=ii+1,innloop
149 znorm=znorm+a(ii,j)*a(ii,j)
150 END DO
151 znorm=sqrt(znorm)
152 zbeta=sqrt(znorm*znorm+a(ii,ii)*a(ii,ii))
153 zbetas=zbeta
154 zbeta=-sign(zbeta,a(ii,ii))
155 tau(ii)=(zbeta-a(ii,ii))/zbeta
156
157 DO j=ii+1,innloop
158 arow(j)=a(ii,j)
159 a(ii,j)=a(ii,j)/(a(ii,ii)-zbeta)
160
161
162
163 adfac=ad_a(ii,j)/(a(ii,ii)-zbeta)
164 ad_zbeta=ad_zbeta+a(ii,j)*adfac
165 ad_a(ii,ii)=ad_a(ii,ii)-a(ii,j)*adfac
166 ad_a(ii,j )=adfac
167 END DO
168
169
170 adfac=ad_tau(ii)/zbeta
171 ad_zbeta=ad_zbeta+adfac-tau(ii)*adfac
172 ad_a(ii,ii)=ad_a(ii,ii)-adfac
173 ad_tau(ii)=0.0_r8
174
175
176 ad_zbeta=-sign(1.0_r8,zbetas)*sign(1.0_r8,a(ii,ii))*ad_zbeta
177 IF (zbetas.NE.0.0_r8) THEN
178
179
180 adfac=ad_zbeta/zbetas
181 ad_znorm=ad_znorm+znorm*adfac
182 ad_a(ii,ii)=ad_a(ii,ii)+a(ii,ii)*adfac
183 ad_zbeta=0.0_r8
184 ELSE
185
186
187 ad_zbeta=0.0_r8
188 END IF
189 IF (znorm.NE.0.0_r8) THEN
190
191
192 ad_znorm=0.5_r8*ad_znorm/znorm
193 ELSE
194
195
196 ad_znorm=0.0_r8
197 END IF
198 DO j=ii+1,innloop
199
200
201 ad_a(ii,j)=ad_a(ii,j)+2.0_r8*arow(j)*ad_znorm
202 END DO
203
204
205 ad_znorm=0.0_r8
206 END DO
207 DO i=1,innloop
208
209
210 ad_tau(i)=0.0_r8
211 END DO
212
213 RETURN
subroutine reclqbs(iflag, kk, innloop, a, tau, y)