/[MITgcm]/MITgcm_contrib/ecco_utils/ecco_v4_release3_optimization/lsopt/dgscale.F
ViewVC logotype

Contents of /MITgcm_contrib/ecco_utils/ecco_v4_release3_optimization/lsopt/dgscale.F

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


Revision 1.1 - (show annotations) (download)
Wed Jan 3 17:13:47 2018 UTC (7 years, 6 months ago) by ou.wang
Branch: MAIN
CVS Tags: HEAD
Check in the optimization used in ECCO v4r3

1
2 subroutine dgscale( nn, gold, xdiff, diag, rmin )
3
4
5 c ==================================================================
6 c SUBROUTINE dgscale
7 c ==================================================================
8 c
9 c o computes new preconditioner and writes it to OPWARMD
10 c
11 c o started: ??? not reproducible
12 c
13 c o Version: 2.1.0, 02-Mar-2000: Patrick Heimbach, MIT/EAPS
14 c
15 c ==================================================================
16 c SUBROUTINE dgscale
17 c ==================================================================
18
19 implicit none
20
21 #include "blas1.h"
22
23 integer nn
24 double precision gold(nn), xdiff(nn), diag(nn)
25
26 integer i
27 double precision r1, rmin, den
28
29
30 c-----------------------------------------
31 c read diagonal
32 c-----------------------------------------
33 call dostore( nn, diag, .false., 3 )
34
35 r1 = 0.
36 do i = 1, nn
37 r1 = r1 + gold(i)*gold(i)*diag(i)
38 end do
39 r1 = 1.0 / r1
40
41 call SSCAL( nn, r1, diag, 1 )
42
43 c-----------------------------------------
44 c update the diagonal
45 c (gg is used as an auxiliary vector)
46 c-----------------------------------------
47
48 den = 0.0
49
50 do i = 1, nn
51 cph(
52 if (diag(i).LE.rmin) then
53 cph print *, 'pathei-lsopt: in dgscale; diag = 0 for i=', i
54 diag(i) = rmin
55 end if
56 cph)
57 den = den + xdiff(i)*xdiff(i) / diag(i)
58 end do
59
60 do i = 1, nn
61 diag(i) = 1./
62 $ (1./diag(i)+gold(i)**2-(xdiff(i)/diag(i))**2/den)
63 if (diag(i).le.rmin) then
64 diag(i) = rmin
65 endif
66 end do
67
68 c-----------------------------------------
69 c write diagonal
70 c-----------------------------------------
71 call dostore( nn, diag, .true., 3 )
72
73 return
74 end

  ViewVC Help
Powered by ViewVC 1.1.22