/[MITgcm]/MITgcm_contrib/darwin2/pkg/monod/solve_tridiagonal_pivot.F
ViewVC logotype

Contents of /MITgcm_contrib/darwin2/pkg/monod/solve_tridiagonal_pivot.F

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

Revision 1.1 - (show annotations) (download)
Thu Aug 23 21:53:00 2012 UTC (12 years, 10 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt64r_20131210, ctrb_darwin2_ckpt64k_20130723, ctrb_darwin2_ckpt65w_20160512, ctrb_darwin2_ckpt64n_20130826, ctrb_darwin2_ckpt65g_20141120, ctrb_darwin2_ckpt64z_20140711, ctrb_darwin2_ckpt65h_20141217, ctrb_darwin2_ckpt66a_20161020, ctrb_darwin2_ckpt63s_20120908, ctrb_darwin2_ckpt65k_20150402, ctrb_darwin2_ckpt66g_20170424, ctrb_darwin2_ckpt64v_20140411, ctrb_darwin2_ckpt65t_20160221, ctrb_darwin2_ckpt64w_20140502, ctrb_darwin2_ckpt64o_20131024, ctrb_darwin2_ckpt66k_20171025, ctrb_darwin2_ckpt65j_20150225, ctrb_darwin2_ckpt66n_20180118, ctrb_darwin2_ckpt65l_20150504, ctrb_darwin2_ckpt65v_20160409, ctrb_darwin2_ckpt65z_20160929, ctrb_darwin2_ckpt65s_20160114, ctrb_darwin2_ckpt66e_20170314, ctrb_darwin2_ckpt64_20121012, ctrb_darwin2_ckpt64p_20131118, ctrb_darwin2_ckpt64m_20130820, ctrb_darwin2_ckpt64y_20140622, ctrb_darwin2_ckpt66d_20170214, ctrb_darwin2_ckpt65d_20140915, ctrb_darwin2_ckpt64q_20131118, ctrb_darwin2_ckpt65n_20150729, ctrb_darwin2_ckpt65p_20151023, ctrb_darwin2_ckpt64b_20121224, ctrb_darwin2_ckpt66h_20170602, ctrb_darwin2_ckpt64d_20130219, ctrb_darwin2_ckpt65m_20150615, ctrb_darwin2_ckpt65q_20151118, ctrb_darwin2_ckpt64i_20130622, ctrb_darwin2_ckpt64p_20131024, ctrb_darwin2_ckpt64g_20130503, ctrb_darwin2_ckpt65f_20141014, ctrb_darwin2_ckpt65_20140718, ctrb_darwin2_ckpt65u_20160315, ctrb_darwin2_ckpt65r_20151221, ctrb_darwin2_ckpt64f_20130405, ctrb_darwin2_ckpt66i_20170718, ctrb_darwin2_ckpt66b_20161219, ctrb_darwin2_ckpt64u_20140308, ctrb_darwin2_ckpt64c_20130120, ctrb_darwin2_ckpt64s_20140105, ctrb_darwin2_ckpt64a_20121116, ctrb_darwin2_ckpt65i_20150123, ctrb_darwin2_ckpt65e_20140929, ctrb_darwin2_ckpt64j_20130704, ctrb_darwin2_ckpt64x_20140524, ctrb_darwin2_ckpt65o_20150914, ctrb_darwin2_ckpt65c_20140830, ctrb_darwin2_ckpt66j_20170815, ctrb_darwin2_ckpt64e_20130305, ctrb_darwin2_ckpt65y_20160801, ctrb_darwin2_ckpt64t_20140202, ctrb_darwin2_ckpt66c_20170121, ctrb_darwin2_ckpt66l_20171025, ctrb_darwin2_ckpt66o_20180209, ctrb_darwin2_ckpt65a_20140728, ctrb_darwin2_ckpt64l_20130806, ctrb_darwin2_ckpt65x_20160612, ctrb_darwin2_ckpt64h_20130528, ctrb_darwin2_ckpt65b_20140812, ctrb_darwin2_ckpt66f_20170407, ctrb_darwin2_ckpt66m_20171213, HEAD 
add files
1 C $Header$
2 C $Name$
3
4 #include "CPP_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: SOLVE_TRIDIAGONAL_PIVOT
8 C !INTERFACE:
9 SUBROUTINE SOLVE_TRIDIAGONAL_PIVOT(
10 U a3d, b3d, c3d,
11 U y3d,
12 I n, myThid )
13 C !DESCRIPTION: \bv
14 C *==========================================================*
15 C | S/R SOLVE_TRIDIAGONAL_PIVOT
16 C | o Solve a tri-diagonal system A*X=Y (dimension n)
17 C | by Gaussian elimination with pivoting
18 C *==========================================================*
19 C | o input arrays are overwritten
20 C | o result is returned in y3d
21 C *==========================================================*
22 C \ev
23
24 C !USES:
25 IMPLICIT NONE
26 C == Global data ==
27 #include "SIZE.h"
28
29 C !INPUT/OUTPUT PARAMETERS:
30 C == Routine Arguments ==
31 C a3d(2:n) :: matrix lower diagnonal
32 C b3d(1:n) :: matrix main diagnonal
33 C c3d(1:n-1) :: matrix upper diagnonal
34 C y3d(1:n) :: Input = Y vector ; Output = X = solution of A*X=Y
35 _RL a3d(2*Nr)
36 _RL b3d(2*Nr)
37 _RL c3d(2*Nr)
38 _RL y3d(2*Nr)
39 INTEGER n,myThid
40 CEOP
41
42 C !LOCAL VARIABLES:
43 C == Local variables ==
44 INTEGER k
45 _RL tmpc, tmpy
46 _RL recVar
47 C d3d(1:Nr-2) :: second upper diagnonal, generated by pivoting
48 _RL d3d(2*Nr)
49
50 c3d(n) = 0.
51
52 c eliminate lower diagonal
53 c only c3d, d3d and y3d are modified and will be used later
54 DO k=1,n-1
55 IF(ABS(b3d(k)).GE.ABS(a3d(k+1)))THEN
56 c scale current row, subtract from next
57 recVar = 1. / b3d(k)
58 c3d(k) = c3d(k)*recVar
59 y3d(k) = y3d(k)*recVar
60 b3d(k+1) = b3d(k+1) - a3d(k+1)*c3d(k)
61 y3d(k+1) = y3d(k+1) - a3d(k+1)*y3d(k)
62 d3d(k) = 0.
63 ELSE
64 c swap rows, then scale current and subtract from next
65 recVar = 1. / a3d(k+1)
66 tmpc = c3d(k)
67 tmpy = y3d(k)
68 c3d(k) = b3d(k+1)*recVar
69 d3d(k) = c3d(k+1)*recVar
70 y3d(k) = y3d(k+1)*recVar
71 b3d(k+1) = tmpc - b3d(k)*c3d(k)
72 c3d(k+1) = - b3d(k)*d3d(k)
73 y3d(k+1) = tmpy - b3d(k)*y3d(k)
74 ENDIF
75 ENDDO
76 recVar = 1. / b3d(n)
77 y3d(n) = y3d(n)*recVar
78
79 c backsubstitute
80 c y3d(n) is already good
81 y3d(n-1) = y3d(n-1) - c3d(n-1)*y3d(n)
82 DO k=n-2,1,-1
83 y3d(k) = y3d(k) - c3d(k)*y3d(k+1) - d3d(k)*y3d(k+2)
84 ENDDO
85
86 RETURN
87 END
88
  ViewVC Help
Powered by ViewVC 1.1.22