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

Annotation 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 - (hide 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 jahn 1.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