-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathline_vector_mod.F90
239 lines (192 loc) · 7.33 KB
/
line_vector_mod.F90
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
!-------------------------------------------------------------------------------
! Copyright (c) 2017, Met Office, on behalf of HMSO and Queen's Printer
! For further details please refer to the file LICENCE.original which you
! should have received as part of this distribution.
!-------------------------------------------------------------------------------
module line_vector_mod
use constants_mod, only : r_def, i_def
use vector_mod, only : abstract_vector_type
use log_mod, only : log_event, log_scratch_space, LOG_LEVEL_ERROR, LOG_LEVEL_INFO
private
type, public, extends(abstract_vector_type) :: line_vector_type
integer(kind=i_def) :: ndata
real(kind=r_def), public, allocatable, dimension(:) :: vdata
contains
procedure, public :: set_scalar => set_scalar_lv
procedure, public :: axpy => axpy_lv
procedure, public :: aypx => aypx_lv
procedure, public :: norm => norm_lv
procedure, public :: dot => dot_lv
procedure, public :: scale => scale_lv
procedure, public :: duplicate => duplicate_lv
procedure, public :: copy => copy_lv
procedure, public :: p => p_lv
! private procedures
procedure, private :: set_scalar_lv
procedure, private :: axpy_lv
procedure, private :: aypx_lv
procedure, private :: norm_lv
procedure, private :: dot_lv
procedure, private :: scale_lv
procedure, private :: duplicate_lv
procedure, private :: copy_lv
procedure, private :: p_lv
final :: line_vector_destructor
end type line_vector_type
interface line_vector_type
module procedure line_vector_constructor, line_vector_two_val_constructor
end interface line_vector_type
contains
function line_vector_constructor(ndata) result(self)
implicit none
integer, intent(in) :: ndata
type(line_vector_type) :: self
allocate(self%vdata(ndata))
self%ndata = ndata
end function line_vector_constructor
function line_vector_two_val_constructor(ndata, d1, d2) result(self)
implicit none
integer, intent(in) :: ndata
real(kind=r_def), intent(in) :: d1
real(kind=r_def), intent(in) :: d2
type(line_vector_type) :: self
integer(kind=i_def) :: nlp1, nlp2
allocate(self%vdata(ndata))
self%ndata = ndata
do nlp1 = 1, self%ndata/2
self%vdata(nlp1) = d1
end do
do nlp2 = nlp1, self%ndata
self%vdata(nlp2) = d2
end do
end function line_vector_two_val_constructor
subroutine set_scalar_lv(self, scalar)
implicit none
class(line_vector_type), intent(inout) :: self
real(kind=r_def), intent(in) :: scalar
integer :: nlp
do nlp = 1, self%ndata
self%vdata(nlp) = scalar
end do
end subroutine set_scalar_lv
subroutine axpy_lv(self, alpha, x)
implicit none
class(line_vector_type), intent(inout) :: self
real(kind=r_def), intent(in) :: alpha
class(abstract_vector_type), intent(inout) :: x
integer(kind=i_def) :: nlp
select type(x)
type is(line_vector_type)
! compute y = alpha *x + y
do nlp = 1, self%ndata
self%vdata(nlp) = self%vdata(nlp) + alpha*x%vdata(nlp)
end do
class default
write(log_scratch_space,'(A)') "line_vector:axpy: type of x is not line_vector_type"
call log_event(log_scratch_space, LOG_LEVEL_ERROR)
end select
end subroutine axpy_lv
subroutine aypx_lv(self, alpha, x)
implicit none
class(line_vector_type), intent(inout) :: self
real(kind=r_def), intent(in) :: alpha
class(abstract_vector_type), intent(inout) :: x
integer(kind=i_def) :: nlp
select type(x)
type is(line_vector_type)
! compute y = alpha *y + x
do nlp = 1, self%ndata
self%vdata(nlp) = alpha*self%vdata(nlp) + x%vdata(nlp)
end do
class default
write(log_scratch_space,'(A)') "line_vector:aypx: type of x is not line_vector_type"
call log_event(log_scratch_space, LOG_LEVEL_ERROR)
end select
end subroutine aypx_lv
function norm_lv(self) result(normal)
implicit none
class(line_vector_type), intent(in) :: self
real(kind=r_def) :: normal
integer(kind=i_def) :: nlp
normal = 0.0_r_def
do nlp = 1, self%ndata
normal = normal + ( self%vdata(nlp) * self%vdata(nlp) )
end do
normal = sqrt(normal)
end function norm_lv
function dot_lv(self, x) result(dot_prod)
implicit none
class(line_vector_type), intent(in) :: self
class(abstract_vector_type), intent(in) :: x
real(kind=r_def) :: dot_prod
integer(kind=i_def) :: nlp
select type(x)
type is(line_vector_type)
dot_prod = 0.0_r_def
do nlp = 1, self%ndata
dot_prod = dot_prod + ( self%vdata(nlp) * x%vdata(nlp) )
end do
class default
write(log_scratch_space,'(A)') "line_vector:dot: type of x is not line_vector_type"
call log_event(log_scratch_space, LOG_LEVEL_ERROR)
end select
end function dot_lv
subroutine scale_lv(self, scalar)
implicit none
class(line_vector_type), intent(inout) :: self
real(kind=r_def), intent(in) :: scalar
self%vdata(:) = self%vdata(:) * scalar
end subroutine scale_lv
subroutine duplicate_lv(self, vec)
implicit none
class(line_vector_type), intent(in) :: self
class(abstract_vector_type), allocatable, intent(inout) :: vec
allocate(line_vector_type::vec)
select type(vec)
type is (line_vector_type)
if(.not.(allocated(vec%vdata))) then
vec%ndata = self%ndata
allocate(vec%vdata(vec%ndata))
end if
class default
write(log_scratch_space,'(A)') "line_vector:duplicate: type of vec is not line_vector_type"
call log_event(log_scratch_space, LOG_LEVEL_ERROR)
end select
end subroutine duplicate_lv
subroutine copy_lv(self, source)
implicit none
class(line_vector_type), intent(out) :: self
class(abstract_vector_type), intent(in) :: source
select type(source)
type is(line_vector_type)
if(.not.(allocated(self%vdata)))then
self%ndata = source%ndata
allocate(self%vdata(self%ndata))
end if
self%vdata(:) = source%vdata(:)
class default
write(log_scratch_space,'(A)') "line_vector:copy: type of source is not line_vector_type"
call log_event(log_scratch_space, LOG_LEVEL_ERROR)
end select
end subroutine copy_lv
subroutine p_lv(self)
implicit none
class(line_vector_type), intent(inout) :: self
select type(self)
type is(line_vector_type)
write(log_scratch_space,'(A,E16.8)') "line_vector:p:1st element:",self%vdata(1)
! call log_event(log_scratch_space, LOG_LEVEL_INFO)
write(*,*) trim(log_scratch_space)
class default
write(log_scratch_space,'(A)') "line_vector:p: type of vector is not line_vector_type"
call log_event(log_scratch_space, LOG_LEVEL_ERROR)
end select
end subroutine p_lv
subroutine line_vector_destructor(self)
implicit none
type(line_vector_type), intent(inout) :: self
if(allocated(self%vdata)) then
deallocate(self%vdata)
end if
end subroutine line_vector_destructor
end module line_vector_mod