GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-04 18:54:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merged in jdh8/eigen (pull request PR-16)

Browse Source
This commit is contained in:
Gael Guennebaud 2012-09-08 12:16:49 +02:00
commit 24f371bdb4
16 changed files with 1112 additions and 291 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
Eigen/src/Core/MathFunctions.h
View File

@ -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<typename Scalar, bool IsInteger>
struct atanh2_default_impl
{
typedef Scalar retval;
typedef typename NumTraits<Scalar>::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<RealScalar>::epsilon()))
return RealScalar(0.5) * log((y + x) / (y - x));
else
return z + z*z*z / RealScalar(3);
}
};
template<typename Scalar>
struct atanh2_default_impl<Scalar, true>
{
static inline Scalar run(const Scalar&, const Scalar&)
{
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
return Scalar(0);
}
};
template<typename Scalar>
struct atanh2_impl : atanh2_default_impl<Scalar, NumTraits<Scalar>::IsInteger> {};
template<typename Scalar>
struct atanh2_retval
{
typedef Scalar type;
};
template<typename Scalar>
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 *
****************************************************************************/

30
blas/GeneralRank1Update.h Normal file
View File

@ -0,0 +1,30 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
//
// 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<typename Scalar, typename Index, bool ConjRhs>
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<Scalar,Dynamic,1> PlainVector;
internal::conj_if<ConjRhs> cj;
for (Index i=0; i<cols; ++i)
Map<PlainVector>(mat+stride*i,rows) += alpha * cj(v[i]) * Map<const PlainVector>(u,rows);
}
};
} // end namespace internal
#endif // EIGEN_GENERAL_RANK1UPDATE_H

57
blas/Rank2Update.h Normal file
View File

@ -0,0 +1,57 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
//
// 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<typename Scalar, typename Index, int UpLo>
struct rank2_update_selector;
template<typename Scalar, typename Index>
struct rank2_update_selector<Scalar,Index,Upper>
{
static void run(Index size, Scalar* mat, Index stride, const Scalar* _u, const Scalar* _v, Scalar alpha)
{
typedef Matrix<Scalar,Dynamic,1> PlainVector;
Map<const PlainVector> u(_u, size), v(_v, size);
for (Index i=0; i<size; ++i)
{
Map<PlainVector>(mat+stride*i, i+1) +=
conj(alpha) * conj(_u[i]) * v.head(i+1)
+ alpha * conj(_v[i]) * u.head(i+1);
}
}
};
template<typename Scalar, typename Index>
struct rank2_update_selector<Scalar,Index,Lower>
{
static void run(Index size, Scalar* mat, Index stride, const Scalar* _u, const Scalar* _v, Scalar alpha)
{
typedef Matrix<Scalar,Dynamic,1> PlainVector;
Map<const PlainVector> u(_u, size), v(_v, size);
for (Index i=0; i<size; ++i)
{
Map<PlainVector>(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

4
blas/common.h
View File

@ -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;

14
blas/double.cpp
View File

@ -2,6 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
//
// 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<double>().cwiseProduct(vector(y,*n).cast<double>())).sum();
else if(*incx>0 && *incy>0) return (vector(x,*n,*incx).cast<double>().cwiseProduct(vector(y,*n,*incy).cast<double>())).sum();
else if(*incx<0 && *incy>0) return (vector(x,*n,-*incx).reverse().cast<double>().cwiseProduct(vector(y,*n,*incy).cast<double>())).sum();
else if(*incx>0 && *incy<0) return (vector(x,*n,*incx).cast<double>().cwiseProduct(vector(y,*n,-*incy).reverse().cast<double>())).sum();
else if(*incx<0 && *incy<0) return (vector(x,*n,-*incx).reverse().cast<double>().cwiseProduct(vector(y,*n,-*incy).reverse().cast<double>())).sum();
else return 0;
}

57
blas/level2_cplx_impl.h
View File

@ -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<Scalar,int,ColMajor,Upper,false,Conj>::run);
func[LO] = (selfadjoint_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run);
init = true;
}
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* a = reinterpret_cast<Scalar*>(pa);
RealScalar alpha = *reinterpret_cast<RealScalar*>(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<Lower>().rankUpdate(vector(x_cpy,*n), alpha);
// else if(UPLO(*uplo)==UP) matrix(a,*n,*n,*lda).selfadjointView<Upper>().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<Scalar,int,Upper>::run);
func[LO] = (internal::rank2_update_selector<Scalar,int,Lower>::run);
init = true;
}
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* y = reinterpret_cast<Scalar*>(py);
Scalar* a = reinterpret_cast<Scalar*>(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<Lower>().rankUpdate(vector(x_cpy,*n),vector(y_cpy,*n),alpha);
else if(UPLO(*uplo)==UP) matrix(a,*n,*n,*lda).selfadjointView<Upper>().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<Scalar,int,false>::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<Scalar,int,Conj>::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha);
if(x_cpy!=x) delete[] x_cpy;
if(y_cpy!=y) delete[] y_cpy;

46
blas/level2_impl.h
View File

@ -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<Scalar*>(px);
Scalar* y = reinterpret_cast<Scalar*>(py);
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar alpha = *reinterpret_cast<Scalar*>(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<std::max(1,*m)) info = 9;
if(info)
return xerbla_(SCALAR_SUFFIX_UP"GER ",&info,6);
if(alpha==Scalar(0))
return 1;
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();
if(x_cpy!=x) delete[] x_cpy;
if(y_cpy!=y) delete[] y_cpy;
return 1;
}

89
blas/level2_real_impl.h
View File

@ -68,6 +68,20 @@ int EIGEN_BLAS_FUNC(syr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px,
// init = true;
// }
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<Scalar,int,ColMajor,Upper,false,Conj>::run);
func[LO] = (selfadjoint_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run);
init = true;
}
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* c = reinterpret_cast<Scalar*>(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<Scalar,Dynamic,Dynamic> 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<Lower>().rankUpdate(vector(x_cpy,*n), alpha);
// else if(UPLO(*uplo)==UP) matrix(c,*n,*n,*ldc).selfadjointView<Upper>().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<Scalar,int,Upper>::run);
func[LO] = (internal::rank2_update_selector<Scalar,int,Lower>::run);
init = true;
}
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* y = reinterpret_cast<Scalar*>(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<Lower>().rankUpdate(vector(x_cpy,*n), vector(y_cpy,*n), alpha);
else if(UPLO(*uplo)==UP) matrix(c,*n,*n,*ldc).selfadjointView<Upper>().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<Scalar*>(px);
Scalar* y = reinterpret_cast<Scalar*>(py);
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar alpha = *reinterpret_cast<Scalar*>(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<std::max(1,*m)) info = 9;
if(info)
return xerbla_(SCALAR_SUFFIX_UP"GER ",&info,6);
if(alpha==Scalar(0))
return 1;
Scalar* x_cpy = get_compact_vector(x,*m,*incx);
Scalar* y_cpy = get_compact_vector(y,*n,*incy);
internal::general_rank1_update<Scalar,int,false>::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;
}

2
blas/level3_impl.h
View File

@ -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<Scalar*>(pa);
Scalar* c = reinterpret_cast<Scalar*>(pc);

3
blas/single.cpp
View File

@ -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); }

476
blas/testing/dblat1.f
View File

@ -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

430
blas/testing/sblat1.f
View File

@ -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

2
doc/I02_HiPerformance.dox
View File

@ -79,7 +79,7 @@ temp = m2 * m3;
m1 += temp.adjoint(); \endcode</td>
<td>\code
m1.noalias() += m3.adjoint()
* * m2.adjoint(); \endcode</td>
* * m2.adjoint(); \endcode</td>
<td>This is because the product expression has the EvalBeforeNesting bit which
enforces the evaluation of the product by the Tranpose expression.</td>
</tr>

107
doc/I10_Assertions.dox
View File

@ -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 <stdexcept>
#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<bool condition>
struct static_assertion {};
template<>
struct static_assertion<true>
{
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<bool(CONDITION)>::MSG. If the condition evaluates into \c false, your compiler displays a lot of messages explaining there is no MSG in static_assert<false>. Nevertheless, this is \a not in what we are interested. As you can see, all members of static_assert<true> are ALL_CAPS_AND_THEY_ARE_SHOUTING.
\warning
When using this macro, MSG should be a member of static_assertion<true>, 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.
- <b>EIGEN_STATIC_ASSERT_SAME_XPR_KIND(Derived1, Derived2)</b> - 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 <b>may not</b> 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.
- <b>EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(TYPE, SIZE)</b> - passes if \a TYPE must be a vector of the given size.
- <b>EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(TYPE, ROWS, COLS)</b> - 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 <tt>eigen_assert</tt>'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.
*/
}

17
unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h
View File

@ -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<MatrixType>::compute(const MatrixType& A)
return result;
}
/** \brief Compute atanh (inverse hyperbolic tangent) for \f$ y / x \f$. */
template <typename MatrixType>
typename MatrixType::Scalar MatrixLogarithmAtomic<MatrixType>::atanh2(Scalar y, Scalar x)
{
using std::abs;
using std::sqrt;
Scalar z = y / x;
if (abs(z) > sqrt(NumTraits<Scalar>::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 <typename MatrixType>
void MatrixLogarithmAtomic<MatrixType>::compute2x2(const MatrixType& A, MatrixType& result)
@ -131,7 +116,7 @@ void MatrixLogarithmAtomic<MatrixType>::compute2x2(const MatrixType& A, MatrixTy
// computation in previous branch is inaccurate if A(1,1) \approx A(0,0)
int unwindingNumber = static_cast<int>(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;
}
}

22
unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
View File

@ -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<MatrixType,PlainObject>::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<MatrixType,PlainObject>::getFractionalExponent()
}
}
template<typename MatrixType, typename PlainObject>
std::complex<typename MatrixType::RealScalar>
MatrixPower<MatrixType,PlainObject>::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<RealScalar>::epsilon()))
return RealScalar(0.5) * log((x + y) / (x - y));
else
return z + z*z*z / RealScalar(3);
}
template<typename MatrixType, typename PlainObject>
void MatrixPower<MatrixType,PlainObject>::compute2x2(RealScalar p)
{
@ -337,7 +319,7 @@ void MatrixPower<MatrixType,PlainObject>::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));
}