/[MITgcm]/MITgcm/model/src/solve_pentadiagonal.F
ViewVC logotype

Diff of /MITgcm/model/src/solve_pentadiagonal.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.3 by jmc, Tue Mar 16 00:08:27 2010 UTC revision 1.4 by gforget, Tue Aug 10 17:58:30 2010 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6    C o Switch to code that has the k-loop inside the
7    C   ij-loops, which matters in adjoint mode.
8    #ifdef ALLOW_AUTODIFF
9    #define ALLOW_SOLVERS_KLOOPINSIDE
10    #endif
11    
12  CBOP  CBOP
13  C     !ROUTINE: SOLVE_PENTADIAGONAL  C     !ROUTINE: SOLVE_PENTADIAGONAL
14  C     !INTERFACE:  C     !INTERFACE:
# Line 57  C     errCode :: > 0 if singular matrix Line 63  C     errCode :: > 0 if singular matrix
63  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
64  C     == Local variables ==  C     == Local variables ==
65        INTEGER i,j,k        INTEGER i,j,k
66    #ifdef ALLOW_SOLVERS_KLOOPINSIDE
67          _RL y5d_m1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
68          _RL a5d_prime(Nr), b5d_prime(Nr)
69          _RL c5d_prime(Nr), d5d_prime(Nr), e5d_prime(Nr)
70          _RL y5d_prime(Nr), y5d_update(Nr), tmpval
71    #endif
72  CEOP  CEOP
73    
74    #ifndef ALLOW_SOLVERS_KLOOPINSIDE
75    
76        errCode = 0        errCode = 0
77    
78        DO k=1,Nr        DO k=1,Nr
# Line 149  C-      end if k= .. ; end of k loop Line 163  C-      end if k= .. ; end of k loop
163         ENDIF         ENDIF
164        ENDDO        ENDDO
165    
166    #else  /* ALLOW_SOLVERS_KLOOPINSIDE */
167    
168          errCode = 0
169    
170    C--   Temporary array
171          DO j=jMin,jMax
172          DO i=iMin,iMax
173          DO k=1,Nr
174             y5d_m1(i,j,k) = y5d(i,j,k,bi,bj)
175          ENDDO
176          ENDDO
177          ENDDO
178    
179    C--   Main loop
180          DO j=jMin,jMax
181          DO i=iMin,iMax
182    
183          DO k=1,Nr
184            a5d_prime(k) = 0. _d 0
185            b5d_prime(k) = 0. _d 0
186            c5d_prime(k) = 0. _d 0
187            d5d_prime(k) = 0. _d 0
188            e5d_prime(k) = 0. _d 0
189            y5d_prime(k) = 0. _d 0
190            y5d_update(k) = 0. _d 0
191          ENDDO
192    
193          DO k=1,Nr
194    C--   forward sweep (starting from top)
195    
196              IF (k.EQ.1) THEN
197    c just copy terms
198               a5d_prime(k) = 0. _d 0
199               b5d_prime(k) = 0. _d 0
200               c5d_prime(k) = c5d(i,j,k)
201               d5d_prime(k) = d5d(i,j,k)
202               e5d_prime(k) = e5d(i,j,k)
203               y5d_prime(k) = y5d_m1(i,j,k)
204              ELSEIF (k.EQ.2) THEN
205    c subtract one term
206               a5d_prime(k) = 0. _d 0
207               b5d_prime(k) = 0. _d 0
208               c5d_prime(k) = c5d(i,j,k)
209         &      -b5d(i,j,k)*d5d_prime(k-1)
210               d5d_prime(k) = d5d(i,j,k)
211         &      -b5d(i,j,k)*e5d_prime(k-1)
212               e5d_prime(k) = e5d(i,j,k)
213               y5d_prime(k) = y5d_m1(i,j,k)
214         &      -b5d(i,j,k)*y5d_prime(k-1)
215              ELSE
216    c subtract two terms
217               a5d_prime(k) = 0. _d 0
218               b5d_prime(k) = 0. _d 0
219               c5d_prime(k) = c5d(i,j,k)
220         &      -a5d(i,j,k)*e5d_prime(k-2)
221         &      -(b5d(i,j,k)-a5d(i,j,k)*d5d_prime(k-2))*d5d_prime(k-1)
222               d5d_prime(k) = d5d(i,j,k)
223         &      -(b5d(i,j,k)-a5d(i,j,k)*d5d_prime(k-2))*e5d_prime(k-1)
224               e5d_prime(k) = e5d(i,j,k)
225               y5d_prime(k) = y5d_m1(i,j,k)
226         &      -a5d(i,j,k)*y5d_prime(k-2)
227         &      -(b5d(i,j,k)-a5d(i,j,k)*d5d_prime(k-2))*y5d_prime(k-1)
228              ENDIF
229    
230    c normalization
231              tmpval=c5d_prime(k)
232              IF ( tmpval.NE.0. _d 0 ) THEN
233               a5d_prime(k) = a5d_prime(k) / tmpval
234               b5d_prime(k) = b5d_prime(k) / tmpval
235               c5d_prime(k) = 1. _d 0
236               d5d_prime(k) = d5d_prime(k) / tmpval
237               e5d_prime(k) = e5d_prime(k) / tmpval
238               y5d_prime(k) = y5d_prime(k) / tmpval
239              ELSE
240               a5d_prime(k) = 0. _d 0
241               b5d_prime(k) = 0. _d 0
242               c5d_prime(k) = 0. _d 0
243               d5d_prime(k) = 0. _d 0
244               e5d_prime(k) = 0. _d 0
245               y5d_prime(k) = 0. _d 0
246               errCode = 1
247              ENDIF
248    
249          ENDDO
250    
251    C--   Backward sweep (starting from bottom)
252          DO k=Nr,1,-1
253           IF (k.EQ.Nr) THEN
254              y5d_update(k) =   y5d_prime(k)
255           ELSEIF (k.EQ.Nr-1) THEN
256              y5d_update(k) =   y5d_prime(k)
257         &     - y5d_update(k+1)*d5d_prime(k)
258           ELSE
259              y5d_update(k) =   y5d_prime(k)
260         &     - y5d_update(k+1)*d5d_prime(k)
261         &     - y5d_update(k+2)*e5d_prime(k)
262           ENDIF
263          ENDDO
264    
265    C--   Update array
266          DO k=1,Nr
267             y5d(i,j,k,bi,bj)=y5d_update(k)
268          ENDDO
269    
270          ENDDO
271          ENDDO
272    
273    #endif  /* ALLOW_SOLVERS_KLOOPINSIDE */
274    
275        RETURN        RETURN
276        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22