diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 05e913f2f..5b57c2ff2 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -519,6 +519,53 @@ inline EIGEN_MATHFUNC_RETVAL(atan2, Scalar) atan2(const Scalar& x, const Scalar& return EIGEN_MATHFUNC_IMPL(atan2, Scalar)::run(x, y); } +/**************************************************************************** +* Implementation of atanh2 * +****************************************************************************/ + +template +struct atanh2_default_impl +{ + typedef Scalar retval; + typedef typename NumTraits::Real RealScalar; + static inline Scalar run(const Scalar& x, const Scalar& y) + { + using std::abs; + using std::log; + using std::sqrt; + Scalar z = x / y; + if (abs(z) > sqrt(NumTraits::epsilon())) + return RealScalar(0.5) * log((y + x) / (y - x)); + else + return z + z*z*z / RealScalar(3); + } +}; + +template +struct atanh2_default_impl +{ + static inline Scalar run(const Scalar&, const Scalar&) + { + EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar) + return Scalar(0); + } +}; + +template +struct atanh2_impl : atanh2_default_impl::IsInteger> {}; + +template +struct atanh2_retval +{ + typedef Scalar type; +}; + +template +inline EIGEN_MATHFUNC_RETVAL(atanh2, Scalar) atanh2(const Scalar& x, const Scalar& y) +{ + return EIGEN_MATHFUNC_IMPL(atanh2, Scalar)::run(x, y); +} + /**************************************************************************** * Implementation of pow * ****************************************************************************/ diff --git a/blas/GeneralRank1Update.h b/blas/GeneralRank1Update.h new file mode 100644 index 000000000..a3301ed92 --- /dev/null +++ b/blas/GeneralRank1Update.h @@ -0,0 +1,30 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Chen-Pang He +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_GENERAL_RANK1UPDATE_H +#define EIGEN_GENERAL_RANK1UPDATE_H + +namespace internal { + +/* Optimized matrix += alpha * uv' */ +template +struct general_rank1_update +{ + static void run(Index rows, Index cols, Scalar* mat, Index stride, const Scalar* u, const Scalar* v, Scalar alpha) + { + typedef Matrix PlainVector; + internal::conj_if cj; + for (Index i=0; i(mat+stride*i,rows) += alpha * cj(v[i]) * Map(u,rows); + } +}; + +} // end namespace internal + +#endif // EIGEN_GENERAL_RANK1UPDATE_H diff --git a/blas/Rank2Update.h b/blas/Rank2Update.h new file mode 100644 index 000000000..e7a5eeaba --- /dev/null +++ b/blas/Rank2Update.h @@ -0,0 +1,57 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Chen-Pang He +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_RANK2UPDATE_H +#define EIGEN_RANK2UPDATE_H + +namespace internal { + +/* Optimized selfadjoint matrix += alpha * uv' + conj(alpha)*vu' + * This is the low-level version of SelfadjointRank2Update.h + */ +template +struct rank2_update_selector; + +template +struct rank2_update_selector +{ + static void run(Index size, Scalar* mat, Index stride, const Scalar* _u, const Scalar* _v, Scalar alpha) + { + typedef Matrix PlainVector; + Map u(_u, size), v(_v, size); + + for (Index i=0; i(mat+stride*i, i+1) += + conj(alpha) * conj(_u[i]) * v.head(i+1) + + alpha * conj(_v[i]) * u.head(i+1); + } + } +}; + +template +struct rank2_update_selector +{ + static void run(Index size, Scalar* mat, Index stride, const Scalar* _u, const Scalar* _v, Scalar alpha) + { + typedef Matrix PlainVector; + Map u(_u, size), v(_v, size); + + for (Index i=0; i(mat+(stride+1)*i, size-i) += + conj(alpha) * conj(_u[i]) * v.tail(size-i) + + alpha * conj(_v[i]) * u.tail(size-i); + } + } +}; + +} // end namespace internal + +#endif // EIGEN_RANK2UPDATE_H diff --git a/blas/common.h b/blas/common.h index b598c4e45..26b4ed5a3 100644 --- a/blas/common.h +++ b/blas/common.h @@ -48,7 +48,7 @@ : ((X)=='L' || (X)=='l') ? LO \ : INVALID) -#define DIAG(X) ( ((X)=='N' || (X)=='N') ? NUNIT \ +#define DIAG(X) ( ((X)=='N' || (X)=='n') ? NUNIT \ : ((X)=='U' || (X)=='u') ? UNIT \ : INVALID) @@ -74,6 +74,8 @@ inline bool check_uplo(const char* uplo) namespace Eigen { #include "BandTriangularSolver.h" +#include "GeneralRank1Update.h" +#include "Rank2Update.h" } using namespace Eigen; diff --git a/blas/double.cpp b/blas/double.cpp index cad2f63ec..8fd0709ba 100644 --- a/blas/double.cpp +++ b/blas/double.cpp @@ -2,6 +2,7 @@ // for linear algebra. // // Copyright (C) 2009 Gael Guennebaud +// Copyright (C) 2012 Chen-Pang He // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -17,3 +18,16 @@ #include "level2_impl.h" #include "level2_real_impl.h" #include "level3_impl.h" + +double BLASFUNC(dsdot)(int* n, float* x, int* incx, float* y, int* incy) +{ + if(*n<=0) return 0; + + if(*incx==1 && *incy==1) return (vector(x,*n).cast().cwiseProduct(vector(y,*n).cast())).sum(); + else if(*incx>0 && *incy>0) return (vector(x,*n,*incx).cast().cwiseProduct(vector(y,*n,*incy).cast())).sum(); + else if(*incx<0 && *incy>0) return (vector(x,*n,-*incx).reverse().cast().cwiseProduct(vector(y,*n,*incy).cast())).sum(); + else if(*incx>0 && *incy<0) return (vector(x,*n,*incx).cast().cwiseProduct(vector(y,*n,-*incy).reverse().cast())).sum(); + else if(*incx<0 && *incy<0) return (vector(x,*n,-*incx).reverse().cast().cwiseProduct(vector(y,*n,-*incy).reverse().cast())).sum(); + else return 0; +} + diff --git a/blas/level2_cplx_impl.h b/blas/level2_cplx_impl.h index 7878f2a16..477f6d649 100644 --- a/blas/level2_cplx_impl.h +++ b/blas/level2_cplx_impl.h @@ -117,6 +117,21 @@ int EIGEN_BLAS_FUNC(hemv)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa */ int EIGEN_BLAS_FUNC(her)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pa, int *lda) { + typedef void (*functype)(int, Scalar*, int, const Scalar*, Scalar); + static functype func[2]; + + static bool init = false; + if(!init) + { + for(int k=0; k<2; ++k) + func[k] = 0; + + func[UP] = (selfadjoint_rank1_update::run); + func[LO] = (selfadjoint_rank1_update::run); + + init = true; + } + Scalar* x = reinterpret_cast(px); Scalar* a = reinterpret_cast(pa); RealScalar alpha = *reinterpret_cast(palpha); @@ -134,16 +149,11 @@ int EIGEN_BLAS_FUNC(her)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, Scalar* x_cpy = get_compact_vector(x, *n, *incx); - // TODO perform direct calls to underlying implementation -// if(UPLO(*uplo)==LO) matrix(a,*n,*n,*lda).selfadjointView().rankUpdate(vector(x_cpy,*n), alpha); -// else if(UPLO(*uplo)==UP) matrix(a,*n,*n,*lda).selfadjointView().rankUpdate(vector(x_cpy,*n), alpha); + int code = UPLO(*uplo); + if(code>=2 || func[code]==0) + return 0; - if(UPLO(*uplo)==LO) - for(int j=0;j<*n;++j) - matrix(a,*n,*n,*lda).col(j).tail(*n-j) += alpha * internal::conj(x_cpy[j]) * vector(x_cpy+j,*n-j); - else - for(int j=0;j<*n;++j) - matrix(a,*n,*n,*lda).col(j).head(j+1) += alpha * internal::conj(x_cpy[j]) * vector(x_cpy,j+1); + func[code](*n, a, *lda, x_cpy, alpha); matrix(a,*n,*n,*lda).diagonal().imag().setZero(); @@ -161,6 +171,21 @@ int EIGEN_BLAS_FUNC(her)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, */ int EIGEN_BLAS_FUNC(her2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pa, int *lda) { + typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, Scalar); + static functype func[2]; + + static bool init = false; + if(!init) + { + for(int k=0; k<2; ++k) + func[k] = 0; + + func[UP] = (internal::rank2_update_selector::run); + func[LO] = (internal::rank2_update_selector::run); + + init = true; + } + Scalar* x = reinterpret_cast(px); Scalar* y = reinterpret_cast(py); Scalar* a = reinterpret_cast(pa); @@ -181,9 +206,11 @@ int EIGEN_BLAS_FUNC(her2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px Scalar* x_cpy = get_compact_vector(x, *n, *incx); Scalar* y_cpy = get_compact_vector(y, *n, *incy); - // TODO perform direct calls to underlying implementation - if(UPLO(*uplo)==LO) matrix(a,*n,*n,*lda).selfadjointView().rankUpdate(vector(x_cpy,*n),vector(y_cpy,*n),alpha); - else if(UPLO(*uplo)==UP) matrix(a,*n,*n,*lda).selfadjointView().rankUpdate(vector(x_cpy,*n),vector(y_cpy,*n),alpha); + int code = UPLO(*uplo); + if(code>=2 || func[code]==0) + return 0; + + func[code](*n, a, *lda, x_cpy, y_cpy, alpha); matrix(a,*n,*n,*lda).diagonal().imag().setZero(); @@ -222,8 +249,7 @@ int EIGEN_BLAS_FUNC(geru)(int *m, int *n, RealScalar *palpha, RealScalar *px, in Scalar* x_cpy = get_compact_vector(x,*m,*incx); Scalar* y_cpy = get_compact_vector(y,*n,*incy); - // TODO perform direct calls to underlying implementation - matrix(a,*m,*n,*lda) += alpha * vector(x_cpy,*m) * vector(y_cpy,*n).transpose(); + internal::general_rank1_update::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha); if(x_cpy!=x) delete[] x_cpy; if(y_cpy!=y) delete[] y_cpy; @@ -260,8 +286,7 @@ int EIGEN_BLAS_FUNC(gerc)(int *m, int *n, RealScalar *palpha, RealScalar *px, in Scalar* x_cpy = get_compact_vector(x,*m,*incx); Scalar* y_cpy = get_compact_vector(y,*n,*incy); - // TODO perform direct calls to underlying implementation - matrix(a,*m,*n,*lda) += alpha * vector(x_cpy,*m) * vector(y_cpy,*n).adjoint(); + internal::general_rank1_update::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha); if(x_cpy!=x) delete[] x_cpy; if(y_cpy!=y) delete[] y_cpy; diff --git a/blas/level2_impl.h b/blas/level2_impl.h index 7099cf96d..f1f7371ee 100644 --- a/blas/level2_impl.h +++ b/blas/level2_impl.h @@ -49,7 +49,8 @@ int EIGEN_BLAS_FUNC(gemv)(char *opa, int *m, int *n, RealScalar *palpha, RealSca int actual_m = *m; int actual_n = *n; - if(OP(*opa)!=NOTR) + int code = OP(*opa); + if(code!=NOTR) std::swap(actual_m,actual_n); Scalar* actual_b = get_compact_vector(b,actual_n,*incb); @@ -61,7 +62,9 @@ int EIGEN_BLAS_FUNC(gemv)(char *opa, int *m, int *n, RealScalar *palpha, RealSca else vector(actual_c, actual_m) *= beta; } - int code = OP(*opa); + if(code>=4 || func[code]==0) + return 0; + func[code](actual_m, actual_n, a, *lda, actual_b, 1, actual_c, 1, alpha); if(actual_b!=b) delete[] actual_b; @@ -416,42 +419,3 @@ int EIGEN_BLAS_FUNC(tbsv)(char *uplo, char *op, char *diag, int *n, int *k, Real // return 1; // } -/** DGER performs the rank 1 operation - * - * A := alpha*x*y' + A, - * - * where alpha is a scalar, x is an m element vector, y is an n element - * vector and A is an m by n matrix. - */ -int EIGEN_BLAS_FUNC(ger)(int *m, int *n, Scalar *palpha, Scalar *px, int *incx, Scalar *py, int *incy, Scalar *pa, int *lda) -{ - Scalar* x = reinterpret_cast(px); - Scalar* y = reinterpret_cast(py); - Scalar* a = reinterpret_cast(pa); - Scalar alpha = *reinterpret_cast(palpha); - - int info = 0; - if(*m<0) info = 1; - else if(*n<0) info = 2; - else if(*incx==0) info = 5; - else if(*incy==0) info = 7; - else if(*lda::run); + func[LO] = (selfadjoint_rank1_update::run); + + init = true; + } Scalar* x = reinterpret_cast(px); Scalar* c = reinterpret_cast(pc); @@ -86,18 +100,11 @@ int EIGEN_BLAS_FUNC(syr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, // if the increment is not 1, let's copy it to a temporary vector to enable vectorization Scalar* x_cpy = get_compact_vector(x,*n,*incx); - Matrix m2(matrix(c,*n,*n,*ldc)); - - // TODO check why this is not accurate enough for lapack tests -// if(UPLO(*uplo)==LO) matrix(c,*n,*n,*ldc).selfadjointView().rankUpdate(vector(x_cpy,*n), alpha); -// else if(UPLO(*uplo)==UP) matrix(c,*n,*n,*ldc).selfadjointView().rankUpdate(vector(x_cpy,*n), alpha); + int code = UPLO(*uplo); + if(code>=2 || func[code]==0) + return 0; - if(UPLO(*uplo)==LO) - for(int j=0;j<*n;++j) - matrix(c,*n,*n,*ldc).col(j).tail(*n-j) += alpha * x_cpy[j] * vector(x_cpy+j,*n-j); - else - for(int j=0;j<*n;++j) - matrix(c,*n,*n,*ldc).col(j).head(j+1) += alpha * x_cpy[j] * vector(x_cpy,j+1); + func[code](*n, c, *ldc, x_cpy, alpha); if(x_cpy!=x) delete[] x_cpy; @@ -121,6 +128,20 @@ int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px // // init = true; // } + typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, Scalar); + static functype func[2]; + + static bool init = false; + if(!init) + { + for(int k=0; k<2; ++k) + func[k] = 0; + + func[UP] = (internal::rank2_update_selector::run); + func[LO] = (internal::rank2_update_selector::run); + + init = true; + } Scalar* x = reinterpret_cast(px); Scalar* y = reinterpret_cast(py); @@ -141,10 +162,12 @@ int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px Scalar* x_cpy = get_compact_vector(x,*n,*incx); Scalar* y_cpy = get_compact_vector(y,*n,*incy); + + int code = UPLO(*uplo); + if(code>=2 || func[code]==0) + return 0; - // TODO perform direct calls to underlying implementation - if(UPLO(*uplo)==LO) matrix(c,*n,*n,*ldc).selfadjointView().rankUpdate(vector(x_cpy,*n), vector(y_cpy,*n), alpha); - else if(UPLO(*uplo)==UP) matrix(c,*n,*n,*ldc).selfadjointView().rankUpdate(vector(x_cpy,*n), vector(y_cpy,*n), alpha); + func[code](*n, c, *ldc, x_cpy, y_cpy, alpha); if(x_cpy!=x) delete[] x_cpy; if(y_cpy!=y) delete[] y_cpy; @@ -208,3 +231,41 @@ int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px // return 1; // } +/** DGER performs the rank 1 operation + * + * A := alpha*x*y' + A, + * + * where alpha is a scalar, x is an m element vector, y is an n element + * vector and A is an m by n matrix. + */ +int EIGEN_BLAS_FUNC(ger)(int *m, int *n, Scalar *palpha, Scalar *px, int *incx, Scalar *py, int *incy, Scalar *pa, int *lda) +{ + Scalar* x = reinterpret_cast(px); + Scalar* y = reinterpret_cast(py); + Scalar* a = reinterpret_cast(pa); + Scalar alpha = *reinterpret_cast(palpha); + + int info = 0; + if(*m<0) info = 1; + else if(*n<0) info = 2; + else if(*incx==0) info = 5; + else if(*incy==0) info = 7; + else if(*lda::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha); + + if(x_cpy!=x) delete[] x_cpy; + if(y_cpy!=y) delete[] y_cpy; + + return 1; +} + + diff --git a/blas/level3_impl.h b/blas/level3_impl.h index 2371f25c3..84c9f4f2b 100644 --- a/blas/level3_impl.h +++ b/blas/level3_impl.h @@ -305,6 +305,7 @@ int EIGEN_BLAS_FUNC(symm)(char *side, char *uplo, int *m, int *n, RealScalar *pa int EIGEN_BLAS_FUNC(syrk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc) { // std::cerr << "in syrk " << *uplo << " " << *op << " " << *n << " " << *k << " " << *palpha << " " << *lda << " " << *pbeta << " " << *ldc << "\n"; + #if !ISCOMPLEX typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar); static functype func[8]; @@ -324,6 +325,7 @@ int EIGEN_BLAS_FUNC(syrk)(char *uplo, char *op, int *n, int *k, RealScalar *palp init = true; } + #endif Scalar* a = reinterpret_cast(pa); Scalar* c = reinterpret_cast(pc); diff --git a/blas/single.cpp b/blas/single.cpp index 1b7775aed..836e3eee2 100644 --- a/blas/single.cpp +++ b/blas/single.cpp @@ -17,3 +17,6 @@ #include "level2_impl.h" #include "level2_real_impl.h" #include "level3_impl.h" + +float BLASFUNC(sdsdot)(int* n, float* alpha, float* x, int* incx, float* y, int* incy) +{ return *alpha + BLASFUNC(dsdot)(n, x, incx, y, incy); } diff --git a/blas/testing/dblat1.f b/blas/testing/dblat1.f index 5a45d69f4..30691f9bf 100644 --- a/blas/testing/dblat1.f +++ b/blas/testing/dblat1.f @@ -1,12 +1,54 @@ +*> \brief \b DBLAT1 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* PROGRAM DBLAT1 +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> Test program for the DOUBLE PRECISION Level 1 BLAS. +*> +*> Based upon the original BLAS test routine together with: +*> F06EAF Example Program Text +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date April 2012 +* +*> \ingroup double_blas_testing +* +* ===================================================================== PROGRAM DBLAT1 -* Test program for the DOUBLE PRECISION Level 1 BLAS. -* Based upon the original BLAS test routine together with: -* F06EAF Example Program Text +* +* -- Reference BLAS test routine (version 3.4.1) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* April 2012 +* +* ===================================================================== +* * .. Parameters .. INTEGER NOUT PARAMETER (NOUT=6) * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. DOUBLE PRECISION SFAC @@ -14,31 +56,30 @@ * .. External Subroutines .. EXTERNAL CHECK0, CHECK1, CHECK2, CHECK3, HEADER * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. DATA SFAC/9.765625D-4/ * .. Executable Statements .. WRITE (NOUT,99999) - DO 20 IC = 1, 10 + DO 20 IC = 1, 13 ICASE = IC CALL HEADER * -* .. Initialize PASS, INCX, INCY, and MODE for a new case. .. -* .. the value 9999 for INCX, INCY or MODE will appear in the .. +* .. Initialize PASS, INCX, and INCY for a new case. .. +* .. the value 9999 for INCX or INCY will appear in the .. * .. detailed output, if any, for cases that do not involve .. * .. these parameters .. * PASS = .TRUE. INCX = 9999 INCY = 9999 - MODE = 9999 - IF (ICASE.EQ.3) THEN + IF (ICASE.EQ.3 .OR. ICASE.EQ.11) THEN CALL CHECK0(SFAC) ELSE IF (ICASE.EQ.7 .OR. ICASE.EQ.8 .OR. ICASE.EQ.9 .OR. + ICASE.EQ.10) THEN CALL CHECK1(SFAC) ELSE IF (ICASE.EQ.1 .OR. ICASE.EQ.2 .OR. ICASE.EQ.5 .OR. - + ICASE.EQ.6) THEN + + ICASE.EQ.6 .OR. ICASE.EQ.12 .OR. ICASE.EQ.13) THEN CALL CHECK2(SFAC) ELSE IF (ICASE.EQ.4) THEN CALL CHECK3(SFAC) @@ -56,12 +97,12 @@ INTEGER NOUT PARAMETER (NOUT=6) * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Arrays .. - CHARACTER*6 L(10) + CHARACTER*6 L(13) * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. DATA L(1)/' DDOT '/ DATA L(2)/'DAXPY '/ @@ -73,6 +114,9 @@ DATA L(8)/'DASUM '/ DATA L(9)/'DSCAL '/ DATA L(10)/'IDAMAX'/ + DATA L(11)/'DROTMG'/ + DATA L(12)/'DROTM '/ + DATA L(13)/'DSDOT '/ * .. Executable Statements .. WRITE (NOUT,99999) ICASE, L(ICASE) RETURN @@ -86,18 +130,18 @@ * .. Scalar Arguments .. DOUBLE PRECISION SFAC * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. - DOUBLE PRECISION D12, SA, SB, SC, SS - INTEGER K + DOUBLE PRECISION SA, SB, SC, SS, D12 + INTEGER I, K * .. Local Arrays .. DOUBLE PRECISION DA1(8), DATRUE(8), DB1(8), DBTRUE(8), DC1(8), - + DS1(8) + $ DS1(8), DAB(4,9), DTEMP(9), DTRUE(9,9) * .. External Subroutines .. - EXTERNAL DROTG, STEST1 + EXTERNAL DROTG, DROTMG, STEST1 * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. DATA DA1/0.3D0, 0.4D0, -0.3D0, -0.4D0, -0.3D0, 0.0D0, + 0.0D0, 1.0D0/ @@ -111,7 +155,52 @@ + 0.0D0, 1.0D0, 1.0D0/ DATA DBTRUE/0.0D0, 0.6D0, 0.0D0, -0.6D0, 0.0D0, + 0.0D0, 1.0D0, 0.0D0/ - DATA D12/4096.0D0/ +* INPUT FOR MODIFIED GIVENS + DATA DAB/ .1D0,.3D0,1.2D0,.2D0, + A .7D0, .2D0, .6D0, 4.2D0, + B 0.D0,0.D0,0.D0,0.D0, + C 4.D0, -1.D0, 2.D0, 4.D0, + D 6.D-10, 2.D-2, 1.D5, 10.D0, + E 4.D10, 2.D-2, 1.D-5, 10.D0, + F 2.D-10, 4.D-2, 1.D5, 10.D0, + G 2.D10, 4.D-2, 1.D-5, 10.D0, + H 4.D0, -2.D0, 8.D0, 4.D0 / +* TRUE RESULTS FOR MODIFIED GIVENS + DATA DTRUE/0.D0,0.D0, 1.3D0, .2D0, 0.D0,0.D0,0.D0, .5D0, 0.D0, + A 0.D0,0.D0, 4.5D0, 4.2D0, 1.D0, .5D0, 0.D0,0.D0,0.D0, + B 0.D0,0.D0,0.D0,0.D0, -2.D0, 0.D0,0.D0,0.D0,0.D0, + C 0.D0,0.D0,0.D0, 4.D0, -1.D0, 0.D0,0.D0,0.D0,0.D0, + D 0.D0, 15.D-3, 0.D0, 10.D0, -1.D0, 0.D0, -1.D-4, + E 0.D0, 1.D0, + F 0.D0,0.D0, 6144.D-5, 10.D0, -1.D0, 4096.D0, -1.D6, + G 0.D0, 1.D0, + H 0.D0,0.D0,15.D0,10.D0,-1.D0, 5.D-5, 0.D0,1.D0,0.D0, + I 0.D0,0.D0, 15.D0, 10.D0, -1. D0, 5.D5, -4096.D0, + J 1.D0, 4096.D-6, + K 0.D0,0.D0, 7.D0, 4.D0, 0.D0,0.D0, -.5D0, -.25D0, 0.D0/ +* 4096 = 2 ** 12 + DATA D12 /4096.D0/ + DTRUE(1,1) = 12.D0 / 130.D0 + DTRUE(2,1) = 36.D0 / 130.D0 + DTRUE(7,1) = -1.D0 / 6.D0 + DTRUE(1,2) = 14.D0 / 75.D0 + DTRUE(2,2) = 49.D0 / 75.D0 + DTRUE(9,2) = 1.D0 / 7.D0 + DTRUE(1,5) = 45.D-11 * (D12 * D12) + DTRUE(3,5) = 4.D5 / (3.D0 * D12) + DTRUE(6,5) = 1.D0 / D12 + DTRUE(8,5) = 1.D4 / (3.D0 * D12) + DTRUE(1,6) = 4.D10 / (1.5D0 * D12 * D12) + DTRUE(2,6) = 2.D-2 / 1.5D0 + DTRUE(8,6) = 5.D-7 * D12 + DTRUE(1,7) = 4.D0 / 150.D0 + DTRUE(2,7) = (2.D-10 / 1.5D0) * (D12 * D12) + DTRUE(7,7) = -DTRUE(6,5) + DTRUE(9,7) = 1.D4 / D12 + DTRUE(1,8) = DTRUE(1,7) + DTRUE(2,8) = 2.D10 / (1.5D0 * D12 * D12) + DTRUE(1,9) = 32.D0 / 7.D0 + DTRUE(2,9) = -16.D0 / 7.D0 * .. Executable Statements .. * * Compute true values which cannot be prestored @@ -134,6 +223,15 @@ CALL STEST1(SB,DBTRUE(K),DBTRUE(K),SFAC) CALL STEST1(SC,DC1(K),DC1(K),SFAC) CALL STEST1(SS,DS1(K),DS1(K),SFAC) + ELSEIF (ICASE.EQ.11) THEN +* .. DROTMG .. + DO I=1,4 + DTEMP(I)= DAB(I,K) + DTEMP(I+4) = 0.0 + END DO + DTEMP(9) = 0.0 + CALL DROTMG(DTEMP(1),DTEMP(2),DTEMP(3),DTEMP(4),DTEMP(5)) + CALL STEST(9,DTEMP,DTRUE(1,K),DTRUE(1,K),SFAC) ELSE WRITE (NOUT,*) ' Shouldn''t be here in CHECK0' STOP @@ -148,7 +246,7 @@ * .. Scalar Arguments .. DOUBLE PRECISION SFAC * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. INTEGER I, LEN, NP1 @@ -165,7 +263,7 @@ * .. Intrinsic Functions .. INTRINSIC MAX * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. DATA SA/0.3D0, -1.0D0, 0.0D0, 1.0D0, 0.3D0, 0.3D0, + 0.3D0, 0.3D0, 0.3D0, 0.3D0/ @@ -212,11 +310,11 @@ IF (ICASE.EQ.7) THEN * .. DNRM2 .. STEMP(1) = DTRUE1(NP1) - CALL STEST1(DNRM2(N,SX,INCX),STEMP,STEMP,SFAC) + CALL STEST1(DNRM2(N,SX,INCX),STEMP(1),STEMP,SFAC) ELSE IF (ICASE.EQ.8) THEN * .. DASUM .. STEMP(1) = DTRUE3(NP1) - CALL STEST1(DASUM(N,SX,INCX),STEMP,STEMP,SFAC) + CALL STEST1(DASUM(N,SX,INCX),STEMP(1),STEMP,SFAC) ELSE IF (ICASE.EQ.9) THEN * .. DSCAL .. CALL DSCAL(N,SA((INCX-1)*5+NP1),SX,INCX) @@ -242,27 +340,39 @@ * .. Scalar Arguments .. DOUBLE PRECISION SFAC * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. - DOUBLE PRECISION SA, SC, SS - INTEGER I, J, KI, KN, KSIZE, LENX, LENY, MX, MY + DOUBLE PRECISION SA + INTEGER I, J, KI, KN, KNI, KPAR, KSIZE, LENX, LENY, + $ MX, MY * .. Local Arrays .. DOUBLE PRECISION DT10X(7,4,4), DT10Y(7,4,4), DT7(4,4), - + DT8(7,4,4), DT9X(7,4,4), DT9Y(7,4,4), DX1(7), - + DY1(7), SSIZE1(4), SSIZE2(14,2), STX(7), STY(7), - + SX(7), SY(7) + $ DT8(7,4,4), DX1(7), + $ DY1(7), SSIZE1(4), SSIZE2(14,2), SSIZE(7), + $ STX(7), STY(7), SX(7), SY(7), + $ DPAR(5,4), DT19X(7,4,16),DT19XA(7,4,4), + $ DT19XB(7,4,4), DT19XC(7,4,4),DT19XD(7,4,4), + $ DT19Y(7,4,16), DT19YA(7,4,4),DT19YB(7,4,4), + $ DT19YC(7,4,4), DT19YD(7,4,4), DTEMP(5) INTEGER INCXS(4), INCYS(4), LENS(4,2), NS(4) * .. External Functions .. - DOUBLE PRECISION DDOT - EXTERNAL DDOT + DOUBLE PRECISION DDOT, DSDOT + EXTERNAL DDOT, DSDOT * .. External Subroutines .. - EXTERNAL DAXPY, DCOPY, DSWAP, STEST, STEST1 + EXTERNAL DAXPY, DCOPY, DROTM, DSWAP, STEST, STEST1 * .. Intrinsic Functions .. INTRINSIC ABS, MIN * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. + EQUIVALENCE (DT19X(1,1,1),DT19XA(1,1,1)),(DT19X(1,1,5), + A DT19XB(1,1,1)),(DT19X(1,1,9),DT19XC(1,1,1)), + B (DT19X(1,1,13),DT19XD(1,1,1)) + EQUIVALENCE (DT19Y(1,1,1),DT19YA(1,1,1)),(DT19Y(1,1,5), + A DT19YB(1,1,1)),(DT19Y(1,1,9),DT19YC(1,1,1)), + B (DT19Y(1,1,13),DT19YD(1,1,1)) + DATA SA/0.3D0/ DATA INCXS/1, 2, -2, -1/ DATA INCYS/1, -2, 1, -2/ @@ -272,7 +382,6 @@ + -0.4D0/ DATA DY1/0.5D0, -0.9D0, 0.3D0, 0.7D0, -0.6D0, 0.2D0, + 0.8D0/ - DATA SC, SS/0.8D0, 0.6D0/ DATA DT7/0.0D0, 0.30D0, 0.21D0, 0.62D0, 0.0D0, + 0.30D0, -0.07D0, 0.85D0, 0.0D0, 0.30D0, -0.79D0, + -0.74D0, 0.0D0, 0.30D0, 0.33D0, 1.27D0/ @@ -295,44 +404,6 @@ + 0.0D0, 0.68D0, -0.9D0, 0.33D0, 0.0D0, 0.0D0, + 0.0D0, 0.0D0, 0.68D0, -0.9D0, 0.33D0, 0.7D0, + -0.75D0, 0.2D0, 1.04D0/ - DATA DT9X/0.6D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, - + 0.0D0, 0.78D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, - + 0.0D0, 0.0D0, 0.78D0, -0.46D0, 0.0D0, 0.0D0, - + 0.0D0, 0.0D0, 0.0D0, 0.78D0, -0.46D0, -0.22D0, - + 1.06D0, 0.0D0, 0.0D0, 0.0D0, 0.6D0, 0.0D0, - + 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.78D0, - + 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, - + 0.66D0, 0.1D0, -0.1D0, 0.0D0, 0.0D0, 0.0D0, - + 0.0D0, 0.96D0, 0.1D0, -0.76D0, 0.8D0, 0.90D0, - + -0.3D0, -0.02D0, 0.6D0, 0.0D0, 0.0D0, 0.0D0, - + 0.0D0, 0.0D0, 0.0D0, 0.78D0, 0.0D0, 0.0D0, - + 0.0D0, 0.0D0, 0.0D0, 0.0D0, -0.06D0, 0.1D0, - + -0.1D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.90D0, - + 0.1D0, -0.22D0, 0.8D0, 0.18D0, -0.3D0, -0.02D0, - + 0.6D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, - + 0.78D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, - + 0.0D0, 0.78D0, 0.26D0, 0.0D0, 0.0D0, 0.0D0, - + 0.0D0, 0.0D0, 0.78D0, 0.26D0, -0.76D0, 1.12D0, - + 0.0D0, 0.0D0, 0.0D0/ - DATA DT9Y/0.5D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, - + 0.0D0, 0.04D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, - + 0.0D0, 0.0D0, 0.04D0, -0.78D0, 0.0D0, 0.0D0, - + 0.0D0, 0.0D0, 0.0D0, 0.04D0, -0.78D0, 0.54D0, - + 0.08D0, 0.0D0, 0.0D0, 0.0D0, 0.5D0, 0.0D0, - + 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.04D0, - + 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.7D0, - + -0.9D0, -0.12D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, - + 0.64D0, -0.9D0, -0.30D0, 0.7D0, -0.18D0, 0.2D0, - + 0.28D0, 0.5D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, - + 0.0D0, 0.0D0, 0.04D0, 0.0D0, 0.0D0, 0.0D0, - + 0.0D0, 0.0D0, 0.0D0, 0.7D0, -1.08D0, 0.0D0, - + 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.64D0, -1.26D0, - + 0.54D0, 0.20D0, 0.0D0, 0.0D0, 0.0D0, 0.5D0, - + 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, - + 0.04D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, - + 0.0D0, 0.04D0, -0.9D0, 0.18D0, 0.0D0, 0.0D0, - + 0.0D0, 0.0D0, 0.04D0, -0.9D0, 0.18D0, 0.7D0, - + -0.18D0, 0.2D0, 0.16D0/ DATA DT10X/0.6D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, + 0.0D0, 0.5D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, + 0.0D0, 0.5D0, -0.9D0, 0.0D0, 0.0D0, 0.0D0, @@ -375,6 +446,150 @@ + 0.0D0, 1.17D0, 1.17D0, 1.17D0, 1.17D0, 1.17D0, + 1.17D0, 1.17D0, 1.17D0, 1.17D0, 1.17D0, 1.17D0, + 1.17D0, 1.17D0, 1.17D0/ +* +* FOR DROTM +* + DATA DPAR/-2.D0, 0.D0,0.D0,0.D0,0.D0, + A -1.D0, 2.D0, -3.D0, -4.D0, 5.D0, + B 0.D0, 0.D0, 2.D0, -3.D0, 0.D0, + C 1.D0, 5.D0, 2.D0, 0.D0, -4.D0/ +* TRUE X RESULTS F0R ROTATIONS DROTM + DATA DT19XA/.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + A .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + B .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + C .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + D .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + E -.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + F -.9D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + G 3.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + H .6D0, .1D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + I -.8D0, 3.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + J -.9D0, 2.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + K 3.5D0, -.4D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + L .6D0, .1D0, -.5D0, .8D0, 0.D0,0.D0,0.D0, + M -.8D0, 3.8D0, -2.2D0, -1.2D0, 0.D0,0.D0,0.D0, + N -.9D0, 2.8D0, -1.4D0, -1.3D0, 0.D0,0.D0,0.D0, + O 3.5D0, -.4D0, -2.2D0, 4.7D0, 0.D0,0.D0,0.D0/ +* + DATA DT19XB/.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + A .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + B .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + C .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + D .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + E -.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + F -.9D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + G 3.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + H .6D0, .1D0, -.5D0, 0.D0,0.D0,0.D0,0.D0, + I 0.D0, .1D0, -3.0D0, 0.D0,0.D0,0.D0,0.D0, + J -.3D0, .1D0, -2.0D0, 0.D0,0.D0,0.D0,0.D0, + K 3.3D0, .1D0, -2.0D0, 0.D0,0.D0,0.D0,0.D0, + L .6D0, .1D0, -.5D0, .8D0, .9D0, -.3D0, -.4D0, + M -2.0D0, .1D0, 1.4D0, .8D0, .6D0, -.3D0, -2.8D0, + N -1.8D0, .1D0, 1.3D0, .8D0, 0.D0, -.3D0, -1.9D0, + O 3.8D0, .1D0, -3.1D0, .8D0, 4.8D0, -.3D0, -1.5D0 / +* + DATA DT19XC/.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + A .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + B .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + C .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + D .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + E -.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + F -.9D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + G 3.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + H .6D0, .1D0, -.5D0, 0.D0,0.D0,0.D0,0.D0, + I 4.8D0, .1D0, -3.0D0, 0.D0,0.D0,0.D0,0.D0, + J 3.3D0, .1D0, -2.0D0, 0.D0,0.D0,0.D0,0.D0, + K 2.1D0, .1D0, -2.0D0, 0.D0,0.D0,0.D0,0.D0, + L .6D0, .1D0, -.5D0, .8D0, .9D0, -.3D0, -.4D0, + M -1.6D0, .1D0, -2.2D0, .8D0, 5.4D0, -.3D0, -2.8D0, + N -1.5D0, .1D0, -1.4D0, .8D0, 3.6D0, -.3D0, -1.9D0, + O 3.7D0, .1D0, -2.2D0, .8D0, 3.6D0, -.3D0, -1.5D0 / +* + DATA DT19XD/.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + A .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + B .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + C .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + D .6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + E -.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + F -.9D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + G 3.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + H .6D0, .1D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + I -.8D0, -1.0D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + J -.9D0, -.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + K 3.5D0, .8D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + L .6D0, .1D0, -.5D0, .8D0, 0.D0,0.D0,0.D0, + M -.8D0, -1.0D0, 1.4D0, -1.6D0, 0.D0,0.D0,0.D0, + N -.9D0, -.8D0, 1.3D0, -1.6D0, 0.D0,0.D0,0.D0, + O 3.5D0, .8D0, -3.1D0, 4.8D0, 0.D0,0.D0,0.D0/ +* TRUE Y RESULTS FOR ROTATIONS DROTM + DATA DT19YA/.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + A .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + B .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + C .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + D .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + E .7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + F 1.7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + G -2.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + H .5D0, -.9D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + I .7D0, -4.8D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + J 1.7D0, -.7D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + K -2.6D0, 3.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + L .5D0, -.9D0, .3D0, .7D0, 0.D0,0.D0,0.D0, + M .7D0, -4.8D0, 3.0D0, 1.1D0, 0.D0,0.D0,0.D0, + N 1.7D0, -.7D0, -.7D0, 2.3D0, 0.D0,0.D0,0.D0, + O -2.6D0, 3.5D0, -.7D0, -3.6D0, 0.D0,0.D0,0.D0/ +* + DATA DT19YB/.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + A .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + B .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + C .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + D .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + E .7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + F 1.7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + G -2.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + H .5D0, -.9D0, .3D0, 0.D0,0.D0,0.D0,0.D0, + I 4.0D0, -.9D0, -.3D0, 0.D0,0.D0,0.D0,0.D0, + J -.5D0, -.9D0, 1.5D0, 0.D0,0.D0,0.D0,0.D0, + K -1.5D0, -.9D0, -1.8D0, 0.D0,0.D0,0.D0,0.D0, + L .5D0, -.9D0, .3D0, .7D0, -.6D0, .2D0, .8D0, + M 3.7D0, -.9D0, -1.2D0, .7D0, -1.5D0, .2D0, 2.2D0, + N -.3D0, -.9D0, 2.1D0, .7D0, -1.6D0, .2D0, 2.0D0, + O -1.6D0, -.9D0, -2.1D0, .7D0, 2.9D0, .2D0, -3.8D0 / +* + DATA DT19YC/.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + A .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + B .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + C .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + D .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + E .7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + F 1.7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + G -2.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + H .5D0, -.9D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + I 4.0D0, -6.3D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + J -.5D0, .3D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + K -1.5D0, 3.0D0, 0.D0,0.D0,0.D0,0.D0,0.D0, + L .5D0, -.9D0, .3D0, .7D0, 0.D0,0.D0,0.D0, + M 3.7D0, -7.2D0, 3.0D0, 1.7D0, 0.D0,0.D0,0.D0, + N -.3D0, .9D0, -.7D0, 1.9D0, 0.D0,0.D0,0.D0, + O -1.6D0, 2.7D0, -.7D0, -3.4D0, 0.D0,0.D0,0.D0/ +* + DATA DT19YD/.5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + A .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + B .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + C .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + D .5D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + E .7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + F 1.7D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + G -2.6D0, 0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, + H .5D0, -.9D0, .3D0, 0.D0,0.D0,0.D0,0.D0, + I .7D0, -.9D0, 1.2D0, 0.D0,0.D0,0.D0,0.D0, + J 1.7D0, -.9D0, .5D0, 0.D0,0.D0,0.D0,0.D0, + K -2.6D0, -.9D0, -1.3D0, 0.D0,0.D0,0.D0,0.D0, + L .5D0, -.9D0, .3D0, .7D0, -.6D0, .2D0, .8D0, + M .7D0, -.9D0, 1.2D0, .7D0, -1.5D0, .2D0, 1.6D0, + N 1.7D0, -.9D0, .5D0, .7D0, -1.6D0, .2D0, 2.4D0, + O -2.6D0, -.9D0, -1.3D0, .7D0, 2.9D0, .2D0, -4.0D0 / +* * .. Executable Statements .. * DO 120 KI = 1, 4 @@ -421,6 +636,39 @@ 80 CONTINUE CALL STEST(LENX,SX,STX,SSIZE2(1,1),1.0D0) CALL STEST(LENY,SY,STY,SSIZE2(1,1),1.0D0) + ELSE IF (ICASE.EQ.12) THEN +* .. DROTM .. + KNI=KN+4*(KI-1) + DO KPAR=1,4 + DO I=1,7 + SX(I) = DX1(I) + SY(I) = DY1(I) + STX(I)= DT19X(I,KPAR,KNI) + STY(I)= DT19Y(I,KPAR,KNI) + END DO +* + DO I=1,5 + DTEMP(I) = DPAR(I,KPAR) + END DO +* + DO I=1,LENX + SSIZE(I)=STX(I) + END DO +* SEE REMARK ABOVE ABOUT DT11X(1,2,7) +* AND DT11X(5,3,8). + IF ((KPAR .EQ. 2) .AND. (KNI .EQ. 7)) + $ SSIZE(1) = 2.4D0 + IF ((KPAR .EQ. 3) .AND. (KNI .EQ. 8)) + $ SSIZE(5) = 1.8D0 +* + CALL DROTM(N,SX,INCX,SY,INCY,DTEMP) + CALL STEST(LENX,SX,STX,SSIZE,SFAC) + CALL STEST(LENY,SY,STY,STY,SFAC) + END DO + ELSE IF (ICASE.EQ.13) THEN +* .. DSDOT .. + CALL TESTDSDOT(REAL(DSDOT(N,REAL(SX),INCX,REAL(SY),INCY)), + $ REAL(DT7(KN,KI)),REAL(SSIZE1(KN)), .3125E-1) ELSE WRITE (NOUT,*) ' Shouldn''t be here in CHECK2' STOP @@ -436,10 +684,10 @@ * .. Scalar Arguments .. DOUBLE PRECISION SFAC * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. - DOUBLE PRECISION SA, SC, SS + DOUBLE PRECISION SC, SS INTEGER I, K, KI, KN, KSIZE, LENX, LENY, MX, MY * .. Local Arrays .. DOUBLE PRECISION COPYX(5), COPYY(5), DT9X(7,4,4), DT9Y(7,4,4), @@ -454,9 +702,8 @@ * .. Intrinsic Functions .. INTRINSIC ABS, MIN * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. - DATA SA/0.3D0/ DATA INCXS/1, 2, -2, -1/ DATA INCYS/1, -2, 1, -2/ DATA LENS/1, 1, 2, 4, 1, 1, 3, 7/ @@ -647,14 +894,15 @@ * * .. Parameters .. INTEGER NOUT - PARAMETER (NOUT=6) + DOUBLE PRECISION ZERO + PARAMETER (NOUT=6, ZERO=0.0D0) * .. Scalar Arguments .. DOUBLE PRECISION SFAC INTEGER LEN * .. Array Arguments .. DOUBLE PRECISION SCOMP(LEN), SSIZE(LEN), STRUE(LEN) * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. DOUBLE PRECISION SD @@ -665,12 +913,12 @@ * .. Intrinsic Functions .. INTRINSIC ABS * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Executable Statements .. * DO 40 I = 1, LEN SD = SCOMP(I) - STRUE(I) - IF (SDIFF(ABS(SSIZE(I))+ABS(SFAC*SD),ABS(SSIZE(I))).EQ.0.0D0) + IF (ABS(SFAC*SD) .LE. ABS(SSIZE(I))*EPSILON(ZERO)) + GO TO 40 * * HERE SCOMP(I) IS NOT CLOSE TO STRUE(I). @@ -680,16 +928,64 @@ PASS = .FALSE. WRITE (NOUT,99999) WRITE (NOUT,99998) - 20 WRITE (NOUT,99997) ICASE, N, INCX, INCY, MODE, I, SCOMP(I), + 20 WRITE (NOUT,99997) ICASE, N, INCX, INCY, I, SCOMP(I), + STRUE(I), SD, SSIZE(I) 40 CONTINUE RETURN * 99999 FORMAT (' FAIL') -99998 FORMAT (/' CASE N INCX INCY MODE I ', +99998 FORMAT (/' CASE N INCX INCY I ', + ' COMP(I) TRUE(I) DIFFERENCE', + ' SIZE(I)',/1X) -99997 FORMAT (1X,I4,I3,3I5,I3,2D36.8,2D12.4) +99997 FORMAT (1X,I4,I3,2I5,I3,2D36.8,2D12.4) + END + SUBROUTINE TESTDSDOT(SCOMP,STRUE,SSIZE,SFAC) +* ********************************* STEST ************************** +* +* THIS SUBR COMPARES ARRAYS SCOMP() AND STRUE() OF LENGTH LEN TO +* SEE IF THE TERM BY TERM DIFFERENCES, MULTIPLIED BY SFAC, ARE +* NEGLIGIBLE. +* +* C. L. LAWSON, JPL, 1974 DEC 10 +* +* .. Parameters .. + INTEGER NOUT + REAL ZERO + PARAMETER (NOUT=6, ZERO=0.0E0) +* .. Scalar Arguments .. + REAL SFAC, SCOMP, SSIZE, STRUE +* .. Scalars in Common .. + INTEGER ICASE, INCX, INCY, N + LOGICAL PASS +* .. Local Scalars .. + REAL SD +* .. Intrinsic Functions .. + INTRINSIC ABS +* .. Common blocks .. + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS +* .. Executable Statements .. +* + SD = SCOMP - STRUE + IF (ABS(SFAC*SD) .LE. ABS(SSIZE) * EPSILON(ZERO)) + + GO TO 40 +* +* HERE SCOMP(I) IS NOT CLOSE TO STRUE(I). +* + IF ( .NOT. PASS) GO TO 20 +* PRINT FAIL MESSAGE AND HEADER. + PASS = .FALSE. + WRITE (NOUT,99999) + WRITE (NOUT,99998) + 20 WRITE (NOUT,99997) ICASE, N, INCX, INCY, SCOMP, + + STRUE, SD, SSIZE + 40 CONTINUE + RETURN +* +99999 FORMAT (' FAIL') +99998 FORMAT (/' CASE N INCX INCY ', + + ' COMP(I) TRUE(I) DIFFERENCE', + + ' SIZE(I)',/1X) +99997 FORMAT (1X,I4,I3,1I5,I3,2E36.8,2E12.4) END SUBROUTINE STEST1(SCOMP1,STRUE1,SSIZE,SFAC) * ************************* STEST1 ***************************** @@ -739,12 +1035,12 @@ * .. Scalar Arguments .. INTEGER ICOMP, ITRUE * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. INTEGER ID * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Executable Statements .. * IF (ICOMP.EQ.ITRUE) GO TO 40 @@ -757,13 +1053,13 @@ WRITE (NOUT,99999) WRITE (NOUT,99998) 20 ID = ICOMP - ITRUE - WRITE (NOUT,99997) ICASE, N, INCX, INCY, MODE, ICOMP, ITRUE, ID + WRITE (NOUT,99997) ICASE, N, INCX, INCY, ICOMP, ITRUE, ID 40 CONTINUE RETURN * 99999 FORMAT (' FAIL') -99998 FORMAT (/' CASE N INCX INCY MODE ', +99998 FORMAT (/' CASE N INCX INCY ', + ' COMP TRUE DIFFERENCE', + /1X) -99997 FORMAT (1X,I4,I3,3I5,2I36,I12) +99997 FORMAT (1X,I4,I3,2I5,2I36,I12) END diff --git a/blas/testing/sblat1.f b/blas/testing/sblat1.f index a982d1852..6657c2693 100644 --- a/blas/testing/sblat1.f +++ b/blas/testing/sblat1.f @@ -1,12 +1,54 @@ +*> \brief \b SBLAT1 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* PROGRAM SBLAT1 +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> Test program for the REAL Level 1 BLAS. +*> +*> Based upon the original BLAS test routine together with: +*> F06EAF Example Program Text +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date April 2012 +* +*> \ingroup single_blas_testing +* +* ===================================================================== PROGRAM SBLAT1 -* Test program for the REAL Level 1 BLAS. -* Based upon the original BLAS test routine together with: -* F06EAF Example Program Text +* +* -- Reference BLAS test routine (version 3.4.1) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* April 2012 +* +* ===================================================================== +* * .. Parameters .. INTEGER NOUT PARAMETER (NOUT=6) * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. REAL SFAC @@ -14,31 +56,30 @@ * .. External Subroutines .. EXTERNAL CHECK0, CHECK1, CHECK2, CHECK3, HEADER * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. DATA SFAC/9.765625E-4/ * .. Executable Statements .. WRITE (NOUT,99999) - DO 20 IC = 1, 10 + DO 20 IC = 1, 13 ICASE = IC CALL HEADER * -* .. Initialize PASS, INCX, INCY, and MODE for a new case. .. -* .. the value 9999 for INCX, INCY or MODE will appear in the .. +* .. Initialize PASS, INCX, and INCY for a new case. .. +* .. the value 9999 for INCX or INCY will appear in the .. * .. detailed output, if any, for cases that do not involve .. * .. these parameters .. * PASS = .TRUE. INCX = 9999 INCY = 9999 - MODE = 9999 - IF (ICASE.EQ.3) THEN + IF (ICASE.EQ.3 .OR. ICASE.EQ.11) THEN CALL CHECK0(SFAC) ELSE IF (ICASE.EQ.7 .OR. ICASE.EQ.8 .OR. ICASE.EQ.9 .OR. + ICASE.EQ.10) THEN CALL CHECK1(SFAC) ELSE IF (ICASE.EQ.1 .OR. ICASE.EQ.2 .OR. ICASE.EQ.5 .OR. - + ICASE.EQ.6) THEN + + ICASE.EQ.6 .OR. ICASE.EQ.12 .OR. ICASE.EQ.13) THEN CALL CHECK2(SFAC) ELSE IF (ICASE.EQ.4) THEN CALL CHECK3(SFAC) @@ -56,12 +97,12 @@ INTEGER NOUT PARAMETER (NOUT=6) * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Arrays .. - CHARACTER*6 L(10) + CHARACTER*6 L(13) * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. DATA L(1)/' SDOT '/ DATA L(2)/'SAXPY '/ @@ -73,6 +114,9 @@ DATA L(8)/'SASUM '/ DATA L(9)/'SSCAL '/ DATA L(10)/'ISAMAX'/ + DATA L(11)/'SROTMG'/ + DATA L(12)/'SROTM '/ + DATA L(13)/'SDSDOT'/ * .. Executable Statements .. WRITE (NOUT,99999) ICASE, L(ICASE) RETURN @@ -86,18 +130,18 @@ * .. Scalar Arguments .. REAL SFAC * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. REAL D12, SA, SB, SC, SS - INTEGER K + INTEGER I, K * .. Local Arrays .. REAL DA1(8), DATRUE(8), DB1(8), DBTRUE(8), DC1(8), - + DS1(8) + + DS1(8), DAB(4,9), DTEMP(9), DTRUE(9,9) * .. External Subroutines .. - EXTERNAL SROTG, STEST1 + EXTERNAL SROTG, SROTMG, STEST1 * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. DATA DA1/0.3E0, 0.4E0, -0.3E0, -0.4E0, -0.3E0, 0.0E0, + 0.0E0, 1.0E0/ @@ -111,7 +155,52 @@ + 0.0E0, 1.0E0, 1.0E0/ DATA DBTRUE/0.0E0, 0.6E0, 0.0E0, -0.6E0, 0.0E0, + 0.0E0, 1.0E0, 0.0E0/ - DATA D12/4096.0E0/ +* INPUT FOR MODIFIED GIVENS + DATA DAB/ .1E0,.3E0,1.2E0,.2E0, + A .7E0, .2E0, .6E0, 4.2E0, + B 0.E0,0.E0,0.E0,0.E0, + C 4.E0, -1.E0, 2.E0, 4.E0, + D 6.E-10, 2.E-2, 1.E5, 10.E0, + E 4.E10, 2.E-2, 1.E-5, 10.E0, + F 2.E-10, 4.E-2, 1.E5, 10.E0, + G 2.E10, 4.E-2, 1.E-5, 10.E0, + H 4.E0, -2.E0, 8.E0, 4.E0 / +* TRUE RESULTS FOR MODIFIED GIVENS + DATA DTRUE/0.E0,0.E0, 1.3E0, .2E0, 0.E0,0.E0,0.E0, .5E0, 0.E0, + A 0.E0,0.E0, 4.5E0, 4.2E0, 1.E0, .5E0, 0.E0,0.E0,0.E0, + B 0.E0,0.E0,0.E0,0.E0, -2.E0, 0.E0,0.E0,0.E0,0.E0, + C 0.E0,0.E0,0.E0, 4.E0, -1.E0, 0.E0,0.E0,0.E0,0.E0, + D 0.E0, 15.E-3, 0.E0, 10.E0, -1.E0, 0.E0, -1.E-4, + E 0.E0, 1.E0, + F 0.E0,0.E0, 6144.E-5, 10.E0, -1.E0, 4096.E0, -1.E6, + G 0.E0, 1.E0, + H 0.E0,0.E0,15.E0,10.E0,-1.E0, 5.E-5, 0.E0,1.E0,0.E0, + I 0.E0,0.E0, 15.E0, 10.E0, -1. E0, 5.E5, -4096.E0, + J 1.E0, 4096.E-6, + K 0.E0,0.E0, 7.E0, 4.E0, 0.E0,0.E0, -.5E0, -.25E0, 0.E0/ +* 4096 = 2 ** 12 + DATA D12 /4096.E0/ + DTRUE(1,1) = 12.E0 / 130.E0 + DTRUE(2,1) = 36.E0 / 130.E0 + DTRUE(7,1) = -1.E0 / 6.E0 + DTRUE(1,2) = 14.E0 / 75.E0 + DTRUE(2,2) = 49.E0 / 75.E0 + DTRUE(9,2) = 1.E0 / 7.E0 + DTRUE(1,5) = 45.E-11 * (D12 * D12) + DTRUE(3,5) = 4.E5 / (3.E0 * D12) + DTRUE(6,5) = 1.E0 / D12 + DTRUE(8,5) = 1.E4 / (3.E0 * D12) + DTRUE(1,6) = 4.E10 / (1.5E0 * D12 * D12) + DTRUE(2,6) = 2.E-2 / 1.5E0 + DTRUE(8,6) = 5.E-7 * D12 + DTRUE(1,7) = 4.E0 / 150.E0 + DTRUE(2,7) = (2.E-10 / 1.5E0) * (D12 * D12) + DTRUE(7,7) = -DTRUE(6,5) + DTRUE(9,7) = 1.E4 / D12 + DTRUE(1,8) = DTRUE(1,7) + DTRUE(2,8) = 2.E10 / (1.5E0 * D12 * D12) + DTRUE(1,9) = 32.E0 / 7.E0 + DTRUE(2,9) = -16.E0 / 7.E0 * .. Executable Statements .. * * Compute true values which cannot be prestored @@ -134,6 +223,15 @@ CALL STEST1(SB,DBTRUE(K),DBTRUE(K),SFAC) CALL STEST1(SC,DC1(K),DC1(K),SFAC) CALL STEST1(SS,DS1(K),DS1(K),SFAC) + ELSEIF (ICASE.EQ.11) THEN +* .. SROTMG .. + DO I=1,4 + DTEMP(I)= DAB(I,K) + DTEMP(I+4) = 0.0 + END DO + DTEMP(9) = 0.0 + CALL SROTMG(DTEMP(1),DTEMP(2),DTEMP(3),DTEMP(4),DTEMP(5)) + CALL STEST(9,DTEMP,DTRUE(1,K),DTRUE(1,K),SFAC) ELSE WRITE (NOUT,*) ' Shouldn''t be here in CHECK0' STOP @@ -148,7 +246,7 @@ * .. Scalar Arguments .. REAL SFAC * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. INTEGER I, LEN, NP1 @@ -165,7 +263,7 @@ * .. Intrinsic Functions .. INTRINSIC MAX * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. DATA SA/0.3E0, -1.0E0, 0.0E0, 1.0E0, 0.3E0, 0.3E0, + 0.3E0, 0.3E0, 0.3E0, 0.3E0/ @@ -212,11 +310,11 @@ IF (ICASE.EQ.7) THEN * .. SNRM2 .. STEMP(1) = DTRUE1(NP1) - CALL STEST1(SNRM2(N,SX,INCX),STEMP,STEMP,SFAC) + CALL STEST1(SNRM2(N,SX,INCX),STEMP(1),STEMP,SFAC) ELSE IF (ICASE.EQ.8) THEN * .. SASUM .. STEMP(1) = DTRUE3(NP1) - CALL STEST1(SASUM(N,SX,INCX),STEMP,STEMP,SFAC) + CALL STEST1(SASUM(N,SX,INCX),STEMP(1),STEMP,SFAC) ELSE IF (ICASE.EQ.9) THEN * .. SSCAL .. CALL SSCAL(N,SA((INCX-1)*5+NP1),SX,INCX) @@ -242,27 +340,40 @@ * .. Scalar Arguments .. REAL SFAC * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. - REAL SA, SC, SS - INTEGER I, J, KI, KN, KSIZE, LENX, LENY, MX, MY + REAL SA + INTEGER I, J, KI, KN, KNI, KPAR, KSIZE, LENX, LENY, + $ MX, MY * .. Local Arrays .. REAL DT10X(7,4,4), DT10Y(7,4,4), DT7(4,4), - + DT8(7,4,4), DT9X(7,4,4), DT9Y(7,4,4), DX1(7), - + DY1(7), SSIZE1(4), SSIZE2(14,2), STX(7), STY(7), - + SX(7), SY(7) + $ DT8(7,4,4), DX1(7), + $ DY1(7), SSIZE1(4), SSIZE2(14,2), SSIZE3(4), + $ SSIZE(7), STX(7), STY(7), SX(7), SY(7), + $ DPAR(5,4), DT19X(7,4,16),DT19XA(7,4,4), + $ DT19XB(7,4,4), DT19XC(7,4,4),DT19XD(7,4,4), + $ DT19Y(7,4,16), DT19YA(7,4,4),DT19YB(7,4,4), + $ DT19YC(7,4,4), DT19YD(7,4,4), DTEMP(5), + $ ST7B(4,4) INTEGER INCXS(4), INCYS(4), LENS(4,2), NS(4) * .. External Functions .. - REAL SDOT - EXTERNAL SDOT + REAL SDOT, SDSDOT + EXTERNAL SDOT, SDSDOT * .. External Subroutines .. - EXTERNAL SAXPY, SCOPY, SSWAP, STEST, STEST1 + EXTERNAL SAXPY, SCOPY, SROTM, SSWAP, STEST, STEST1 * .. Intrinsic Functions .. INTRINSIC ABS, MIN * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. + EQUIVALENCE (DT19X(1,1,1),DT19XA(1,1,1)),(DT19X(1,1,5), + A DT19XB(1,1,1)),(DT19X(1,1,9),DT19XC(1,1,1)), + B (DT19X(1,1,13),DT19XD(1,1,1)) + EQUIVALENCE (DT19Y(1,1,1),DT19YA(1,1,1)),(DT19Y(1,1,5), + A DT19YB(1,1,1)),(DT19Y(1,1,9),DT19YC(1,1,1)), + B (DT19Y(1,1,13),DT19YD(1,1,1)) + DATA SA/0.3E0/ DATA INCXS/1, 2, -2, -1/ DATA INCYS/1, -2, 1, -2/ @@ -272,10 +383,11 @@ + -0.4E0/ DATA DY1/0.5E0, -0.9E0, 0.3E0, 0.7E0, -0.6E0, 0.2E0, + 0.8E0/ - DATA SC, SS/0.8E0, 0.6E0/ DATA DT7/0.0E0, 0.30E0, 0.21E0, 0.62E0, 0.0E0, + 0.30E0, -0.07E0, 0.85E0, 0.0E0, 0.30E0, -0.79E0, + -0.74E0, 0.0E0, 0.30E0, 0.33E0, 1.27E0/ + DATA ST7B/ .1, .4, .31, .72, .1, .4, .03, .95, + + .1, .4, -.69, -.64, .1, .4, .43, 1.37/ DATA DT8/0.5E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, + 0.0E0, 0.68E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, + 0.0E0, 0.0E0, 0.68E0, -0.87E0, 0.0E0, 0.0E0, @@ -295,44 +407,6 @@ + 0.0E0, 0.68E0, -0.9E0, 0.33E0, 0.0E0, 0.0E0, + 0.0E0, 0.0E0, 0.68E0, -0.9E0, 0.33E0, 0.7E0, + -0.75E0, 0.2E0, 1.04E0/ - DATA DT9X/0.6E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, - + 0.0E0, 0.78E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, - + 0.0E0, 0.0E0, 0.78E0, -0.46E0, 0.0E0, 0.0E0, - + 0.0E0, 0.0E0, 0.0E0, 0.78E0, -0.46E0, -0.22E0, - + 1.06E0, 0.0E0, 0.0E0, 0.0E0, 0.6E0, 0.0E0, - + 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.78E0, - + 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, - + 0.66E0, 0.1E0, -0.1E0, 0.0E0, 0.0E0, 0.0E0, - + 0.0E0, 0.96E0, 0.1E0, -0.76E0, 0.8E0, 0.90E0, - + -0.3E0, -0.02E0, 0.6E0, 0.0E0, 0.0E0, 0.0E0, - + 0.0E0, 0.0E0, 0.0E0, 0.78E0, 0.0E0, 0.0E0, - + 0.0E0, 0.0E0, 0.0E0, 0.0E0, -0.06E0, 0.1E0, - + -0.1E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.90E0, - + 0.1E0, -0.22E0, 0.8E0, 0.18E0, -0.3E0, -0.02E0, - + 0.6E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, - + 0.78E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, - + 0.0E0, 0.78E0, 0.26E0, 0.0E0, 0.0E0, 0.0E0, - + 0.0E0, 0.0E0, 0.78E0, 0.26E0, -0.76E0, 1.12E0, - + 0.0E0, 0.0E0, 0.0E0/ - DATA DT9Y/0.5E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, - + 0.0E0, 0.04E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, - + 0.0E0, 0.0E0, 0.04E0, -0.78E0, 0.0E0, 0.0E0, - + 0.0E0, 0.0E0, 0.0E0, 0.04E0, -0.78E0, 0.54E0, - + 0.08E0, 0.0E0, 0.0E0, 0.0E0, 0.5E0, 0.0E0, - + 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.04E0, - + 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.7E0, - + -0.9E0, -0.12E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, - + 0.64E0, -0.9E0, -0.30E0, 0.7E0, -0.18E0, 0.2E0, - + 0.28E0, 0.5E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, - + 0.0E0, 0.0E0, 0.04E0, 0.0E0, 0.0E0, 0.0E0, - + 0.0E0, 0.0E0, 0.0E0, 0.7E0, -1.08E0, 0.0E0, - + 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.64E0, -1.26E0, - + 0.54E0, 0.20E0, 0.0E0, 0.0E0, 0.0E0, 0.5E0, - + 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, - + 0.04E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, - + 0.0E0, 0.04E0, -0.9E0, 0.18E0, 0.0E0, 0.0E0, - + 0.0E0, 0.0E0, 0.04E0, -0.9E0, 0.18E0, 0.7E0, - + -0.18E0, 0.2E0, 0.16E0/ DATA DT10X/0.6E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, + 0.0E0, 0.5E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, + 0.0E0, 0.5E0, -0.9E0, 0.0E0, 0.0E0, 0.0E0, @@ -375,6 +449,151 @@ + 0.0E0, 1.17E0, 1.17E0, 1.17E0, 1.17E0, 1.17E0, + 1.17E0, 1.17E0, 1.17E0, 1.17E0, 1.17E0, 1.17E0, + 1.17E0, 1.17E0, 1.17E0/ + DATA SSIZE3/ .1, .4, 1.7, 3.3 / +* +* FOR DROTM +* + DATA DPAR/-2.E0, 0.E0,0.E0,0.E0,0.E0, + A -1.E0, 2.E0, -3.E0, -4.E0, 5.E0, + B 0.E0, 0.E0, 2.E0, -3.E0, 0.E0, + C 1.E0, 5.E0, 2.E0, 0.E0, -4.E0/ +* TRUE X RESULTS F0R ROTATIONS DROTM + DATA DT19XA/.6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + A .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + B .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + C .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + D .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + E -.8E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + F -.9E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + G 3.5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + H .6E0, .1E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + I -.8E0, 3.8E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + J -.9E0, 2.8E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + K 3.5E0, -.4E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + L .6E0, .1E0, -.5E0, .8E0, 0.E0,0.E0,0.E0, + M -.8E0, 3.8E0, -2.2E0, -1.2E0, 0.E0,0.E0,0.E0, + N -.9E0, 2.8E0, -1.4E0, -1.3E0, 0.E0,0.E0,0.E0, + O 3.5E0, -.4E0, -2.2E0, 4.7E0, 0.E0,0.E0,0.E0/ +* + DATA DT19XB/.6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + A .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + B .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + C .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + D .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + E -.8E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + F -.9E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + G 3.5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + H .6E0, .1E0, -.5E0, 0.E0,0.E0,0.E0,0.E0, + I 0.E0, .1E0, -3.0E0, 0.E0,0.E0,0.E0,0.E0, + J -.3E0, .1E0, -2.0E0, 0.E0,0.E0,0.E0,0.E0, + K 3.3E0, .1E0, -2.0E0, 0.E0,0.E0,0.E0,0.E0, + L .6E0, .1E0, -.5E0, .8E0, .9E0, -.3E0, -.4E0, + M -2.0E0, .1E0, 1.4E0, .8E0, .6E0, -.3E0, -2.8E0, + N -1.8E0, .1E0, 1.3E0, .8E0, 0.E0, -.3E0, -1.9E0, + O 3.8E0, .1E0, -3.1E0, .8E0, 4.8E0, -.3E0, -1.5E0 / +* + DATA DT19XC/.6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + A .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + B .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + C .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + D .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + E -.8E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + F -.9E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + G 3.5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + H .6E0, .1E0, -.5E0, 0.E0,0.E0,0.E0,0.E0, + I 4.8E0, .1E0, -3.0E0, 0.E0,0.E0,0.E0,0.E0, + J 3.3E0, .1E0, -2.0E0, 0.E0,0.E0,0.E0,0.E0, + K 2.1E0, .1E0, -2.0E0, 0.E0,0.E0,0.E0,0.E0, + L .6E0, .1E0, -.5E0, .8E0, .9E0, -.3E0, -.4E0, + M -1.6E0, .1E0, -2.2E0, .8E0, 5.4E0, -.3E0, -2.8E0, + N -1.5E0, .1E0, -1.4E0, .8E0, 3.6E0, -.3E0, -1.9E0, + O 3.7E0, .1E0, -2.2E0, .8E0, 3.6E0, -.3E0, -1.5E0 / +* + DATA DT19XD/.6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + A .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + B .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + C .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + D .6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + E -.8E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + F -.9E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + G 3.5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + H .6E0, .1E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + I -.8E0, -1.0E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + J -.9E0, -.8E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + K 3.5E0, .8E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + L .6E0, .1E0, -.5E0, .8E0, 0.E0,0.E0,0.E0, + M -.8E0, -1.0E0, 1.4E0, -1.6E0, 0.E0,0.E0,0.E0, + N -.9E0, -.8E0, 1.3E0, -1.6E0, 0.E0,0.E0,0.E0, + O 3.5E0, .8E0, -3.1E0, 4.8E0, 0.E0,0.E0,0.E0/ +* TRUE Y RESULTS FOR ROTATIONS DROTM + DATA DT19YA/.5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + A .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + B .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + C .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + D .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + E .7E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + F 1.7E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + G -2.6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + H .5E0, -.9E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + I .7E0, -4.8E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + J 1.7E0, -.7E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + K -2.6E0, 3.5E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + L .5E0, -.9E0, .3E0, .7E0, 0.E0,0.E0,0.E0, + M .7E0, -4.8E0, 3.0E0, 1.1E0, 0.E0,0.E0,0.E0, + N 1.7E0, -.7E0, -.7E0, 2.3E0, 0.E0,0.E0,0.E0, + O -2.6E0, 3.5E0, -.7E0, -3.6E0, 0.E0,0.E0,0.E0/ +* + DATA DT19YB/.5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + A .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + B .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + C .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + D .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + E .7E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + F 1.7E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + G -2.6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + H .5E0, -.9E0, .3E0, 0.E0,0.E0,0.E0,0.E0, + I 4.0E0, -.9E0, -.3E0, 0.E0,0.E0,0.E0,0.E0, + J -.5E0, -.9E0, 1.5E0, 0.E0,0.E0,0.E0,0.E0, + K -1.5E0, -.9E0, -1.8E0, 0.E0,0.E0,0.E0,0.E0, + L .5E0, -.9E0, .3E0, .7E0, -.6E0, .2E0, .8E0, + M 3.7E0, -.9E0, -1.2E0, .7E0, -1.5E0, .2E0, 2.2E0, + N -.3E0, -.9E0, 2.1E0, .7E0, -1.6E0, .2E0, 2.0E0, + O -1.6E0, -.9E0, -2.1E0, .7E0, 2.9E0, .2E0, -3.8E0 / +* + DATA DT19YC/.5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + A .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + B .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + C .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + D .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + E .7E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + F 1.7E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + G -2.6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + H .5E0, -.9E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + I 4.0E0, -6.3E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + J -.5E0, .3E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + K -1.5E0, 3.0E0, 0.E0,0.E0,0.E0,0.E0,0.E0, + L .5E0, -.9E0, .3E0, .7E0, 0.E0,0.E0,0.E0, + M 3.7E0, -7.2E0, 3.0E0, 1.7E0, 0.E0,0.E0,0.E0, + N -.3E0, .9E0, -.7E0, 1.9E0, 0.E0,0.E0,0.E0, + O -1.6E0, 2.7E0, -.7E0, -3.4E0, 0.E0,0.E0,0.E0/ +* + DATA DT19YD/.5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + A .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + B .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + C .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + D .5E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + E .7E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + F 1.7E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + G -2.6E0, 0.E0,0.E0,0.E0,0.E0,0.E0,0.E0, + H .5E0, -.9E0, .3E0, 0.E0,0.E0,0.E0,0.E0, + I .7E0, -.9E0, 1.2E0, 0.E0,0.E0,0.E0,0.E0, + J 1.7E0, -.9E0, .5E0, 0.E0,0.E0,0.E0,0.E0, + K -2.6E0, -.9E0, -1.3E0, 0.E0,0.E0,0.E0,0.E0, + L .5E0, -.9E0, .3E0, .7E0, -.6E0, .2E0, .8E0, + M .7E0, -.9E0, 1.2E0, .7E0, -1.5E0, .2E0, 1.6E0, + N 1.7E0, -.9E0, .5E0, .7E0, -1.6E0, .2E0, 2.4E0, + O -2.6E0, -.9E0, -1.3E0, .7E0, 2.9E0, .2E0, -4.0E0 / +* * .. Executable Statements .. * DO 120 KI = 1, 4 @@ -421,6 +640,39 @@ 80 CONTINUE CALL STEST(LENX,SX,STX,SSIZE2(1,1),1.0E0) CALL STEST(LENY,SY,STY,SSIZE2(1,1),1.0E0) + ELSEIF (ICASE.EQ.12) THEN +* .. SROTM .. + KNI=KN+4*(KI-1) + DO KPAR=1,4 + DO I=1,7 + SX(I) = DX1(I) + SY(I) = DY1(I) + STX(I)= DT19X(I,KPAR,KNI) + STY(I)= DT19Y(I,KPAR,KNI) + END DO +* + DO I=1,5 + DTEMP(I) = DPAR(I,KPAR) + END DO +* + DO I=1,LENX + SSIZE(I)=STX(I) + END DO +* SEE REMARK ABOVE ABOUT DT11X(1,2,7) +* AND DT11X(5,3,8). + IF ((KPAR .EQ. 2) .AND. (KNI .EQ. 7)) + $ SSIZE(1) = 2.4E0 + IF ((KPAR .EQ. 3) .AND. (KNI .EQ. 8)) + $ SSIZE(5) = 1.8E0 +* + CALL SROTM(N,SX,INCX,SY,INCY,DTEMP) + CALL STEST(LENX,SX,STX,SSIZE,SFAC) + CALL STEST(LENY,SY,STY,STY,SFAC) + END DO + ELSEIF (ICASE.EQ.13) THEN +* .. SDSROT .. + CALL STEST1 (SDSDOT(N,.1,SX,INCX,SY,INCY), + $ ST7B(KN,KI),SSIZE3(KN),SFAC) ELSE WRITE (NOUT,*) ' Shouldn''t be here in CHECK2' STOP @@ -436,10 +688,10 @@ * .. Scalar Arguments .. REAL SFAC * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. - REAL SA, SC, SS + REAL SC, SS INTEGER I, K, KI, KN, KSIZE, LENX, LENY, MX, MY * .. Local Arrays .. REAL COPYX(5), COPYY(5), DT9X(7,4,4), DT9Y(7,4,4), @@ -454,9 +706,8 @@ * .. Intrinsic Functions .. INTRINSIC ABS, MIN * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Data statements .. - DATA SA/0.3E0/ DATA INCXS/1, 2, -2, -1/ DATA INCYS/1, -2, 1, -2/ DATA LENS/1, 1, 2, 4, 1, 1, 3, 7/ @@ -647,14 +898,15 @@ * * .. Parameters .. INTEGER NOUT - PARAMETER (NOUT=6) + REAL ZERO + PARAMETER (NOUT=6, ZERO=0.0E0) * .. Scalar Arguments .. REAL SFAC INTEGER LEN * .. Array Arguments .. REAL SCOMP(LEN), SSIZE(LEN), STRUE(LEN) * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. REAL SD @@ -665,12 +917,12 @@ * .. Intrinsic Functions .. INTRINSIC ABS * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Executable Statements .. * DO 40 I = 1, LEN SD = SCOMP(I) - STRUE(I) - IF (SDIFF(ABS(SSIZE(I))+ABS(SFAC*SD),ABS(SSIZE(I))).EQ.0.0E0) + IF (ABS(SFAC*SD) .LE. ABS(SSIZE(I))*EPSILON(ZERO)) + GO TO 40 * * HERE SCOMP(I) IS NOT CLOSE TO STRUE(I). @@ -680,16 +932,16 @@ PASS = .FALSE. WRITE (NOUT,99999) WRITE (NOUT,99998) - 20 WRITE (NOUT,99997) ICASE, N, INCX, INCY, MODE, I, SCOMP(I), + 20 WRITE (NOUT,99997) ICASE, N, INCX, INCY, I, SCOMP(I), + STRUE(I), SD, SSIZE(I) 40 CONTINUE RETURN * 99999 FORMAT (' FAIL') -99998 FORMAT (/' CASE N INCX INCY MODE I ', +99998 FORMAT (/' CASE N INCX INCY I ', + ' COMP(I) TRUE(I) DIFFERENCE', + ' SIZE(I)',/1X) -99997 FORMAT (1X,I4,I3,3I5,I3,2E36.8,2E12.4) +99997 FORMAT (1X,I4,I3,2I5,I3,2E36.8,2E12.4) END SUBROUTINE STEST1(SCOMP1,STRUE1,SSIZE,SFAC) * ************************* STEST1 ***************************** @@ -739,12 +991,12 @@ * .. Scalar Arguments .. INTEGER ICOMP, ITRUE * .. Scalars in Common .. - INTEGER ICASE, INCX, INCY, MODE, N + INTEGER ICASE, INCX, INCY, N LOGICAL PASS * .. Local Scalars .. INTEGER ID * .. Common blocks .. - COMMON /COMBLA/ICASE, N, INCX, INCY, MODE, PASS + COMMON /COMBLA/ICASE, N, INCX, INCY, PASS * .. Executable Statements .. * IF (ICOMP.EQ.ITRUE) GO TO 40 @@ -757,13 +1009,13 @@ WRITE (NOUT,99999) WRITE (NOUT,99998) 20 ID = ICOMP - ITRUE - WRITE (NOUT,99997) ICASE, N, INCX, INCY, MODE, ICOMP, ITRUE, ID + WRITE (NOUT,99997) ICASE, N, INCX, INCY, ICOMP, ITRUE, ID 40 CONTINUE RETURN * 99999 FORMAT (' FAIL') -99998 FORMAT (/' CASE N INCX INCY MODE ', +99998 FORMAT (/' CASE N INCX INCY ', + ' COMP TRUE DIFFERENCE', + /1X) -99997 FORMAT (1X,I4,I3,3I5,2I36,I12) +99997 FORMAT (1X,I4,I3,2I5,2I36,I12) END diff --git a/doc/I02_HiPerformance.dox b/doc/I02_HiPerformance.dox index d7a02fb5c..ab6cdfd44 100644 --- a/doc/I02_HiPerformance.dox +++ b/doc/I02_HiPerformance.dox @@ -79,7 +79,7 @@ temp = m2 * m3; m1 += temp.adjoint(); \endcode \code m1.noalias() += m3.adjoint() -* * m2.adjoint(); \endcode +* * m2.adjoint(); \endcode This is because the product expression has the EvalBeforeNesting bit which enforces the evaluation of the product by the Tranpose expression. diff --git a/doc/I10_Assertions.dox b/doc/I10_Assertions.dox index d5697fcee..e5bcbe536 100644 --- a/doc/I10_Assertions.dox +++ b/doc/I10_Assertions.dox @@ -2,12 +2,113 @@ namespace Eigen { /** \page TopicAssertions Assertions +\b Table \b of \b contents + - \ref PlainAssert + - \ref RedefineAssert + - \ref DisableAssert + - \ref StaticAssert + - \ref DerivedStaticAssert + - \ref DisableStaticAssert -TODO: write this dox page! +\section PlainAssert Assertions -Is linked from the tutorial on matrix arithmetic. +The macro eigen_assert is defined to be \c eigen_plain_assert by default. We use eigen_plain_assert instead of \c assert to work around a known bug for GCC <= 4.3. Basically, eigen_plain_assert \a is \c assert. -\sa Section \ref TopicPreprocessorDirectivesAssertions on page \ref TopicPreprocessorDirectives. +\subsection RedefineAssert Redefining assertions + +Both eigen_assert and eigen_plain_assert are defined in Macros.h. Defining eigen_assert indirectly gives you a chance to change its behavior. You can redefine this macro if you want to do something else such as throwing an exception, and fall back to its default behavior with eigen_plain_assert. The code below tells Eigen to throw an std::runtime_error: + +\code +#include +#undef eigen_assert +#define eigen_assert(x) \ + if (!x) { throw (std::runtime_error("Put your message here")); } +\endcode + +\subsection DisableAssert Disabling assertions + +Assertions cost run time and can be turned off. You can suppress eigen_assert by defining \c EIGEN_NO_DEBUG \b before including Eigen headers. \c EIGEN_NO_DEBUG is undefined by default unless \c NDEBUG is defined. + +\section StaticAssert Static assertions + +Static assertions are not standardized until C++11. However, in the Eigen library, there are many conditions can and should be detectedat compile time. For instance, we use static assertions to prevent the code below from compiling. + +\code +Matrix3d() + Matrix4d(); // adding matrices of different sizes +Matrix4cd() * Vector3cd(); // invalid product known at compile time +\endcode + +Static assertions are defined in StaticAssert.h. If there is native static_assert, we use it. Otherwise, we have implemented an assertion macro that can show a limited range of messages. + +One can easily come up with static assertions without messages, such as: + +\code +#define STATIC_ASSERT(x) \ + switch(0) { case 0: case x:; } +\endcode + +However, the example above obviously cannot tell why the assertion failed. Therefore, we define a \c struct in namespace Eigen::internal to handle available messages. + +\code +template +struct static_assertion {}; + +template<> +struct static_assertion +{ + enum { + YOU_TRIED_CALLING_A_VECTOR_METHOD_ON_A_MATRIX, + YOU_MIXED_VECTORS_OF_DIFFERENT_SIZES, + // see StaticAssert.h for all enums. + }; +}; +\endcode + +And then, we define EIGEN_STATIC_ASSERT(CONDITION,MSG) to access Eigen::internal::static_assertion::MSG. If the condition evaluates into \c false, your compiler displays a lot of messages explaining there is no MSG in static_assert. Nevertheless, this is \a not in what we are interested. As you can see, all members of static_assert are ALL_CAPS_AND_THEY_ARE_SHOUTING. + +\warning +When using this macro, MSG should be a member of static_assertion, or the static assertion \b always fails. +Currently, it can only be used in function scope. + +\subsection DerivedStaticAssert Derived static assertions + +There are other macros derived from EIGEN_STATIC_ASSERT to enhance readability. Their names are self-explanatory. + +- \b EIGEN_STATIC_ASSERT_FIXED_SIZE(TYPE) - passes if \a TYPE is fixed size. +- \b EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(TYPE) - passes if \a TYPE is dynamic size. +- \b EIGEN_STATIC_ASSERT_LVALUE(Derived) - failes if \a Derived is read-only. +- \b EIGEN_STATIC_ASSERT_ARRAYXPR(Derived) - passes if \a Derived is an array expression. +- EIGEN_STATIC_ASSERT_SAME_XPR_KIND(Derived1, Derived2) - failes if the two expressions are an array one and a matrix one. + +Because Eigen handles both fixed-size and dynamic-size expressions, some conditions cannot be clearly determined at compile time. We classify them into strict assertions and permissive assertions. + +\subsubsection StrictAssertions Strict assertions + +These assertions fail if the condition may not be met. For example, MatrixXd may not be a vector, so it fails EIGEN_STATIC_ASSERT_VECTOR_ONLY. + +- \b EIGEN_STATIC_ASSERT_VECTOR_ONLY(TYPE) - passes if \a TYPE must be a vector type. +- EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(TYPE, SIZE) - passes if \a TYPE must be a vector of the given size. +- EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(TYPE, ROWS, COLS) - passes if \a TYPE must be a matrix with given rows and columns. + +\subsubsection PermissiveAssertions Permissive assertions + +These assertions fail if the condition \b cannot be met. For example, MatrixXd and Matrix4d may have the same size, so they pass EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE. + +- \b EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(TYPE0,TYPE1) - fails if the two vector expression types must have different sizes. +- \b EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(TYPE0,TYPE1) - fails if the two matrix expression types must have different sizes. +- \b EIGEN_STATIC_ASSERT_SIZE_1x1(TYPE) - fails if \a TYPE cannot be an 1x1 expression. + +See StaticAssert.h for details such as what messages they throw. + +\subsection DisableStaticAssert Disabling static assertions + +If \c EIGEN_NO_STATIC_ASSERT is defined, static assertions turn into eigen_assert's, working like: + +\code +#define EIGEN_STATIC_ASSERT(CONDITION,MSG) eigen_assert((CONDITION) && #MSG); +\endcode + +This saves compile time but consumes more run time. \c EIGEN_NO_STATIC_ASSERT is undefined by default. */ } diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h index 18bcf3d0d..e1e5b770c 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h @@ -51,7 +51,6 @@ private: void compute2x2(const MatrixType& A, MatrixType& result); void computeBig(const MatrixType& A, MatrixType& result); - static Scalar atanh2(Scalar y, Scalar x); int getPadeDegree(float normTminusI); int getPadeDegree(double normTminusI); int getPadeDegree(long double normTminusI); @@ -93,20 +92,6 @@ MatrixType MatrixLogarithmAtomic::compute(const MatrixType& A) return result; } -/** \brief Compute atanh (inverse hyperbolic tangent) for \f$ y / x \f$. */ -template -typename MatrixType::Scalar MatrixLogarithmAtomic::atanh2(Scalar y, Scalar x) -{ - using std::abs; - using std::sqrt; - - Scalar z = y / x; - if (abs(z) > sqrt(NumTraits::epsilon())) - return Scalar(0.5) * log((x + y) / (x - y)); - else - return z + z*z*z / Scalar(3); -} - /** \brief Compute logarithm of 2x2 triangular matrix. */ template void MatrixLogarithmAtomic::compute2x2(const MatrixType& A, MatrixType& result) @@ -131,7 +116,7 @@ void MatrixLogarithmAtomic::compute2x2(const MatrixType& A, MatrixTy // computation in previous branch is inaccurate if A(1,1) \approx A(0,0) int unwindingNumber = static_cast(ceil((imag(logA11 - logA00) - M_PI) / (2*M_PI))); Scalar y = A(1,1) - A(0,0), x = A(1,1) + A(0,0); - result(0,1) = A(0,1) * (Scalar(2) * atanh2(y,x) + Scalar(0,2*M_PI*unwindingNumber)) / y; + result(0,1) = A(0,1) * (Scalar(2) * internal::atanh2(y,x) + Scalar(0,2*M_PI*unwindingNumber)) / y; } } diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index 2a46d2cc0..7238501ed 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -111,9 +111,6 @@ class MatrixPower */ void getFractionalExponent(); - /** \brief Compute atanh (inverse hyperbolic tangent) for \f$ y / x \f$. */ - static ComplexScalar atanh2(const ComplexScalar& y, const ComplexScalar& x); - /** \brief Compute power of 2x2 triangular matrix. */ void compute2x2(RealScalar p); @@ -223,7 +220,7 @@ void MatrixPower::computeChainProduct(ResultType& result int cost = computeCost(p); if (m_pInt < RealScalar(0)) { - if (p * m_dimb <= cost * m_dimA) { + if (p * m_dimb <= cost * m_dimA && m_dimA > 2) { partialPivLuSolve(result, p); return; } else { @@ -296,21 +293,6 @@ void MatrixPower::getFractionalExponent() } } -template -std::complex -MatrixPower::atanh2(const ComplexScalar& y, const ComplexScalar& x) -{ - using std::abs; - using std::log; - using std::sqrt; - const ComplexScalar z = y / x; - - if (abs(z) > sqrt(NumTraits::epsilon())) - return RealScalar(0.5) * log((x + y) / (x - y)); - else - return z + z*z*z / RealScalar(3); -} - template void MatrixPower::compute2x2(RealScalar p) { @@ -337,7 +319,7 @@ void MatrixPower::compute2x2(RealScalar p) } else { // computation in previous branch is inaccurate if abs(m_T(j,j)) \approx abs(m_T(i,i)) unwindingNumber = ceil((imag(m_logTdiag[j] - m_logTdiag[i]) - M_PI) / (2 * M_PI)); - w = atanh2(m_T(j,j) - m_T(i,i), m_T(j,j) + m_T(i,i)) + ComplexScalar(0, M_PI * unwindingNumber); + w = internal::atanh2(m_T(j,j) - m_T(i,i), m_T(j,j) + m_T(i,i)) + ComplexScalar(0, M_PI * unwindingNumber); m_fT(i,j) = m_T(i,j) * RealScalar(2) * exp(RealScalar(0.5) * p * (m_logTdiag[j] + m_logTdiag[i])) * sinh(p * w) / (m_T(j,j) - m_T(i,i)); }