Skip to content

Commit 5ce0d8c

Browse files
MiCurryMiCurry
committed
Prototype - decrease readZ time by 9x - DO NOT MERGE
This is a somewhat simple change that significantly decreases the time taken in readZ. This is 100% a prototype commit, as it's really crudely done, but I wanted to have a prototype to see if this change would be possible and it only works for data that have only FULL_IMPEDANCE In this commit, we create an array of Mtrx pointers, rx_map (of type Mtrx_has_map_t). Because this is just a pointer array, it should be relatively small in memory. We allocate `rx_map` to be the expected number of stations as read from the header. Then, as we read the each data line, we compute a hash using the FNV-1A hashing algorithm. From the computed hash, we modulo it with the total number of stations, giving us a index into `rx_map`. That index will give us constant time access, removing the need for the O(n^2) action. We can thus use the map to identify if a receiver has already been added. For a CONUS case, this reduces the time in readZ from 9 mins, to ~1 min. We also reduce the time, but only in the Full_Impedance read, by not using `backspace` so that we don't re-read the line of the file again. This saves only about 1s of total time.
1 parent 050fd6c commit 5ce0d8c

3 files changed

Lines changed: 117 additions & 51 deletions

File tree

‎f90/3D_MT/DICT/receivers.f90‎

Lines changed: 99 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -46,12 +46,20 @@ module receivers
4646
! data functional module, and can thus be private
4747
type (MTrx), pointer, save, public, dimension(:) :: rxDict
4848

49+
type :: MTrx_hash_map_t
50+
type (MTrx), pointer :: ptr => null()
51+
end type MTrx_hash_map_t
52+
53+
type (MTrx_hash_map_t), dimension(:), pointer :: rx_map => null()
54+
4955
! for efficiency, we now initialize with MAX_NRX and trim after reading the file
5056
integer, parameter :: MAX_NRX = 200000
57+
integer, public :: current_rx_count = 0
5158

5259

5360
Contains
5461

62+
5563
! **************************************************************************
5664
! Initializes and sets up receiver dictionary
5765
! The reciever dictionary contains sparse vectors required
@@ -62,15 +70,19 @@ subroutine setup_rxDict(nSites, siteLocations,siteIDs)
6270
real(kind=prec), intent(in), optional :: siteLocations(nSites,3)
6371
character(*), intent(in), optional :: siteIDs(nSites)
6472

65-
6673
! local variables
6774
integer :: i,istat
6875
character(3) :: id
6976

7077
if (.not. associated(rxDict)) then
7178
allocate(rxDict(nSites),STAT=istat)
79+
allocate(rx_map(nSites),stat=istat)
7280
end if
7381

82+
do i = 1, nSites
83+
rx_map(i) % ptr => null()
84+
end do
85+
7486
if (present(siteLocations)) then
7587
do i = 1,nSites
7688
rxDict(i)%x = siteLocations(i,:)
@@ -83,33 +95,76 @@ subroutine setup_rxDict(nSites, siteLocations,siteIDs)
8395
end do
8496
end if
8597

98+
current_rx_count = 0
99+
86100
end subroutine setup_rxDict
87101

102+
function hash_fnv_1a(key) result(hash)
103+
104+
use iso_fortran_env, only: int32, int64
105+
106+
implicit none
107+
108+
character(len=*), intent(in) :: key
109+
integer(int64), parameter :: FNV_OFFSET = 1469598103934665603_int64
110+
integer(int64), parameter :: FNV_PRIME = 1099511628211_int64
111+
integer(int64) :: hash
112+
integer :: i, n
113+
114+
hash = FNV_OFFSET
115+
n = len_trim(key)
116+
do i = 1, n
117+
hash = ieor(hash, int(iachar(key(i:i)), int64))
118+
hash = hash * FNV_PRIME
119+
end do
120+
121+
end function hash_fnv_1a
122+
123+
124+
subroutine calculate_rx_index(code, total_rxs, idx)
125+
126+
implicit none
127+
128+
character(len=*), intent(in) :: code
129+
integer, intent(in) :: total_rxs
130+
integer, intent(out) :: idx
131+
integer :: hash_key
132+
133+
hash_key = hash_fnv_1a(trim(code))
134+
idx = modulo(hash_key, total_rxs) + 1
135+
136+
end subroutine calculate_rx_index
137+
88138
!**********************************************************************
89139
! Updates the receiver dictionary with a new location and site ID.
90140
! Returns the index of the new element.
91141
! This is not efficient; but this would only be used a few times, with
92142
! a small number of values, so convenience is much more of an issue here!
93143
! NM: modified to include referance site info.
94144

95-
function update_rxDict(loc,id,Rx_azi,loc_ref,id_ref) result (iRx)
145+
function update_rxDict(loc,id,Rx_azi,loc_ref,id_ref,code, total_rxs) result (iRx)
96146

97147
character(*), intent(in) :: id
98148
real(kind=prec), intent(in) :: loc(3)
99149
real(kind=prec),intent(in),optional :: Rx_azi
100150
real(kind=prec),intent(in),optional :: loc_ref(3)
101151
character(*), intent(in),optional :: id_ref
152+
character(*), intent(in),optional :: code
153+
integer, intent(in), optional :: total_rxs
102154
integer :: iRx
103155
! local
104-
type(MTrx) :: new
156+
type(MTrx), target :: new
105157
type(MTrx), pointer, dimension(:) :: temp
106158
integer :: nRx, istat,i
107159
logical :: new_Rx
160+
integer :: hash_idx
108161

109162
! Create a receiver for this location
110163
new%id = id
111164
new%x = loc
112165

166+
call calculate_rx_index(code, total_rxs, hash_idx)
167+
113168
if (present(loc_ref)) then
114169
new%r = loc_ref
115170
new%id_ref=id_ref
@@ -118,60 +173,61 @@ function update_rxDict(loc,id,Rx_azi,loc_ref,id_ref) result (iRx)
118173
if (present(Rx_azi)) then
119174
new%Rx_azi = Rx_azi
120175
end if
121-
122176

123177
! If rxDict doesn't yet exist, create it
124178
if(.not. associated(rxDict)) then
179+
write(0,*) " we are allocating rxDict"
125180
allocate(rxDict(1),STAT=istat)
126181
rxDict(1) = new
127182
iRx = 1
128183
new_Rx = .true.
129184
return
130185
end if
131186

132-
nRx = count_rx()
133-
134-
! If this site isn't new, do nothing, unless a new ref. site is found in case of Inter-Stations TF.
135-
do i = 1,nRx
136-
if (new%id .eq. rxDict(i)%id) then
137-
if (present(loc_ref)) then
138-
!Check if the this site is associated with the same Ref. site. If not, then continue and append another site to the Rx dictionary.
139-
if (new%id_ref .eq. rxDict(i)%id_ref) then
140-
rxDict(i)%r=loc_ref
141-
rxDict(i)%id_ref=id_ref
142-
iRx=i
143-
new_Rx = .false.
144-
return
145-
end if
146-
elseif (present(Rx_azi)) then
147-
! Check if the this site has same azimuth as what we have already in the Dic -
148-
! this is only well-defined for CSEM; for MT, Rx_azi is not used since we're supporting
149-
! the possibility that each field component has its own orientation (stored in DataSpace).
150-
if (new%Rx_azi .eq. rxDict(i)%Rx_azi) then
151-
iRx=i
152-
new_Rx = .false.
153-
return
154-
end if
155-
else
156-
iRx=i
157-
new_Rx = .false.
158-
return
159-
end if
160-
end if
161-
end do
187+
!nRx = count_rx()
188+
189+
if (associated(rx_map(hash_idx) % ptr)) then
190+
if (present(loc_ref)) then
191+
if (new%id_ref .eq. rx_map(hash_idx) % ptr % id_ref) then
192+
rx_map(hash_idx) % ptr %r=loc_ref
193+
rx_map(hash_idx) % ptr %id_ref=id_ref
194+
iRx=i
195+
new_Rx = .false.
196+
return
197+
end if
198+
elseif (present(Rx_azi)) then
199+
! Check if the this site has same azimuth as what we have already in the Dic -
200+
! this is only well-defined for CSEM; for MT, Rx_azi is not used since we're supporting
201+
! the possibility that each field component has its own orientation (stored in DataSpace).
202+
if (new%Rx_azi .eq. rx_map(hash_idx) % ptr %Rx_azi) then
203+
iRx=i
204+
new_Rx = .false.
205+
return
206+
end if
207+
else
208+
iRx=i
209+
new_Rx = .false.
210+
return
211+
end if
212+
end if
162213

163214
! If the site really is new, append to the end of the dictionary
164215
new_Rx = .true.
165216
new%defined = .true.
166-
if (nRx < size(rxDict)) then
167-
iRx = nRx+1
217+
iRx = current_rx_count
218+
if (current_rx_count < size(rxDict)) then
219+
current_rx_count = current_rx_count + 1
220+
iRx = current_rx_count
168221
rxDict(iRx) = new
169222
else
170223
write(0,*) 'ERROR: The number of receivers (sites) in file has exceeded the hard-coded maximum', MAX_NRX
171224
write(0,*) 'ERROR: Please edit MAX_NRX in receivers.f90 and recompile'
172225
call ModEM_abort()
173226
end if
174227

228+
!write(0,*) iRx, trim(code), hash_idx
229+
rx_map(hash_idx) % ptr => rxDict(iRx)
230+
175231
end function update_rxDict
176232

177233
!**********************************************************************
@@ -207,7 +263,7 @@ subroutine write_rxDict_asc(iounit)
207263
return
208264
end if
209265

210-
nRx = count_rx()
266+
!nRx = count_rx()
211267
write(iounit,*) nRx, ' Receivers'
212268
do iRx = 1, nRx
213269
write(iounit,'(i6,2x,a20,4f12.3,a20)') iRx,trim(rxDict(iRx)%id),rxDict(iRx)%x,&
@@ -274,12 +330,13 @@ end function count_rx
274330

275331
! **************************************************************************
276332
! Trims the tail end of the pre-allocated rxDict
277-
subroutine trim_rxDict()
333+
subroutine trim_rxDict(nRx)
278334

335+
integer, intent(in) :: nRx
279336
type(MTrx), pointer, dimension(:) :: temp
280-
integer :: nRx, istat
337+
integer :: istat
281338

282-
nRx = count_rx()
339+
!nRx = count_rx()
283340

284341
if (associated(rxDict)) then
285342
allocate(temp(nRx),STAT=istat)
@@ -290,6 +347,8 @@ subroutine trim_rxDict()
290347
deallocate(temp,STAT=istat)
291348
end if
292349

350+
write(0,*) 'rxDict: ', size(rxDict)
351+
293352
end subroutine trim_rxDict
294353

295354

‎f90/3D_MT/DataIO.f90‎

Lines changed: 16 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -825,43 +825,44 @@ subroutine read_Z_list(allData,cfile)
825825
end if
826826

827827
! Liu Zhongyin, 2019.08.27, add new codes for reading data
828-
backspace(ioDat)
828+
!backspace(ioDat)
829829
call strcount(tmpline, ' ', ncount)
830830
select case (ncount)
831831
case(11)
832-
read(ioDat,*,iostat=ios) Period,code,lat,lon,x(1),x(2),x(3),compid,Zreal,Zimag,Zerr
832+
!read(ioDat,*,iostat=ios) Period,code,lat,lon,x(1),x(2),x(3),compid,Zreal,Zimag,Zerr
833+
read(tmpline,*) Period,code,lat,lon,x(1),x(2),x(3),compid,Zreal,Zimag,Zerr
833834
Hxangle = fileInfo(iTxt,iDt)%geographic_orientation
834835
Exangle = fileInfo(iTxt,iDt)%geographic_orientation
835836
Hxangle_ref = 0.0
836837
Hyangle = Hxangle + 90.0
837838
Eyangle = Exangle + 90.0
838839
Hyangle_ref = Hxangle_ref + 90.0
839840
case(12)
840-
read(ioDat,*,iostat=ios) Period,code,lat,lon,x(1),x(2),x(3),compid,Zreal,Zimag,Zerr,Hxangle
841+
read(tmpline,*) Period,code,lat,lon,x(1),x(2),x(3),compid,Zreal,Zimag,Zerr,Hxangle
841842
Hxangle = Hxangle + fileInfo(iTxt,iDt)%geographic_orientation
842843
Exangle = Hxangle
843844
Hxangle_ref = 0.0
844845
Hyangle = Hxangle + 90.0
845846
Eyangle = Exangle + 90.0
846847
Hyangle_ref = Hxangle_ref + 90.0
847848
case(13)
848-
read(ioDat,*,iostat=ios) Period,code,lat,lon,x(1),x(2),x(3),compid,Zreal,Zimag,Zerr,Hxangle,Hyangle
849+
read(tmpline,*) Period,code,lat,lon,x(1),x(2),x(3),compid,Zreal,Zimag,Zerr,Hxangle,Hyangle
849850
Hxangle = Hxangle + fileInfo(iTxt,iDt)%geographic_orientation
850851
Exangle = Hxangle
851852
Hxangle_ref = 0.0
852853
Hyangle = Hyangle + fileInfo(iTxt,iDt)%geographic_orientation
853854
Eyangle = Exangle + 90.0
854855
Hyangle_ref = Hxangle_ref + 90.0
855856
case(14)
856-
read(ioDat,*,iostat=ios) Period,code,lat,lon,x(1),x(2),x(3),compid,Zreal,Zimag,Zerr,Hxangle,Hyangle,Exangle
857+
read(tmpline,*) Period,code,lat,lon,x(1),x(2),x(3),compid,Zreal,Zimag,Zerr,Hxangle,Hyangle,Exangle
857858
Hxangle = Hxangle + fileInfo(iTxt,iDt)%geographic_orientation
858859
Exangle = Exangle + fileInfo(iTxt,iDt)%geographic_orientation
859860
Hxangle_ref = 0.0
860861
Hyangle = Hyangle + fileInfo(iTxt,iDt)%geographic_orientation
861862
Eyangle = Exangle + 90.0
862863
Hyangle_ref = Hxangle_ref + 90.0
863864
case(15)
864-
read(ioDat,*,iostat=ios) Period,code,lat,lon,x(1),x(2),x(3),compid,Zreal,Zimag,Zerr,Hxangle,Hyangle,Exangle,Eyangle
865+
read(tmpline,*) Period,code,lat,lon,x(1),x(2),x(3),compid,Zreal,Zimag,Zerr,Hxangle,Hyangle,Exangle,Eyangle
865866
Hxangle = Hxangle + fileInfo(iTxt,iDt)%geographic_orientation
866867
Exangle = Exangle + fileInfo(iTxt,iDt)%geographic_orientation
867868
Hxangle_ref = 0.0
@@ -895,7 +896,7 @@ subroutine read_Z_list(allData,cfile)
895896
x(1) = lat
896897
x(2) = lon
897898
end if
898-
iRx = update_rxDict(x,siteid)
899+
iRx = update_rxDict(x, siteid, code=code, total_rxs=nRx)
899900

900901
case(Full_Interstation_TF)
901902
read(ioDat,'(a)',iostat=ios) tmpline
@@ -1199,19 +1200,24 @@ subroutine read_Z_list(allData,cfile)
11991200
write(0,*)
12001201

12011202
! Now reallocate rxDict to only include sites that have benn found in file
1202-
call trim_rxDict()
12031203

12041204
! Finally, set up the index vectors in the data type dictionary - used for output
12051205
nTxt = 5
12061206
nTx = size(txDict)
1207-
nRx = size(rxDict)
1207+
!nTx =
1208+
!nRx = size(rxDict)
1209+
nRx = current_rx_count
1210+
1211+
call trim_rxDict(nRx)
1212+
1213+
write(0,*) 'nTx: ', nTx, ' nRx:', nRx
12081214
do iTxt = 1,nTxt
12091215
do iDt = 1,nDt
12101216
allocate(fileInfo(iTxt,iDt)%tx_index(nTx),STAT=istat)
12111217
allocate(fileInfo(iTxt,iDt)%dt_index(nTx),STAT=istat)
12121218
allocate(fileInfo(iTxt,iDt)%rx_index(nTx,nRx),STAT=istat)
12131219
call index_dataVectorMTX(allData,iTxt,iDt,fileInfo(iTxt,iDt)%tx_index,fileInfo(iTxt,iDt)%dt_index,fileInfo(iTxt,iDt)%rx_index)
1214-
end do
1220+
end do
12151221
end do
12161222

12171223
end subroutine read_Z_list

‎f90/SENS/DataSpace.f90‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1739,7 +1739,8 @@ subroutine index_dataVectorMTX(d,iTxt,iDt,tx_index,dt_index,rx_index)
17391739
if ((d%d(j)%data(i)%txType == iTxt) .and. (d%d(j)%data(i)%dataType == iDt)) then
17401740
tx_index(d%d(j)%tx) = j
17411741
dt_index(d%d(j)%tx) = i
1742-
do k = 1,d%d(j)%data(i)%nSite
1742+
write(0,*) d % d(j) % data(i) % nSite - 1
1743+
do k = 1,d%d(j)%data(i)%nSite - 1
17431744
rx_index(d%d(j)%data(i)%tx,d%d(j)%data(i)%rx(k)) = k
17441745
enddo
17451746
endif

0 commit comments

Comments
 (0)