Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
On Thursday, 7th July from 1 to 3 pm there will be a maintenance with a short downtime of GitLab.
Open sidebar
elpa
elpa
Commits
f0bdb04a
Commit
f0bdb04a
authored
Feb 18, 2021
by
Andreas Marek
Browse files
Remove obsolete code in merge_systems
parent
8db88f22
Changes
1
Pipelines
1
Hide whitespace changes
Inline
Side-by-side
src/solve_tridi/merge_systems_template.F90
View file @
f0bdb04a
...
...
@@ -548,9 +548,11 @@
call
obj
%
timer
%
start
(
"OpenMP parallel"
//
PRECISION_SUFFIX
)
!$OMP PARALLEL DO PRIVATE(i) SHARED(na1, my_proc, n_procs, &
!$OMP d1,dbase, ddiff, z, ev_scale, obj) &
!$OMP DEFAULT(NONE)
!$omp PARALLEL DO &
!$omp default(none) &
!$omp private(i) &
!$omp SHARED(na1, my_proc, n_procs, &
!$OMP d1,dbase, ddiff, z, ev_scale, obj)
#endif
DO
i
=
my_proc
+1
,
na1
,
n_procs
! work distributed over all processors
...
...
@@ -915,345 +917,5 @@
return
#if 0
contains
subroutine
add_tmp_
&
&
PRECISION
&
&(
obj
,
d1
,
dbase
,
ddiff
,
z
,
ev_scale_value
,
na1
,
i
)
use
precision
use
elpa_abstract_impl
implicit
none
class
(
elpa_abstract_impl_t
),
intent
(
inout
)
::
obj
integer
(
kind
=
ik
),
intent
(
in
)
::
na1
,
i
real
(
kind
=
REAL_DATATYPE
),
intent
(
in
)
::
d1
(:),
dbase
(:),
ddiff
(:),
z
(:)
real
(
kind
=
REAL_DATATYPE
),
intent
(
inout
)
::
ev_scale_value
real
(
kind
=
REAL_DATATYPE
)
::
tmp
(
1
:
na1
)
! tmp(1:na1) = z(1:na1) / delta(1:na1,i) ! original code
! tmp(1:na1) = z(1:na1) / (d1(1:na1)-d(i))! bad results
! All we want to calculate is tmp = (d1(1:na1)-dbase(i))+ddiff(i)
! in exactly this order, but we want to prevent compiler optimization
tmp
(
1
:
na1
)
=
d1
(
1
:
na1
)
-
dbase
(
i
)
call
v_add_s_
&
&
PRECISION
&
&(
obj
,
tmp
(
1
:
na1
),
na1
,
ddiff
(
i
))
tmp
(
1
:
na1
)
=
z
(
1
:
na1
)
/
tmp
(
1
:
na1
)
ev_scale_value
=
1.0_rk
/
sqrt
(
dot_product
(
tmp
(
1
:
na1
),
tmp
(
1
:
na1
)))
end
subroutine
add_tmp_
&
&
PRECISION
#endif
#if 0
subroutine
resort_ev_
&
&
PRECISION
&
&(
obj
,
idx_ev
,
nLength
)
use
precision
use
elpa_abstract_impl
implicit
none
class
(
elpa_abstract_impl_t
),
intent
(
inout
)
::
obj
integer
(
kind
=
ik
),
intent
(
in
)
::
nLength
integer
(
kind
=
ik
)
::
idx_ev
(
nLength
)
integer
(
kind
=
ik
)
::
i
,
nc
,
pc1
,
pc2
,
lc1
,
lc2
,
l_cols_out
real
(
kind
=
REAL_DATATYPE
),
allocatable
::
qtmp
(:,:)
integer
(
kind
=
ik
)
::
istat
character
(
200
)
::
errorMessage
if
(
l_rows
==
0
)
return
! My processor column has no work to do
! Resorts eigenvectors so that q_new(:,i) = q_old(:,idx_ev(i))
l_cols_out
=
COUNT
(
p_col_out
(
1
:
na
)
==
my_pcol
)
allocate
(
qtmp
(
l_rows
,
l_cols_out
),
stat
=
istat
,
errmsg
=
errorMessage
)
check_allocate
(
"resort_ev: qtmp"
,
istat
,
errorMessage
)
nc
=
0
do
i
=
1
,
na
pc1
=
p_col
(
idx_ev
(
i
))
lc1
=
l_col
(
idx_ev
(
i
))
pc2
=
p_col_out
(
i
)
if
(
pc2
<
0
)
cycle
! This column is not needed in output
if
(
pc2
==
my_pcol
)
nc
=
nc
+1
! Counter for output columns
if
(
pc1
==
my_pcol
)
then
if
(
pc2
==
my_pcol
)
then
! send and recieve column are local
qtmp
(
1
:
l_rows
,
nc
)
=
q
(
l_rqs
:
l_rqe
,
lc1
)
else
#ifdef WITH_MPI
call
obj
%
timer
%
start
(
"mpi_communication"
)
call
mpi_send
(
q
(
l_rqs
,
lc1
),
int
(
l_rows
,
kind
=
MPI_KIND
),
MPI_REAL_PRECISION
,
pc2
,
int
(
mod
(
i
,
4096
),
kind
=
MPI_KIND
),
&
int
(
mpi_comm_cols
,
kind
=
MPI_KIND
),
mpierr
)
call
obj
%
timer
%
stop
(
"mpi_communication"
)
#endif /* WITH_MPI */
endif
else
if
(
pc2
==
my_pcol
)
then
#ifdef WITH_MPI
call
obj
%
timer
%
start
(
"mpi_communication"
)
call
mpi_recv
(
qtmp
(
1
,
nc
),
int
(
l_rows
,
kind
=
MPI_KIND
),
MPI_REAL_PRECISION
,
pc1
,
int
(
mod
(
i
,
4096
),
kind
=
MPI_KIND
),
&
int
(
mpi_comm_cols
,
kind
=
MPI_KIND
),
MPI_STATUS_IGNORE
,
mpierr
)
call
obj
%
timer
%
stop
(
"mpi_communication"
)
#else /* WITH_MPI */
qtmp
(
1
:
l_rows
,
nc
)
=
q
(
l_rqs
:
l_rqe
,
lc1
)
#endif /* WITH_MPI */
endif
enddo
! Insert qtmp into (output) q
nc
=
0
do
i
=
1
,
na
pc2
=
p_col_out
(
i
)
lc2
=
l_col_out
(
i
)
if
(
pc2
==
my_pcol
)
then
nc
=
nc
+1
q
(
l_rqs
:
l_rqe
,
lc2
)
=
qtmp
(
1
:
l_rows
,
nc
)
endif
enddo
deallocate
(
qtmp
,
stat
=
istat
,
errmsg
=
errorMessage
)
check_deallocate
(
"resort_ev: qtmp"
,
istat
,
errorMessage
)
end
subroutine
resort_ev_
&
&
PRECISION
#endif
#if 0
subroutine
transform_columns_
&
&
PRECISION
&
&(
obj
,
col1
,
col2
)
use
precision
use
elpa_abstract_impl
implicit
none
class
(
elpa_abstract_impl_t
),
intent
(
inout
)
::
obj
integer
(
kind
=
ik
)
::
col1
,
col2
integer
(
kind
=
ik
)
::
pc1
,
pc2
,
lc1
,
lc2
if
(
l_rows
==
0
)
return
! My processor column has no work to do
pc1
=
p_col
(
col1
)
lc1
=
l_col
(
col1
)
pc2
=
p_col
(
col2
)
lc2
=
l_col
(
col2
)
if
(
pc1
==
my_pcol
)
then
if
(
pc2
==
my_pcol
)
then
! both columns are local
tmp
(
1
:
l_rows
)
=
q
(
l_rqs
:
l_rqe
,
lc1
)
*
qtrans
(
1
,
1
)
+
q
(
l_rqs
:
l_rqe
,
lc2
)
*
qtrans
(
2
,
1
)
q
(
l_rqs
:
l_rqe
,
lc2
)
=
q
(
l_rqs
:
l_rqe
,
lc1
)
*
qtrans
(
1
,
2
)
+
q
(
l_rqs
:
l_rqe
,
lc2
)
*
qtrans
(
2
,
2
)
q
(
l_rqs
:
l_rqe
,
lc1
)
=
tmp
(
1
:
l_rows
)
else
#ifdef WITH_MPI
call
obj
%
timer
%
start
(
"mpi_communication"
)
call
mpi_sendrecv
(
q
(
l_rqs
,
lc1
),
int
(
l_rows
,
kind
=
MPI_KIND
),
MPI_REAL_PRECISION
,
pc2
,
1_MPI_KIND
,
&
tmp
,
int
(
l_rows
,
kind
=
MPI_KIND
),
MPI_REAL_PRECISION
,
pc2
,
1_MPI_KIND
,
&
int
(
mpi_comm_cols
,
kind
=
MPI_KIND
),
MPI_STATUS_IGNORE
,
mpierr
)
call
obj
%
timer
%
stop
(
"mpi_communication"
)
#else /* WITH_MPI */
tmp
(
1
:
l_rows
)
=
q
(
l_rqs
:
l_rqe
,
lc1
)
#endif /* WITH_MPI */
q
(
l_rqs
:
l_rqe
,
lc1
)
=
q
(
l_rqs
:
l_rqe
,
lc1
)
*
qtrans
(
1
,
1
)
+
tmp
(
1
:
l_rows
)
*
qtrans
(
2
,
1
)
endif
else
if
(
pc2
==
my_pcol
)
then
#ifdef WITH_MPI
call
obj
%
timer
%
start
(
"mpi_communication"
)
call
mpi_sendrecv
(
q
(
l_rqs
,
lc2
),
int
(
l_rows
,
kind
=
MPI_KIND
),
MPI_REAL_PRECISION
,
pc1
,
1_MPI_KIND
,
&
tmp
,
int
(
l_rows
,
kind
=
MPI_KIND
),
MPI_REAL_PRECISION
,
pc1
,
1_MPI_KIND
,
&
int
(
mpi_comm_cols
,
kind
=
MPI_KIND
),
MPI_STATUS_IGNORE
,
mpierr
)
call
obj
%
timer
%
stop
(
"mpi_communication"
)
#else /* WITH_MPI */
tmp
(
1
:
l_rows
)
=
q
(
l_rqs
:
l_rqe
,
lc2
)
#endif /* WITH_MPI */
q
(
l_rqs
:
l_rqe
,
lc2
)
=
tmp
(
1
:
l_rows
)
*
qtrans
(
1
,
2
)
+
q
(
l_rqs
:
l_rqe
,
lc2
)
*
qtrans
(
2
,
2
)
endif
end
subroutine
transform_columns_
&
&
PRECISION
#endif
#if 0
subroutine
global_gather_
&
&
PRECISION
&
&(
obj
,
z
,
n
)
! This routine sums up z over all processors.
! It should only be used for gathering distributed results,
! i.e. z(i) should be nonzero on exactly 1 processor column,
! otherways the results may be numerically different on different columns
use
precision
use
elpa_abstract_impl
implicit
none
class
(
elpa_abstract_impl_t
),
intent
(
inout
)
::
obj
integer
(
kind
=
ik
)
::
n
real
(
kind
=
REAL_DATATYPE
)
::
z
(
n
)
real
(
kind
=
REAL_DATATYPE
)
::
tmp
(
n
)
if
(
npc_n
==
1
.and.
np_rows
==
1
)
return
! nothing to do
! Do an mpi_allreduce over processor rows
#ifdef WITH_MPI
call
obj
%
timer
%
start
(
"mpi_communication"
)
call
mpi_allreduce
(
z
,
tmp
,
int
(
n
,
kind
=
MPI_KIND
),
MPI_REAL_PRECISION
,
MPI_SUM
,
int
(
mpi_comm_rows
,
kind
=
MPI_KIND
),
mpierr
)
call
obj
%
timer
%
stop
(
"mpi_communication"
)
#else /* WITH_MPI */
tmp
=
z
#endif /* WITH_MPI */
! If only 1 processor column, we are done
if
(
npc_n
==
1
)
then
z
(:)
=
tmp
(:)
return
endif
! If all processor columns are involved, we can use mpi_allreduce
if
(
npc_n
==
np_cols
)
then
#ifdef WITH_MPI
call
obj
%
timer
%
start
(
"mpi_communication"
)
call
mpi_allreduce
(
tmp
,
z
,
int
(
n
,
kind
=
MPI_KIND
),
MPI_REAL_PRECISION
,
MPI_SUM
,
int
(
mpi_comm_cols
,
kind
=
MPI_KIND
),
mpierr
)
call
obj
%
timer
%
stop
(
"mpi_communication"
)
#else /* WITH_MPI */
tmp
=
z
#endif /* WITH_MPI */
return
endif
! Do a ring send over processor columns
z
(:)
=
0
do
np
=
1
,
npc_n
z
(:)
=
z
(:)
+
tmp
(:)
#ifdef WITH_MPI
call
obj
%
timer
%
start
(
"mpi_communication"
)
call
MPI_Sendrecv_replace
(
z
,
int
(
n
,
kind
=
MPI_KIND
),
MPI_REAL_PRECISION
,
int
(
np_next
,
kind
=
MPI_KIND
),
1111_MPI_KIND
,
&
int
(
np_prev
,
kind
=
MPI_KIND
),
1111_MPI_KIND
,
&
int
(
mpi_comm_cols
,
kind
=
MPI_KIND
),
MPI_STATUS_IGNORE
,
mpierr
)
call
obj
%
timer
%
stop
(
"mpi_communication"
)
#endif /* WITH_MPI */
enddo
end
subroutine
global_gather_
&
&
PRECISION
#endif
#if 0
subroutine
global_product_
&
&
PRECISION
&
&(
obj
,
z
,
n
)
! This routine calculates the global product of z.
use
precision
use
elpa_abstract_impl
implicit
none
class
(
elpa_abstract_impl_t
),
intent
(
inout
)
::
obj
integer
(
kind
=
ik
)
::
n
real
(
kind
=
REAL_DATATYPE
)
::
z
(
n
)
real
(
kind
=
REAL_DATATYPE
)
::
tmp
(
n
)
if
(
npc_n
==
1
.and.
np_rows
==
1
)
return
! nothing to do
! Do an mpi_allreduce over processor rows
#ifdef WITH_MPI
call
obj
%
timer
%
start
(
"mpi_communication"
)
call
mpi_allreduce
(
z
,
tmp
,
int
(
n
,
kind
=
MPI_KIND
),
MPI_REAL_PRECISION
,
MPI_PROD
,
int
(
mpi_comm_rows
,
kind
=
MPI_KIND
),
mpierr
)
call
obj
%
timer
%
stop
(
"mpi_communication"
)
#else /* WITH_MPI */
tmp
=
z
#endif /* WITH_MPI */
! If only 1 processor column, we are done
if
(
npc_n
==
1
)
then
z
(:)
=
tmp
(:)
return
endif
! If all processor columns are involved, we can use mpi_allreduce
if
(
npc_n
==
np_cols
)
then
#ifdef WITH_MPI
call
obj
%
timer
%
start
(
"mpi_communication"
)
call
mpi_allreduce
(
tmp
,
z
,
int
(
n
,
kind
=
MPI_KIND
),
MPI_REAL_PRECISION
,
MPI_PROD
,
int
(
mpi_comm_cols
,
kind
=
MPI_KIND
),
mpierr
)
call
obj
%
timer
%
stop
(
"mpi_communication"
)
#else /* WITH_MPI */
z
=
tmp
#endif /* WITH_MPI */
return
endif
! We send all vectors to the first proc, do the product there
! and redistribute the result.
if
(
my_pcol
==
npc_0
)
then
z
(
1
:
n
)
=
tmp
(
1
:
n
)
do
np
=
npc_0
+1
,
npc_0
+
npc_n
-1
#ifdef WITH_MPI
call
obj
%
timer
%
start
(
"mpi_communication"
)
call
mpi_recv
(
tmp
,
int
(
n
,
kind
=
MPI_KIND
),
MPI_REAL_PRECISION
,
int
(
np
,
kind
=
MPI_KIND
),
1111_MPI_KIND
,
&
int
(
mpi_comm_cols
,
kind
=
MPI_KIND
),
MPI_STATUS_IGNORE
,
mpierr
)
call
obj
%
timer
%
stop
(
"mpi_communication"
)
#else /* WITH_MPI */
tmp
(
1
:
n
)
=
z
(
1
:
n
)
#endif /* WITH_MPI */
z
(
1
:
n
)
=
z
(
1
:
n
)
*
tmp
(
1
:
n
)
enddo
do
np
=
npc_0
+1
,
npc_0
+
npc_n
-1
#ifdef WITH_MPI
call
obj
%
timer
%
start
(
"mpi_communication"
)
call
mpi_send
(
z
,
int
(
n
,
kind
=
MPI_KIND
),
MPI_REAL_PRECISION
,
int
(
np
,
kind
=
MPI_KIND
),
1111_MPI_KIND
,
&
int
(
mpi_comm_cols
,
kind
=
MPI_KIND
),
mpierr
)
call
obj
%
timer
%
stop
(
"mpi_communication"
)
#endif /* WITH_MPI */
enddo
else
#ifdef WITH_MPI
call
obj
%
timer
%
start
(
"mpi_communication"
)
call
mpi_send
(
tmp
,
int
(
n
,
kind
=
MPI_KIND
),
MPI_REAL_PRECISION
,
int
(
npc_0
,
kind
=
MPI_KIND
),
1111_MPI_KIND
,
&
int
(
mpi_comm_cols
,
kind
=
MPI_KIND
),
mpierr
)
call
mpi_recv
(
z
,
int
(
n
,
kind
=
MPI_KIND
),
MPI_REAL_PRECISION
,
int
(
npc_0
,
kind
=
MPI_KIND
),
1111_MPI_KIND
,
&
int
(
mpi_comm_cols
,
kind
=
MPI_KIND
),
MPI_STATUS_IGNORE
,
mpierr
)
call
obj
%
timer
%
stop
(
"mpi_communication"
)
#else /* WITH_MPI */
z
(
1
:
n
)
=
tmp
(
1
:
n
)
#endif /* WITH_MPI */
endif
end
subroutine
global_product_
&
&
PRECISION
#endif
#if 0
subroutine
check_monotony_
&
&
PRECISION
&
&(
obj
,
n
,
d
,
text
,
wantDebug
,
success
)
! This is a test routine for checking if the eigenvalues are monotonically increasing.
! It is for debug purposes only, an error should never be triggered!
use
precision
use
elpa_abstract_impl
implicit
none
class
(
elpa_abstract_impl_t
),
intent
(
inout
)
::
obj
integer
(
kind
=
ik
)
::
n
real
(
kind
=
REAL_DATATYPE
)
::
d
(
n
)
character
*
(
*
)
::
text
integer
(
kind
=
ik
)
::
i
logical
,
intent
(
in
)
::
wantDebug
logical
,
intent
(
out
)
::
success
success
=
.true.
do
i
=
1
,
n
-1
if
(
d
(
i
+1
)
<
d
(
i
))
then
if
(
wantDebug
)
write
(
error_unit
,
'(a,a,i8,2g25.17)'
)
'ELPA1_check_monotony: Monotony error on '
,
text
,
i
,
d
(
i
),
d
(
i
+1
)
success
=
.false.
return
endif
enddo
end
subroutine
check_monotony_
&
&
PRECISION
#endif
end
subroutine
merge_systems_
&
&
PRECISION
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment