diff options
author | Carlos Hernandez <chernand@google.com> | 2014-08-05 17:53:32 -0700 |
---|---|---|
committer | Carlos Hernandez <chernand@google.com> | 2014-08-05 17:54:05 -0700 |
commit | 7faaa9f3f0df9d23790277834d426c3d992ac3ba (patch) | |
tree | b788ae3b96daf9f5a79d8ec434e1e9edd56b3a72 /blas | |
parent | 810535bb0c575a003b32076e5352ab8fd3f23a1c (diff) | |
download | eigen-7faaa9f3f0df9d23790277834d426c3d992ac3ba.tar.gz |
Update Eigen to the latest stable release, 3.2.2android-wear-5.1.1_r1android-wear-5.1.0_r1android-wear-5.0.0_r1android-l-preview_r2android-cts-5.1_r9android-cts-5.1_r8android-cts-5.1_r7android-cts-5.1_r6android-cts-5.1_r5android-cts-5.1_r4android-cts-5.1_r3android-cts-5.1_r28android-cts-5.1_r27android-cts-5.1_r26android-cts-5.1_r25android-cts-5.1_r24android-cts-5.1_r23android-cts-5.1_r22android-cts-5.1_r21android-cts-5.1_r20android-cts-5.1_r2android-cts-5.1_r19android-cts-5.1_r18android-cts-5.1_r17android-cts-5.1_r16android-cts-5.1_r15android-cts-5.1_r14android-cts-5.1_r13android-cts-5.1_r10android-cts-5.1_r1android-cts-5.0_r9android-cts-5.0_r8android-cts-5.0_r7android-cts-5.0_r6android-cts-5.0_r5android-cts-5.0_r4android-cts-5.0_r3android-5.1.1_r9android-5.1.1_r8android-5.1.1_r7android-5.1.1_r6android-5.1.1_r5android-5.1.1_r4android-5.1.1_r38android-5.1.1_r37android-5.1.1_r36android-5.1.1_r35android-5.1.1_r34android-5.1.1_r33android-5.1.1_r30android-5.1.1_r3android-5.1.1_r29android-5.1.1_r28android-5.1.1_r26android-5.1.1_r25android-5.1.1_r24android-5.1.1_r23android-5.1.1_r22android-5.1.1_r20android-5.1.1_r2android-5.1.1_r19android-5.1.1_r18android-5.1.1_r17android-5.1.1_r16android-5.1.1_r15android-5.1.1_r14android-5.1.1_r13android-5.1.1_r12android-5.1.1_r10android-5.1.1_r1android-5.1.0_r5android-5.1.0_r4android-5.1.0_r3android-5.1.0_r1android-5.0.2_r3android-5.0.2_r1android-5.0.1_r1android-5.0.0_r7android-5.0.0_r6android-5.0.0_r5.1android-5.0.0_r5android-5.0.0_r4android-5.0.0_r3android-5.0.0_r2android-5.0.0_r1lollipop-wear-releaselollipop-releaselollipop-mr1-wfc-releaselollipop-mr1-releaselollipop-mr1-fi-releaselollipop-mr1-devlollipop-mr1-cts-releaselollipop-devlollipop-cts-releasel-preview
./Eigen/src/Core/util/NonMPL2.h is left untouched, so that
usage of non MPL2 code is disabled.
Change-Id: I86fc9257b3c30d0ca15b268d4ef07bf038bba7ca
Diffstat (limited to 'blas')
34 files changed, 1549 insertions, 4622 deletions
diff --git a/blas/CMakeLists.txt b/blas/CMakeLists.txt index 453d5874c..a9bc05137 100644 --- a/blas/CMakeLists.txt +++ b/blas/CMakeLists.txt @@ -7,6 +7,9 @@ workaround_9220(Fortran EIGEN_Fortran_COMPILER_WORKS) if(EIGEN_Fortran_COMPILER_WORKS) enable_language(Fortran OPTIONAL) + if(NOT CMAKE_Fortran_COMPILER) + set(EIGEN_Fortran_COMPILER_WORKS OFF) + endif() endif() add_custom_target(blas) @@ -18,10 +21,10 @@ if(EIGEN_Fortran_COMPILER_WORKS) set(EigenBlas_SRCS ${EigenBlas_SRCS} complexdots.f srotm.f srotmg.f drotm.f drotmg.f - lsame.f chpr2.f dspmv.f dtpsv.f ssbmv.f sspr.f stpmv.f - zhpr2.f chbmv.f chpr.f ctpmv.f dspr2.f sspmv.f stpsv.f - zhbmv.f zhpr.f ztpmv.f chpmv.f ctpsv.f dsbmv.f dspr.f dtpmv.f sspr2.f - zhpmv.f ztpsv.f + lsame.f dspmv.f ssbmv.f + chbmv.f sspmv.f + zhbmv.f chpmv.f dsbmv.f + zhpmv.f dtbmv.f stbmv.f ctbmv.f ztbmv.f ) else() diff --git a/blas/GeneralRank1Update.h b/blas/GeneralRank1Update.h new file mode 100644 index 000000000..07d388c88 --- /dev/null +++ b/blas/GeneralRank1Update.h @@ -0,0 +1,44 @@ +// 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, int StorageOrder, bool ConjLhs, bool ConjRhs> +struct general_rank1_update; + +template<typename Scalar, typename Index, bool ConjLhs, bool ConjRhs> +struct general_rank1_update<Scalar,Index,ColMajor,ConjLhs,ConjRhs> +{ + static void run(Index rows, Index cols, Scalar* mat, Index stride, const Scalar* u, const Scalar* v, Scalar alpha) + { + typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap; + typedef typename conj_expr_if<ConjLhs,OtherMap>::type ConjRhsType; + conj_if<ConjRhs> cj; + + for (Index i=0; i<cols; ++i) + Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i,rows) += alpha * cj(v[i]) * ConjRhsType(OtherMap(u,rows)); + } +}; + +template<typename Scalar, typename Index, bool ConjLhs, bool ConjRhs> +struct general_rank1_update<Scalar,Index,RowMajor,ConjLhs,ConjRhs> +{ + static void run(Index rows, Index cols, Scalar* mat, Index stride, const Scalar* u, const Scalar* v, Scalar alpha) + { + general_rank1_update<Scalar,Index,ColMajor,ConjRhs,ConjRhs>::run(rows,cols,mat,stride,u,v,alpha); + } +}; + +} // end namespace internal + +#endif // EIGEN_GENERAL_RANK1UPDATE_H diff --git a/blas/PackedSelfadjointProduct.h b/blas/PackedSelfadjointProduct.h new file mode 100644 index 000000000..07327a246 --- /dev/null +++ b/blas/PackedSelfadjointProduct.h @@ -0,0 +1,53 @@ +// 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_SELFADJOINT_PACKED_PRODUCT_H +#define EIGEN_SELFADJOINT_PACKED_PRODUCT_H + +namespace internal { + +/* Optimized matrix += alpha * uv' + * The matrix is in packed form. + */ +template<typename Scalar, typename Index, int StorageOrder, int UpLo, bool ConjLhs, bool ConjRhs> +struct selfadjoint_packed_rank1_update; + +template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs> +struct selfadjoint_packed_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRhs> +{ + typedef typename NumTraits<Scalar>::Real RealScalar; + static void run(Index size, Scalar* mat, const Scalar* vec, RealScalar alpha) + { + typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap; + typedef typename conj_expr_if<ConjLhs,OtherMap>::type ConjRhsType; + conj_if<ConjRhs> cj; + + for (Index i=0; i<size; ++i) + { + Map<Matrix<Scalar,Dynamic,1> >(mat, UpLo==Lower ? size-i : (i+1)) += alpha * cj(vec[i]) * ConjRhsType(OtherMap(vec+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1))); + //FIXME This should be handled outside. + mat[UpLo==Lower ? 0 : i] = numext::real(mat[UpLo==Lower ? 0 : i]); + mat += UpLo==Lower ? size-i : (i+1); + } + } +}; + +template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs> +struct selfadjoint_packed_rank1_update<Scalar,Index,RowMajor,UpLo,ConjLhs,ConjRhs> +{ + typedef typename NumTraits<Scalar>::Real RealScalar; + static void run(Index size, Scalar* mat, const Scalar* vec, RealScalar alpha) + { + selfadjoint_packed_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,vec,alpha); + } +}; + +} // end namespace internal + +#endif // EIGEN_SELFADJOINT_PACKED_PRODUCT_H diff --git a/blas/PackedTriangularMatrixVector.h b/blas/PackedTriangularMatrixVector.h new file mode 100644 index 000000000..e9886d56f --- /dev/null +++ b/blas/PackedTriangularMatrixVector.h @@ -0,0 +1,79 @@ +// 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_PACKED_TRIANGULAR_MATRIX_VECTOR_H +#define EIGEN_PACKED_TRIANGULAR_MATRIX_VECTOR_H + +namespace internal { + +template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs, int StorageOrder> +struct packed_triangular_matrix_vector_product; + +template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs> +struct packed_triangular_matrix_vector_product<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,ColMajor> +{ + typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar; + enum { + IsLower = (Mode & Lower) ==Lower, + HasUnitDiag = (Mode & UnitDiag)==UnitDiag, + HasZeroDiag = (Mode & ZeroDiag)==ZeroDiag + }; + static void run(Index size, const LhsScalar* lhs, const RhsScalar* rhs, ResScalar* res, ResScalar alpha) + { + internal::conj_if<ConjRhs> cj; + typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap; + typedef typename conj_expr_if<ConjLhs,LhsMap>::type ConjLhsType; + typedef Map<Matrix<ResScalar,Dynamic,1> > ResMap; + + for (Index i=0; i<size; ++i) + { + Index s = IsLower&&(HasUnitDiag||HasZeroDiag) ? 1 : 0; + Index r = IsLower ? size-i: i+1; + if (EIGEN_IMPLIES(HasUnitDiag||HasZeroDiag, (--r)>0)) + ResMap(res+(IsLower ? s+i : 0),r) += alpha * cj(rhs[i]) * ConjLhsType(LhsMap(lhs+s,r)); + if (HasUnitDiag) + res[i] += alpha * cj(rhs[i]); + lhs += IsLower ? size-i: i+1; + } + }; +}; + +template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs> +struct packed_triangular_matrix_vector_product<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,RowMajor> +{ + typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar; + enum { + IsLower = (Mode & Lower) ==Lower, + HasUnitDiag = (Mode & UnitDiag)==UnitDiag, + HasZeroDiag = (Mode & ZeroDiag)==ZeroDiag + }; + static void run(Index size, const LhsScalar* lhs, const RhsScalar* rhs, ResScalar* res, ResScalar alpha) + { + internal::conj_if<ConjRhs> cj; + typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap; + typedef typename conj_expr_if<ConjLhs,LhsMap>::type ConjLhsType; + typedef Map<const Matrix<RhsScalar,Dynamic,1> > RhsMap; + typedef typename conj_expr_if<ConjRhs,RhsMap>::type ConjRhsType; + + for (Index i=0; i<size; ++i) + { + Index s = !IsLower&&(HasUnitDiag||HasZeroDiag) ? 1 : 0; + Index r = IsLower ? i+1 : size-i; + if (EIGEN_IMPLIES(HasUnitDiag||HasZeroDiag, (--r)>0)) + res[i] += alpha * (ConjLhsType(LhsMap(lhs+s,r)).cwiseProduct(ConjRhsType(RhsMap(rhs+(IsLower ? 0 : s+i),r)))).sum(); + if (HasUnitDiag) + res[i] += alpha * cj(rhs[i]); + lhs += IsLower ? i+1 : size-i; + } + }; +}; + +} // end namespace internal + +#endif // EIGEN_PACKED_TRIANGULAR_MATRIX_VECTOR_H diff --git a/blas/PackedTriangularSolverVector.h b/blas/PackedTriangularSolverVector.h new file mode 100644 index 000000000..5c0bb4bd6 --- /dev/null +++ b/blas/PackedTriangularSolverVector.h @@ -0,0 +1,88 @@ +// 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_PACKED_TRIANGULAR_SOLVER_VECTOR_H +#define EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H + +namespace internal { + +template<typename LhsScalar, typename RhsScalar, typename Index, int Side, int Mode, bool Conjugate, int StorageOrder> +struct packed_triangular_solve_vector; + +// forward and backward substitution, row-major, rhs is a vector +template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate> +struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, RowMajor> +{ + enum { + IsLower = (Mode&Lower)==Lower + }; + static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) + { + internal::conj_if<Conjugate> cj; + typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap; + typedef typename conj_expr_if<Conjugate,LhsMap>::type ConjLhsType; + + lhs += IsLower ? 0 : (size*(size+1)>>1)-1; + for(Index pi=0; pi<size; ++pi) + { + Index i = IsLower ? pi : size-pi-1; + Index s = IsLower ? 0 : 1; + if (pi>0) + rhs[i] -= (ConjLhsType(LhsMap(lhs+s,pi)) + .cwiseProduct(Map<const Matrix<RhsScalar,Dynamic,1> >(rhs+(IsLower ? 0 : i+1),pi))).sum(); + if (!(Mode & UnitDiag)) + rhs[i] /= cj(lhs[IsLower ? i : 0]); + IsLower ? lhs += pi+1 : lhs -= pi+2; + } + } +}; + +// forward and backward substitution, column-major, rhs is a vector +template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate> +struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, ColMajor> +{ + enum { + IsLower = (Mode&Lower)==Lower + }; + static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) + { + internal::conj_if<Conjugate> cj; + typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap; + typedef typename conj_expr_if<Conjugate,LhsMap>::type ConjLhsType; + + lhs += IsLower ? 0 : size*(size-1)>>1; + for(Index pi=0; pi<size; ++pi) + { + Index i = IsLower ? pi : size-pi-1; + Index r = size - pi - 1; + if (!(Mode & UnitDiag)) + rhs[i] /= cj(lhs[IsLower ? 0 : i]); + if (r>0) + Map<Matrix<RhsScalar,Dynamic,1> >(rhs+(IsLower? i+1 : 0),r) -= + rhs[i] * ConjLhsType(LhsMap(lhs+(IsLower? 1 : 0),r)); + IsLower ? lhs += size-pi : lhs -= r; + } + } +}; + +template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder> +struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheRight, Mode, Conjugate, StorageOrder> +{ + static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) + { + packed_triangular_solve_vector<LhsScalar,RhsScalar,Index,OnTheLeft, + ((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag), + Conjugate,StorageOrder==RowMajor?ColMajor:RowMajor + >::run(size, lhs, rhs); + } +}; + +} // end namespace internal + +#endif // EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H diff --git a/blas/README.txt b/blas/README.txt index 07a8bd92a..63a5203b9 100644 --- a/blas/README.txt +++ b/blas/README.txt @@ -1,9 +1,6 @@ This directory contains a BLAS library built on top of Eigen. -This is currently a work in progress which is far to be ready for use, -but feel free to contribute to it if you wish. - This module is not built by default. In order to compile it, you need to type 'make blas' from within your build dir. diff --git a/blas/Rank2Update.h b/blas/Rank2Update.h new file mode 100644 index 000000000..138d70fd6 --- /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 <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 +{ + static void run(Index size, Scalar* mat, Index stride, const Scalar* u, const Scalar* v, Scalar alpha) + { + typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap; + for (Index i=0; i<size; ++i) + { + Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1)) += + numext::conj(alpha) * numext::conj(u[i]) * OtherMap(v+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1)) + + alpha * numext::conj(v[i]) * OtherMap(u+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1)); + } + } +}; + +/* Optimized selfadjoint matrix += alpha * uv' + conj(alpha)*vu' + * The matrix is in packed form. + */ +template<typename Scalar, typename Index, int UpLo> +struct packed_rank2_update_selector +{ + static void run(Index size, Scalar* mat, const Scalar* u, const Scalar* v, Scalar alpha) + { + typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap; + Index offset = 0; + for (Index i=0; i<size; ++i) + { + Map<Matrix<Scalar,Dynamic,1> >(mat+offset, UpLo==Lower ? size-i : (i+1)) += + numext::conj(alpha) * numext::conj(u[i]) * OtherMap(v+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1)) + + alpha * numext::conj(v[i]) * OtherMap(u+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1)); + //FIXME This should be handled outside. + mat[offset+(UpLo==Lower ? 0 : i)] = numext::real(mat[offset+(UpLo==Lower ? 0 : i)]); + offset += UpLo==Lower ? size-i : (i+1); + } + } +}; + +} // end namespace internal + +#endif // EIGEN_RANK2UPDATE_H diff --git a/blas/chpr.f b/blas/chpr.f deleted file mode 100644 index 11bd5c6ee..000000000 --- a/blas/chpr.f +++ /dev/null @@ -1,220 +0,0 @@ - SUBROUTINE CHPR(UPLO,N,ALPHA,X,INCX,AP) -* .. Scalar Arguments .. - REAL ALPHA - INTEGER INCX,N - CHARACTER UPLO -* .. -* .. Array Arguments .. - COMPLEX AP(*),X(*) -* .. -* -* Purpose -* ======= -* -* CHPR performs the hermitian rank 1 operation -* -* A := alpha*x*conjg( x' ) + A, -* -* where alpha is a real scalar, x is an n element vector and A is an -* n by n hermitian matrix, supplied in packed form. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the upper or lower -* triangular part of the matrix A is supplied in the packed -* array AP as follows: -* -* UPLO = 'U' or 'u' The upper triangular part of A is -* supplied in AP. -* -* UPLO = 'L' or 'l' The lower triangular part of A is -* supplied in AP. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* ALPHA - REAL . -* On entry, ALPHA specifies the scalar alpha. -* Unchanged on exit. -* -* X - COMPLEX array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element vector x. -* Unchanged on exit. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* AP - COMPLEX array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular part of the hermitian matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) -* and a( 2, 2 ) respectively, and so on. On exit, the array -* AP is overwritten by the upper triangular part of the -* updated matrix. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular part of the hermitian matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) -* and a( 3, 1 ) respectively, and so on. On exit, the array -* AP is overwritten by the lower triangular part of the -* updated matrix. -* Note that the imaginary parts of the diagonal elements need -* not be set, they are assumed to be zero, and on exit they -* are set to zero. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - COMPLEX ZERO - PARAMETER (ZERO= (0.0E+0,0.0E+0)) -* .. -* .. Local Scalars .. - COMPLEX TEMP - INTEGER I,INFO,IX,J,JX,K,KK,KX -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC CONJG,REAL -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (N.LT.0) THEN - INFO = 2 - ELSE IF (INCX.EQ.0) THEN - INFO = 5 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('CHPR ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF ((N.EQ.0) .OR. (ALPHA.EQ.REAL(ZERO))) RETURN -* -* Set the start point in X if the increment is not unity. -* - IF (INCX.LE.0) THEN - KX = 1 - (N-1)*INCX - ELSE IF (INCX.NE.1) THEN - KX = 1 - END IF -* -* Start the operations. In this version the elements of the array AP -* are accessed sequentially with one pass through AP. -* - KK = 1 - IF (LSAME(UPLO,'U')) THEN -* -* Form A when upper triangle is stored in AP. -* - IF (INCX.EQ.1) THEN - DO 20 J = 1,N - IF (X(J).NE.ZERO) THEN - TEMP = ALPHA*CONJG(X(J)) - K = KK - DO 10 I = 1,J - 1 - AP(K) = AP(K) + X(I)*TEMP - K = K + 1 - 10 CONTINUE - AP(KK+J-1) = REAL(AP(KK+J-1)) + REAL(X(J)*TEMP) - ELSE - AP(KK+J-1) = REAL(AP(KK+J-1)) - END IF - KK = KK + J - 20 CONTINUE - ELSE - JX = KX - DO 40 J = 1,N - IF (X(JX).NE.ZERO) THEN - TEMP = ALPHA*CONJG(X(JX)) - IX = KX - DO 30 K = KK,KK + J - 2 - AP(K) = AP(K) + X(IX)*TEMP - IX = IX + INCX - 30 CONTINUE - AP(KK+J-1) = REAL(AP(KK+J-1)) + REAL(X(JX)*TEMP) - ELSE - AP(KK+J-1) = REAL(AP(KK+J-1)) - END IF - JX = JX + INCX - KK = KK + J - 40 CONTINUE - END IF - ELSE -* -* Form A when lower triangle is stored in AP. -* - IF (INCX.EQ.1) THEN - DO 60 J = 1,N - IF (X(J).NE.ZERO) THEN - TEMP = ALPHA*CONJG(X(J)) - AP(KK) = REAL(AP(KK)) + REAL(TEMP*X(J)) - K = KK + 1 - DO 50 I = J + 1,N - AP(K) = AP(K) + X(I)*TEMP - K = K + 1 - 50 CONTINUE - ELSE - AP(KK) = REAL(AP(KK)) - END IF - KK = KK + N - J + 1 - 60 CONTINUE - ELSE - JX = KX - DO 80 J = 1,N - IF (X(JX).NE.ZERO) THEN - TEMP = ALPHA*CONJG(X(JX)) - AP(KK) = REAL(AP(KK)) + REAL(TEMP*X(JX)) - IX = JX - DO 70 K = KK + 1,KK + N - J - IX = IX + INCX - AP(K) = AP(K) + X(IX)*TEMP - 70 CONTINUE - ELSE - AP(KK) = REAL(AP(KK)) - END IF - JX = JX + INCX - KK = KK + N - J + 1 - 80 CONTINUE - END IF - END IF -* - RETURN -* -* End of CHPR . -* - END diff --git a/blas/chpr2.f b/blas/chpr2.f deleted file mode 100644 index a0020ef3e..000000000 --- a/blas/chpr2.f +++ /dev/null @@ -1,255 +0,0 @@ - SUBROUTINE CHPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP) -* .. Scalar Arguments .. - COMPLEX ALPHA - INTEGER INCX,INCY,N - CHARACTER UPLO -* .. -* .. Array Arguments .. - COMPLEX AP(*),X(*),Y(*) -* .. -* -* Purpose -* ======= -* -* CHPR2 performs the hermitian rank 2 operation -* -* A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, -* -* where alpha is a scalar, x and y are n element vectors and A is an -* n by n hermitian matrix, supplied in packed form. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the upper or lower -* triangular part of the matrix A is supplied in the packed -* array AP as follows: -* -* UPLO = 'U' or 'u' The upper triangular part of A is -* supplied in AP. -* -* UPLO = 'L' or 'l' The lower triangular part of A is -* supplied in AP. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* ALPHA - COMPLEX . -* On entry, ALPHA specifies the scalar alpha. -* Unchanged on exit. -* -* X - COMPLEX array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element vector x. -* Unchanged on exit. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* Y - COMPLEX array of dimension at least -* ( 1 + ( n - 1 )*abs( INCY ) ). -* Before entry, the incremented array Y must contain the n -* element vector y. -* Unchanged on exit. -* -* INCY - INTEGER. -* On entry, INCY specifies the increment for the elements of -* Y. INCY must not be zero. -* Unchanged on exit. -* -* AP - COMPLEX array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular part of the hermitian matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) -* and a( 2, 2 ) respectively, and so on. On exit, the array -* AP is overwritten by the upper triangular part of the -* updated matrix. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular part of the hermitian matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) -* and a( 3, 1 ) respectively, and so on. On exit, the array -* AP is overwritten by the lower triangular part of the -* updated matrix. -* Note that the imaginary parts of the diagonal elements need -* not be set, they are assumed to be zero, and on exit they -* are set to zero. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - COMPLEX ZERO - PARAMETER (ZERO= (0.0E+0,0.0E+0)) -* .. -* .. Local Scalars .. - COMPLEX TEMP1,TEMP2 - INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC CONJG,REAL -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (N.LT.0) THEN - INFO = 2 - ELSE IF (INCX.EQ.0) THEN - INFO = 5 - ELSE IF (INCY.EQ.0) THEN - INFO = 7 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('CHPR2 ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN -* -* Set up the start points in X and Y if the increments are not both -* unity. -* - IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN - IF (INCX.GT.0) THEN - KX = 1 - ELSE - KX = 1 - (N-1)*INCX - END IF - IF (INCY.GT.0) THEN - KY = 1 - ELSE - KY = 1 - (N-1)*INCY - END IF - JX = KX - JY = KY - END IF -* -* Start the operations. In this version the elements of the array AP -* are accessed sequentially with one pass through AP. -* - KK = 1 - IF (LSAME(UPLO,'U')) THEN -* -* Form A when upper triangle is stored in AP. -* - IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN - DO 20 J = 1,N - IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN - TEMP1 = ALPHA*CONJG(Y(J)) - TEMP2 = CONJG(ALPHA*X(J)) - K = KK - DO 10 I = 1,J - 1 - AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2 - K = K + 1 - 10 CONTINUE - AP(KK+J-1) = REAL(AP(KK+J-1)) + - + REAL(X(J)*TEMP1+Y(J)*TEMP2) - ELSE - AP(KK+J-1) = REAL(AP(KK+J-1)) - END IF - KK = KK + J - 20 CONTINUE - ELSE - DO 40 J = 1,N - IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN - TEMP1 = ALPHA*CONJG(Y(JY)) - TEMP2 = CONJG(ALPHA*X(JX)) - IX = KX - IY = KY - DO 30 K = KK,KK + J - 2 - AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2 - IX = IX + INCX - IY = IY + INCY - 30 CONTINUE - AP(KK+J-1) = REAL(AP(KK+J-1)) + - + REAL(X(JX)*TEMP1+Y(JY)*TEMP2) - ELSE - AP(KK+J-1) = REAL(AP(KK+J-1)) - END IF - JX = JX + INCX - JY = JY + INCY - KK = KK + J - 40 CONTINUE - END IF - ELSE -* -* Form A when lower triangle is stored in AP. -* - IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN - DO 60 J = 1,N - IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN - TEMP1 = ALPHA*CONJG(Y(J)) - TEMP2 = CONJG(ALPHA*X(J)) - AP(KK) = REAL(AP(KK)) + - + REAL(X(J)*TEMP1+Y(J)*TEMP2) - K = KK + 1 - DO 50 I = J + 1,N - AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2 - K = K + 1 - 50 CONTINUE - ELSE - AP(KK) = REAL(AP(KK)) - END IF - KK = KK + N - J + 1 - 60 CONTINUE - ELSE - DO 80 J = 1,N - IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN - TEMP1 = ALPHA*CONJG(Y(JY)) - TEMP2 = CONJG(ALPHA*X(JX)) - AP(KK) = REAL(AP(KK)) + - + REAL(X(JX)*TEMP1+Y(JY)*TEMP2) - IX = JX - IY = JY - DO 70 K = KK + 1,KK + N - J - IX = IX + INCX - IY = IY + INCY - AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2 - 70 CONTINUE - ELSE - AP(KK) = REAL(AP(KK)) - END IF - JX = JX + INCX - JY = JY + INCY - KK = KK + N - J + 1 - 80 CONTINUE - END IF - END IF -* - RETURN -* -* End of CHPR2 . -* - END diff --git a/blas/common.h b/blas/common.h index b598c4e45..2bf642c6b 100644 --- a/blas/common.h +++ b/blas/common.h @@ -10,6 +10,9 @@ #ifndef EIGEN_BLAS_COMMON_H #define EIGEN_BLAS_COMMON_H +#include <Eigen/Core> +#include <Eigen/Jacobi> + #include <iostream> #include <complex> @@ -48,7 +51,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) @@ -68,12 +71,14 @@ inline bool check_uplo(const char* uplo) return UPLO(*uplo)!=0xff; } -#include <Eigen/Core> -#include <Eigen/Jacobi> - namespace Eigen { #include "BandTriangularSolver.h" +#include "GeneralRank1Update.h" +#include "PackedSelfadjointProduct.h" +#include "PackedTriangularMatrixVector.h" +#include "PackedTriangularSolverVector.h" +#include "Rank2Update.h" } using namespace Eigen; diff --git a/blas/ctpmv.f b/blas/ctpmv.f deleted file mode 100644 index b63742ccb..000000000 --- a/blas/ctpmv.f +++ /dev/null @@ -1,329 +0,0 @@ - SUBROUTINE CTPMV(UPLO,TRANS,DIAG,N,AP,X,INCX) -* .. Scalar Arguments .. - INTEGER INCX,N - CHARACTER DIAG,TRANS,UPLO -* .. -* .. Array Arguments .. - COMPLEX AP(*),X(*) -* .. -* -* Purpose -* ======= -* -* CTPMV performs one of the matrix-vector operations -* -* x := A*x, or x := A'*x, or x := conjg( A' )*x, -* -* where x is an n element vector and A is an n by n unit, or non-unit, -* upper or lower triangular matrix, supplied in packed form. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the matrix is an upper or -* lower triangular matrix as follows: -* -* UPLO = 'U' or 'u' A is an upper triangular matrix. -* -* UPLO = 'L' or 'l' A is a lower triangular matrix. -* -* Unchanged on exit. -* -* TRANS - CHARACTER*1. -* On entry, TRANS specifies the operation to be performed as -* follows: -* -* TRANS = 'N' or 'n' x := A*x. -* -* TRANS = 'T' or 't' x := A'*x. -* -* TRANS = 'C' or 'c' x := conjg( A' )*x. -* -* Unchanged on exit. -* -* DIAG - CHARACTER*1. -* On entry, DIAG specifies whether or not A is unit -* triangular as follows: -* -* DIAG = 'U' or 'u' A is assumed to be unit triangular. -* -* DIAG = 'N' or 'n' A is not assumed to be unit -* triangular. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* AP - COMPLEX array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) -* respectively, and so on. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) -* respectively, and so on. -* Note that when DIAG = 'U' or 'u', the diagonal elements of -* A are not referenced, but are assumed to be unity. -* Unchanged on exit. -* -* X - COMPLEX array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element vector x. On exit, X is overwritten with the -* tranformed vector x. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - COMPLEX ZERO - PARAMETER (ZERO= (0.0E+0,0.0E+0)) -* .. -* .. Local Scalars .. - COMPLEX TEMP - INTEGER I,INFO,IX,J,JX,K,KK,KX - LOGICAL NOCONJ,NOUNIT -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC CONJG -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND. - + .NOT.LSAME(TRANS,'C')) THEN - INFO = 2 - ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN - INFO = 3 - ELSE IF (N.LT.0) THEN - INFO = 4 - ELSE IF (INCX.EQ.0) THEN - INFO = 7 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('CTPMV ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF (N.EQ.0) RETURN -* - NOCONJ = LSAME(TRANS,'T') - NOUNIT = LSAME(DIAG,'N') -* -* Set up the start point in X if the increment is not unity. This -* will be ( N - 1 )*INCX too small for descending loops. -* - IF (INCX.LE.0) THEN - KX = 1 - (N-1)*INCX - ELSE IF (INCX.NE.1) THEN - KX = 1 - END IF -* -* Start the operations. In this version the elements of AP are -* accessed sequentially with one pass through AP. -* - IF (LSAME(TRANS,'N')) THEN -* -* Form x:= A*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = 1 - IF (INCX.EQ.1) THEN - DO 20 J = 1,N - IF (X(J).NE.ZERO) THEN - TEMP = X(J) - K = KK - DO 10 I = 1,J - 1 - X(I) = X(I) + TEMP*AP(K) - K = K + 1 - 10 CONTINUE - IF (NOUNIT) X(J) = X(J)*AP(KK+J-1) - END IF - KK = KK + J - 20 CONTINUE - ELSE - JX = KX - DO 40 J = 1,N - IF (X(JX).NE.ZERO) THEN - TEMP = X(JX) - IX = KX - DO 30 K = KK,KK + J - 2 - X(IX) = X(IX) + TEMP*AP(K) - IX = IX + INCX - 30 CONTINUE - IF (NOUNIT) X(JX) = X(JX)*AP(KK+J-1) - END IF - JX = JX + INCX - KK = KK + J - 40 CONTINUE - END IF - ELSE - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 60 J = N,1,-1 - IF (X(J).NE.ZERO) THEN - TEMP = X(J) - K = KK - DO 50 I = N,J + 1,-1 - X(I) = X(I) + TEMP*AP(K) - K = K - 1 - 50 CONTINUE - IF (NOUNIT) X(J) = X(J)*AP(KK-N+J) - END IF - KK = KK - (N-J+1) - 60 CONTINUE - ELSE - KX = KX + (N-1)*INCX - JX = KX - DO 80 J = N,1,-1 - IF (X(JX).NE.ZERO) THEN - TEMP = X(JX) - IX = KX - DO 70 K = KK,KK - (N- (J+1)),-1 - X(IX) = X(IX) + TEMP*AP(K) - IX = IX - INCX - 70 CONTINUE - IF (NOUNIT) X(JX) = X(JX)*AP(KK-N+J) - END IF - JX = JX - INCX - KK = KK - (N-J+1) - 80 CONTINUE - END IF - END IF - ELSE -* -* Form x := A'*x or x := conjg( A' )*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 110 J = N,1,-1 - TEMP = X(J) - K = KK - 1 - IF (NOCONJ) THEN - IF (NOUNIT) TEMP = TEMP*AP(KK) - DO 90 I = J - 1,1,-1 - TEMP = TEMP + AP(K)*X(I) - K = K - 1 - 90 CONTINUE - ELSE - IF (NOUNIT) TEMP = TEMP*CONJG(AP(KK)) - DO 100 I = J - 1,1,-1 - TEMP = TEMP + CONJG(AP(K))*X(I) - K = K - 1 - 100 CONTINUE - END IF - X(J) = TEMP - KK = KK - J - 110 CONTINUE - ELSE - JX = KX + (N-1)*INCX - DO 140 J = N,1,-1 - TEMP = X(JX) - IX = JX - IF (NOCONJ) THEN - IF (NOUNIT) TEMP = TEMP*AP(KK) - DO 120 K = KK - 1,KK - J + 1,-1 - IX = IX - INCX - TEMP = TEMP + AP(K)*X(IX) - 120 CONTINUE - ELSE - IF (NOUNIT) TEMP = TEMP*CONJG(AP(KK)) - DO 130 K = KK - 1,KK - J + 1,-1 - IX = IX - INCX - TEMP = TEMP + CONJG(AP(K))*X(IX) - 130 CONTINUE - END IF - X(JX) = TEMP - JX = JX - INCX - KK = KK - J - 140 CONTINUE - END IF - ELSE - KK = 1 - IF (INCX.EQ.1) THEN - DO 170 J = 1,N - TEMP = X(J) - K = KK + 1 - IF (NOCONJ) THEN - IF (NOUNIT) TEMP = TEMP*AP(KK) - DO 150 I = J + 1,N - TEMP = TEMP + AP(K)*X(I) - K = K + 1 - 150 CONTINUE - ELSE - IF (NOUNIT) TEMP = TEMP*CONJG(AP(KK)) - DO 160 I = J + 1,N - TEMP = TEMP + CONJG(AP(K))*X(I) - K = K + 1 - 160 CONTINUE - END IF - X(J) = TEMP - KK = KK + (N-J+1) - 170 CONTINUE - ELSE - JX = KX - DO 200 J = 1,N - TEMP = X(JX) - IX = JX - IF (NOCONJ) THEN - IF (NOUNIT) TEMP = TEMP*AP(KK) - DO 180 K = KK + 1,KK + N - J - IX = IX + INCX - TEMP = TEMP + AP(K)*X(IX) - 180 CONTINUE - ELSE - IF (NOUNIT) TEMP = TEMP*CONJG(AP(KK)) - DO 190 K = KK + 1,KK + N - J - IX = IX + INCX - TEMP = TEMP + CONJG(AP(K))*X(IX) - 190 CONTINUE - END IF - X(JX) = TEMP - JX = JX + INCX - KK = KK + (N-J+1) - 200 CONTINUE - END IF - END IF - END IF -* - RETURN -* -* End of CTPMV . -* - END diff --git a/blas/ctpsv.f b/blas/ctpsv.f deleted file mode 100644 index 1804797ea..000000000 --- a/blas/ctpsv.f +++ /dev/null @@ -1,332 +0,0 @@ - SUBROUTINE CTPSV(UPLO,TRANS,DIAG,N,AP,X,INCX) -* .. Scalar Arguments .. - INTEGER INCX,N - CHARACTER DIAG,TRANS,UPLO -* .. -* .. Array Arguments .. - COMPLEX AP(*),X(*) -* .. -* -* Purpose -* ======= -* -* CTPSV solves one of the systems of equations -* -* A*x = b, or A'*x = b, or conjg( A' )*x = b, -* -* where b and x are n element vectors and A is an n by n unit, or -* non-unit, upper or lower triangular matrix, supplied in packed form. -* -* No test for singularity or near-singularity is included in this -* routine. Such tests must be performed before calling this routine. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the matrix is an upper or -* lower triangular matrix as follows: -* -* UPLO = 'U' or 'u' A is an upper triangular matrix. -* -* UPLO = 'L' or 'l' A is a lower triangular matrix. -* -* Unchanged on exit. -* -* TRANS - CHARACTER*1. -* On entry, TRANS specifies the equations to be solved as -* follows: -* -* TRANS = 'N' or 'n' A*x = b. -* -* TRANS = 'T' or 't' A'*x = b. -* -* TRANS = 'C' or 'c' conjg( A' )*x = b. -* -* Unchanged on exit. -* -* DIAG - CHARACTER*1. -* On entry, DIAG specifies whether or not A is unit -* triangular as follows: -* -* DIAG = 'U' or 'u' A is assumed to be unit triangular. -* -* DIAG = 'N' or 'n' A is not assumed to be unit -* triangular. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* AP - COMPLEX array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) -* respectively, and so on. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) -* respectively, and so on. -* Note that when DIAG = 'U' or 'u', the diagonal elements of -* A are not referenced, but are assumed to be unity. -* Unchanged on exit. -* -* X - COMPLEX array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element right-hand side vector b. On exit, X is overwritten -* with the solution vector x. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - COMPLEX ZERO - PARAMETER (ZERO= (0.0E+0,0.0E+0)) -* .. -* .. Local Scalars .. - COMPLEX TEMP - INTEGER I,INFO,IX,J,JX,K,KK,KX - LOGICAL NOCONJ,NOUNIT -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC CONJG -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND. - + .NOT.LSAME(TRANS,'C')) THEN - INFO = 2 - ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN - INFO = 3 - ELSE IF (N.LT.0) THEN - INFO = 4 - ELSE IF (INCX.EQ.0) THEN - INFO = 7 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('CTPSV ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF (N.EQ.0) RETURN -* - NOCONJ = LSAME(TRANS,'T') - NOUNIT = LSAME(DIAG,'N') -* -* Set up the start point in X if the increment is not unity. This -* will be ( N - 1 )*INCX too small for descending loops. -* - IF (INCX.LE.0) THEN - KX = 1 - (N-1)*INCX - ELSE IF (INCX.NE.1) THEN - KX = 1 - END IF -* -* Start the operations. In this version the elements of AP are -* accessed sequentially with one pass through AP. -* - IF (LSAME(TRANS,'N')) THEN -* -* Form x := inv( A )*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 20 J = N,1,-1 - IF (X(J).NE.ZERO) THEN - IF (NOUNIT) X(J) = X(J)/AP(KK) - TEMP = X(J) - K = KK - 1 - DO 10 I = J - 1,1,-1 - X(I) = X(I) - TEMP*AP(K) - K = K - 1 - 10 CONTINUE - END IF - KK = KK - J - 20 CONTINUE - ELSE - JX = KX + (N-1)*INCX - DO 40 J = N,1,-1 - IF (X(JX).NE.ZERO) THEN - IF (NOUNIT) X(JX) = X(JX)/AP(KK) - TEMP = X(JX) - IX = JX - DO 30 K = KK - 1,KK - J + 1,-1 - IX = IX - INCX - X(IX) = X(IX) - TEMP*AP(K) - 30 CONTINUE - END IF - JX = JX - INCX - KK = KK - J - 40 CONTINUE - END IF - ELSE - KK = 1 - IF (INCX.EQ.1) THEN - DO 60 J = 1,N - IF (X(J).NE.ZERO) THEN - IF (NOUNIT) X(J) = X(J)/AP(KK) - TEMP = X(J) - K = KK + 1 - DO 50 I = J + 1,N - X(I) = X(I) - TEMP*AP(K) - K = K + 1 - 50 CONTINUE - END IF - KK = KK + (N-J+1) - 60 CONTINUE - ELSE - JX = KX - DO 80 J = 1,N - IF (X(JX).NE.ZERO) THEN - IF (NOUNIT) X(JX) = X(JX)/AP(KK) - TEMP = X(JX) - IX = JX - DO 70 K = KK + 1,KK + N - J - IX = IX + INCX - X(IX) = X(IX) - TEMP*AP(K) - 70 CONTINUE - END IF - JX = JX + INCX - KK = KK + (N-J+1) - 80 CONTINUE - END IF - END IF - ELSE -* -* Form x := inv( A' )*x or x := inv( conjg( A' ) )*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = 1 - IF (INCX.EQ.1) THEN - DO 110 J = 1,N - TEMP = X(J) - K = KK - IF (NOCONJ) THEN - DO 90 I = 1,J - 1 - TEMP = TEMP - AP(K)*X(I) - K = K + 1 - 90 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK+J-1) - ELSE - DO 100 I = 1,J - 1 - TEMP = TEMP - CONJG(AP(K))*X(I) - K = K + 1 - 100 CONTINUE - IF (NOUNIT) TEMP = TEMP/CONJG(AP(KK+J-1)) - END IF - X(J) = TEMP - KK = KK + J - 110 CONTINUE - ELSE - JX = KX - DO 140 J = 1,N - TEMP = X(JX) - IX = KX - IF (NOCONJ) THEN - DO 120 K = KK,KK + J - 2 - TEMP = TEMP - AP(K)*X(IX) - IX = IX + INCX - 120 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK+J-1) - ELSE - DO 130 K = KK,KK + J - 2 - TEMP = TEMP - CONJG(AP(K))*X(IX) - IX = IX + INCX - 130 CONTINUE - IF (NOUNIT) TEMP = TEMP/CONJG(AP(KK+J-1)) - END IF - X(JX) = TEMP - JX = JX + INCX - KK = KK + J - 140 CONTINUE - END IF - ELSE - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 170 J = N,1,-1 - TEMP = X(J) - K = KK - IF (NOCONJ) THEN - DO 150 I = N,J + 1,-1 - TEMP = TEMP - AP(K)*X(I) - K = K - 1 - 150 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK-N+J) - ELSE - DO 160 I = N,J + 1,-1 - TEMP = TEMP - CONJG(AP(K))*X(I) - K = K - 1 - 160 CONTINUE - IF (NOUNIT) TEMP = TEMP/CONJG(AP(KK-N+J)) - END IF - X(J) = TEMP - KK = KK - (N-J+1) - 170 CONTINUE - ELSE - KX = KX + (N-1)*INCX - JX = KX - DO 200 J = N,1,-1 - TEMP = X(JX) - IX = KX - IF (NOCONJ) THEN - DO 180 K = KK,KK - (N- (J+1)),-1 - TEMP = TEMP - AP(K)*X(IX) - IX = IX - INCX - 180 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK-N+J) - ELSE - DO 190 K = KK,KK - (N- (J+1)),-1 - TEMP = TEMP - CONJG(AP(K))*X(IX) - IX = IX - INCX - 190 CONTINUE - IF (NOUNIT) TEMP = TEMP/CONJG(AP(KK-N+J)) - END IF - X(JX) = TEMP - JX = JX - INCX - KK = KK - (N-J+1) - 200 CONTINUE - END IF - END IF - END IF -* - RETURN -* -* End of CTPSV . -* - END 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 <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; +} + diff --git a/blas/dspr.f b/blas/dspr.f deleted file mode 100644 index 538e4f76b..000000000 --- a/blas/dspr.f +++ /dev/null @@ -1,202 +0,0 @@ - SUBROUTINE DSPR(UPLO,N,ALPHA,X,INCX,AP) -* .. Scalar Arguments .. - DOUBLE PRECISION ALPHA - INTEGER INCX,N - CHARACTER UPLO -* .. -* .. Array Arguments .. - DOUBLE PRECISION AP(*),X(*) -* .. -* -* Purpose -* ======= -* -* DSPR performs the symmetric rank 1 operation -* -* A := alpha*x*x' + A, -* -* where alpha is a real scalar, x is an n element vector and A is an -* n by n symmetric matrix, supplied in packed form. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the upper or lower -* triangular part of the matrix A is supplied in the packed -* array AP as follows: -* -* UPLO = 'U' or 'u' The upper triangular part of A is -* supplied in AP. -* -* UPLO = 'L' or 'l' The lower triangular part of A is -* supplied in AP. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* ALPHA - DOUBLE PRECISION. -* On entry, ALPHA specifies the scalar alpha. -* Unchanged on exit. -* -* X - DOUBLE PRECISION array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element vector x. -* Unchanged on exit. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* AP - DOUBLE PRECISION array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular part of the symmetric matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) -* and a( 2, 2 ) respectively, and so on. On exit, the array -* AP is overwritten by the upper triangular part of the -* updated matrix. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular part of the symmetric matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) -* and a( 3, 1 ) respectively, and so on. On exit, the array -* AP is overwritten by the lower triangular part of the -* updated matrix. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - DOUBLE PRECISION ZERO - PARAMETER (ZERO=0.0D+0) -* .. -* .. Local Scalars .. - DOUBLE PRECISION TEMP - INTEGER I,INFO,IX,J,JX,K,KK,KX -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (N.LT.0) THEN - INFO = 2 - ELSE IF (INCX.EQ.0) THEN - INFO = 5 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('DSPR ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN -* -* Set the start point in X if the increment is not unity. -* - IF (INCX.LE.0) THEN - KX = 1 - (N-1)*INCX - ELSE IF (INCX.NE.1) THEN - KX = 1 - END IF -* -* Start the operations. In this version the elements of the array AP -* are accessed sequentially with one pass through AP. -* - KK = 1 - IF (LSAME(UPLO,'U')) THEN -* -* Form A when upper triangle is stored in AP. -* - IF (INCX.EQ.1) THEN - DO 20 J = 1,N - IF (X(J).NE.ZERO) THEN - TEMP = ALPHA*X(J) - K = KK - DO 10 I = 1,J - AP(K) = AP(K) + X(I)*TEMP - K = K + 1 - 10 CONTINUE - END IF - KK = KK + J - 20 CONTINUE - ELSE - JX = KX - DO 40 J = 1,N - IF (X(JX).NE.ZERO) THEN - TEMP = ALPHA*X(JX) - IX = KX - DO 30 K = KK,KK + J - 1 - AP(K) = AP(K) + X(IX)*TEMP - IX = IX + INCX - 30 CONTINUE - END IF - JX = JX + INCX - KK = KK + J - 40 CONTINUE - END IF - ELSE -* -* Form A when lower triangle is stored in AP. -* - IF (INCX.EQ.1) THEN - DO 60 J = 1,N - IF (X(J).NE.ZERO) THEN - TEMP = ALPHA*X(J) - K = KK - DO 50 I = J,N - AP(K) = AP(K) + X(I)*TEMP - K = K + 1 - 50 CONTINUE - END IF - KK = KK + N - J + 1 - 60 CONTINUE - ELSE - JX = KX - DO 80 J = 1,N - IF (X(JX).NE.ZERO) THEN - TEMP = ALPHA*X(JX) - IX = JX - DO 70 K = KK,KK + N - J - AP(K) = AP(K) + X(IX)*TEMP - IX = IX + INCX - 70 CONTINUE - END IF - JX = JX + INCX - KK = KK + N - J + 1 - 80 CONTINUE - END IF - END IF -* - RETURN -* -* End of DSPR . -* - END diff --git a/blas/dspr2.f b/blas/dspr2.f deleted file mode 100644 index 6f6b54a8c..000000000 --- a/blas/dspr2.f +++ /dev/null @@ -1,233 +0,0 @@ - SUBROUTINE DSPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP) -* .. Scalar Arguments .. - DOUBLE PRECISION ALPHA - INTEGER INCX,INCY,N - CHARACTER UPLO -* .. -* .. Array Arguments .. - DOUBLE PRECISION AP(*),X(*),Y(*) -* .. -* -* Purpose -* ======= -* -* DSPR2 performs the symmetric rank 2 operation -* -* A := alpha*x*y' + alpha*y*x' + A, -* -* where alpha is a scalar, x and y are n element vectors and A is an -* n by n symmetric matrix, supplied in packed form. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the upper or lower -* triangular part of the matrix A is supplied in the packed -* array AP as follows: -* -* UPLO = 'U' or 'u' The upper triangular part of A is -* supplied in AP. -* -* UPLO = 'L' or 'l' The lower triangular part of A is -* supplied in AP. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* ALPHA - DOUBLE PRECISION. -* On entry, ALPHA specifies the scalar alpha. -* Unchanged on exit. -* -* X - DOUBLE PRECISION array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element vector x. -* Unchanged on exit. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* Y - DOUBLE PRECISION array of dimension at least -* ( 1 + ( n - 1 )*abs( INCY ) ). -* Before entry, the incremented array Y must contain the n -* element vector y. -* Unchanged on exit. -* -* INCY - INTEGER. -* On entry, INCY specifies the increment for the elements of -* Y. INCY must not be zero. -* Unchanged on exit. -* -* AP - DOUBLE PRECISION array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular part of the symmetric matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) -* and a( 2, 2 ) respectively, and so on. On exit, the array -* AP is overwritten by the upper triangular part of the -* updated matrix. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular part of the symmetric matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) -* and a( 3, 1 ) respectively, and so on. On exit, the array -* AP is overwritten by the lower triangular part of the -* updated matrix. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - DOUBLE PRECISION ZERO - PARAMETER (ZERO=0.0D+0) -* .. -* .. Local Scalars .. - DOUBLE PRECISION TEMP1,TEMP2 - INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (N.LT.0) THEN - INFO = 2 - ELSE IF (INCX.EQ.0) THEN - INFO = 5 - ELSE IF (INCY.EQ.0) THEN - INFO = 7 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('DSPR2 ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN -* -* Set up the start points in X and Y if the increments are not both -* unity. -* - IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN - IF (INCX.GT.0) THEN - KX = 1 - ELSE - KX = 1 - (N-1)*INCX - END IF - IF (INCY.GT.0) THEN - KY = 1 - ELSE - KY = 1 - (N-1)*INCY - END IF - JX = KX - JY = KY - END IF -* -* Start the operations. In this version the elements of the array AP -* are accessed sequentially with one pass through AP. -* - KK = 1 - IF (LSAME(UPLO,'U')) THEN -* -* Form A when upper triangle is stored in AP. -* - IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN - DO 20 J = 1,N - IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN - TEMP1 = ALPHA*Y(J) - TEMP2 = ALPHA*X(J) - K = KK - DO 10 I = 1,J - AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2 - K = K + 1 - 10 CONTINUE - END IF - KK = KK + J - 20 CONTINUE - ELSE - DO 40 J = 1,N - IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN - TEMP1 = ALPHA*Y(JY) - TEMP2 = ALPHA*X(JX) - IX = KX - IY = KY - DO 30 K = KK,KK + J - 1 - AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2 - IX = IX + INCX - IY = IY + INCY - 30 CONTINUE - END IF - JX = JX + INCX - JY = JY + INCY - KK = KK + J - 40 CONTINUE - END IF - ELSE -* -* Form A when lower triangle is stored in AP. -* - IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN - DO 60 J = 1,N - IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN - TEMP1 = ALPHA*Y(J) - TEMP2 = ALPHA*X(J) - K = KK - DO 50 I = J,N - AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2 - K = K + 1 - 50 CONTINUE - END IF - KK = KK + N - J + 1 - 60 CONTINUE - ELSE - DO 80 J = 1,N - IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN - TEMP1 = ALPHA*Y(JY) - TEMP2 = ALPHA*X(JX) - IX = JX - IY = JY - DO 70 K = KK,KK + N - J - AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2 - IX = IX + INCX - IY = IY + INCY - 70 CONTINUE - END IF - JX = JX + INCX - JY = JY + INCY - KK = KK + N - J + 1 - 80 CONTINUE - END IF - END IF -* - RETURN -* -* End of DSPR2 . -* - END diff --git a/blas/dtpmv.f b/blas/dtpmv.f deleted file mode 100644 index c5bc112dc..000000000 --- a/blas/dtpmv.f +++ /dev/null @@ -1,293 +0,0 @@ - SUBROUTINE DTPMV(UPLO,TRANS,DIAG,N,AP,X,INCX) -* .. Scalar Arguments .. - INTEGER INCX,N - CHARACTER DIAG,TRANS,UPLO -* .. -* .. Array Arguments .. - DOUBLE PRECISION AP(*),X(*) -* .. -* -* Purpose -* ======= -* -* DTPMV performs one of the matrix-vector operations -* -* x := A*x, or x := A'*x, -* -* where x is an n element vector and A is an n by n unit, or non-unit, -* upper or lower triangular matrix, supplied in packed form. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the matrix is an upper or -* lower triangular matrix as follows: -* -* UPLO = 'U' or 'u' A is an upper triangular matrix. -* -* UPLO = 'L' or 'l' A is a lower triangular matrix. -* -* Unchanged on exit. -* -* TRANS - CHARACTER*1. -* On entry, TRANS specifies the operation to be performed as -* follows: -* -* TRANS = 'N' or 'n' x := A*x. -* -* TRANS = 'T' or 't' x := A'*x. -* -* TRANS = 'C' or 'c' x := A'*x. -* -* Unchanged on exit. -* -* DIAG - CHARACTER*1. -* On entry, DIAG specifies whether or not A is unit -* triangular as follows: -* -* DIAG = 'U' or 'u' A is assumed to be unit triangular. -* -* DIAG = 'N' or 'n' A is not assumed to be unit -* triangular. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* AP - DOUBLE PRECISION array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) -* respectively, and so on. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) -* respectively, and so on. -* Note that when DIAG = 'U' or 'u', the diagonal elements of -* A are not referenced, but are assumed to be unity. -* Unchanged on exit. -* -* X - DOUBLE PRECISION array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element vector x. On exit, X is overwritten with the -* tranformed vector x. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - DOUBLE PRECISION ZERO - PARAMETER (ZERO=0.0D+0) -* .. -* .. Local Scalars .. - DOUBLE PRECISION TEMP - INTEGER I,INFO,IX,J,JX,K,KK,KX - LOGICAL NOUNIT -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND. - + .NOT.LSAME(TRANS,'C')) THEN - INFO = 2 - ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN - INFO = 3 - ELSE IF (N.LT.0) THEN - INFO = 4 - ELSE IF (INCX.EQ.0) THEN - INFO = 7 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('DTPMV ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF (N.EQ.0) RETURN -* - NOUNIT = LSAME(DIAG,'N') -* -* Set up the start point in X if the increment is not unity. This -* will be ( N - 1 )*INCX too small for descending loops. -* - IF (INCX.LE.0) THEN - KX = 1 - (N-1)*INCX - ELSE IF (INCX.NE.1) THEN - KX = 1 - END IF -* -* Start the operations. In this version the elements of AP are -* accessed sequentially with one pass through AP. -* - IF (LSAME(TRANS,'N')) THEN -* -* Form x:= A*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = 1 - IF (INCX.EQ.1) THEN - DO 20 J = 1,N - IF (X(J).NE.ZERO) THEN - TEMP = X(J) - K = KK - DO 10 I = 1,J - 1 - X(I) = X(I) + TEMP*AP(K) - K = K + 1 - 10 CONTINUE - IF (NOUNIT) X(J) = X(J)*AP(KK+J-1) - END IF - KK = KK + J - 20 CONTINUE - ELSE - JX = KX - DO 40 J = 1,N - IF (X(JX).NE.ZERO) THEN - TEMP = X(JX) - IX = KX - DO 30 K = KK,KK + J - 2 - X(IX) = X(IX) + TEMP*AP(K) - IX = IX + INCX - 30 CONTINUE - IF (NOUNIT) X(JX) = X(JX)*AP(KK+J-1) - END IF - JX = JX + INCX - KK = KK + J - 40 CONTINUE - END IF - ELSE - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 60 J = N,1,-1 - IF (X(J).NE.ZERO) THEN - TEMP = X(J) - K = KK - DO 50 I = N,J + 1,-1 - X(I) = X(I) + TEMP*AP(K) - K = K - 1 - 50 CONTINUE - IF (NOUNIT) X(J) = X(J)*AP(KK-N+J) - END IF - KK = KK - (N-J+1) - 60 CONTINUE - ELSE - KX = KX + (N-1)*INCX - JX = KX - DO 80 J = N,1,-1 - IF (X(JX).NE.ZERO) THEN - TEMP = X(JX) - IX = KX - DO 70 K = KK,KK - (N- (J+1)),-1 - X(IX) = X(IX) + TEMP*AP(K) - IX = IX - INCX - 70 CONTINUE - IF (NOUNIT) X(JX) = X(JX)*AP(KK-N+J) - END IF - JX = JX - INCX - KK = KK - (N-J+1) - 80 CONTINUE - END IF - END IF - ELSE -* -* Form x := A'*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 100 J = N,1,-1 - TEMP = X(J) - IF (NOUNIT) TEMP = TEMP*AP(KK) - K = KK - 1 - DO 90 I = J - 1,1,-1 - TEMP = TEMP + AP(K)*X(I) - K = K - 1 - 90 CONTINUE - X(J) = TEMP - KK = KK - J - 100 CONTINUE - ELSE - JX = KX + (N-1)*INCX - DO 120 J = N,1,-1 - TEMP = X(JX) - IX = JX - IF (NOUNIT) TEMP = TEMP*AP(KK) - DO 110 K = KK - 1,KK - J + 1,-1 - IX = IX - INCX - TEMP = TEMP + AP(K)*X(IX) - 110 CONTINUE - X(JX) = TEMP - JX = JX - INCX - KK = KK - J - 120 CONTINUE - END IF - ELSE - KK = 1 - IF (INCX.EQ.1) THEN - DO 140 J = 1,N - TEMP = X(J) - IF (NOUNIT) TEMP = TEMP*AP(KK) - K = KK + 1 - DO 130 I = J + 1,N - TEMP = TEMP + AP(K)*X(I) - K = K + 1 - 130 CONTINUE - X(J) = TEMP - KK = KK + (N-J+1) - 140 CONTINUE - ELSE - JX = KX - DO 160 J = 1,N - TEMP = X(JX) - IX = JX - IF (NOUNIT) TEMP = TEMP*AP(KK) - DO 150 K = KK + 1,KK + N - J - IX = IX + INCX - TEMP = TEMP + AP(K)*X(IX) - 150 CONTINUE - X(JX) = TEMP - JX = JX + INCX - KK = KK + (N-J+1) - 160 CONTINUE - END IF - END IF - END IF -* - RETURN -* -* End of DTPMV . -* - END diff --git a/blas/dtpsv.f b/blas/dtpsv.f deleted file mode 100644 index c7e58d32f..000000000 --- a/blas/dtpsv.f +++ /dev/null @@ -1,296 +0,0 @@ - SUBROUTINE DTPSV(UPLO,TRANS,DIAG,N,AP,X,INCX) -* .. Scalar Arguments .. - INTEGER INCX,N - CHARACTER DIAG,TRANS,UPLO -* .. -* .. Array Arguments .. - DOUBLE PRECISION AP(*),X(*) -* .. -* -* Purpose -* ======= -* -* DTPSV solves one of the systems of equations -* -* A*x = b, or A'*x = b, -* -* where b and x are n element vectors and A is an n by n unit, or -* non-unit, upper or lower triangular matrix, supplied in packed form. -* -* No test for singularity or near-singularity is included in this -* routine. Such tests must be performed before calling this routine. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the matrix is an upper or -* lower triangular matrix as follows: -* -* UPLO = 'U' or 'u' A is an upper triangular matrix. -* -* UPLO = 'L' or 'l' A is a lower triangular matrix. -* -* Unchanged on exit. -* -* TRANS - CHARACTER*1. -* On entry, TRANS specifies the equations to be solved as -* follows: -* -* TRANS = 'N' or 'n' A*x = b. -* -* TRANS = 'T' or 't' A'*x = b. -* -* TRANS = 'C' or 'c' A'*x = b. -* -* Unchanged on exit. -* -* DIAG - CHARACTER*1. -* On entry, DIAG specifies whether or not A is unit -* triangular as follows: -* -* DIAG = 'U' or 'u' A is assumed to be unit triangular. -* -* DIAG = 'N' or 'n' A is not assumed to be unit -* triangular. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* AP - DOUBLE PRECISION array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) -* respectively, and so on. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) -* respectively, and so on. -* Note that when DIAG = 'U' or 'u', the diagonal elements of -* A are not referenced, but are assumed to be unity. -* Unchanged on exit. -* -* X - DOUBLE PRECISION array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element right-hand side vector b. On exit, X is overwritten -* with the solution vector x. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - DOUBLE PRECISION ZERO - PARAMETER (ZERO=0.0D+0) -* .. -* .. Local Scalars .. - DOUBLE PRECISION TEMP - INTEGER I,INFO,IX,J,JX,K,KK,KX - LOGICAL NOUNIT -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND. - + .NOT.LSAME(TRANS,'C')) THEN - INFO = 2 - ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN - INFO = 3 - ELSE IF (N.LT.0) THEN - INFO = 4 - ELSE IF (INCX.EQ.0) THEN - INFO = 7 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('DTPSV ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF (N.EQ.0) RETURN -* - NOUNIT = LSAME(DIAG,'N') -* -* Set up the start point in X if the increment is not unity. This -* will be ( N - 1 )*INCX too small for descending loops. -* - IF (INCX.LE.0) THEN - KX = 1 - (N-1)*INCX - ELSE IF (INCX.NE.1) THEN - KX = 1 - END IF -* -* Start the operations. In this version the elements of AP are -* accessed sequentially with one pass through AP. -* - IF (LSAME(TRANS,'N')) THEN -* -* Form x := inv( A )*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 20 J = N,1,-1 - IF (X(J).NE.ZERO) THEN - IF (NOUNIT) X(J) = X(J)/AP(KK) - TEMP = X(J) - K = KK - 1 - DO 10 I = J - 1,1,-1 - X(I) = X(I) - TEMP*AP(K) - K = K - 1 - 10 CONTINUE - END IF - KK = KK - J - 20 CONTINUE - ELSE - JX = KX + (N-1)*INCX - DO 40 J = N,1,-1 - IF (X(JX).NE.ZERO) THEN - IF (NOUNIT) X(JX) = X(JX)/AP(KK) - TEMP = X(JX) - IX = JX - DO 30 K = KK - 1,KK - J + 1,-1 - IX = IX - INCX - X(IX) = X(IX) - TEMP*AP(K) - 30 CONTINUE - END IF - JX = JX - INCX - KK = KK - J - 40 CONTINUE - END IF - ELSE - KK = 1 - IF (INCX.EQ.1) THEN - DO 60 J = 1,N - IF (X(J).NE.ZERO) THEN - IF (NOUNIT) X(J) = X(J)/AP(KK) - TEMP = X(J) - K = KK + 1 - DO 50 I = J + 1,N - X(I) = X(I) - TEMP*AP(K) - K = K + 1 - 50 CONTINUE - END IF - KK = KK + (N-J+1) - 60 CONTINUE - ELSE - JX = KX - DO 80 J = 1,N - IF (X(JX).NE.ZERO) THEN - IF (NOUNIT) X(JX) = X(JX)/AP(KK) - TEMP = X(JX) - IX = JX - DO 70 K = KK + 1,KK + N - J - IX = IX + INCX - X(IX) = X(IX) - TEMP*AP(K) - 70 CONTINUE - END IF - JX = JX + INCX - KK = KK + (N-J+1) - 80 CONTINUE - END IF - END IF - ELSE -* -* Form x := inv( A' )*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = 1 - IF (INCX.EQ.1) THEN - DO 100 J = 1,N - TEMP = X(J) - K = KK - DO 90 I = 1,J - 1 - TEMP = TEMP - AP(K)*X(I) - K = K + 1 - 90 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK+J-1) - X(J) = TEMP - KK = KK + J - 100 CONTINUE - ELSE - JX = KX - DO 120 J = 1,N - TEMP = X(JX) - IX = KX - DO 110 K = KK,KK + J - 2 - TEMP = TEMP - AP(K)*X(IX) - IX = IX + INCX - 110 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK+J-1) - X(JX) = TEMP - JX = JX + INCX - KK = KK + J - 120 CONTINUE - END IF - ELSE - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 140 J = N,1,-1 - TEMP = X(J) - K = KK - DO 130 I = N,J + 1,-1 - TEMP = TEMP - AP(K)*X(I) - K = K - 1 - 130 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK-N+J) - X(J) = TEMP - KK = KK - (N-J+1) - 140 CONTINUE - ELSE - KX = KX + (N-1)*INCX - JX = KX - DO 160 J = N,1,-1 - TEMP = X(JX) - IX = KX - DO 150 K = KK,KK - (N- (J+1)),-1 - TEMP = TEMP - AP(K)*X(IX) - IX = IX - INCX - 150 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK-N+J) - X(JX) = TEMP - JX = JX - INCX - KK = KK - (N-J+1) - 160 CONTINUE - END IF - END IF - END IF -* - RETURN -* -* End of DTPSV . -* - END diff --git a/blas/level1_cplx_impl.h b/blas/level1_cplx_impl.h index 8d6b92829..283b9f827 100644 --- a/blas/level1_cplx_impl.h +++ b/blas/level1_cplx_impl.h @@ -12,7 +12,7 @@ struct scalar_norm1_op { typedef RealScalar result_type; EIGEN_EMPTY_STRUCT_CTOR(scalar_norm1_op) - inline RealScalar operator() (const Scalar& a) const { return internal::norm1(a); } + inline RealScalar operator() (const Scalar& a) const { return numext::norm1(a); } }; namespace Eigen { namespace internal { diff --git a/blas/level1_impl.h b/blas/level1_impl.h index 95ea220af..b08c2f6be 100644 --- a/blas/level1_impl.h +++ b/blas/level1_impl.h @@ -75,6 +75,9 @@ int EIGEN_CAT(EIGEN_CAT(i,SCALAR_SUFFIX),amin_)(int *n, RealScalar *px, int *inc int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar *ps) { + using std::sqrt; + using std::abs; + Scalar& a = *reinterpret_cast<Scalar*>(pa); Scalar& b = *reinterpret_cast<Scalar*>(pb); RealScalar* c = pc; @@ -82,8 +85,8 @@ int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealSc #if !ISCOMPLEX Scalar r,z; - Scalar aa = internal::abs(a); - Scalar ab = internal::abs(b); + Scalar aa = abs(a); + Scalar ab = abs(b); if((aa+ab)==Scalar(0)) { *c = 1; @@ -93,7 +96,7 @@ int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealSc } else { - r = internal::sqrt(a*a + b*b); + r = sqrt(a*a + b*b); Scalar amax = aa>ab ? a : b; r = amax>0 ? r : -r; *c = a/r; @@ -108,7 +111,7 @@ int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealSc #else Scalar alpha; RealScalar norm,scale; - if(internal::abs(a)==RealScalar(0)) + if(abs(a)==RealScalar(0)) { *c = RealScalar(0); *s = Scalar(1); @@ -116,11 +119,11 @@ int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealSc } else { - scale = internal::abs(a) + internal::abs(b); - norm = scale*internal::sqrt((internal::abs2(a/scale))+ (internal::abs2(b/scale))); - alpha = a/internal::abs(a); - *c = internal::abs(a)/norm; - *s = alpha*internal::conj(b)/norm; + scale = abs(a) + abs(b); + norm = scale*sqrt((numext::abs2(a/scale)) + (numext::abs2(b/scale))); + alpha = a/abs(a); + *c = abs(a)/norm; + *s = alpha*numext::conj(b)/norm; a = alpha*norm; } #endif diff --git a/blas/level2_cplx_impl.h b/blas/level2_cplx_impl.h index 7878f2a16..b850b6cd1 100644 --- a/blas/level2_cplx_impl.h +++ b/blas/level2_cplx_impl.h @@ -18,6 +18,21 @@ */ int EIGEN_BLAS_FUNC(hemv)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy) { + typedef void (*functype)(int, const Scalar*, int, const Scalar*, int, 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::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run); + func[LO] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run); + + init = true; + } + Scalar* a = reinterpret_cast<Scalar*>(pa); Scalar* x = reinterpret_cast<Scalar*>(px); Scalar* y = reinterpret_cast<Scalar*>(py); @@ -48,9 +63,11 @@ int EIGEN_BLAS_FUNC(hemv)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa if(alpha!=Scalar(0)) { - // TODO performs a direct call to the underlying implementation function - if(UPLO(*uplo)==UP) vector(actual_y,*n).noalias() += matrix(a,*n,*n,*lda).selfadjointView<Upper>() * (alpha * vector(actual_x,*n)); - else if(UPLO(*uplo)==LO) vector(actual_y,*n).noalias() += matrix(a,*n,*n,*lda).selfadjointView<Lower>() * (alpha * vector(actual_x,*n)); + int code = UPLO(*uplo); + if(code>=2 || func[code]==0) + return 0; + + func[code](*n, a, *lda, actual_x, 1, actual_y, alpha); } if(actual_x!=x) delete[] actual_x; @@ -91,10 +108,49 @@ int EIGEN_BLAS_FUNC(hemv)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa * where alpha is a real scalar, x is an n element vector and A is an * n by n hermitian matrix, supplied in packed form. */ -// int EIGEN_BLAS_FUNC(hpr)(char *uplo, int *n, RealScalar *alpha, RealScalar *x, int *incx, RealScalar *ap) -// { -// return 1; -// } +int EIGEN_BLAS_FUNC(hpr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pap) +{ + typedef void (*functype)(int, Scalar*, const Scalar*, RealScalar); + static functype func[2]; + + static bool init = false; + if(!init) + { + for(int k=0; k<2; ++k) + func[k] = 0; + + func[UP] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run); + func[LO] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run); + + init = true; + } + + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar* ap = reinterpret_cast<Scalar*>(pap); + RealScalar alpha = *palpha; + + int info = 0; + if(UPLO(*uplo)==INVALID) info = 1; + else if(*n<0) info = 2; + else if(*incx==0) info = 5; + if(info) + return xerbla_(SCALAR_SUFFIX_UP"HPR ",&info,6); + + if(alpha==Scalar(0)) + return 1; + + Scalar* x_cpy = get_compact_vector(x, *n, *incx); + + int code = UPLO(*uplo); + if(code>=2 || func[code]==0) + return 0; + + func[code](*n, ap, x_cpy, alpha); + + if(x_cpy!=x) delete[] x_cpy; + + return 1; +} /** ZHPR2 performs the hermitian rank 2 operation * @@ -103,10 +159,53 @@ int EIGEN_BLAS_FUNC(hemv)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa * where alpha is a scalar, x and y are n element vectors and A is an * n by n hermitian matrix, supplied in packed form. */ -// int EIGEN_BLAS_FUNC(hpr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *x, int *incx, RealScalar *y, int *incy, RealScalar *ap) -// { -// return 1; -// } +int EIGEN_BLAS_FUNC(hpr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pap) +{ + typedef void (*functype)(int, Scalar*, 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::packed_rank2_update_selector<Scalar,int,Upper>::run); + func[LO] = (internal::packed_rank2_update_selector<Scalar,int,Lower>::run); + + init = true; + } + + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar* y = reinterpret_cast<Scalar*>(py); + Scalar* ap = reinterpret_cast<Scalar*>(pap); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + + int info = 0; + if(UPLO(*uplo)==INVALID) info = 1; + else if(*n<0) info = 2; + else if(*incx==0) info = 5; + else if(*incy==0) info = 7; + if(info) + return xerbla_(SCALAR_SUFFIX_UP"HPR2 ",&info,6); + + if(alpha==Scalar(0)) + return 1; + + 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; + + func[code](*n, ap, x_cpy, y_cpy, alpha); + + if(x_cpy!=x) delete[] x_cpy; + if(y_cpy!=y) delete[] y_cpy; + + return 1; +} /** ZHER performs the hermitian rank 1 operation * @@ -117,6 +216,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*, const Scalar*, const 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 +248,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, x_cpy, alpha); matrix(a,*n,*n,*lda).diagonal().imag().setZero(); @@ -161,6 +270,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 +305,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 +348,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,ColMajor,false,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 +385,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,ColMajor,false,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; diff --git a/blas/level2_impl.h b/blas/level2_impl.h index 7099cf96d..5f3941975 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; @@ -127,7 +130,7 @@ int EIGEN_BLAS_FUNC(trsv)(char *uplo, char *opa, char *diag, int *n, RealScalar int EIGEN_BLAS_FUNC(trmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb) { - typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar); + typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, const Scalar&); static functype func[16]; static bool init = false; @@ -184,7 +187,7 @@ int EIGEN_BLAS_FUNC(trmv)(char *uplo, char *opa, char *diag, int *n, RealScalar copy_back(res.data(),b,*n,*incb); if(actual_b!=b) delete[] actual_b; - return 0; + return 1; } /** GBMV performs one of the matrix-vector operations @@ -396,10 +399,66 @@ int EIGEN_BLAS_FUNC(tbsv)(char *uplo, char *op, char *diag, int *n, int *k, Real * where x is an n element vector and A is an n by n unit, or non-unit, * upper or lower triangular matrix, supplied in packed form. */ -// int EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *trans, char *diag, int *n, RealScalar *ap, RealScalar *x, int *incx) -// { -// return 1; -// } +int EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx) +{ + typedef void (*functype)(int, const Scalar*, const Scalar*, Scalar*, Scalar); + static functype func[16]; + + static bool init = false; + if(!init) + { + for(int k=0; k<16; ++k) + func[k] = 0; + + func[NOTR | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,ColMajor>::run); + func[TR | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,RowMajor>::run); + func[ADJ | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|0, Scalar,Conj, Scalar,false,RowMajor>::run); + + func[NOTR | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,ColMajor>::run); + func[TR | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,RowMajor>::run); + func[ADJ | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|0, Scalar,Conj, Scalar,false,RowMajor>::run); + + func[NOTR | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run); + func[TR | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run); + func[ADJ | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run); + + func[NOTR | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run); + func[TR | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run); + func[ADJ | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run); + + init = true; + } + + Scalar* ap = reinterpret_cast<Scalar*>(pap); + Scalar* x = reinterpret_cast<Scalar*>(px); + + int info = 0; + if(UPLO(*uplo)==INVALID) info = 1; + else if(OP(*opa)==INVALID) info = 2; + else if(DIAG(*diag)==INVALID) info = 3; + else if(*n<0) info = 4; + else if(*incx==0) info = 7; + if(info) + return xerbla_(SCALAR_SUFFIX_UP"TPMV ",&info,6); + + if(*n==0) + return 1; + + Scalar* actual_x = get_compact_vector(x,*n,*incx); + Matrix<Scalar,Dynamic,1> res(*n); + res.setZero(); + + int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); + if(code>=16 || func[code]==0) + return 0; + + func[code](*n, ap, actual_x, res.data(), Scalar(1)); + + copy_back(res.data(),x,*n,*incx); + if(actual_x!=x) delete[] actual_x; + + return 1; +} /** DTPSV solves one of the systems of equations * @@ -411,47 +470,55 @@ int EIGEN_BLAS_FUNC(tbsv)(char *uplo, char *op, char *diag, int *n, int *k, Real * No test for singularity or near-singularity is included in this * routine. Such tests must be performed before calling this routine. */ -// int EIGEN_BLAS_FUNC(tpsv)(char *uplo, char *trans, char *diag, int *n, RealScalar *ap, RealScalar *x, int *incx) -// { -// 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) +int EIGEN_BLAS_FUNC(tpsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx) { + typedef void (*functype)(int, const Scalar*, Scalar*); + static functype func[16]; + + static bool init = false; + if(!init) + { + for(int k=0; k<16; ++k) + func[k] = 0; + + func[NOTR | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,ColMajor>::run); + func[TR | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,RowMajor>::run); + func[ADJ | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, Conj, RowMajor>::run); + + func[NOTR | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,ColMajor>::run); + func[TR | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,RowMajor>::run); + func[ADJ | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, Conj, RowMajor>::run); + + func[NOTR | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run); + func[TR | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run); + func[ADJ | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run); + + func[NOTR | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run); + func[TR | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run); + func[ADJ | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run); + + init = true; + } + + Scalar* ap = reinterpret_cast<Scalar*>(pap); 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(UPLO(*uplo)==INVALID) info = 1; + else if(OP(*opa)==INVALID) info = 2; + else if(DIAG(*diag)==INVALID) info = 3; + else if(*n<0) info = 4; + else if(*incx==0) info = 7; if(info) - return xerbla_(SCALAR_SUFFIX_UP"GER ",&info,6); + return xerbla_(SCALAR_SUFFIX_UP"TPSV ",&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); + Scalar* actual_x = get_compact_vector(x,*n,*incx); - // TODO perform direct calls to underlying implementation - matrix(a,*m,*n,*lda) += alpha * vector(x_cpy,*m) * vector(y_cpy,*n).adjoint(); + int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); + func[code](*n, ap, actual_x); - if(x_cpy!=x) delete[] x_cpy; - if(y_cpy!=y) delete[] y_cpy; + if(actual_x!=x) delete[] copy_back(actual_x,x,*n,*incx); return 1; } - diff --git a/blas/level2_real_impl.h b/blas/level2_real_impl.h index cd8332973..8d56eaaa1 100644 --- a/blas/level2_real_impl.h +++ b/blas/level2_real_impl.h @@ -12,6 +12,21 @@ // y = alpha*A*x + beta*y int EIGEN_BLAS_FUNC(symv) (char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy) { + typedef void (*functype)(int, const Scalar*, int, const Scalar*, int, 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::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run); + func[LO] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run); + + init = true; + } + Scalar* a = reinterpret_cast<Scalar*>(pa); Scalar* x = reinterpret_cast<Scalar*>(px); Scalar* y = reinterpret_cast<Scalar*>(py); @@ -40,9 +55,11 @@ int EIGEN_BLAS_FUNC(symv) (char *uplo, int *n, RealScalar *palpha, RealScalar *p else vector(actual_y, *n) *= beta; } - // TODO performs a direct call to the underlying implementation function - if(UPLO(*uplo)==UP) vector(actual_y,*n).noalias() += matrix(a,*n,*n,*lda).selfadjointView<Upper>() * (alpha * vector(actual_x,*n)); - else if(UPLO(*uplo)==LO) vector(actual_y,*n).noalias() += matrix(a,*n,*n,*lda).selfadjointView<Lower>() * (alpha * vector(actual_x,*n)); + int code = UPLO(*uplo); + if(code>=2 || func[code]==0) + return 0; + + func[code](*n, a, *lda, actual_x, 1, actual_y, alpha); if(actual_x!=x) delete[] actual_x; if(actual_y!=y) delete[] copy_back(actual_y,y,*n,*incy); @@ -68,6 +85,20 @@ int EIGEN_BLAS_FUNC(syr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, // init = true; // } + typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, const 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 +117,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, x_cpy, alpha); if(x_cpy!=x) delete[] x_cpy; @@ -121,6 +145,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 +179,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; @@ -191,10 +231,49 @@ int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px * where alpha is a real scalar, x is an n element vector and A is an * n by n symmetric matrix, supplied in packed form. */ -// int EIGEN_BLAS_FUNC(spr)(char *uplo, int *n, Scalar *alpha, Scalar *x, int *incx, Scalar *ap) -// { -// return 1; -// } +int EIGEN_BLAS_FUNC(spr)(char *uplo, int *n, Scalar *palpha, Scalar *px, int *incx, Scalar *pap) +{ + typedef void (*functype)(int, 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::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,false>::run); + func[LO] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,false>::run); + + init = true; + } + + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar* ap = reinterpret_cast<Scalar*>(pap); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + + int info = 0; + if(UPLO(*uplo)==INVALID) info = 1; + else if(*n<0) info = 2; + else if(*incx==0) info = 5; + if(info) + return xerbla_(SCALAR_SUFFIX_UP"SPR ",&info,6); + + if(alpha==Scalar(0)) + return 1; + + Scalar* x_cpy = get_compact_vector(x, *n, *incx); + + int code = UPLO(*uplo); + if(code>=2 || func[code]==0) + return 0; + + func[code](*n, ap, x_cpy, alpha); + + if(x_cpy!=x) delete[] x_cpy; + + return 1; +} /** DSPR2 performs the symmetric rank 2 operation * @@ -203,8 +282,89 @@ int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px * where alpha is a scalar, x and y are n element vectors and A is an * n by n symmetric matrix, supplied in packed form. */ -// int EIGEN_BLAS_FUNC(spr2)(char *uplo, int *n, RealScalar *alpha, RealScalar *x, int *incx, RealScalar *y, int *incy, RealScalar *ap) -// { -// return 1; -// } +int EIGEN_BLAS_FUNC(spr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pap) +{ + typedef void (*functype)(int, Scalar*, 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::packed_rank2_update_selector<Scalar,int,Upper>::run); + func[LO] = (internal::packed_rank2_update_selector<Scalar,int,Lower>::run); + + init = true; + } + + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar* y = reinterpret_cast<Scalar*>(py); + Scalar* ap = reinterpret_cast<Scalar*>(pap); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + + int info = 0; + if(UPLO(*uplo)==INVALID) info = 1; + else if(*n<0) info = 2; + else if(*incx==0) info = 5; + else if(*incy==0) info = 7; + if(info) + return xerbla_(SCALAR_SUFFIX_UP"SPR2 ",&info,6); + + if(alpha==Scalar(0)) + return 1; + + 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; + + func[code](*n, ap, x_cpy, y_cpy, alpha); + + if(x_cpy!=x) delete[] x_cpy; + if(y_cpy!=y) delete[] y_cpy; + + 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,ColMajor,false,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; +} + diff --git a/blas/level3_impl.h b/blas/level3_impl.h index 2371f25c3..07dbc22ff 100644 --- a/blas/level3_impl.h +++ b/blas/level3_impl.h @@ -152,7 +152,7 @@ int EIGEN_BLAS_FUNC(trsm)(char *side, char *uplo, char *opa, char *diag, int *m, int EIGEN_BLAS_FUNC(trmm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb) { // std::cerr << "in trmm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << " " << *n << " " << *lda << " " << *ldb << " " << *palpha << "\n"; - typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar, internal::level3_blocking<Scalar,Scalar>&); + typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&, internal::level3_blocking<Scalar,Scalar>&); static functype func[32]; static bool init = false; if(!init) @@ -305,7 +305,8 @@ 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"; - typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar); + #if !ISCOMPLEX + typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&); static functype func[8]; static bool init = false; @@ -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); @@ -498,7 +500,7 @@ int EIGEN_BLAS_FUNC(hemm)(char *side, char *uplo, int *m, int *n, RealScalar *pa // c = alpha*conj(a')*a + beta*c for op = 'C'or'c' int EIGEN_BLAS_FUNC(herk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc) { - typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar); + typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&); static functype func[8]; static bool init = false; @@ -606,24 +608,24 @@ int EIGEN_BLAS_FUNC(her2k)(char *uplo, char *op, int *n, int *k, RealScalar *pal if(UPLO(*uplo)==UP) { matrix(c, *n, *n, *ldc).triangularView<Upper>() - += alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint() - + internal::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint(); + += alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint() + + numext::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint(); } else if(UPLO(*uplo)==LO) matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha*matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint() - + internal::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint(); + + numext::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint(); } else if(OP(*op)==ADJ) { if(UPLO(*uplo)==UP) matrix(c, *n, *n, *ldc).triangularView<Upper>() - += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb) - + internal::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda); + += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb) + + numext::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda); else if(UPLO(*uplo)==LO) matrix(c, *n, *n, *ldc).triangularView<Lower>() - += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb) - + internal::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda); + += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb) + + numext::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda); } return 1; 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/sspr.f b/blas/sspr.f deleted file mode 100644 index bae92612e..000000000 --- a/blas/sspr.f +++ /dev/null @@ -1,202 +0,0 @@ - SUBROUTINE SSPR(UPLO,N,ALPHA,X,INCX,AP) -* .. Scalar Arguments .. - REAL ALPHA - INTEGER INCX,N - CHARACTER UPLO -* .. -* .. Array Arguments .. - REAL AP(*),X(*) -* .. -* -* Purpose -* ======= -* -* SSPR performs the symmetric rank 1 operation -* -* A := alpha*x*x' + A, -* -* where alpha is a real scalar, x is an n element vector and A is an -* n by n symmetric matrix, supplied in packed form. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the upper or lower -* triangular part of the matrix A is supplied in the packed -* array AP as follows: -* -* UPLO = 'U' or 'u' The upper triangular part of A is -* supplied in AP. -* -* UPLO = 'L' or 'l' The lower triangular part of A is -* supplied in AP. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* ALPHA - REAL . -* On entry, ALPHA specifies the scalar alpha. -* Unchanged on exit. -* -* X - REAL array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element vector x. -* Unchanged on exit. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* AP - REAL array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular part of the symmetric matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) -* and a( 2, 2 ) respectively, and so on. On exit, the array -* AP is overwritten by the upper triangular part of the -* updated matrix. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular part of the symmetric matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) -* and a( 3, 1 ) respectively, and so on. On exit, the array -* AP is overwritten by the lower triangular part of the -* updated matrix. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - REAL ZERO - PARAMETER (ZERO=0.0E+0) -* .. -* .. Local Scalars .. - REAL TEMP - INTEGER I,INFO,IX,J,JX,K,KK,KX -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (N.LT.0) THEN - INFO = 2 - ELSE IF (INCX.EQ.0) THEN - INFO = 5 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('SSPR ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN -* -* Set the start point in X if the increment is not unity. -* - IF (INCX.LE.0) THEN - KX = 1 - (N-1)*INCX - ELSE IF (INCX.NE.1) THEN - KX = 1 - END IF -* -* Start the operations. In this version the elements of the array AP -* are accessed sequentially with one pass through AP. -* - KK = 1 - IF (LSAME(UPLO,'U')) THEN -* -* Form A when upper triangle is stored in AP. -* - IF (INCX.EQ.1) THEN - DO 20 J = 1,N - IF (X(J).NE.ZERO) THEN - TEMP = ALPHA*X(J) - K = KK - DO 10 I = 1,J - AP(K) = AP(K) + X(I)*TEMP - K = K + 1 - 10 CONTINUE - END IF - KK = KK + J - 20 CONTINUE - ELSE - JX = KX - DO 40 J = 1,N - IF (X(JX).NE.ZERO) THEN - TEMP = ALPHA*X(JX) - IX = KX - DO 30 K = KK,KK + J - 1 - AP(K) = AP(K) + X(IX)*TEMP - IX = IX + INCX - 30 CONTINUE - END IF - JX = JX + INCX - KK = KK + J - 40 CONTINUE - END IF - ELSE -* -* Form A when lower triangle is stored in AP. -* - IF (INCX.EQ.1) THEN - DO 60 J = 1,N - IF (X(J).NE.ZERO) THEN - TEMP = ALPHA*X(J) - K = KK - DO 50 I = J,N - AP(K) = AP(K) + X(I)*TEMP - K = K + 1 - 50 CONTINUE - END IF - KK = KK + N - J + 1 - 60 CONTINUE - ELSE - JX = KX - DO 80 J = 1,N - IF (X(JX).NE.ZERO) THEN - TEMP = ALPHA*X(JX) - IX = JX - DO 70 K = KK,KK + N - J - AP(K) = AP(K) + X(IX)*TEMP - IX = IX + INCX - 70 CONTINUE - END IF - JX = JX + INCX - KK = KK + N - J + 1 - 80 CONTINUE - END IF - END IF -* - RETURN -* -* End of SSPR . -* - END diff --git a/blas/sspr2.f b/blas/sspr2.f deleted file mode 100644 index cd27c734b..000000000 --- a/blas/sspr2.f +++ /dev/null @@ -1,233 +0,0 @@ - SUBROUTINE SSPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP) -* .. Scalar Arguments .. - REAL ALPHA - INTEGER INCX,INCY,N - CHARACTER UPLO -* .. -* .. Array Arguments .. - REAL AP(*),X(*),Y(*) -* .. -* -* Purpose -* ======= -* -* SSPR2 performs the symmetric rank 2 operation -* -* A := alpha*x*y' + alpha*y*x' + A, -* -* where alpha is a scalar, x and y are n element vectors and A is an -* n by n symmetric matrix, supplied in packed form. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the upper or lower -* triangular part of the matrix A is supplied in the packed -* array AP as follows: -* -* UPLO = 'U' or 'u' The upper triangular part of A is -* supplied in AP. -* -* UPLO = 'L' or 'l' The lower triangular part of A is -* supplied in AP. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* ALPHA - REAL . -* On entry, ALPHA specifies the scalar alpha. -* Unchanged on exit. -* -* X - REAL array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element vector x. -* Unchanged on exit. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* Y - REAL array of dimension at least -* ( 1 + ( n - 1 )*abs( INCY ) ). -* Before entry, the incremented array Y must contain the n -* element vector y. -* Unchanged on exit. -* -* INCY - INTEGER. -* On entry, INCY specifies the increment for the elements of -* Y. INCY must not be zero. -* Unchanged on exit. -* -* AP - REAL array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular part of the symmetric matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) -* and a( 2, 2 ) respectively, and so on. On exit, the array -* AP is overwritten by the upper triangular part of the -* updated matrix. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular part of the symmetric matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) -* and a( 3, 1 ) respectively, and so on. On exit, the array -* AP is overwritten by the lower triangular part of the -* updated matrix. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - REAL ZERO - PARAMETER (ZERO=0.0E+0) -* .. -* .. Local Scalars .. - REAL TEMP1,TEMP2 - INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (N.LT.0) THEN - INFO = 2 - ELSE IF (INCX.EQ.0) THEN - INFO = 5 - ELSE IF (INCY.EQ.0) THEN - INFO = 7 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('SSPR2 ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN -* -* Set up the start points in X and Y if the increments are not both -* unity. -* - IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN - IF (INCX.GT.0) THEN - KX = 1 - ELSE - KX = 1 - (N-1)*INCX - END IF - IF (INCY.GT.0) THEN - KY = 1 - ELSE - KY = 1 - (N-1)*INCY - END IF - JX = KX - JY = KY - END IF -* -* Start the operations. In this version the elements of the array AP -* are accessed sequentially with one pass through AP. -* - KK = 1 - IF (LSAME(UPLO,'U')) THEN -* -* Form A when upper triangle is stored in AP. -* - IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN - DO 20 J = 1,N - IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN - TEMP1 = ALPHA*Y(J) - TEMP2 = ALPHA*X(J) - K = KK - DO 10 I = 1,J - AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2 - K = K + 1 - 10 CONTINUE - END IF - KK = KK + J - 20 CONTINUE - ELSE - DO 40 J = 1,N - IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN - TEMP1 = ALPHA*Y(JY) - TEMP2 = ALPHA*X(JX) - IX = KX - IY = KY - DO 30 K = KK,KK + J - 1 - AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2 - IX = IX + INCX - IY = IY + INCY - 30 CONTINUE - END IF - JX = JX + INCX - JY = JY + INCY - KK = KK + J - 40 CONTINUE - END IF - ELSE -* -* Form A when lower triangle is stored in AP. -* - IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN - DO 60 J = 1,N - IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN - TEMP1 = ALPHA*Y(J) - TEMP2 = ALPHA*X(J) - K = KK - DO 50 I = J,N - AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2 - K = K + 1 - 50 CONTINUE - END IF - KK = KK + N - J + 1 - 60 CONTINUE - ELSE - DO 80 J = 1,N - IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN - TEMP1 = ALPHA*Y(JY) - TEMP2 = ALPHA*X(JX) - IX = JX - IY = JY - DO 70 K = KK,KK + N - J - AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2 - IX = IX + INCX - IY = IY + INCY - 70 CONTINUE - END IF - JX = JX + INCX - JY = JY + INCY - KK = KK + N - J + 1 - 80 CONTINUE - END IF - END IF -* - RETURN -* -* End of SSPR2 . -* - END diff --git a/blas/stpmv.f b/blas/stpmv.f deleted file mode 100644 index 71ea49a36..000000000 --- a/blas/stpmv.f +++ /dev/null @@ -1,293 +0,0 @@ - SUBROUTINE STPMV(UPLO,TRANS,DIAG,N,AP,X,INCX) -* .. Scalar Arguments .. - INTEGER INCX,N - CHARACTER DIAG,TRANS,UPLO -* .. -* .. Array Arguments .. - REAL AP(*),X(*) -* .. -* -* Purpose -* ======= -* -* STPMV performs one of the matrix-vector operations -* -* x := A*x, or x := A'*x, -* -* where x is an n element vector and A is an n by n unit, or non-unit, -* upper or lower triangular matrix, supplied in packed form. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the matrix is an upper or -* lower triangular matrix as follows: -* -* UPLO = 'U' or 'u' A is an upper triangular matrix. -* -* UPLO = 'L' or 'l' A is a lower triangular matrix. -* -* Unchanged on exit. -* -* TRANS - CHARACTER*1. -* On entry, TRANS specifies the operation to be performed as -* follows: -* -* TRANS = 'N' or 'n' x := A*x. -* -* TRANS = 'T' or 't' x := A'*x. -* -* TRANS = 'C' or 'c' x := A'*x. -* -* Unchanged on exit. -* -* DIAG - CHARACTER*1. -* On entry, DIAG specifies whether or not A is unit -* triangular as follows: -* -* DIAG = 'U' or 'u' A is assumed to be unit triangular. -* -* DIAG = 'N' or 'n' A is not assumed to be unit -* triangular. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* AP - REAL array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) -* respectively, and so on. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) -* respectively, and so on. -* Note that when DIAG = 'U' or 'u', the diagonal elements of -* A are not referenced, but are assumed to be unity. -* Unchanged on exit. -* -* X - REAL array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element vector x. On exit, X is overwritten with the -* tranformed vector x. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - REAL ZERO - PARAMETER (ZERO=0.0E+0) -* .. -* .. Local Scalars .. - REAL TEMP - INTEGER I,INFO,IX,J,JX,K,KK,KX - LOGICAL NOUNIT -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND. - + .NOT.LSAME(TRANS,'C')) THEN - INFO = 2 - ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN - INFO = 3 - ELSE IF (N.LT.0) THEN - INFO = 4 - ELSE IF (INCX.EQ.0) THEN - INFO = 7 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('STPMV ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF (N.EQ.0) RETURN -* - NOUNIT = LSAME(DIAG,'N') -* -* Set up the start point in X if the increment is not unity. This -* will be ( N - 1 )*INCX too small for descending loops. -* - IF (INCX.LE.0) THEN - KX = 1 - (N-1)*INCX - ELSE IF (INCX.NE.1) THEN - KX = 1 - END IF -* -* Start the operations. In this version the elements of AP are -* accessed sequentially with one pass through AP. -* - IF (LSAME(TRANS,'N')) THEN -* -* Form x:= A*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = 1 - IF (INCX.EQ.1) THEN - DO 20 J = 1,N - IF (X(J).NE.ZERO) THEN - TEMP = X(J) - K = KK - DO 10 I = 1,J - 1 - X(I) = X(I) + TEMP*AP(K) - K = K + 1 - 10 CONTINUE - IF (NOUNIT) X(J) = X(J)*AP(KK+J-1) - END IF - KK = KK + J - 20 CONTINUE - ELSE - JX = KX - DO 40 J = 1,N - IF (X(JX).NE.ZERO) THEN - TEMP = X(JX) - IX = KX - DO 30 K = KK,KK + J - 2 - X(IX) = X(IX) + TEMP*AP(K) - IX = IX + INCX - 30 CONTINUE - IF (NOUNIT) X(JX) = X(JX)*AP(KK+J-1) - END IF - JX = JX + INCX - KK = KK + J - 40 CONTINUE - END IF - ELSE - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 60 J = N,1,-1 - IF (X(J).NE.ZERO) THEN - TEMP = X(J) - K = KK - DO 50 I = N,J + 1,-1 - X(I) = X(I) + TEMP*AP(K) - K = K - 1 - 50 CONTINUE - IF (NOUNIT) X(J) = X(J)*AP(KK-N+J) - END IF - KK = KK - (N-J+1) - 60 CONTINUE - ELSE - KX = KX + (N-1)*INCX - JX = KX - DO 80 J = N,1,-1 - IF (X(JX).NE.ZERO) THEN - TEMP = X(JX) - IX = KX - DO 70 K = KK,KK - (N- (J+1)),-1 - X(IX) = X(IX) + TEMP*AP(K) - IX = IX - INCX - 70 CONTINUE - IF (NOUNIT) X(JX) = X(JX)*AP(KK-N+J) - END IF - JX = JX - INCX - KK = KK - (N-J+1) - 80 CONTINUE - END IF - END IF - ELSE -* -* Form x := A'*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 100 J = N,1,-1 - TEMP = X(J) - IF (NOUNIT) TEMP = TEMP*AP(KK) - K = KK - 1 - DO 90 I = J - 1,1,-1 - TEMP = TEMP + AP(K)*X(I) - K = K - 1 - 90 CONTINUE - X(J) = TEMP - KK = KK - J - 100 CONTINUE - ELSE - JX = KX + (N-1)*INCX - DO 120 J = N,1,-1 - TEMP = X(JX) - IX = JX - IF (NOUNIT) TEMP = TEMP*AP(KK) - DO 110 K = KK - 1,KK - J + 1,-1 - IX = IX - INCX - TEMP = TEMP + AP(K)*X(IX) - 110 CONTINUE - X(JX) = TEMP - JX = JX - INCX - KK = KK - J - 120 CONTINUE - END IF - ELSE - KK = 1 - IF (INCX.EQ.1) THEN - DO 140 J = 1,N - TEMP = X(J) - IF (NOUNIT) TEMP = TEMP*AP(KK) - K = KK + 1 - DO 130 I = J + 1,N - TEMP = TEMP + AP(K)*X(I) - K = K + 1 - 130 CONTINUE - X(J) = TEMP - KK = KK + (N-J+1) - 140 CONTINUE - ELSE - JX = KX - DO 160 J = 1,N - TEMP = X(JX) - IX = JX - IF (NOUNIT) TEMP = TEMP*AP(KK) - DO 150 K = KK + 1,KK + N - J - IX = IX + INCX - TEMP = TEMP + AP(K)*X(IX) - 150 CONTINUE - X(JX) = TEMP - JX = JX + INCX - KK = KK + (N-J+1) - 160 CONTINUE - END IF - END IF - END IF -* - RETURN -* -* End of STPMV . -* - END diff --git a/blas/stpsv.f b/blas/stpsv.f deleted file mode 100644 index 7d95efbde..000000000 --- a/blas/stpsv.f +++ /dev/null @@ -1,296 +0,0 @@ - SUBROUTINE STPSV(UPLO,TRANS,DIAG,N,AP,X,INCX) -* .. Scalar Arguments .. - INTEGER INCX,N - CHARACTER DIAG,TRANS,UPLO -* .. -* .. Array Arguments .. - REAL AP(*),X(*) -* .. -* -* Purpose -* ======= -* -* STPSV solves one of the systems of equations -* -* A*x = b, or A'*x = b, -* -* where b and x are n element vectors and A is an n by n unit, or -* non-unit, upper or lower triangular matrix, supplied in packed form. -* -* No test for singularity or near-singularity is included in this -* routine. Such tests must be performed before calling this routine. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the matrix is an upper or -* lower triangular matrix as follows: -* -* UPLO = 'U' or 'u' A is an upper triangular matrix. -* -* UPLO = 'L' or 'l' A is a lower triangular matrix. -* -* Unchanged on exit. -* -* TRANS - CHARACTER*1. -* On entry, TRANS specifies the equations to be solved as -* follows: -* -* TRANS = 'N' or 'n' A*x = b. -* -* TRANS = 'T' or 't' A'*x = b. -* -* TRANS = 'C' or 'c' A'*x = b. -* -* Unchanged on exit. -* -* DIAG - CHARACTER*1. -* On entry, DIAG specifies whether or not A is unit -* triangular as follows: -* -* DIAG = 'U' or 'u' A is assumed to be unit triangular. -* -* DIAG = 'N' or 'n' A is not assumed to be unit -* triangular. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* AP - REAL array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) -* respectively, and so on. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) -* respectively, and so on. -* Note that when DIAG = 'U' or 'u', the diagonal elements of -* A are not referenced, but are assumed to be unity. -* Unchanged on exit. -* -* X - REAL array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element right-hand side vector b. On exit, X is overwritten -* with the solution vector x. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - REAL ZERO - PARAMETER (ZERO=0.0E+0) -* .. -* .. Local Scalars .. - REAL TEMP - INTEGER I,INFO,IX,J,JX,K,KK,KX - LOGICAL NOUNIT -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND. - + .NOT.LSAME(TRANS,'C')) THEN - INFO = 2 - ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN - INFO = 3 - ELSE IF (N.LT.0) THEN - INFO = 4 - ELSE IF (INCX.EQ.0) THEN - INFO = 7 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('STPSV ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF (N.EQ.0) RETURN -* - NOUNIT = LSAME(DIAG,'N') -* -* Set up the start point in X if the increment is not unity. This -* will be ( N - 1 )*INCX too small for descending loops. -* - IF (INCX.LE.0) THEN - KX = 1 - (N-1)*INCX - ELSE IF (INCX.NE.1) THEN - KX = 1 - END IF -* -* Start the operations. In this version the elements of AP are -* accessed sequentially with one pass through AP. -* - IF (LSAME(TRANS,'N')) THEN -* -* Form x := inv( A )*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 20 J = N,1,-1 - IF (X(J).NE.ZERO) THEN - IF (NOUNIT) X(J) = X(J)/AP(KK) - TEMP = X(J) - K = KK - 1 - DO 10 I = J - 1,1,-1 - X(I) = X(I) - TEMP*AP(K) - K = K - 1 - 10 CONTINUE - END IF - KK = KK - J - 20 CONTINUE - ELSE - JX = KX + (N-1)*INCX - DO 40 J = N,1,-1 - IF (X(JX).NE.ZERO) THEN - IF (NOUNIT) X(JX) = X(JX)/AP(KK) - TEMP = X(JX) - IX = JX - DO 30 K = KK - 1,KK - J + 1,-1 - IX = IX - INCX - X(IX) = X(IX) - TEMP*AP(K) - 30 CONTINUE - END IF - JX = JX - INCX - KK = KK - J - 40 CONTINUE - END IF - ELSE - KK = 1 - IF (INCX.EQ.1) THEN - DO 60 J = 1,N - IF (X(J).NE.ZERO) THEN - IF (NOUNIT) X(J) = X(J)/AP(KK) - TEMP = X(J) - K = KK + 1 - DO 50 I = J + 1,N - X(I) = X(I) - TEMP*AP(K) - K = K + 1 - 50 CONTINUE - END IF - KK = KK + (N-J+1) - 60 CONTINUE - ELSE - JX = KX - DO 80 J = 1,N - IF (X(JX).NE.ZERO) THEN - IF (NOUNIT) X(JX) = X(JX)/AP(KK) - TEMP = X(JX) - IX = JX - DO 70 K = KK + 1,KK + N - J - IX = IX + INCX - X(IX) = X(IX) - TEMP*AP(K) - 70 CONTINUE - END IF - JX = JX + INCX - KK = KK + (N-J+1) - 80 CONTINUE - END IF - END IF - ELSE -* -* Form x := inv( A' )*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = 1 - IF (INCX.EQ.1) THEN - DO 100 J = 1,N - TEMP = X(J) - K = KK - DO 90 I = 1,J - 1 - TEMP = TEMP - AP(K)*X(I) - K = K + 1 - 90 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK+J-1) - X(J) = TEMP - KK = KK + J - 100 CONTINUE - ELSE - JX = KX - DO 120 J = 1,N - TEMP = X(JX) - IX = KX - DO 110 K = KK,KK + J - 2 - TEMP = TEMP - AP(K)*X(IX) - IX = IX + INCX - 110 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK+J-1) - X(JX) = TEMP - JX = JX + INCX - KK = KK + J - 120 CONTINUE - END IF - ELSE - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 140 J = N,1,-1 - TEMP = X(J) - K = KK - DO 130 I = N,J + 1,-1 - TEMP = TEMP - AP(K)*X(I) - K = K - 1 - 130 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK-N+J) - X(J) = TEMP - KK = KK - (N-J+1) - 140 CONTINUE - ELSE - KX = KX + (N-1)*INCX - JX = KX - DO 160 J = N,1,-1 - TEMP = X(JX) - IX = KX - DO 150 K = KK,KK - (N- (J+1)),-1 - TEMP = TEMP - AP(K)*X(IX) - IX = IX - INCX - 150 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK-N+J) - X(JX) = TEMP - JX = JX - INCX - KK = KK - (N-J+1) - 160 CONTINUE - END IF - END IF - END IF -* - RETURN -* -* End of STPSV . -* - END 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,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,3I5,I3,2D36.8,2D12.4) +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/blas/zhpr.f b/blas/zhpr.f deleted file mode 100644 index 40efbc7d5..000000000 --- a/blas/zhpr.f +++ /dev/null @@ -1,220 +0,0 @@ - SUBROUTINE ZHPR(UPLO,N,ALPHA,X,INCX,AP) -* .. Scalar Arguments .. - DOUBLE PRECISION ALPHA - INTEGER INCX,N - CHARACTER UPLO -* .. -* .. Array Arguments .. - DOUBLE COMPLEX AP(*),X(*) -* .. -* -* Purpose -* ======= -* -* ZHPR performs the hermitian rank 1 operation -* -* A := alpha*x*conjg( x' ) + A, -* -* where alpha is a real scalar, x is an n element vector and A is an -* n by n hermitian matrix, supplied in packed form. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the upper or lower -* triangular part of the matrix A is supplied in the packed -* array AP as follows: -* -* UPLO = 'U' or 'u' The upper triangular part of A is -* supplied in AP. -* -* UPLO = 'L' or 'l' The lower triangular part of A is -* supplied in AP. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* ALPHA - DOUBLE PRECISION. -* On entry, ALPHA specifies the scalar alpha. -* Unchanged on exit. -* -* X - COMPLEX*16 array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element vector x. -* Unchanged on exit. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* AP - COMPLEX*16 array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular part of the hermitian matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) -* and a( 2, 2 ) respectively, and so on. On exit, the array -* AP is overwritten by the upper triangular part of the -* updated matrix. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular part of the hermitian matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) -* and a( 3, 1 ) respectively, and so on. On exit, the array -* AP is overwritten by the lower triangular part of the -* updated matrix. -* Note that the imaginary parts of the diagonal elements need -* not be set, they are assumed to be zero, and on exit they -* are set to zero. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - DOUBLE COMPLEX ZERO - PARAMETER (ZERO= (0.0D+0,0.0D+0)) -* .. -* .. Local Scalars .. - DOUBLE COMPLEX TEMP - INTEGER I,INFO,IX,J,JX,K,KK,KX -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC DBLE,DCONJG -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (N.LT.0) THEN - INFO = 2 - ELSE IF (INCX.EQ.0) THEN - INFO = 5 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('ZHPR ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF ((N.EQ.0) .OR. (ALPHA.EQ.DBLE(ZERO))) RETURN -* -* Set the start point in X if the increment is not unity. -* - IF (INCX.LE.0) THEN - KX = 1 - (N-1)*INCX - ELSE IF (INCX.NE.1) THEN - KX = 1 - END IF -* -* Start the operations. In this version the elements of the array AP -* are accessed sequentially with one pass through AP. -* - KK = 1 - IF (LSAME(UPLO,'U')) THEN -* -* Form A when upper triangle is stored in AP. -* - IF (INCX.EQ.1) THEN - DO 20 J = 1,N - IF (X(J).NE.ZERO) THEN - TEMP = ALPHA*DCONJG(X(J)) - K = KK - DO 10 I = 1,J - 1 - AP(K) = AP(K) + X(I)*TEMP - K = K + 1 - 10 CONTINUE - AP(KK+J-1) = DBLE(AP(KK+J-1)) + DBLE(X(J)*TEMP) - ELSE - AP(KK+J-1) = DBLE(AP(KK+J-1)) - END IF - KK = KK + J - 20 CONTINUE - ELSE - JX = KX - DO 40 J = 1,N - IF (X(JX).NE.ZERO) THEN - TEMP = ALPHA*DCONJG(X(JX)) - IX = KX - DO 30 K = KK,KK + J - 2 - AP(K) = AP(K) + X(IX)*TEMP - IX = IX + INCX - 30 CONTINUE - AP(KK+J-1) = DBLE(AP(KK+J-1)) + DBLE(X(JX)*TEMP) - ELSE - AP(KK+J-1) = DBLE(AP(KK+J-1)) - END IF - JX = JX + INCX - KK = KK + J - 40 CONTINUE - END IF - ELSE -* -* Form A when lower triangle is stored in AP. -* - IF (INCX.EQ.1) THEN - DO 60 J = 1,N - IF (X(J).NE.ZERO) THEN - TEMP = ALPHA*DCONJG(X(J)) - AP(KK) = DBLE(AP(KK)) + DBLE(TEMP*X(J)) - K = KK + 1 - DO 50 I = J + 1,N - AP(K) = AP(K) + X(I)*TEMP - K = K + 1 - 50 CONTINUE - ELSE - AP(KK) = DBLE(AP(KK)) - END IF - KK = KK + N - J + 1 - 60 CONTINUE - ELSE - JX = KX - DO 80 J = 1,N - IF (X(JX).NE.ZERO) THEN - TEMP = ALPHA*DCONJG(X(JX)) - AP(KK) = DBLE(AP(KK)) + DBLE(TEMP*X(JX)) - IX = JX - DO 70 K = KK + 1,KK + N - J - IX = IX + INCX - AP(K) = AP(K) + X(IX)*TEMP - 70 CONTINUE - ELSE - AP(KK) = DBLE(AP(KK)) - END IF - JX = JX + INCX - KK = KK + N - J + 1 - 80 CONTINUE - END IF - END IF -* - RETURN -* -* End of ZHPR . -* - END diff --git a/blas/zhpr2.f b/blas/zhpr2.f deleted file mode 100644 index 99977462e..000000000 --- a/blas/zhpr2.f +++ /dev/null @@ -1,255 +0,0 @@ - SUBROUTINE ZHPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP) -* .. Scalar Arguments .. - DOUBLE COMPLEX ALPHA - INTEGER INCX,INCY,N - CHARACTER UPLO -* .. -* .. Array Arguments .. - DOUBLE COMPLEX AP(*),X(*),Y(*) -* .. -* -* Purpose -* ======= -* -* ZHPR2 performs the hermitian rank 2 operation -* -* A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, -* -* where alpha is a scalar, x and y are n element vectors and A is an -* n by n hermitian matrix, supplied in packed form. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the upper or lower -* triangular part of the matrix A is supplied in the packed -* array AP as follows: -* -* UPLO = 'U' or 'u' The upper triangular part of A is -* supplied in AP. -* -* UPLO = 'L' or 'l' The lower triangular part of A is -* supplied in AP. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* ALPHA - COMPLEX*16 . -* On entry, ALPHA specifies the scalar alpha. -* Unchanged on exit. -* -* X - COMPLEX*16 array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element vector x. -* Unchanged on exit. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* Y - COMPLEX*16 array of dimension at least -* ( 1 + ( n - 1 )*abs( INCY ) ). -* Before entry, the incremented array Y must contain the n -* element vector y. -* Unchanged on exit. -* -* INCY - INTEGER. -* On entry, INCY specifies the increment for the elements of -* Y. INCY must not be zero. -* Unchanged on exit. -* -* AP - COMPLEX*16 array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular part of the hermitian matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) -* and a( 2, 2 ) respectively, and so on. On exit, the array -* AP is overwritten by the upper triangular part of the -* updated matrix. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular part of the hermitian matrix -* packed sequentially, column by column, so that AP( 1 ) -* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) -* and a( 3, 1 ) respectively, and so on. On exit, the array -* AP is overwritten by the lower triangular part of the -* updated matrix. -* Note that the imaginary parts of the diagonal elements need -* not be set, they are assumed to be zero, and on exit they -* are set to zero. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - DOUBLE COMPLEX ZERO - PARAMETER (ZERO= (0.0D+0,0.0D+0)) -* .. -* .. Local Scalars .. - DOUBLE COMPLEX TEMP1,TEMP2 - INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC DBLE,DCONJG -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (N.LT.0) THEN - INFO = 2 - ELSE IF (INCX.EQ.0) THEN - INFO = 5 - ELSE IF (INCY.EQ.0) THEN - INFO = 7 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('ZHPR2 ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN -* -* Set up the start points in X and Y if the increments are not both -* unity. -* - IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN - IF (INCX.GT.0) THEN - KX = 1 - ELSE - KX = 1 - (N-1)*INCX - END IF - IF (INCY.GT.0) THEN - KY = 1 - ELSE - KY = 1 - (N-1)*INCY - END IF - JX = KX - JY = KY - END IF -* -* Start the operations. In this version the elements of the array AP -* are accessed sequentially with one pass through AP. -* - KK = 1 - IF (LSAME(UPLO,'U')) THEN -* -* Form A when upper triangle is stored in AP. -* - IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN - DO 20 J = 1,N - IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN - TEMP1 = ALPHA*DCONJG(Y(J)) - TEMP2 = DCONJG(ALPHA*X(J)) - K = KK - DO 10 I = 1,J - 1 - AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2 - K = K + 1 - 10 CONTINUE - AP(KK+J-1) = DBLE(AP(KK+J-1)) + - + DBLE(X(J)*TEMP1+Y(J)*TEMP2) - ELSE - AP(KK+J-1) = DBLE(AP(KK+J-1)) - END IF - KK = KK + J - 20 CONTINUE - ELSE - DO 40 J = 1,N - IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN - TEMP1 = ALPHA*DCONJG(Y(JY)) - TEMP2 = DCONJG(ALPHA*X(JX)) - IX = KX - IY = KY - DO 30 K = KK,KK + J - 2 - AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2 - IX = IX + INCX - IY = IY + INCY - 30 CONTINUE - AP(KK+J-1) = DBLE(AP(KK+J-1)) + - + DBLE(X(JX)*TEMP1+Y(JY)*TEMP2) - ELSE - AP(KK+J-1) = DBLE(AP(KK+J-1)) - END IF - JX = JX + INCX - JY = JY + INCY - KK = KK + J - 40 CONTINUE - END IF - ELSE -* -* Form A when lower triangle is stored in AP. -* - IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN - DO 60 J = 1,N - IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN - TEMP1 = ALPHA*DCONJG(Y(J)) - TEMP2 = DCONJG(ALPHA*X(J)) - AP(KK) = DBLE(AP(KK)) + - + DBLE(X(J)*TEMP1+Y(J)*TEMP2) - K = KK + 1 - DO 50 I = J + 1,N - AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2 - K = K + 1 - 50 CONTINUE - ELSE - AP(KK) = DBLE(AP(KK)) - END IF - KK = KK + N - J + 1 - 60 CONTINUE - ELSE - DO 80 J = 1,N - IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN - TEMP1 = ALPHA*DCONJG(Y(JY)) - TEMP2 = DCONJG(ALPHA*X(JX)) - AP(KK) = DBLE(AP(KK)) + - + DBLE(X(JX)*TEMP1+Y(JY)*TEMP2) - IX = JX - IY = JY - DO 70 K = KK + 1,KK + N - J - IX = IX + INCX - IY = IY + INCY - AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2 - 70 CONTINUE - ELSE - AP(KK) = DBLE(AP(KK)) - END IF - JX = JX + INCX - JY = JY + INCY - KK = KK + N - J + 1 - 80 CONTINUE - END IF - END IF -* - RETURN -* -* End of ZHPR2 . -* - END diff --git a/blas/ztpmv.f b/blas/ztpmv.f deleted file mode 100644 index 5a7b3b8b7..000000000 --- a/blas/ztpmv.f +++ /dev/null @@ -1,329 +0,0 @@ - SUBROUTINE ZTPMV(UPLO,TRANS,DIAG,N,AP,X,INCX) -* .. Scalar Arguments .. - INTEGER INCX,N - CHARACTER DIAG,TRANS,UPLO -* .. -* .. Array Arguments .. - DOUBLE COMPLEX AP(*),X(*) -* .. -* -* Purpose -* ======= -* -* ZTPMV performs one of the matrix-vector operations -* -* x := A*x, or x := A'*x, or x := conjg( A' )*x, -* -* where x is an n element vector and A is an n by n unit, or non-unit, -* upper or lower triangular matrix, supplied in packed form. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the matrix is an upper or -* lower triangular matrix as follows: -* -* UPLO = 'U' or 'u' A is an upper triangular matrix. -* -* UPLO = 'L' or 'l' A is a lower triangular matrix. -* -* Unchanged on exit. -* -* TRANS - CHARACTER*1. -* On entry, TRANS specifies the operation to be performed as -* follows: -* -* TRANS = 'N' or 'n' x := A*x. -* -* TRANS = 'T' or 't' x := A'*x. -* -* TRANS = 'C' or 'c' x := conjg( A' )*x. -* -* Unchanged on exit. -* -* DIAG - CHARACTER*1. -* On entry, DIAG specifies whether or not A is unit -* triangular as follows: -* -* DIAG = 'U' or 'u' A is assumed to be unit triangular. -* -* DIAG = 'N' or 'n' A is not assumed to be unit -* triangular. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* AP - COMPLEX*16 array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) -* respectively, and so on. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) -* respectively, and so on. -* Note that when DIAG = 'U' or 'u', the diagonal elements of -* A are not referenced, but are assumed to be unity. -* Unchanged on exit. -* -* X - COMPLEX*16 array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element vector x. On exit, X is overwritten with the -* tranformed vector x. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - DOUBLE COMPLEX ZERO - PARAMETER (ZERO= (0.0D+0,0.0D+0)) -* .. -* .. Local Scalars .. - DOUBLE COMPLEX TEMP - INTEGER I,INFO,IX,J,JX,K,KK,KX - LOGICAL NOCONJ,NOUNIT -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC DCONJG -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND. - + .NOT.LSAME(TRANS,'C')) THEN - INFO = 2 - ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN - INFO = 3 - ELSE IF (N.LT.0) THEN - INFO = 4 - ELSE IF (INCX.EQ.0) THEN - INFO = 7 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('ZTPMV ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF (N.EQ.0) RETURN -* - NOCONJ = LSAME(TRANS,'T') - NOUNIT = LSAME(DIAG,'N') -* -* Set up the start point in X if the increment is not unity. This -* will be ( N - 1 )*INCX too small for descending loops. -* - IF (INCX.LE.0) THEN - KX = 1 - (N-1)*INCX - ELSE IF (INCX.NE.1) THEN - KX = 1 - END IF -* -* Start the operations. In this version the elements of AP are -* accessed sequentially with one pass through AP. -* - IF (LSAME(TRANS,'N')) THEN -* -* Form x:= A*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = 1 - IF (INCX.EQ.1) THEN - DO 20 J = 1,N - IF (X(J).NE.ZERO) THEN - TEMP = X(J) - K = KK - DO 10 I = 1,J - 1 - X(I) = X(I) + TEMP*AP(K) - K = K + 1 - 10 CONTINUE - IF (NOUNIT) X(J) = X(J)*AP(KK+J-1) - END IF - KK = KK + J - 20 CONTINUE - ELSE - JX = KX - DO 40 J = 1,N - IF (X(JX).NE.ZERO) THEN - TEMP = X(JX) - IX = KX - DO 30 K = KK,KK + J - 2 - X(IX) = X(IX) + TEMP*AP(K) - IX = IX + INCX - 30 CONTINUE - IF (NOUNIT) X(JX) = X(JX)*AP(KK+J-1) - END IF - JX = JX + INCX - KK = KK + J - 40 CONTINUE - END IF - ELSE - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 60 J = N,1,-1 - IF (X(J).NE.ZERO) THEN - TEMP = X(J) - K = KK - DO 50 I = N,J + 1,-1 - X(I) = X(I) + TEMP*AP(K) - K = K - 1 - 50 CONTINUE - IF (NOUNIT) X(J) = X(J)*AP(KK-N+J) - END IF - KK = KK - (N-J+1) - 60 CONTINUE - ELSE - KX = KX + (N-1)*INCX - JX = KX - DO 80 J = N,1,-1 - IF (X(JX).NE.ZERO) THEN - TEMP = X(JX) - IX = KX - DO 70 K = KK,KK - (N- (J+1)),-1 - X(IX) = X(IX) + TEMP*AP(K) - IX = IX - INCX - 70 CONTINUE - IF (NOUNIT) X(JX) = X(JX)*AP(KK-N+J) - END IF - JX = JX - INCX - KK = KK - (N-J+1) - 80 CONTINUE - END IF - END IF - ELSE -* -* Form x := A'*x or x := conjg( A' )*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 110 J = N,1,-1 - TEMP = X(J) - K = KK - 1 - IF (NOCONJ) THEN - IF (NOUNIT) TEMP = TEMP*AP(KK) - DO 90 I = J - 1,1,-1 - TEMP = TEMP + AP(K)*X(I) - K = K - 1 - 90 CONTINUE - ELSE - IF (NOUNIT) TEMP = TEMP*DCONJG(AP(KK)) - DO 100 I = J - 1,1,-1 - TEMP = TEMP + DCONJG(AP(K))*X(I) - K = K - 1 - 100 CONTINUE - END IF - X(J) = TEMP - KK = KK - J - 110 CONTINUE - ELSE - JX = KX + (N-1)*INCX - DO 140 J = N,1,-1 - TEMP = X(JX) - IX = JX - IF (NOCONJ) THEN - IF (NOUNIT) TEMP = TEMP*AP(KK) - DO 120 K = KK - 1,KK - J + 1,-1 - IX = IX - INCX - TEMP = TEMP + AP(K)*X(IX) - 120 CONTINUE - ELSE - IF (NOUNIT) TEMP = TEMP*DCONJG(AP(KK)) - DO 130 K = KK - 1,KK - J + 1,-1 - IX = IX - INCX - TEMP = TEMP + DCONJG(AP(K))*X(IX) - 130 CONTINUE - END IF - X(JX) = TEMP - JX = JX - INCX - KK = KK - J - 140 CONTINUE - END IF - ELSE - KK = 1 - IF (INCX.EQ.1) THEN - DO 170 J = 1,N - TEMP = X(J) - K = KK + 1 - IF (NOCONJ) THEN - IF (NOUNIT) TEMP = TEMP*AP(KK) - DO 150 I = J + 1,N - TEMP = TEMP + AP(K)*X(I) - K = K + 1 - 150 CONTINUE - ELSE - IF (NOUNIT) TEMP = TEMP*DCONJG(AP(KK)) - DO 160 I = J + 1,N - TEMP = TEMP + DCONJG(AP(K))*X(I) - K = K + 1 - 160 CONTINUE - END IF - X(J) = TEMP - KK = KK + (N-J+1) - 170 CONTINUE - ELSE - JX = KX - DO 200 J = 1,N - TEMP = X(JX) - IX = JX - IF (NOCONJ) THEN - IF (NOUNIT) TEMP = TEMP*AP(KK) - DO 180 K = KK + 1,KK + N - J - IX = IX + INCX - TEMP = TEMP + AP(K)*X(IX) - 180 CONTINUE - ELSE - IF (NOUNIT) TEMP = TEMP*DCONJG(AP(KK)) - DO 190 K = KK + 1,KK + N - J - IX = IX + INCX - TEMP = TEMP + DCONJG(AP(K))*X(IX) - 190 CONTINUE - END IF - X(JX) = TEMP - JX = JX + INCX - KK = KK + (N-J+1) - 200 CONTINUE - END IF - END IF - END IF -* - RETURN -* -* End of ZTPMV . -* - END diff --git a/blas/ztpsv.f b/blas/ztpsv.f deleted file mode 100644 index b56e1d8c4..000000000 --- a/blas/ztpsv.f +++ /dev/null @@ -1,332 +0,0 @@ - SUBROUTINE ZTPSV(UPLO,TRANS,DIAG,N,AP,X,INCX) -* .. Scalar Arguments .. - INTEGER INCX,N - CHARACTER DIAG,TRANS,UPLO -* .. -* .. Array Arguments .. - DOUBLE COMPLEX AP(*),X(*) -* .. -* -* Purpose -* ======= -* -* ZTPSV solves one of the systems of equations -* -* A*x = b, or A'*x = b, or conjg( A' )*x = b, -* -* where b and x are n element vectors and A is an n by n unit, or -* non-unit, upper or lower triangular matrix, supplied in packed form. -* -* No test for singularity or near-singularity is included in this -* routine. Such tests must be performed before calling this routine. -* -* Arguments -* ========== -* -* UPLO - CHARACTER*1. -* On entry, UPLO specifies whether the matrix is an upper or -* lower triangular matrix as follows: -* -* UPLO = 'U' or 'u' A is an upper triangular matrix. -* -* UPLO = 'L' or 'l' A is a lower triangular matrix. -* -* Unchanged on exit. -* -* TRANS - CHARACTER*1. -* On entry, TRANS specifies the equations to be solved as -* follows: -* -* TRANS = 'N' or 'n' A*x = b. -* -* TRANS = 'T' or 't' A'*x = b. -* -* TRANS = 'C' or 'c' conjg( A' )*x = b. -* -* Unchanged on exit. -* -* DIAG - CHARACTER*1. -* On entry, DIAG specifies whether or not A is unit -* triangular as follows: -* -* DIAG = 'U' or 'u' A is assumed to be unit triangular. -* -* DIAG = 'N' or 'n' A is not assumed to be unit -* triangular. -* -* Unchanged on exit. -* -* N - INTEGER. -* On entry, N specifies the order of the matrix A. -* N must be at least zero. -* Unchanged on exit. -* -* AP - COMPLEX*16 array of DIMENSION at least -* ( ( n*( n + 1 ) )/2 ). -* Before entry with UPLO = 'U' or 'u', the array AP must -* contain the upper triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) -* respectively, and so on. -* Before entry with UPLO = 'L' or 'l', the array AP must -* contain the lower triangular matrix packed sequentially, -* column by column, so that AP( 1 ) contains a( 1, 1 ), -* AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) -* respectively, and so on. -* Note that when DIAG = 'U' or 'u', the diagonal elements of -* A are not referenced, but are assumed to be unity. -* Unchanged on exit. -* -* X - COMPLEX*16 array of dimension at least -* ( 1 + ( n - 1 )*abs( INCX ) ). -* Before entry, the incremented array X must contain the n -* element right-hand side vector b. On exit, X is overwritten -* with the solution vector x. -* -* INCX - INTEGER. -* On entry, INCX specifies the increment for the elements of -* X. INCX must not be zero. -* Unchanged on exit. -* -* Further Details -* =============== -* -* Level 2 Blas routine. -* -* -- Written on 22-October-1986. -* Jack Dongarra, Argonne National Lab. -* Jeremy Du Croz, Nag Central Office. -* Sven Hammarling, Nag Central Office. -* Richard Hanson, Sandia National Labs. -* -* ===================================================================== -* -* .. Parameters .. - DOUBLE COMPLEX ZERO - PARAMETER (ZERO= (0.0D+0,0.0D+0)) -* .. -* .. Local Scalars .. - DOUBLE COMPLEX TEMP - INTEGER I,INFO,IX,J,JX,K,KK,KX - LOGICAL NOCONJ,NOUNIT -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL XERBLA -* .. -* .. Intrinsic Functions .. - INTRINSIC DCONJG -* .. -* -* Test the input parameters. -* - INFO = 0 - IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN - INFO = 1 - ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND. - + .NOT.LSAME(TRANS,'C')) THEN - INFO = 2 - ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN - INFO = 3 - ELSE IF (N.LT.0) THEN - INFO = 4 - ELSE IF (INCX.EQ.0) THEN - INFO = 7 - END IF - IF (INFO.NE.0) THEN - CALL XERBLA('ZTPSV ',INFO) - RETURN - END IF -* -* Quick return if possible. -* - IF (N.EQ.0) RETURN -* - NOCONJ = LSAME(TRANS,'T') - NOUNIT = LSAME(DIAG,'N') -* -* Set up the start point in X if the increment is not unity. This -* will be ( N - 1 )*INCX too small for descending loops. -* - IF (INCX.LE.0) THEN - KX = 1 - (N-1)*INCX - ELSE IF (INCX.NE.1) THEN - KX = 1 - END IF -* -* Start the operations. In this version the elements of AP are -* accessed sequentially with one pass through AP. -* - IF (LSAME(TRANS,'N')) THEN -* -* Form x := inv( A )*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 20 J = N,1,-1 - IF (X(J).NE.ZERO) THEN - IF (NOUNIT) X(J) = X(J)/AP(KK) - TEMP = X(J) - K = KK - 1 - DO 10 I = J - 1,1,-1 - X(I) = X(I) - TEMP*AP(K) - K = K - 1 - 10 CONTINUE - END IF - KK = KK - J - 20 CONTINUE - ELSE - JX = KX + (N-1)*INCX - DO 40 J = N,1,-1 - IF (X(JX).NE.ZERO) THEN - IF (NOUNIT) X(JX) = X(JX)/AP(KK) - TEMP = X(JX) - IX = JX - DO 30 K = KK - 1,KK - J + 1,-1 - IX = IX - INCX - X(IX) = X(IX) - TEMP*AP(K) - 30 CONTINUE - END IF - JX = JX - INCX - KK = KK - J - 40 CONTINUE - END IF - ELSE - KK = 1 - IF (INCX.EQ.1) THEN - DO 60 J = 1,N - IF (X(J).NE.ZERO) THEN - IF (NOUNIT) X(J) = X(J)/AP(KK) - TEMP = X(J) - K = KK + 1 - DO 50 I = J + 1,N - X(I) = X(I) - TEMP*AP(K) - K = K + 1 - 50 CONTINUE - END IF - KK = KK + (N-J+1) - 60 CONTINUE - ELSE - JX = KX - DO 80 J = 1,N - IF (X(JX).NE.ZERO) THEN - IF (NOUNIT) X(JX) = X(JX)/AP(KK) - TEMP = X(JX) - IX = JX - DO 70 K = KK + 1,KK + N - J - IX = IX + INCX - X(IX) = X(IX) - TEMP*AP(K) - 70 CONTINUE - END IF - JX = JX + INCX - KK = KK + (N-J+1) - 80 CONTINUE - END IF - END IF - ELSE -* -* Form x := inv( A' )*x or x := inv( conjg( A' ) )*x. -* - IF (LSAME(UPLO,'U')) THEN - KK = 1 - IF (INCX.EQ.1) THEN - DO 110 J = 1,N - TEMP = X(J) - K = KK - IF (NOCONJ) THEN - DO 90 I = 1,J - 1 - TEMP = TEMP - AP(K)*X(I) - K = K + 1 - 90 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK+J-1) - ELSE - DO 100 I = 1,J - 1 - TEMP = TEMP - DCONJG(AP(K))*X(I) - K = K + 1 - 100 CONTINUE - IF (NOUNIT) TEMP = TEMP/DCONJG(AP(KK+J-1)) - END IF - X(J) = TEMP - KK = KK + J - 110 CONTINUE - ELSE - JX = KX - DO 140 J = 1,N - TEMP = X(JX) - IX = KX - IF (NOCONJ) THEN - DO 120 K = KK,KK + J - 2 - TEMP = TEMP - AP(K)*X(IX) - IX = IX + INCX - 120 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK+J-1) - ELSE - DO 130 K = KK,KK + J - 2 - TEMP = TEMP - DCONJG(AP(K))*X(IX) - IX = IX + INCX - 130 CONTINUE - IF (NOUNIT) TEMP = TEMP/DCONJG(AP(KK+J-1)) - END IF - X(JX) = TEMP - JX = JX + INCX - KK = KK + J - 140 CONTINUE - END IF - ELSE - KK = (N* (N+1))/2 - IF (INCX.EQ.1) THEN - DO 170 J = N,1,-1 - TEMP = X(J) - K = KK - IF (NOCONJ) THEN - DO 150 I = N,J + 1,-1 - TEMP = TEMP - AP(K)*X(I) - K = K - 1 - 150 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK-N+J) - ELSE - DO 160 I = N,J + 1,-1 - TEMP = TEMP - DCONJG(AP(K))*X(I) - K = K - 1 - 160 CONTINUE - IF (NOUNIT) TEMP = TEMP/DCONJG(AP(KK-N+J)) - END IF - X(J) = TEMP - KK = KK - (N-J+1) - 170 CONTINUE - ELSE - KX = KX + (N-1)*INCX - JX = KX - DO 200 J = N,1,-1 - TEMP = X(JX) - IX = KX - IF (NOCONJ) THEN - DO 180 K = KK,KK - (N- (J+1)),-1 - TEMP = TEMP - AP(K)*X(IX) - IX = IX - INCX - 180 CONTINUE - IF (NOUNIT) TEMP = TEMP/AP(KK-N+J) - ELSE - DO 190 K = KK,KK - (N- (J+1)),-1 - TEMP = TEMP - DCONJG(AP(K))*X(IX) - IX = IX - INCX - 190 CONTINUE - IF (NOUNIT) TEMP = TEMP/DCONJG(AP(KK-N+J)) - END IF - X(JX) = TEMP - JX = JX - INCX - KK = KK - (N-J+1) - 200 CONTINUE - END IF - END IF - END IF -* - RETURN -* -* End of ZTPSV . -* - END |