diff --git a/CMakeLists.txt b/CMakeLists.txt
index 73ab55f3d2..b36fec1540 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -16,9 +16,6 @@ set(CMAKE_MODULE_PATH "${LAPACK_SOURCE_DIR}/CMAKE" ${CMAKE_MODULE_PATH})
# Export all symbols on Windows when building shared libraries
SET(CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS TRUE)
-# Tell CMake that our Fortran sources are written in fixed format.
-set(CMAKE_Fortran_FORMAT FIXED)
-
# Set a default build type if none was specified
if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
message(STATUS "Setting build type to 'Release' as none was specified.")
diff --git a/DOCS/Doxyfile b/DOCS/Doxyfile
index f1d9c3d400..e6053c43bc 100644
--- a/DOCS/Doxyfile
+++ b/DOCS/Doxyfile
@@ -886,6 +886,7 @@ INPUT_ENCODING = UTF-8
FILE_PATTERNS = *.c \
*.f \
+ *.f90 \
*.h
# The RECURSIVE tag can be used to specify whether or not subdirectories should
diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt
index bd3be3c6e8..dd37f0ecdc 100644
--- a/SRC/CMakeLists.txt
+++ b/SRC/CMakeLists.txt
@@ -41,6 +41,7 @@ set(ALLAUX ilaenv.f ilaenv2stage.f ieeeck.f lsamen.f iparmq.f iparam2stage.F
../INSTALL/slamch.f)
set(SCLAUX
+ la_constants.f90
sbdsdc.f
sbdsqr.f sdisna.f slabad.f slacpy.f sladiv.f slae2.f slaebz.f
slaed0.f slaed1.f slaed2.f slaed3.f slaed4.f slaed5.f slaed6.f
@@ -59,6 +60,7 @@ set(SCLAUX
${SECOND_SRC})
set(DZLAUX
+ la_constants.f90
dbdsdc.f
dbdsqr.f ddisna.f dlabad.f dlacpy.f dladiv.f dlae2.f dlaebz.f
dlaed0.f dlaed1.f dlaed2.f dlaed3.f dlaed4.f dlaed5.f dlaed6.f
@@ -217,7 +219,7 @@ set(CLASRC
cposv.f cposvx.f cpotf2.f cpotrf2.f cpotri.f cpstrf.f cpstf2.f
cppcon.f cppequ.f cpprfs.f cppsv.f cppsvx.f cpptrf.f cpptri.f cpptrs.f
cptcon.f cpteqr.f cptrfs.f cptsv.f cptsvx.f cpttrf.f cpttrs.f cptts2.f
- crot.f cspcon.f cspmv.f cspr.f csprfs.f cspsv.f
+ crot.f90 cspcon.f cspmv.f cspr.f csprfs.f cspsv.f
cspsvx.f csptrf.f csptri.f csptrs.f csrscl.f cstedc.f
cstegr.f cstein.f csteqr.f csycon.f csymv.f
csyr.f csyrfs.f csysv.f csysvx.f csytf2.f csytrf.f csytri.f
@@ -412,7 +414,7 @@ set(ZLASRC
zposv.f zposvx.f zpotf2.f zpotrf.f zpotrf2.f zpotri.f zpotrs.f zpstrf.f zpstf2.f
zppcon.f zppequ.f zpprfs.f zppsv.f zppsvx.f zpptrf.f zpptri.f zpptrs.f
zptcon.f zpteqr.f zptrfs.f zptsv.f zptsvx.f zpttrf.f zpttrs.f zptts2.f
- zrot.f zspcon.f zspmv.f zspr.f zsprfs.f zspsv.f
+ zrot.f90 zspcon.f zspmv.f zspr.f zsprfs.f zspsv.f
zspsvx.f zsptrf.f zsptri.f zsptrs.f zdrscl.f zstedc.f
zstegr.f zstein.f zsteqr.f zsycon.f zsymv.f
zsyr.f zsyrfs.f zsysv.f zsysvx.f zsytf2.f zsytrf.f zsytri.f
diff --git a/SRC/Makefile b/SRC/Makefile
index 3895fe5963..25e47a365b 100644
--- a/SRC/Makefile
+++ b/SRC/Makefile
@@ -57,9 +57,19 @@
TOPSRCDIR = ..
include $(TOPSRCDIR)/make.inc
-.SUFFIXES: .F .o
-.F.o:
+ALLMOD = la_constants.mod
+
+.SUFFIXES: .f .F .f90 .F90 .o .mod
+%.o: %.f $(ALLMOD)
+ $(FC) $(FFLAGS) -c -o $@ $<
+%.o: %.F $(ALLMOD)
+ $(FC) $(FFLAGS) -c -o $@ $<
+%.o: %.f90 $(ALLMOD)
+ $(FC) $(FFLAGS) -c -o $@ $<
+%.o: %.F90 $(ALLMOD)
$(FC) $(FFLAGS) -c -o $@ $<
+.o.mod:
+ @true
ALLAUX = ilaenv.o ilaenv2stage.o ieeeck.o lsamen.o xerbla.o xerbla_array.o \
iparmq.o iparam2stage.o \
@@ -67,6 +77,7 @@ ALLAUX = ilaenv.o ilaenv2stage.o ieeeck.o lsamen.o xerbla.o xerbla_array.o \
../INSTALL/ilaver.o ../INSTALL/lsame.o ../INSTALL/slamch.o
SCLAUX = \
+ la_constants.o \
sbdsdc.o \
sbdsqr.o sdisna.o slabad.o slacpy.o sladiv.o slae2.o slaebz.o \
slaed0.o slaed1.o slaed2.o slaed3.o slaed4.o slaed5.o slaed6.o \
@@ -85,6 +96,7 @@ SCLAUX = \
../INSTALL/second_$(TIMER).o
DZLAUX = \
+ la_constants.o \
dbdsdc.o \
dbdsqr.o ddisna.o dlabad.o dlacpy.o dladiv.o dlae2.o dlaebz.o \
dlaed0.o dlaed1.o dlaed2.o dlaed3.o dlaed4.o dlaed5.o dlaed6.o \
@@ -608,7 +620,7 @@ endif
.PHONY: clean cleanobj cleanlib
clean: cleanobj cleanlib
cleanobj:
- rm -f *.o DEPRECATED/*.o
+ rm -f *.o *.mod DEPRECATED/*.o DEPRECATED/*.mod
cleanlib:
rm -f $(LAPACKLIB)
@@ -618,3 +630,8 @@ sla_wwaddw.o: sla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<
dla_wwaddw.o: dla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<
cla_wwaddw.o: cla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<
zla_wwaddw.o: zla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<
+
+# Modules
+la_constants.o: la_constants.f90
+ $(FC) $(FFLAGS) -c -o $@ $<
+
diff --git a/SRC/crot.f b/SRC/crot.f
deleted file mode 100644
index 17b825a746..0000000000
--- a/SRC/crot.f
+++ /dev/null
@@ -1,159 +0,0 @@
-*> \brief \b CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
-*
-* =========== DOCUMENTATION ===========
-*
-* Online html documentation available at
-* http://www.netlib.org/lapack/explore-html/
-*
-*> \htmlonly
-*> Download CROT + dependencies
-*>
-*> [TGZ]
-*>
-*> [ZIP]
-*>
-*> [TXT]
-*> \endhtmlonly
-*
-* Definition:
-* ===========
-*
-* SUBROUTINE CROT( N, CX, INCX, CY, INCY, C, S )
-*
-* .. Scalar Arguments ..
-* INTEGER INCX, INCY, N
-* REAL C
-* COMPLEX S
-* ..
-* .. Array Arguments ..
-* COMPLEX CX( * ), CY( * )
-* ..
-*
-*
-*> \par Purpose:
-* =============
-*>
-*> \verbatim
-*>
-*> CROT applies a plane rotation, where the cos (C) is real and the
-*> sin (S) is complex, and the vectors CX and CY are complex.
-*> \endverbatim
-*
-* Arguments:
-* ==========
-*
-*> \param[in] N
-*> \verbatim
-*> N is INTEGER
-*> The number of elements in the vectors CX and CY.
-*> \endverbatim
-*>
-*> \param[in,out] CX
-*> \verbatim
-*> CX is COMPLEX array, dimension (N)
-*> On input, the vector X.
-*> On output, CX is overwritten with C*X + S*Y.
-*> \endverbatim
-*>
-*> \param[in] INCX
-*> \verbatim
-*> INCX is INTEGER
-*> The increment between successive values of CY. INCX <> 0.
-*> \endverbatim
-*>
-*> \param[in,out] CY
-*> \verbatim
-*> CY is COMPLEX array, dimension (N)
-*> On input, the vector Y.
-*> On output, CY is overwritten with -CONJG(S)*X + C*Y.
-*> \endverbatim
-*>
-*> \param[in] INCY
-*> \verbatim
-*> INCY is INTEGER
-*> The increment between successive values of CY. INCX <> 0.
-*> \endverbatim
-*>
-*> \param[in] C
-*> \verbatim
-*> C is REAL
-*> \endverbatim
-*>
-*> \param[in] S
-*> \verbatim
-*> S is COMPLEX
-*> C and S define a rotation
-*> [ C S ]
-*> [ -conjg(S) C ]
-*> where C*C + S*CONJG(S) = 1.0.
-*> \endverbatim
-*
-* Authors:
-* ========
-*
-*> \author Univ. of Tennessee
-*> \author Univ. of California Berkeley
-*> \author Univ. of Colorado Denver
-*> \author NAG Ltd.
-*
-*> \ingroup complexOTHERauxiliary
-*
-* =====================================================================
- SUBROUTINE CROT( N, CX, INCX, CY, INCY, C, S )
-*
-* -- LAPACK auxiliary routine --
-* -- LAPACK is a software package provided by Univ. of Tennessee, --
-* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*
-* .. Scalar Arguments ..
- INTEGER INCX, INCY, N
- REAL C
- COMPLEX S
-* ..
-* .. Array Arguments ..
- COMPLEX CX( * ), CY( * )
-* ..
-*
-* =====================================================================
-*
-* .. Local Scalars ..
- INTEGER I, IX, IY
- COMPLEX STEMP
-* ..
-* .. Intrinsic Functions ..
- INTRINSIC CONJG
-* ..
-* .. Executable Statements ..
-*
- IF( N.LE.0 )
- $ RETURN
- IF( INCX.EQ.1 .AND. INCY.EQ.1 )
- $ GO TO 20
-*
-* Code for unequal increments or equal increments not equal to 1
-*
- IX = 1
- IY = 1
- IF( INCX.LT.0 )
- $ IX = ( -N+1 )*INCX + 1
- IF( INCY.LT.0 )
- $ IY = ( -N+1 )*INCY + 1
- DO 10 I = 1, N
- STEMP = C*CX( IX ) + S*CY( IY )
- CY( IY ) = C*CY( IY ) - CONJG( S )*CX( IX )
- CX( IX ) = STEMP
- IX = IX + INCX
- IY = IY + INCY
- 10 CONTINUE
- RETURN
-*
-* Code for both increments equal to 1
-*
- 20 CONTINUE
- DO 30 I = 1, N
- STEMP = C*CX( I ) + S*CY( I )
- CY( I ) = C*CY( I ) - CONJG( S )*CX( I )
- CX( I ) = STEMP
- 30 CONTINUE
- RETURN
- END
diff --git a/SRC/crot.f90 b/SRC/crot.f90
new file mode 100644
index 0000000000..f69c66ae5e
--- /dev/null
+++ b/SRC/crot.f90
@@ -0,0 +1,162 @@
+!> \brief \b CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
+!
+! =========== DOCUMENTATION ===========
+!
+! Online html documentation available at
+! http://www.netlib.org/lapack/explore-html/
+!
+! Definition:
+! ===========
+!
+! SUBROUTINE CROT( N, X, INCX, Y, INCY, C, S )
+!
+! .. Scalar Arguments ..
+! INTEGER INCX, INCY, N
+! REAL(wp) C
+! COMPLEX(wp) S
+! ..
+! .. Array Arguments ..
+! COMPLEX(wp) X( * ), Y( * )
+! ..
+!
+!
+!> \par Purpose:
+! =============
+!>
+!> \verbatim
+!>
+!> CROT applies a plane rotation to vectors x and y:
+!>
+!> [ c s ] [ x(1) x(2) ... x(n) ]
+!> [ -conjg(s) c ] [ y(1) y(2) ... y(n) ]
+!>
+!> where c is real, s is complex, and c**2 + conjg(s)*s = 1.
+!>
+!> Below, wp=>sp stands for single precision from LA_CONSTANTS module.
+!>
+!> \endverbatim
+!
+! Arguments:
+! ==========
+!
+!> \param[in] N
+!> \verbatim
+!> N is INTEGER
+!> The number of elements in the vectors X and Y.
+!> \endverbatim
+!>
+!> \param[in,out] X
+!> \verbatim
+!> X is COMPLEX(wp) array, dimension (N)
+!> On input, the vector X.
+!> On output, X is overwritten with C*X + S*Y.
+!> \endverbatim
+!>
+!> \param[in] INCX
+!> \verbatim
+!> INCX is INTEGER
+!> The increment between successive values of X. INCX <> 0.
+!> \endverbatim
+!>
+!> \param[in,out] Y
+!> \verbatim
+!> Y is COMPLEX(wp) array, dimension (N)
+!> On input, the vector Y.
+!> On output, Y is overwritten with -CONJG(S)*X + C*Y.
+!> \endverbatim
+!>
+!> \param[in] INCY
+!> \verbatim
+!> INCY is INTEGER
+!> The increment between successive values of Y. INCX <> 0.
+!> \endverbatim
+!>
+!> \param[in] C
+!> \verbatim
+!> C is REAL(wp)
+!> \endverbatim
+!>
+!> \param[in] S
+!> \verbatim
+!> S is COMPLEX(wp)
+!> \endverbatim
+!
+! Authors:
+! ========
+!
+!> \author Edward Anderson, Lockheed Martin
+!
+!> \date August 2016
+!
+!> \ingroup complexOTHERauxiliary
+!
+!> \par Contributors:
+! ==================
+!>
+!> Weslley Pereira, University of Colorado Denver, USA
+!
+!> \par Further Details:
+! =====================
+!>
+!> \verbatim
+!>
+!> Anderson E. (2017)
+!> Algorithm 978: Safe Scaling in the Level 1 BLAS
+!> ACM Trans Math Softw 44:1--28
+!> https://doi.org/10.1145/3061665
+!>
+!> \endverbatim
+!
+! =====================================================================
+subroutine CROT( n, x, incx, y, incy, c, s )
+ use LA_CONSTANTS, only: wp=>sp
+!
+! -- LAPACK auxiliary routine (version 3.10.0) --
+! -- LAPACK is a software package provided by Univ. of Tennessee, --
+! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+! March 2021
+!
+! .. Scalar Arguments ..
+ integer :: incx, incy, n
+ real(wp) :: c
+ complex(wp) :: s
+! ..
+! .. Array Arguments ..
+ complex(wp) :: x(*), y(*)
+! ..
+! .. Local Scalars ..
+ integer :: i, ix, iy
+ complex(wp) :: stmp
+! ..
+!
+! Quick return if possible
+!
+ if( n <= 0 ) return
+!
+ if( incx == 1 .and. incy == 1 ) then
+!
+! Special case for incx = incy = 1
+!
+ do i = 1, n
+ stmp = c*x(i) + s*y(i)
+ y(i) = c*y(i) - conjg(s)*x(i)
+ x(i) = stmp
+ end do
+ else
+!
+! General increment case
+!
+ ix = 1
+ iy = 1
+ if( incx < 0 ) ix = 1 - (n-1)*incx
+ if( incy < 0 ) iy = 1 - (n-1)*incy
+ do i = 1, n
+ stmp = c*x(ix) + s*y(iy)
+ y(iy) = c*y(iy) - conjg(s)*x(ix)
+ x(ix) = stmp
+ ix = ix + incx
+ iy = iy + incy
+ end do
+ end if
+ return
+end subroutine
diff --git a/SRC/la_constants.f90 b/SRC/la_constants.f90
new file mode 100644
index 0000000000..3b1ab5e8c3
--- /dev/null
+++ b/SRC/la_constants.f90
@@ -0,0 +1,129 @@
+!> \brief \b LA_CONSTANTS is a module for the scaling constants for the compiled Fortran single and double precisions
+!
+! =========== DOCUMENTATION ===========
+!
+! Online html documentation available at
+! http://www.netlib.org/lapack/explore-html/
+!
+! Authors:
+! ========
+!
+!> \author Edward Anderson, Lockheed Martin
+!
+!> \date May 2016
+!
+!> \ingroup OTHERauxiliary
+!
+!> \par Contributors:
+! ==================
+!>
+!> Weslley Pereira, University of Colorado Denver, USA
+!> Nick Papior, Technical University of Denmark, DK
+!
+!> \par Further Details:
+! =====================
+!>
+!> \verbatim
+!>
+!> Anderson E. (2017)
+!> Algorithm 978: Safe Scaling in the Level 1 BLAS
+!> ACM Trans Math Softw 44:1--28
+!> https://doi.org/10.1145/3061665
+!>
+!> Blue, James L. (1978)
+!> A Portable Fortran Program to Find the Euclidean Norm of a Vector
+!> ACM Trans Math Softw 4:15--23
+!> https://doi.org/10.1145/355769.355771
+!>
+!> \endverbatim
+!
+module LA_CONSTANTS
+! -- LAPACK auxiliary module --
+! -- LAPACK is a software package provided by Univ. of Tennessee, --
+! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+
+! Standard constants for
+ integer, parameter :: sp = kind(1.e0)
+
+ real(sp), parameter :: szero = 0.0_sp
+ real(sp), parameter :: shalf = 0.5_sp
+ real(sp), parameter :: sone = 1.0_sp
+ real(sp), parameter :: stwo = 2.0_sp
+ real(sp), parameter :: sthree = 3.0_sp
+ real(sp), parameter :: sfour = 4.0_sp
+ real(sp), parameter :: seight = 8.0_sp
+ real(sp), parameter :: sten = 10.0_sp
+ complex(sp), parameter :: czero = ( 0.0_sp, 0.0_sp )
+ complex(sp), parameter :: chalf = ( 0.5_sp, 0.0_sp )
+ complex(sp), parameter :: cone = ( 1.0_sp, 0.0_sp )
+ character*1, parameter :: sprefix = 'S'
+ character*1, parameter :: cprefix = 'C'
+
+! Scaling constants
+ real(sp), parameter :: sulp = epsilon(0._sp)
+ real(sp), parameter :: seps = sulp * 0.5_sp
+ real(sp), parameter :: ssafmin = real(radix(0._sp),sp)**max( &
+ minexponent(0._sp)-1, &
+ 1-maxexponent(0._sp) &
+ )
+ real(sp), parameter :: ssafmax = sone / ssafmin
+ real(sp), parameter :: ssmlnum = ssafmin / sulp
+ real(sp), parameter :: sbignum = ssafmax * sulp
+ real(sp), parameter :: srtmin = sqrt(ssmlnum)
+ real(sp), parameter :: srtmax = sqrt(sbignum)
+
+! Blue's scaling constants
+ real(sp), parameter :: stsml = real(radix(0._sp), sp)**ceiling( &
+ (minexponent(0._sp) - 1) * 0.5_sp)
+ real(sp), parameter :: stbig = real(radix(0._sp), sp)**floor( &
+ (maxexponent(0._sp) - digits(0._sp) + 1) * 0.5_sp)
+! ssml = 1/s, where s was defined in https://doi.org/10.1145/355769.355771
+ real(sp), parameter :: sssml = real(radix(0._sp), sp)**( - floor( &
+ (minexponent(0._sp) - 1) * 0.5_sp))
+! ssml = 1/S, where S was defined in https://doi.org/10.1145/355769.355771
+ real(sp), parameter :: ssbig = real(radix(0._sp), sp)**( - ceiling( &
+ (maxexponent(0._sp) - digits(0._sp) + 1) * 0.5_sp))
+
+! Standard constants for
+ integer, parameter :: dp = kind(1.d0)
+
+ real(dp), parameter :: dzero = 0.0_dp
+ real(dp), parameter :: dhalf = 0.5_dp
+ real(dp), parameter :: done = 1.0_dp
+ real(dp), parameter :: dtwo = 2.0_dp
+ real(dp), parameter :: dthree = 3.0_dp
+ real(dp), parameter :: dfour = 4.0_dp
+ real(dp), parameter :: deight = 8.0_dp
+ real(dp), parameter :: dten = 10.0_dp
+ complex(dp), parameter :: zzero = ( 0.0_dp, 0.0_dp )
+ complex(dp), parameter :: zhalf = ( 0.5_dp, 0.0_dp )
+ complex(dp), parameter :: zone = ( 1.0_dp, 0.0_dp )
+ character*1, parameter :: dprefix = 'D'
+ character*1, parameter :: zprefix = 'Z'
+
+! Scaling constants
+ real(dp), parameter :: dulp = epsilon(0._dp)
+ real(dp), parameter :: deps = dulp * 0.5_dp
+ real(dp), parameter :: dsafmin = real(radix(0._dp),dp)**max( &
+ minexponent(0._dp)-1, &
+ 1-maxexponent(0._dp) &
+ )
+ real(dp), parameter :: dsafmax = done / dsafmin
+ real(dp), parameter :: dsmlnum = dsafmin / dulp
+ real(dp), parameter :: dbignum = dsafmax * dulp
+ real(dp), parameter :: drtmin = sqrt(dsmlnum)
+ real(dp), parameter :: drtmax = sqrt(dbignum)
+
+! Blue's scaling constants
+ real(dp), parameter :: dtsml = real(radix(0._dp), dp)**ceiling( &
+ (minexponent(0._dp) - 1) * 0.5_dp)
+ real(dp), parameter :: dtbig = real(radix(0._dp), dp)**floor( &
+ (maxexponent(0._dp) - digits(0._dp) + 1) * 0.5_dp)
+! ssml = 1/s, where s was defined in https://doi.org/10.1145/355769.355771
+ real(dp), parameter :: dssml = real(radix(0._dp), dp)**( - floor( &
+ (minexponent(0._dp) - 1) * 0.5_dp))
+! ssml = 1/S, where S was defined in https://doi.org/10.1145/355769.355771
+ real(dp), parameter :: dsbig = real(radix(0._dp), dp)**( - ceiling( &
+ (maxexponent(0._dp) - digits(0._dp) + 1) * 0.5_dp))
+
+end module LA_CONSTANTS
diff --git a/SRC/la_constants32.f90 b/SRC/la_constants32.f90
new file mode 100644
index 0000000000..d495709e92
--- /dev/null
+++ b/SRC/la_constants32.f90
@@ -0,0 +1,40 @@
+module LA_CONSTANTS32
+!
+! -- BLAS/LAPACK module --
+! May 06, 2016
+!
+! Standard constants
+!
+ integer, parameter :: wp = 4
+ real(wp), parameter :: zero = 0.0_wp
+ real(wp), parameter :: half = 0.5_wp
+ real(wp), parameter :: one = 1.0_wp
+ real(wp), parameter :: two = 2.0_wp
+ real(wp), parameter :: three = 3.0_wp
+ real(wp), parameter :: four = 4.0_wp
+ real(wp), parameter :: eight = 8.0_wp
+ real(wp), parameter :: ten = 10.0_wp
+ complex(wp), parameter :: czero = ( 0.0_wp, 0.0_wp )
+ complex(wp), parameter :: chalf = ( 0.5_wp, 0.0_wp )
+ complex(wp), parameter :: cone = ( 1.0_wp, 0.0_wp )
+ character*1, parameter :: sprefix = 'S'
+ character*1, parameter :: cprefix = 'C'
+!
+! Model parameters
+!
+ real(wp), parameter :: eps = 0.5960464478E-07_wp
+ real(wp), parameter :: ulp = 0.1192092896E-06_wp
+ real(wp), parameter :: safmin = 0.1175494351E-37_wp
+ real(wp), parameter :: safmax = 0.8507059173E+38_wp
+ real(wp), parameter :: smlnum = 0.9860761315E-31_wp
+ real(wp), parameter :: bignum = 0.1014120480E+32_wp
+ real(wp), parameter :: rtmin = 0.3140184864E-15_wp
+ real(wp), parameter :: rtmax = 0.3184525782E+16_wp
+!
+! Blue's scaling constants
+!
+ real(wp), parameter :: tsml = 0.1084202172E-18_wp
+ real(wp), parameter :: tbig = 0.4503599627E+16_wp
+ real(wp), parameter :: ssml = 0.3777893186E+23_wp
+ real(wp), parameter :: sbig = 0.1323488980E-22_wp
+end module LA_CONSTANTS32
diff --git a/SRC/zrot.f b/SRC/zrot.f
deleted file mode 100644
index 3401b2638a..0000000000
--- a/SRC/zrot.f
+++ /dev/null
@@ -1,159 +0,0 @@
-*> \brief \b ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
-*
-* =========== DOCUMENTATION ===========
-*
-* Online html documentation available at
-* http://www.netlib.org/lapack/explore-html/
-*
-*> \htmlonly
-*> Download ZROT + dependencies
-*>
-*> [TGZ]
-*>
-*> [ZIP]
-*>
-*> [TXT]
-*> \endhtmlonly
-*
-* Definition:
-* ===========
-*
-* SUBROUTINE ZROT( N, CX, INCX, CY, INCY, C, S )
-*
-* .. Scalar Arguments ..
-* INTEGER INCX, INCY, N
-* DOUBLE PRECISION C
-* COMPLEX*16 S
-* ..
-* .. Array Arguments ..
-* COMPLEX*16 CX( * ), CY( * )
-* ..
-*
-*
-*> \par Purpose:
-* =============
-*>
-*> \verbatim
-*>
-*> ZROT applies a plane rotation, where the cos (C) is real and the
-*> sin (S) is complex, and the vectors CX and CY are complex.
-*> \endverbatim
-*
-* Arguments:
-* ==========
-*
-*> \param[in] N
-*> \verbatim
-*> N is INTEGER
-*> The number of elements in the vectors CX and CY.
-*> \endverbatim
-*>
-*> \param[in,out] CX
-*> \verbatim
-*> CX is COMPLEX*16 array, dimension (N)
-*> On input, the vector X.
-*> On output, CX is overwritten with C*X + S*Y.
-*> \endverbatim
-*>
-*> \param[in] INCX
-*> \verbatim
-*> INCX is INTEGER
-*> The increment between successive values of CY. INCX <> 0.
-*> \endverbatim
-*>
-*> \param[in,out] CY
-*> \verbatim
-*> CY is COMPLEX*16 array, dimension (N)
-*> On input, the vector Y.
-*> On output, CY is overwritten with -CONJG(S)*X + C*Y.
-*> \endverbatim
-*>
-*> \param[in] INCY
-*> \verbatim
-*> INCY is INTEGER
-*> The increment between successive values of CY. INCX <> 0.
-*> \endverbatim
-*>
-*> \param[in] C
-*> \verbatim
-*> C is DOUBLE PRECISION
-*> \endverbatim
-*>
-*> \param[in] S
-*> \verbatim
-*> S is COMPLEX*16
-*> C and S define a rotation
-*> [ C S ]
-*> [ -conjg(S) C ]
-*> where C*C + S*CONJG(S) = 1.0.
-*> \endverbatim
-*
-* Authors:
-* ========
-*
-*> \author Univ. of Tennessee
-*> \author Univ. of California Berkeley
-*> \author Univ. of Colorado Denver
-*> \author NAG Ltd.
-*
-*> \ingroup complex16OTHERauxiliary
-*
-* =====================================================================
- SUBROUTINE ZROT( N, CX, INCX, CY, INCY, C, S )
-*
-* -- LAPACK auxiliary routine --
-* -- LAPACK is a software package provided by Univ. of Tennessee, --
-* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*
-* .. Scalar Arguments ..
- INTEGER INCX, INCY, N
- DOUBLE PRECISION C
- COMPLEX*16 S
-* ..
-* .. Array Arguments ..
- COMPLEX*16 CX( * ), CY( * )
-* ..
-*
-* =====================================================================
-*
-* .. Local Scalars ..
- INTEGER I, IX, IY
- COMPLEX*16 STEMP
-* ..
-* .. Intrinsic Functions ..
- INTRINSIC DCONJG
-* ..
-* .. Executable Statements ..
-*
- IF( N.LE.0 )
- $ RETURN
- IF( INCX.EQ.1 .AND. INCY.EQ.1 )
- $ GO TO 20
-*
-* Code for unequal increments or equal increments not equal to 1
-*
- IX = 1
- IY = 1
- IF( INCX.LT.0 )
- $ IX = ( -N+1 )*INCX + 1
- IF( INCY.LT.0 )
- $ IY = ( -N+1 )*INCY + 1
- DO 10 I = 1, N
- STEMP = C*CX( IX ) + S*CY( IY )
- CY( IY ) = C*CY( IY ) - DCONJG( S )*CX( IX )
- CX( IX ) = STEMP
- IX = IX + INCX
- IY = IY + INCY
- 10 CONTINUE
- RETURN
-*
-* Code for both increments equal to 1
-*
- 20 CONTINUE
- DO 30 I = 1, N
- STEMP = C*CX( I ) + S*CY( I )
- CY( I ) = C*CY( I ) - DCONJG( S )*CX( I )
- CX( I ) = STEMP
- 30 CONTINUE
- RETURN
- END
diff --git a/SRC/zrot.f90 b/SRC/zrot.f90
new file mode 100644
index 0000000000..52612e07eb
--- /dev/null
+++ b/SRC/zrot.f90
@@ -0,0 +1,160 @@
+!> \brief \b ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
+!
+! =========== DOCUMENTATION ===========
+!
+! Online html documentation available at
+! http://www.netlib.org/lapack/explore-html/
+!
+! Definition:
+! ===========
+!
+! SUBROUTINE ZROT( N, X, INCX, Y, INCY, C, S )
+!
+! .. Scalar Arguments ..
+! INTEGER INCX, INCY, N
+! REAL(wp) C
+! COMPLEX(wp) S
+! ..
+! .. Array Arguments ..
+! COMPLEX(wp) X( * ), Y( * )
+! ..
+!
+!
+!> \par Purpose:
+! =============
+!>
+!> \verbatim
+!>
+!> ZROT applies a plane rotation to vectors x and y:
+!>
+!> [ c s ] [ x(1) x(2) ... x(n) ]
+!> [ -conjg(s) c ] [ y(1) y(2) ... y(n) ]
+!>
+!> where c is real, s is complex, and c**2 + conjg(s)*s = 1.
+!>
+!> Below, wp=>dp stands for double precision from LA_CONSTANTS module.
+!> \endverbatim
+!
+! Arguments:
+! ==========
+!
+!> \param[in] N
+!> \verbatim
+!> N is INTEGER
+!> The number of elements in the vectors X and Y.
+!> \endverbatim
+!>
+!> \param[in,out] X
+!> \verbatim
+!> X is COMPLEX(wp) array, dimension (N)
+!> On input, the vector X.
+!> On output, X is overwritten with C*X + S*Y.
+!> \endverbatim
+!>
+!> \param[in] INCX
+!> \verbatim
+!> INCX is INTEGER
+!> The increment between successive values of X. INCX <> 0.
+!> \endverbatim
+!>
+!> \param[in,out] Y
+!> \verbatim
+!> Y is COMPLEX(wp) array, dimension (N)
+!> On input, the vector Y.
+!> On output, Y is overwritten with -CONJG(S)*X + C*Y.
+!> \endverbatim
+!>
+!> \param[in] INCY
+!> \verbatim
+!> INCY is INTEGER
+!> The increment between successive values of Y. INCX <> 0.
+!> \endverbatim
+!>
+!> \param[in] C
+!> \verbatim
+!> C is REAL(wp)
+!> \endverbatim
+!>
+!> \param[in] S
+!> \verbatim
+!> S is COMPLEX(wp)
+!> \endverbatim
+!
+! Authors:
+! ========
+!
+!> \author Edward Anderson, Lockheed Martin
+!
+!> \date August 2016
+!
+!> \ingroup complexOTHERauxiliary
+!
+!> \par Contributors:
+! ==================
+!>
+!> Weslley Pereira, University of Colorado Denver, USA
+!
+!> \par Further Details:
+! =====================
+!>
+!> \verbatim
+!>
+!> Anderson E. (2017)
+!> Algorithm 978: Safe Scaling in the Level 1 BLAS
+!> ACM Trans Math Softw 44:1--28
+!> https://doi.org/10.1145/3061665
+!>
+!> \endverbatim
+!
+! =====================================================================
+subroutine ZROT( n, x, incx, y, incy, c, s )
+ use LA_CONSTANTS, only: wp=>dp
+!
+! Updated Level 1 BLAS
+! E. Anderson
+! August 11, 2016
+!
+! .. Scalar Arguments ..
+ integer :: incx, incy, n
+ real(wp) :: c
+ complex(wp) :: s
+! ..
+! .. Array Arguments ..
+ complex(wp) :: x(*), y(*)
+! ..
+! .. Local Scalars ..
+ integer :: i, ix, iy
+ complex(wp) :: stmp
+! ..
+!
+! Quick return if possible
+!
+ if( n <= 0 ) return
+!
+ if( incx == 1 .and. incy == 1 ) then
+!
+! Special case for incx = incy = 1
+!
+ do i = 1, n
+ stmp = c*x(i) + s*y(i)
+ y(i) = c*y(i) - conjg(s)*x(i)
+ x(i) = stmp
+ end do
+ else
+!
+! General increment case
+!
+ ix = 1
+ iy = 1
+ if( incx < 0 ) ix = 1 - (n-1)*incx
+ if( incy < 0 ) iy = 1 - (n-1)*incy
+ do i = 1, n
+ stmp = c*x(ix) + s*y(iy)
+ y(iy) = c*y(iy) - conjg(s)*x(ix)
+ x(ix) = stmp
+ ix = ix + incx
+ iy = iy + incy
+ end do
+ end if
+ return
+end subroutine