-
Notifications
You must be signed in to change notification settings - Fork 2
/
xux.f
314 lines (296 loc) · 9.18 KB
/
xux.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
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
C->>> -----------------------------------------------> ems_ze_pv_c_v <<<
c Zeroes the nonzero entries in pv_c_v.
c
subroutine ems_ze_pv_c_v(pv_c_v, nw_eta_ix)
implicit none
include 'EMSV.INC'
include 'EMSPM.INC'
include 'ICTVR.INC'
include 'RSMICOM.INC'
integer nw_eta_ix(0:nw_eta_l_ix)
double precision pv_c_v(0:n_r)
integer ix_n, r_n
if (nw_eta_l_ix .le. n_r) then
do 10, ix_n = nw_eta_f_ix, nw_eta_l_ix
pv_c_v(nw_eta_ix(ix_n)) = zero
10 continue
else
do 20, r_n = 1, n_r
pv_c_v(r_n) = zero
20 continue
end if
return
end
C->>> -------------------------------------------------> ems_ze_pi_v <<<
c Zeroes the nonzero entries in pi.
c
subroutine ems_ze_pi_v(pi_v, pi_ix)
implicit none
include 'EMSV.INC'
include 'EMSPM.INC'
include 'ICTVR.INC'
integer pi_ix(0:n_r)
double precision pi_v(0:n_r)
integer ix_n, r_n
if (pi_ix(0) .le. n_r) then
do 10, ix_n = 1, pi_ix(0)
pi_v(pi_ix(ix_n)) = zero
10 continue
else
do 20, r_n = 1, n_r
pi_v(r_n) = zero
20 continue
end if
return
end
C->>> --------------------------------------------> ems_g_vec_ix_bar <<<
c Sets up pointers into vec_ix of the nonzeros which it indicates.
c NB Assumes that vec_ix_bar is zeroed on entry.
c
subroutine ems_g_vec_ix_bar(lc_n_r, vec_ix, vec_ix_bar)
implicit none
include 'EMSV.INC'
include 'EMSPM.INC'
include 'RSMICS.INC'
include 'ICTVR.INC'
include 'EMSMSG.INC'
integer lc_n_r, vec_ix(0:lc_n_r), vec_ix_bar(0:lc_n_r)
integer vec_ix_n
if (vec_ix_bar(0) .ne. 0) then
if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9900)
call ems_msg_wr_li(bug_msg_n)
CM IF (emsol_deb .EQ. 1) THEN
C? call ems_dump
CM ENDIF
end if
if (iand(ck_msk, ze_a_ck_bt) .ne. 0)
& call ems_ck_ze_i_a(lc_n_r, vec_ix_bar)
do 10, vec_ix_n = 1, vec_ix(0)
vec_ix_bar(vec_ix(vec_ix_n)) = vec_ix_n
10 continue
vec_ix_bar(0) = 1
return
9900 format('vec_ix_bar has not been zeroed')
end
C->>> -------------------------------------------> ems_ze_vec_ix_bar <<<
c Zeroes the nonzero entries in vec_ix_bar.
c
subroutine ems_ze_vec_ix_bar(lc_n_r, vec_ix, vec_ix_bar)
implicit none
include 'EMSV.INC'
include 'EMSPM.INC'
integer lc_n_r, vec_ix(0:lc_n_r), vec_ix_bar(0:lc_n_r)
integer ix_n, r_n
if (vec_ix(0) .lt. tl_dse_vec*lc_n_r) then
do 10, ix_n = 1, vec_ix(0)
vec_ix_bar(vec_ix(ix_n)) = 0
10 continue
else
do 20, r_n = 1, lc_n_r
vec_ix_bar(r_n) = 0
20 continue
end if
vec_ix_bar(0) = 0
return
end
C->>> -------------------------------------------> ems_ck_vec_ix_bar <<<
c Checks the consistency of vec_ix_bar with vec_v and vec_ix.
c
subroutine ems_ck_vec_ix_bar(lc_n_r, vec_v, vec_ix, vec_ix_bar)
implicit none
include 'EMSV.INC'
include 'EMSMSG.INC'
integer lc_n_r, vec_ix(0:lc_n_r), vec_ix_bar(0:lc_n_r)
double precision vec_v(0:lc_n_r)
integer r_n, ix_n, n_ix
if (vec_ix_bar(0) .eq. 0) go to 7000
do 10, ix_n = 1, vec_ix(0)
r_n = vec_ix(ix_n)
if (vec_ix_bar(r_n) .ne. ix_n) go to 8001
if (vec_v(r_n) .eq. zero) go to 8002
10 continue
n_ix = 0
do 20, r_n = 1, lc_n_r
if (vec_v(r_n) .ne. zero) n_ix = n_ix + 1
20 continue
if (n_ix .ne. vec_ix(0)) go to 8003
7000 continue
return
8001 continue
if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9810)
call ems_msg_wr_li(bug_msg_n)
CM IF (emsol_deb .EQ. 1) THEN
C? call ems_dump
CM ENDIF
go to 7000
8002 continue
if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9820)
call ems_msg_wr_li(bug_msg_n)
CM IF (emsol_deb .EQ. 1) THEN
C? call ems_dump
CM ENDIF
go to 7000
8003 continue
if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9830)
call ems_msg_wr_li(bug_msg_n)
CM IF (emsol_deb .EQ. 1) THEN
C? call ems_dump
CM ENDIF
go to 7000
9810 format('Inconsistency between vec_ix and vec_ix_bar ')
9820 format('Inconsistency between vec_ix and vec_v ')
9830 format('Inconsistency between n_ix and vec_ix(0)')
end
C->>> ----------------------------------------------> ems_cp_nz_v_ix <<<
c Makes a copy of the nonzero values and indices for an array.
c
subroutine ems_cp_nz_v_ix(lc_n_r, fm_ix, fm_v, t_ix, t_v)
implicit none
include 'EMSV.INC'
include 'EMSPM.INC'
include 'RSMICS.INC'
include 'ICTVR.INC'
integer lc_n_r, fm_ix(0:lc_n_r), t_ix(0:lc_n_r)
double precision fm_v(0:lc_n_r), t_v(0:lc_n_r)
integer ix_n
if (iand(ck_msk, ze_a_ck_bt) .ne. 0)
& call ems_ck_ze_rl_a(lc_n_r, t_v)
t_ix(0) = fm_ix(0)
do 10, ix_n = 1, fm_ix(0)
t_ix(ix_n) = fm_ix(ix_n)
t_v(fm_ix(ix_n)) = fm_v(fm_ix(ix_n))
10 continue
return
end
C->>> ----------------------------------------------> ems_ck_ze_rl_a <<<
c Checks that the array supplied is zeroed.
c
subroutine ems_ck_ze_rl_a(lc_n_r, rl_a)
implicit none
include 'EMSV.INC'
include 'EMSMSG.INC'
integer lc_n_r
double precision rl_a(0:lc_n_r)
integer r_n
logical er_fd
integer pr_pass_1
save pr_pass_1
data pr_pass_1/0/
er_fd = .false.
do 10, r_n = 1, lc_n_r
c
c Keep this .ne. zero just for safety.
c
if (rl_a(r_n) .ne. zero) then
er_fd = .true.
if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9000)
& r_n, rl_a(r_n)
call ems_msg_wr_li(bug_msg_n)
c rl_a(r_n) = zero
end if
10 continue
if (er_fd) then
CM IF (emsol_deb .EQ. 1) THEN
C? call ems_dump
CM ENDIF
else if (pr_pass_1 .eq. 0) then
if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9100)lc_n_r
call ems_msg_wr_li(info_msg_n)
pr_pass_1 = 1
end if
return
9000 format('Entry ', i9, ' in zeroed array has value ', g11.4)
9100 format('Real array of dimension ', i9, ' is zeroed ')
end
C->>> -----------------------------------------------> ems_ck_ze_i_a <<<
c Checks that the integer array supplied is zeroed.
c
subroutine ems_ck_ze_i_a(lc_n_r, i_a)
implicit none
include 'EMSV.INC'
include 'EMSMSG.INC'
integer lc_n_r, i_a(0:lc_n_r)
integer r_n
logical er_fd
integer pr_pass_1
save pr_pass_1
data pr_pass_1/0/
er_fd = .false.
do 10, r_n = 1, lc_n_r
if (i_a(r_n) .ne. 0) then
er_fd = .true.
if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9000)
& r_n, i_a(r_n)
call ems_msg_wr_li(bug_msg_n)
i_a(r_n) = 0
end if
10 continue
if (er_fd) then
CM IF (emsol_deb .EQ. 1) THEN
C? call ems_dump
CM ENDIF
else if (pr_pass_1 .eq. 0) then
if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9100)lc_n_r
call ems_msg_wr_li(info_msg_n)
pr_pass_1 = 1
end if
return
9000 format('Entry ', i9, ' in zeroed array has value ', i9)
9100 format('Integer array of dimension ', i9, ' is zeroed ')
end
C->>> ------------------------------------------> ems_g_mtx_vec_prod <<<
subroutine ems_g_mtx_vec_prod(ty, aa, vec, beta, vec_prod, ds, is)
implicit none
include 'EMSV.INC'
include 'EMSPM.INC'
include 'EMSMMGR.INC'
include 'EMSMEM.INC'
include 'EMSP.INC'
include 'ICTVR.INC'
include 'RLCTVR.INC'
integer ty, is(0:is_n_en_m1)
double precision aa, vec(*), beta, vec_prod(*), ds(0:ds_n_en_m1)
integer c_n, r_n, el_n
if (ty .eq. 1) then
if (beta .eq. zero) then
do 5 r_n=1, n_r
vec_prod(r_n) = zero
5 continue
else if (beta .ne. one) then
do 10 r_n=1, n_r
vec_prod(r_n) = vec_prod(r_n)*beta
10 continue
end if
if (aa .ne. zero) then
do 15 c_n=1, n_c
if (vec(c_n) .eq. zero) go to 15
do 14, el_n = is(p_mtx_c_sa+c_n),
& is(p_mtx_c_sa+c_n+1)-1
vec_prod(is(p_mtx_r_ix+el_n)) =
& vec_prod(is(p_mtx_r_ix+el_n)) +
& aa*vec(c_n)*ds(p_mtx_r_v+el_n)
14 continue
15 continue
end if
else
if (beta .eq. zero) then
do 18 c_n=1, n_c
vec_prod(c_n) = zero
18 continue
else if (beta .ne. one) then
do 20 c_n=1, n_c
vec_prod(c_n) = vec_prod(c_n)*beta
20 continue
end if
if (aa .ne. zero) then
do 25 c_n=1, n_c
do 24, el_n = is(p_mtx_c_sa+c_n),
& is(p_mtx_c_sa+c_n+1)-1
vec_prod(c_n) = vec_prod(c_n) +
& aa*vec(is(p_mtx_r_ix+el_n))*ds(p_mtx_r_v+el_n)
24 continue
25 continue
end if
end if
return
end