-
Notifications
You must be signed in to change notification settings - Fork 22
/
siterativeEof.f
196 lines (166 loc) · 5.72 KB
/
siterativeEof.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
program siterativeEof
use ufileformat
c---------------------------------------------------------------------c
c "iterativeEOF", Version 2.0, December, 2002. c
c A Krylov-based code for computing EOF c
c Vincent TOUMAZOU c
c CNES/MERCATOR-Ocean c
c 18, Avenue Edouard Belin c
c F-31401 TOULOUSE Cedex 4 c
c URL: http:https://www.mercator.com.fr c
c URL: http:https://www.cnes.fr c
c---------------------------------------------------------------------c
c$$$ interface
c$$$ subroutine readMatrix(fname,matrix,imax,jmax)
c$$$ character(len=*),intent(in) :: fname
c$$$ real,intent(out) :: matrix(:,:)
c$$$ integer,intent(out) :: imax,jmax
c$$$ end interface
integer maxm, maxn, maxnev, maxncv, ldv, ldu
parameter (maxm = 40000, maxn=250, maxnev=20, maxncv=25,
& ldu = maxm, ldv=maxn, nlines=22)
Real
& meanA(maxm),
& w(ldu), v(ldv,maxncv), u(ldu, maxnev),
& workl(maxncv*(maxncv+8)), workd(3*maxn),
& s(maxncv,2), resid(maxn),
& tol, t, sumSV
real, pointer :: A(:,:)
logical select(maxncv)
integer M, N, nev, ncv
character*72 comment
character(len=50) list(22)
c
c
write(*,*) '********************************************'
write(*,*) '********** starting program **************'
write(*,*) '********************************************'
c
c Numerical parameters
c
open(1,file='Input/numerical.data')
read(1,500) comment
read(1,500) comment
read(1,500) comment
read(1,500) comment
read(1,500) comment
read(1,*) nev
read(1,500) comment
read(1,500) comment
read(1,*) ncv
read(1,500) comment
read(1,500) comment
read(1,*) tol
close(1)
write(*,*) '***'
write(*,*) 'Numerical data read'
write(*,*) 'Number of required modes = ', nev
write(*,*) 'Size of projection = ', ncv
write(*,*) 'Treshold for Lanczos convergence = ', tol
write(*,*) '***'
c
c Reads the matrix of data
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C CAUTION : Call your subroutine for reading C
C the matrix of data. C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C call readMatrix
C & (Add your parameters here, A, m and n
C & must be the outputs of your subroutine
C & and maxm, maxn the leading dimension of A
C & have to be in the inputs of the subroutine)
C Example : call readMatrix(A, maxm, maxn, m, n)
open(77,file = '/u/alvera/Datafiles/Analist/analistTEM.dat')
do i=1,nlines
read(77,'(a50)') list(i)
end do
call readMatrix(list,A, M, N)
write(*,*) '***'
write(*,*) 'Matrix loaded'
write(*,*) 'Size of the matrix = ', m,' by ', n
write(*,*) '***'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C CAUTION : Uncomment the 2 next calls if you C
C want to center your data with the C
C mean on each point with respect C
C to the total period. C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c
c Computes the Row-means of the matrix
c
call smeanByRow
& (A, m, n, meanA, maxm, maxn)
c
c Transforms the matrix such that its raws have a zero mean
c
call smeanToZero
& (A, m, n, meanA, maxm, maxn)
c
c Estimates the total variance
c
call svariExp
& (A, m, n, maxm, maxn, sumSV, w)
c
c Computes the NEV singular values and the modes of interest
c
write(*,*) '********************************************'
write(*,*) '***** Starting the Lanczos Process *******'
write(*,*) '********************************************'
call ssvd_lancz
& (A, m, n, maxm, maxn, nev, ncv, maxnev, maxncv, ldv,
& ldu, w, v, u, workl, workd, s, resid, select, tol)
write(*,*) '********************************************'
write(*,*) '******* End of the Lanczos Process *******'
write(*,*) '********************************************'
c
c Uncomment the next call if you want to normalize your outputs
c
c call snormal_output
c & (s, u, v, m, n, nev, maxncv, ldu, ldv,
c & maxnev, maxm)
c
c Extracts singular values and associated modes, prepares the
c output file for graphical exploitation.
c
c
c Singular values
c _______________
open(1,file='Output/outputEof.vlsng')
do i=nev,1,-1
write(1,501) s(i,1)
enddo
close(1)
c
c Variance detailed
c _______________
open(1,file='Output/outputEof.varEx')
do i=nev,1,-1
write(1,504) 'Mode ', nev+1-i, ' = ',
& 100.0*s(i,1)*s(i,1)/sumSV, ' %'
enddo
close(1)
c
c Spatial modes
c _______________
open(1,file='Output/outputEof.lftvec')
do j=1,m
write(1,503) (u(j,i),i=nev,1,-1)
enddo
close(1)
c
c Temporal modes
c _______________
open(1,file='Output/outputEof.rghvec')
do j=1,n
write(1,503) (v(j,i),i=nev,1,-1)
enddo
close(1)
close(77)
500 format(a72)
501 format(e10.4)
502 format(3e10.4)
503 format(300(e10.4,1X))
504 format(a5,i3,a3,f8.4,a2)
end