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