Commit 031b5b3a authored by Andreas Marek's avatar Andreas Marek

Start to unify real/complex case: elpa1_tools

parent 3e58514a
......@@ -887,11 +887,10 @@ EXTRA_DIST = \
src/elpa1_compute_real_template.X90 \
src/elpa1_merge_systems_real_template.X90 \
src/elpa1_solve_tridi_real_template.X90 \
src/elpa1_tools_real_template.X90 \
src/elpa1_tools_template.X90 \
src/elpa1_trans_ev_template.X90 \
src/elpa1_tridiag_template.X90 \
src/elpa1_compute_complex_template.X90 \
src/elpa1_tools_complex_template.X90 \
src/elpa2_compute_real_template.X90 \
src/elpa2_compute_complex_template.X90 \
src/precision_macros.h \
......
......@@ -58,7 +58,7 @@
#undef REALCASE
#include "elpa1_tridiag_template.X90"
#include "elpa1_trans_ev_template.X90"
#include "elpa1_tools_template.X90"
#undef COMPLEXCASE
#include "elpa1_tools_complex_template.X90"
#define ALREADY_DEFINED 1
......@@ -67,5 +67,7 @@
#undef REALCASE
#include "elpa1_solve_tridi_real_template.X90"
#include "elpa1_merge_systems_real_template.X90"
#include "elpa1_tools_real_template.X90"
#define REALCASE 1
#undef COMPLEXCASE
#include "elpa1_tools_template.X90"
#undef REALCASE
subroutine hh_transform_complex_PRECISION(alpha, xnorm_sq, xf, tau)
! Similar to LAPACK routine ZLARFP, but uses ||x||**2 instead of x(:)
! and returns the factor xf by which x has to be scaled.
! It also hasn't the special handling for numbers < 1.d-300 or > 1.d150
! since this would be expensive for the parallel implementation.
use precision
#ifdef HAVE_DETAILED_TIMINGS
use timings
#else
use timings_dummy
#endif
implicit none
complex(kind=COMPLEX_DATATYPE), intent(inout) :: alpha
real(kind=REAL_DATATYPE), intent(in) :: xnorm_sq
complex(kind=COMPLEX_DATATYPE), intent(out) :: xf, tau
real(kind=REAL_DATATYPE) :: ALPHR, ALPHI, BETA
call timer%start("hh_transform_complex" // PRECISION_SUFFIX)
ALPHR = real( ALPHA, kind=REAL_DATATYPE )
ALPHI = PRECISION_IMAG( ALPHA )
if ( XNORM_SQ==0. .AND. ALPHI==0. ) then
if ( ALPHR>=0. ) then
TAU = 0.
else
TAU = 2.
ALPHA = -ALPHA
endif
XF = 0.
else
BETA = SIGN( SQRT( ALPHR**2 + ALPHI**2 + XNORM_SQ ), ALPHR )
ALPHA = ALPHA + BETA
IF ( BETA<0 ) THEN
BETA = -BETA
TAU = -ALPHA / BETA
ELSE
ALPHR = ALPHI * (ALPHI/real( ALPHA , kind=KIND_PRECISION))
ALPHR = ALPHR + XNORM_SQ/real( ALPHA, kind=KIND_PRECISION )
TAU = PRECISION_CMPLX( ALPHR/BETA, -ALPHI/BETA )
ALPHA = PRECISION_CMPLX( -ALPHR, ALPHI )
END IF
XF = CONST_REAL_1_0/ALPHA
ALPHA = BETA
endif
call timer%stop("hh_transform_complex" // PRECISION_SUFFIX)
end subroutine hh_transform_complex_PRECISION
......@@ -51,7 +51,9 @@
! with their original authors, but shall adhere to the licensing terms
! distributed along with the original code in the file "COPYING".
#endif
#if REALCASE == 1
subroutine M_v_add_s_PRECISION(v,n,s)
use precision
implicit none
......@@ -225,23 +227,69 @@
end subroutine M_solve_secular_equation_PRECISION
!-------------------------------------------------------------------------------
#endif
#if REALCASE == 1
subroutine M_hh_transform_real_PRECISION(alpha, xnorm_sq, xf, tau)
! Similar to LAPACK routine DLARFP, but uses ||x||**2 instead of x(:)
#endif
#if COMPLEXCASE == 1
subroutine hh_transform_complex_PRECISION(alpha, xnorm_sq, xf, tau)
! Similar to LAPACK routine ZLARFP, but uses ||x||**2 instead of x(:)
#endif
! and returns the factor xf by which x has to be scaled.
! It also hasn't the special handling for numbers < 1.d-300 or > 1.d150
! since this would be expensive for the parallel implementation.
use precision
#ifdef HAVE_DETAILED_TIMINGS
use timings
#else
use timings_dummy
#endif
implicit none
real(kind=REAL_DATATYPE), intent(inout) :: alpha
real(kind=REAL_DATATYPE), intent(in) :: xnorm_sq
real(kind=REAL_DATATYPE), intent(out) :: xf, tau
#if REALCASE == 1
real(kind=REAL_DATATYPE), intent(inout) :: alpha
#endif
#if COMPLEXCASE == 1
complex(kind=COMPLEX_DATATYPE), intent(inout) :: alpha
#endif
real(kind=REAL_DATATYPE), intent(in) :: xnorm_sq
#if REALCASE == 1
real(kind=REAL_DATATYPE), intent(out) :: xf, tau
#endif
#if COMPLEXCASE == 1
complex(kind=COMPLEX_DATATYPE), intent(out) :: xf, tau
real(kind=REAL_DATATYPE) :: ALPHR, ALPHI
#endif
real(kind=REAL_DATATYPE) :: BETA
real(kind=REAL_DATATYPE) :: BETA
#if REALCASE == 1
call timer%start("hh_transform_real" // M_PRECISION_SUFFIX )
#endif
#if COMPLEXCASE == 1
call timer%start("hh_transform_complex" // PRECISION_SUFFIX )
#endif
#if COMPLEXCASE == 1
ALPHR = real( ALPHA, kind=REAL_DATATYPE )
ALPHI = PRECISION_IMAG( ALPHA )
#endif
#if REALCASE == 1
if ( XNORM_SQ==0. ) then
#endif
#if COMPLEXCASE == 1
if ( XNORM_SQ==0. .AND. ALPHI==0. ) then
#endif
#if REALCASE == 1
if ( ALPHA>=0. ) then
#endif
#if COMPLEXCASE == 1
if ( ALPHR>=0. ) then
#endif
TAU = 0.
else
TAU = 2.
......@@ -251,17 +299,48 @@
else
#if REALCASE == 1
BETA = SIGN( SQRT( ALPHA**2 + XNORM_SQ ), ALPHA )
#endif
#if COMPLEXCASE == 1
BETA = SIGN( SQRT( ALPHR**2 + ALPHI**2 + XNORM_SQ ), ALPHR )
#endif
ALPHA = ALPHA + BETA
IF ( BETA<0 ) THEN
BETA = -BETA
TAU = -ALPHA / BETA
TAU = -ALPHA / BETA
ELSE
#if REALCASE == 1
ALPHA = XNORM_SQ / ALPHA
#endif
#if COMPLEXCASE == 1
ALPHR = ALPHI * (ALPHI/real( ALPHA , kind=KIND_PRECISION))
ALPHR = ALPHR + XNORM_SQ/real( ALPHA, kind=KIND_PRECISION )
#endif
#if REALCASE == 1
TAU = ALPHA / BETA
ALPHA = -ALPHA
#endif
#if COMPLEXCASE == 1
TAU = PRECISION_CMPLX( ALPHR/BETA, -ALPHI/BETA )
ALPHA = PRECISION_CMPLX( -ALPHR, ALPHI )
#endif
END IF
XF = 1./ALPHA
XF = 1.0/ALPHA
ALPHA = BETA
endif
#if REALCASE == 1
call timer%stop("hh_transform_real" // M_PRECISION_SUFFIX )
#endif
#if COMPLEXCASE == 1
call timer%stop("hh_transform_complex" // PRECISION_SUFFIX )
#endif
#if REALCASE == 1
end subroutine M_hh_transform_real_PRECISION
#endif
#if COMPLEXCASE == 1
end subroutine hh_transform_complex_PRECISION
#endif
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment