1 |
C $Header: /u/gcmpack/MITgcm/model/src/solve_pentadiagonal.F,v 1.12 2016/09/29 13:02:13 jmc Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "CPP_OPTIONS.h" |
5 |
|
6 |
CBOP |
7 |
C !ROUTINE: SOLVE_PENTADIAGONAL |
8 |
C !INTERFACE: |
9 |
SUBROUTINE SOLVE_PENTADIAGONAL( |
10 |
I iMin,iMax, jMin,jMax, |
11 |
U a5d, b5d, c5d, d5d, e5d, |
12 |
U y5d, |
13 |
O errCode, |
14 |
I bi, bj, myThid ) |
15 |
C !DESCRIPTION: \bv |
16 |
C *==========================================================* |
17 |
C | S/R SOLVE_PENTADIAGONAL |
18 |
C | o Solve a penta-diagonal system A*X=Y (dimension Nr) |
19 |
C *==========================================================* |
20 |
C | o Used to solve implicitly vertical advection & diffusion |
21 |
#ifdef SOLVE_DIAGONAL_LOWMEMORY |
22 |
C | o Allows 2 types of call: |
23 |
C | 1) with INPUT errCode < 0 (first call): |
24 |
C | Solve system A*X=Y by LU decomposition and return |
25 |
C | inverse matrix as output (in a5d,b5d,c5d,d5d,e5d) |
26 |
C | 2) with INPUT errCode = 0 (subsequent calls): |
27 |
C | Solve system A*Xp=Yp by using inverse matrix coeff |
28 |
C | from first call output. |
29 |
#endif /* SOLVE_DIAGONAL_LOWMEMORY */ |
30 |
C *==========================================================* |
31 |
C \ev |
32 |
|
33 |
C !USES: |
34 |
IMPLICIT NONE |
35 |
C == Global data == |
36 |
#include "SIZE.h" |
37 |
#include "EEPARAMS.h" |
38 |
|
39 |
C !INPUT/OUTPUT PARAMETERS: |
40 |
C -- INPUT: -- |
41 |
C iMin,iMax,jMin,jMax :: computational domain |
42 |
C a5d :: 2nd lower diagonal of the pentadiagonal matrix |
43 |
C b5d :: 1rst lower diagonal of the pentadiagonal matrix |
44 |
C c5d :: main diagonal of the pentadiagonal matrix |
45 |
C d5d :: 1rst upper diagonal of the pentadiagonal matrix |
46 |
C e5d :: 2nd upper diagonal of the pentadiagonal matrix |
47 |
C y5d :: Y vector (R.H.S.); |
48 |
C bi,bj :: tile indices |
49 |
C myThid :: thread number |
50 |
C -- OUTPUT: -- |
51 |
C y5d :: X = solution of A*X=Y |
52 |
C errCode :: > 0 if singular matrix |
53 |
#ifdef SOLVE_DIAGONAL_LOWMEMORY |
54 |
C a5d,b5d,c5d,d5d,e5d :: inverse matrix coeff to enable to find Xp solution |
55 |
C of A*Xp=Yp with subsequent call to this routine |
56 |
#endif /* SOLVE_DIAGONAL_LOWMEMORY */ |
57 |
INTEGER iMin,iMax,jMin,jMax |
58 |
_RL a5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
59 |
_RL b5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
60 |
_RL c5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
61 |
_RL d5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
62 |
_RL e5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
63 |
_RL y5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
64 |
INTEGER errCode |
65 |
INTEGER bi, bj, myThid |
66 |
|
67 |
C !LOCAL VARIABLES: |
68 |
C == Local variables == |
69 |
INTEGER i,j,k |
70 |
#ifndef SOLVE_DIAGONAL_LOWMEMORY |
71 |
_RL tmpVar, recVar |
72 |
_RL y5d_m1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
73 |
# ifdef SOLVE_DIAGONAL_KINNER |
74 |
_RL c5d_prime(Nr) |
75 |
_RL d5d_prime(Nr) |
76 |
_RL e5d_prime(Nr) |
77 |
_RL y5d_prime(Nr) |
78 |
_RL y5d_update(Nr) |
79 |
# else |
80 |
_RL c5d_prime(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
81 |
_RL d5d_prime(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
82 |
_RL e5d_prime(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
83 |
_RL y5d_prime(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
84 |
# endif |
85 |
#endif /* SOLVE_DIAGONAL_LOWMEMORY */ |
86 |
CEOP |
87 |
|
88 |
#ifdef SOLVE_DIAGONAL_LOWMEMORY |
89 |
|
90 |
IF ( errCode .LT. 0 ) THEN |
91 |
errCode = 0 |
92 |
|
93 |
DO k=1,Nr |
94 |
C-- forward sweep (starting from top) part-1: matrix L.U decomposition |
95 |
IF (k.EQ.1) THEN |
96 |
DO j=jMin,jMax |
97 |
DO i=iMin,iMax |
98 |
IF ( c5d(i,j,k).NE.0. _d 0 ) THEN |
99 |
c5d(i,j,k) = 1. _d 0 / c5d(i,j,k) |
100 |
ELSE |
101 |
c5d(i,j,k) = 0. _d 0 |
102 |
errCode = 1 |
103 |
ENDIF |
104 |
ENDDO |
105 |
ENDDO |
106 |
|
107 |
ELSEIF (k.EQ.2) THEN |
108 |
DO j=jMin,jMax |
109 |
DO i=iMin,iMax |
110 |
C-- [k] <- [k] - b_k/c_k-1 * [k-1] |
111 |
b5d(i,j,k) = b5d(i,j,k)*c5d(i,j,k-1) |
112 |
c5d(i,j,k) = c5d(i,j,k) - b5d(i,j,k)*d5d(i,j,k-1) |
113 |
d5d(i,j,k) = d5d(i,j,k) - b5d(i,j,k)*e5d(i,j,k-1) |
114 |
IF ( c5d(i,j,k).NE.0. _d 0 ) THEN |
115 |
c5d(i,j,k) = 1. _d 0 / c5d(i,j,k) |
116 |
ELSE |
117 |
c5d(i,j,k) = 0. _d 0 |
118 |
errCode = 1 |
119 |
ENDIF |
120 |
ENDDO |
121 |
ENDDO |
122 |
|
123 |
ELSE |
124 |
C-- Middle of forward sweep |
125 |
DO j=jMin,jMax |
126 |
DO i=iMin,iMax |
127 |
C-- [k] <- [k] - a_k/c_k-2 * [k-2] |
128 |
a5d(i,j,k) = a5d(i,j,k)*c5d(i,j,k-2) |
129 |
b5d(i,j,k) = b5d(i,j,k) - a5d(i,j,k)*d5d(i,j,k-2) |
130 |
c5d(i,j,k) = c5d(i,j,k) - a5d(i,j,k)*e5d(i,j,k-2) |
131 |
C-- [k] <- [k] - b_k/c_k-1 * [k-1] |
132 |
b5d(i,j,k) = b5d(i,j,k)*c5d(i,j,k-1) |
133 |
c5d(i,j,k) = c5d(i,j,k) - b5d(i,j,k)*d5d(i,j,k-1) |
134 |
d5d(i,j,k) = d5d(i,j,k) - b5d(i,j,k)*e5d(i,j,k-1) |
135 |
IF ( c5d(i,j,k).NE.0. _d 0 ) THEN |
136 |
c5d(i,j,k) = 1. _d 0 / c5d(i,j,k) |
137 |
ELSE |
138 |
c5d(i,j,k) = 0. _d 0 |
139 |
errCode = 1 |
140 |
ENDIF |
141 |
ENDDO |
142 |
ENDDO |
143 |
C- end if k= .. ; end of k loop |
144 |
ENDIF |
145 |
ENDDO |
146 |
|
147 |
C-- end if errCode < 0 |
148 |
ENDIF |
149 |
|
150 |
DO k=2,Nr |
151 |
C-- forward sweep (starting from top) part-2: process RHS |
152 |
IF (k.EQ.2) THEN |
153 |
DO j=jMin,jMax |
154 |
DO i=iMin,iMax |
155 |
y5d(i,j,k) = y5d(i,j,k) - b5d(i,j,k)*y5d(i,j,k-1) |
156 |
ENDDO |
157 |
ENDDO |
158 |
ELSE |
159 |
DO j=jMin,jMax |
160 |
DO i=iMin,iMax |
161 |
y5d(i,j,k) = y5d(i,j,k) - b5d(i,j,k)*y5d(i,j,k-1) |
162 |
& - a5d(i,j,k)*y5d(i,j,k-2) |
163 |
ENDDO |
164 |
ENDDO |
165 |
C- end if k= .. ; end of k loop |
166 |
ENDIF |
167 |
ENDDO |
168 |
|
169 |
C-- Backward sweep (starting from bottom) |
170 |
DO k=Nr,1,-1 |
171 |
IF (k.EQ.Nr) THEN |
172 |
DO j=jMin,jMax |
173 |
DO i=iMin,iMax |
174 |
y5d(i,j,k) = y5d(i,j,k)*c5d(i,j,k) |
175 |
ENDDO |
176 |
ENDDO |
177 |
ELSEIF (k.EQ.Nr-1) THEN |
178 |
DO j=jMin,jMax |
179 |
DO i=iMin,iMax |
180 |
y5d(i,j,k) = ( y5d(i,j,k) |
181 |
& - d5d(i,j,k)*y5d(i,j,k+1) |
182 |
& )*c5d(i,j,k) |
183 |
ENDDO |
184 |
ENDDO |
185 |
ELSE |
186 |
DO j=jMin,jMax |
187 |
DO i=iMin,iMax |
188 |
y5d(i,j,k) = ( y5d(i,j,k) |
189 |
& - d5d(i,j,k)*y5d(i,j,k+1) |
190 |
& - e5d(i,j,k)*y5d(i,j,k+2) |
191 |
& )*c5d(i,j,k) |
192 |
ENDDO |
193 |
ENDDO |
194 |
C- end if k= .. ; end of k loop |
195 |
ENDIF |
196 |
ENDDO |
197 |
|
198 |
#else /* ndef SOLVE_DIAGONAL_LOWMEMORY */ |
199 |
|
200 |
errCode = 0 |
201 |
|
202 |
#ifdef SOLVE_DIAGONAL_KINNER |
203 |
C-- Temporary array |
204 |
DO k=1,Nr |
205 |
DO j=1-OLy,sNy+OLy |
206 |
DO i=1-OLx,sNx+OLx |
207 |
y5d_m1(i,j,k) = y5d(i,j,k) |
208 |
ENDDO |
209 |
ENDDO |
210 |
ENDDO |
211 |
|
212 |
C-- Main loop |
213 |
DO j=1-OLy,sNy+OLy |
214 |
DO i=1-OLx,sNx+OLx |
215 |
|
216 |
DO k=1,Nr |
217 |
c5d_prime(k) = 0. _d 0 |
218 |
d5d_prime(k) = 0. _d 0 |
219 |
e5d_prime(k) = 0. _d 0 |
220 |
y5d_prime(k) = 0. _d 0 |
221 |
y5d_update(k) = 0. _d 0 |
222 |
ENDDO |
223 |
|
224 |
DO k=1,Nr |
225 |
C-- forward sweep (starting from top) |
226 |
|
227 |
IF (k.EQ.1) THEN |
228 |
c just copy terms |
229 |
c5d_prime(k) = c5d(i,j,k) |
230 |
d5d_prime(k) = d5d(i,j,k) |
231 |
e5d_prime(k) = e5d(i,j,k) |
232 |
y5d_prime(k) = y5d_m1(i,j,k) |
233 |
ELSEIF (k.EQ.2) THEN |
234 |
c subtract one term |
235 |
c5d_prime(k) = c5d(i,j,k) |
236 |
& -b5d(i,j,k)*d5d_prime(k-1) |
237 |
d5d_prime(k) = d5d(i,j,k) |
238 |
& -b5d(i,j,k)*e5d_prime(k-1) |
239 |
e5d_prime(k) = e5d(i,j,k) |
240 |
y5d_prime(k) = y5d_m1(i,j,k) |
241 |
& -b5d(i,j,k)*y5d_prime(k-1) |
242 |
ELSE |
243 |
c subtract two terms |
244 |
c5d_prime(k) = c5d(i,j,k) |
245 |
& -a5d(i,j,k)*e5d_prime(k-2) |
246 |
& -(b5d(i,j,k)-a5d(i,j,k)*d5d_prime(k-2))*d5d_prime(k-1) |
247 |
d5d_prime(k) = d5d(i,j,k) |
248 |
& -(b5d(i,j,k)-a5d(i,j,k)*d5d_prime(k-2))*e5d_prime(k-1) |
249 |
e5d_prime(k) = e5d(i,j,k) |
250 |
y5d_prime(k) = y5d_m1(i,j,k) |
251 |
& -a5d(i,j,k)*y5d_prime(k-2) |
252 |
& -(b5d(i,j,k)-a5d(i,j,k)*d5d_prime(k-2))*y5d_prime(k-1) |
253 |
ENDIF |
254 |
|
255 |
c normalization |
256 |
tmpVar = c5d_prime(k) |
257 |
IF ( tmpVar.NE.0. _d 0 ) THEN |
258 |
recVar = 1. _d 0 / tmpVar |
259 |
d5d_prime(k) = d5d_prime(k)*recVar |
260 |
e5d_prime(k) = e5d_prime(k)*recVar |
261 |
y5d_prime(k) = y5d_prime(k)*recVar |
262 |
ELSE |
263 |
d5d_prime(k) = 0. _d 0 |
264 |
e5d_prime(k) = 0. _d 0 |
265 |
y5d_prime(k) = 0. _d 0 |
266 |
errCode = 1 |
267 |
ENDIF |
268 |
|
269 |
C-- k loop |
270 |
ENDDO |
271 |
|
272 |
C-- Backward sweep (starting from bottom) |
273 |
DO k=Nr,1,-1 |
274 |
IF (k.EQ.Nr) THEN |
275 |
y5d_update(k) = y5d_prime(k) |
276 |
ELSEIF (k.EQ.Nr-1) THEN |
277 |
y5d_update(k) = y5d_prime(k) |
278 |
& - y5d_update(k+1)*d5d_prime(k) |
279 |
ELSE |
280 |
y5d_update(k) = y5d_prime(k) |
281 |
& - y5d_update(k+1)*d5d_prime(k) |
282 |
& - y5d_update(k+2)*e5d_prime(k) |
283 |
ENDIF |
284 |
ENDDO |
285 |
|
286 |
C-- Update array |
287 |
DO k=1,Nr |
288 |
y5d(i,j,k) = y5d_update(k) |
289 |
ENDDO |
290 |
|
291 |
C- end i,j loop |
292 |
ENDDO |
293 |
ENDDO |
294 |
|
295 |
#else /* ndef SOLVE_DIAGONAL_KINNER */ |
296 |
|
297 |
C-- Init. + copy to temp. array |
298 |
DO k=1,Nr |
299 |
DO j=1-OLy,sNy+OLy |
300 |
DO i=1-OLx,sNx+OLx |
301 |
c5d_prime(i,j,k) = 0. _d 0 |
302 |
d5d_prime(i,j,k) = 0. _d 0 |
303 |
e5d_prime(i,j,k) = 0. _d 0 |
304 |
y5d_prime(i,j,k) = 0. _d 0 |
305 |
y5d_m1(i,j,k) = y5d(i,j,k) |
306 |
ENDDO |
307 |
ENDDO |
308 |
ENDDO |
309 |
|
310 |
CADJ loop = sequential |
311 |
DO k=1,Nr |
312 |
C-- Forward sweep (starting from top) |
313 |
|
314 |
DO j=1-OLy,sNy+OLy |
315 |
DO i=1-OLx,sNx+OLx |
316 |
c DO j=jMin,jMax |
317 |
c DO i=iMin,iMax |
318 |
|
319 |
IF (k.EQ.1) THEN |
320 |
C- just copy terms |
321 |
c5d_prime(i,j,k) = c5d(i,j,k) |
322 |
d5d_prime(i,j,k) = d5d(i,j,k) |
323 |
e5d_prime(i,j,k) = e5d(i,j,k) |
324 |
y5d_prime(i,j,k) = y5d_m1(i,j,k) |
325 |
ELSEIF (k.EQ.2) THEN |
326 |
C- subtract one term |
327 |
c5d_prime(i,j,k) = c5d(i,j,k) |
328 |
& -b5d(i,j,k)*d5d_prime(i,j,k-1) |
329 |
d5d_prime(i,j,k) = d5d(i,j,k) |
330 |
& -b5d(i,j,k)*e5d_prime(i,j,k-1) |
331 |
e5d_prime(i,j,k) = e5d(i,j,k) |
332 |
y5d_prime(i,j,k) = y5d_m1(i,j,k) |
333 |
& -b5d(i,j,k)*y5d_prime(i,j,k-1) |
334 |
ELSE |
335 |
C- subtract two terms |
336 |
c5d_prime(i,j,k) = c5d(i,j,k) |
337 |
& -a5d(i,j,k)*e5d_prime(i,j,k-2) |
338 |
& -(b5d(i,j,k)-a5d(i,j,k)*d5d_prime(i,j,k-2)) |
339 |
& *d5d_prime(i,j,k-1) |
340 |
d5d_prime(i,j,k) = d5d(i,j,k) |
341 |
& -(b5d(i,j,k)-a5d(i,j,k)*d5d_prime(i,j,k-2)) |
342 |
& *e5d_prime(i,j,k-1) |
343 |
e5d_prime(i,j,k) = e5d(i,j,k) |
344 |
y5d_prime(i,j,k) = y5d_m1(i,j,k) |
345 |
& -a5d(i,j,k)*y5d_prime(i,j,k-2) |
346 |
& -(b5d(i,j,k)-a5d(i,j,k)*d5d_prime(i,j,k-2)) |
347 |
& *y5d_prime(i,j,k-1) |
348 |
ENDIF |
349 |
|
350 |
C- normalization |
351 |
tmpVar = c5d_prime(i,j,k) |
352 |
IF ( tmpVar.NE.0. _d 0 ) THEN |
353 |
recVar = 1. _d 0 / tmpVar |
354 |
d5d_prime(i,j,k) = d5d_prime(i,j,k)*recVar |
355 |
e5d_prime(i,j,k) = e5d_prime(i,j,k)*recVar |
356 |
y5d_prime(i,j,k) = y5d_prime(i,j,k)*recVar |
357 |
ELSE |
358 |
d5d_prime(i,j,k) = 0. _d 0 |
359 |
e5d_prime(i,j,k) = 0. _d 0 |
360 |
y5d_prime(i,j,k) = 0. _d 0 |
361 |
errCode = 1 |
362 |
ENDIF |
363 |
|
364 |
ENDDO |
365 |
ENDDO |
366 |
|
367 |
C- end k-loop |
368 |
ENDDO |
369 |
|
370 |
CADJ loop = sequential |
371 |
C-- Backward sweep (starting from bottom) |
372 |
DO k=Nr,1,-1 |
373 |
DO j=1-OLy,sNy+OLy |
374 |
DO i=1-OLx,sNx+OLx |
375 |
c DO j=jMin,jMax |
376 |
c DO i=iMin,iMax |
377 |
IF (k.EQ.Nr) THEN |
378 |
y5d(i,j,k) = y5d_prime(i,j,k) |
379 |
ELSEIF (k.EQ.Nr-1) THEN |
380 |
y5d(i,j,k) = y5d_prime(i,j,k) |
381 |
& - y5d(i,j,k+1)*d5d_prime(i,j,k) |
382 |
ELSE |
383 |
y5d(i,j,k) = y5d_prime(i,j,k) |
384 |
& - y5d(i,j,k+1)*d5d_prime(i,j,k) |
385 |
& - y5d(i,j,k+2)*e5d_prime(i,j,k) |
386 |
ENDIF |
387 |
ENDDO |
388 |
ENDDO |
389 |
C- end k-loop |
390 |
ENDDO |
391 |
|
392 |
#endif /* SOLVE_DIAGONAL_KINNER */ |
393 |
|
394 |
#endif /* SOLVE_DIAGONAL_LOWMEMORY */ |
395 |
|
396 |
RETURN |
397 |
END |