ROMS
Loading...
Searching...
No Matches
ad_sqlq.F File Reference
#include "cppdefs.h"
Include dependency graph for ad_sqlq.F:

Go to the source code of this file.

Functions/Subroutines

subroutine ad_sqlq (innloop, a, ad_a, tau, ad_tau, y, ad_y)
 
subroutine reclqbs (iflag, kk, innloop, a, tau, y)
 

Function/Subroutine Documentation

◆ ad_sqlq()

subroutine ad_sqlq ( integer, intent(in) innloop,
real(r8), dimension(innloop,innloop), intent(inout) a,
real(r8), dimension(innloop,innloop), intent(inout) ad_a,
real(r8), dimension(innloop), intent(inout) tau,
real(r8), dimension(innloop), intent(inout) ad_tau,
real(r8), dimension(innloop), intent(inout) y,
real(r8), dimension(innloop), intent(inout) ad_y )

Definition at line 6 of file ad_sqlq.F.

7!
8!git $Id$
9!================================================== Hernan G. Arango ===
10! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
11! Licensed under a MIT/X style license !
12! See License_ROMS.md !
13!=======================================================================
14! !
15! This routine performs the adjoint of a LQ factorization of the !
16! square matrix A. !
17! !
18! NOTE: A, ad_A, TAU, ad_TAU, Y and ad_Y are overwritten on exit. !
19! !
20!=======================================================================
21!
22 USE mod_kinds
23!
24 implicit none
25!
26! Imported variable declarations
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! Local variable declarations.
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! Adjoint of a LQ factorization of a square matrix.
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! Save a copy of a.
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!^ tl_a(ii,ii)=tl_zaii
67!^
68 ad_zaii=ad_zaii+ad_a(ii,ii)
69 ad_a(ii,ii)=0.0_r8
70!
71! Recompute basic state arrays.
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!^ tl_a(ii+i,j+jj)=tl_a(ii+i,j+jj)+tl_y(i)*ztemp+ &
90!^ & y(i)*tl_ztemp
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!^ tl_ztemp=-tl_tau(ii)*a(ii,j+jj)-tau(ii)*tl_a(ii,j+jj)
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!^ tl_y(i)=tl_y(i)+tl_ztemp*a(ii+i,j+jj)+ &
106!^ & ztemp*tl_a(ii+i,j+jj)
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!^ tl_ztemp=tl_a(ii,j+jj)
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!^ tl_y(i)=0.0_r8
119 ad_y(i)=0.0_r8
120 END DO
121 END IF
122!^ tl_a(ii,ii)=0.0_r8
123!^
124 ad_a(ii,ii)=0.0_r8
125!^ tl_zaii=tl_a(ii,ii)
126!^
127 ad_a(ii,ii)=ad_a(ii,ii)+ad_zaii
128 ad_zaii=0.0_r8
129!^ tl_a(ii,ii)=tl_zbeta
130!^
131 ad_zbeta=ad_zbeta+ad_a(ii,ii)
132 ad_a(ii,ii)=0.0_r8
133!
134! Recompute basic state arrays.
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!^ tl_a(ii,j)=tl_a(ii,j)/(a(ii,ii)-zbeta)- &
161!^ & (tl_a(ii,ii)-tl_zbeta)*a(ii,j)/(a(ii,ii)-zbeta)
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!^ tl_tau(ii)=(tl_zbeta-tl_a(ii,ii))/zbeta-tl_zbeta*tau(ii)/zbeta
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!^ tl_zbeta=-SIGN(1.0_r8,zbeta)*SIGN(1.0_r8,a(ii,ii))*tl_zbeta
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!^ tl_zbeta=(tl_znorm*znorm+tl_a(ii,ii)*a(ii,ii))/zbetas
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!^ tl_zbeta=0.0_r8
186!^
187 ad_zbeta=0.0_r8
188 END IF
189 IF (znorm.NE.0.0_r8) THEN
190!^ tl_znorm=0.5_r8*tl_znorm/znorm
191!^
192 ad_znorm=0.5_r8*ad_znorm/znorm
193 ELSE
194!^ tl_znorm=0.0_r8
195!^
196 ad_znorm=0.0_r8
197 END IF
198 DO j=ii+1,innloop
199!^ tl_znorm=tl_znorm+2.0_r8*a(ii,j)*tl_a(ii,j)
200!^
201 ad_a(ii,j)=ad_a(ii,j)+2.0_r8*arow(j)*ad_znorm
202 END DO
203!^ tl_znorm=0.0_r8
204!^
205 ad_znorm=0.0_r8
206 END DO
207 DO i=1,innloop
208!^ tl_tau(i)=0.0_r8
209!^
210 ad_tau(i)=0.0_r8
211 END DO
212!
213 RETURN
subroutine reclqbs(iflag, kk, innloop, a, tau, y)
Definition ad_sqlq.F:217

References reclqbs().

Here is the call graph for this function:

◆ reclqbs()

subroutine reclqbs ( integer, intent(in) iflag,
integer, intent(in) kk,
integer, intent(in) innloop,
real(r8), dimension(innloop,innloop), intent(inout) a,
real(r8), dimension(innloop), intent(inout) tau,
real(r8), dimension(innloop), intent(inout) y )

Definition at line 216 of file ad_sqlq.F.

217!
218!***********************************************************************
219! !
220! This routine recomputes the intermediate arrays created during !
221! the LQ factorization of A. !
222! !
223!***********************************************************************
224!
225 USE mod_kinds
226!
227 implicit none
228!
229! Imported variable declarations
230!
231 integer, intent(in) :: innLoop, iflag, kk
232 real(r8), dimension(innLoop,innLoop), intent(inout) :: a
233 real(r8), dimension(innLoop), intent(inout) :: tau, y
234!
235! Local variable declarations.
236!
237 integer :: i, j, ii, jj
238 real(r8) :: znorm, zbeta, zaii, ztemp
239!
240!-----------------------------------------------------------------------
241! Recompute intermiediate arrays for the LQ factorization.
242!-----------------------------------------------------------------------
243!
244 DO i=1,innloop
245 tau(i)=0.0_r8
246 END DO
247!
248 DO ii=1,kk
249!
250! Generate elementary reflector H(ii) to annihilate A(ii,ii+1:innLoop).
251!
252 znorm=0.0_r8
253 DO j=ii+1,innloop
254 znorm=znorm+a(ii,j)*a(ii,j)
255 END DO
256 znorm=sqrt(znorm)
257 zbeta=sqrt(znorm*znorm+a(ii,ii)*a(ii,ii))
258 zbeta=-sign(zbeta,a(ii,ii))
259 tau(ii)=(zbeta-a(ii,ii))/zbeta
260!
261 IF ((iflag.eq.2).and.(ii.eq.kk)) RETURN
262!
263 DO j=ii+1,innloop
264 a(ii,j)=a(ii,j)/(a(ii,ii)-zbeta)
265 END DO
266 a(ii,ii)=zbeta
267!
268! Apply H(ii) to A(ii+1:innLoop,ii:innLoop) from the right.
269!
270 zaii = a(ii,ii)
271 a(ii,ii)=1.0_r8
272
273 IF (tau(ii).ne.0.0_r8) THEN
274 DO i=1,innloop
275 y(i)=0.0_r8
276 END DO
277 jj=ii-1
278 DO j=1,innloop-jj
279 ztemp=a(ii,j+jj)
280 DO i=1,innloop-ii
281 y(i)=y(i)+ztemp*a(ii+i,j+jj)
282 END DO
283 END DO
284!
285 IF ((iflag.eq.1).and.(ii.eq.kk)) RETURN
286!
287 DO j=1,innloop-jj
288 ztemp=-tau(ii)*a(ii,j+jj)
289 DO i=1,innloop-ii
290 a(ii+i,j+jj)=a(ii+i,j+jj)+y(i)*ztemp
291 END DO
292 END DO
293 END IF
294 a(ii,ii)=zaii
295 END DO
296!
297 RETURN

Referenced by ad_sqlq().

Here is the caller graph for this function: