ground/gcs/src/libs/eigen/blas/Rank2Update.h

   1 // This file is part of Eigen, a lightweight C++ template library
   2 // for linear algebra.
   3 //
   4 // Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
   5 //
   6 // This Source Code Form is subject to the terms of the Mozilla
   7 // Public License v. 2.0. If a copy of the MPL was not distributed
   8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
   9
  10 #ifndef EIGEN_RANK2UPDATE_H
  11 #define EIGEN_RANK2UPDATE_H
  12
  13 namespace internal {
  14
  15 /* Optimized selfadjoint matrix += alpha * uv' + conj(alpha)*vu'
  16  * This is the low-level version of SelfadjointRank2Update.h
  17  */
  18 template<typename Scalar, typename Index, int UpLo>
  19 struct rank2_update_selector
  20 {
  21   static void run(Index size, Scalar* mat, Index stride, const Scalar* u, const Scalar* v, Scalar alpha)
  22   {
  23     typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap;
  24     for (Index i=0; i<size; ++i)
  25     {
  26       Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1)) +=
  27       numext::conj(alpha) * numext::conj(u[i]) * OtherMap(v+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1))
  28                 + alpha * numext::conj(v[i]) * OtherMap(u+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1));
  29     }
  30   }
  31 };
  32
  33 /* Optimized selfadjoint matrix += alpha * uv' + conj(alpha)*vu'
  34  * The matrix is in packed form.
  35  */
  36 template<typename Scalar, typename Index, int UpLo>
  37 struct packed_rank2_update_selector
  38 {
  39   static void run(Index size, Scalar* mat, const Scalar* u, const Scalar* v, Scalar alpha)
  40   {
  41     typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap;
  42     Index offset = 0;
  43     for (Index i=0; i<size; ++i)
  44     {
  45       Map<Matrix<Scalar,Dynamic,1> >(mat+offset, UpLo==Lower ? size-i : (i+1)) +=
  46       numext::conj(alpha) * numext::conj(u[i]) * OtherMap(v+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1))
  47                 + alpha * numext::conj(v[i]) * OtherMap(u+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1));
  48       //FIXME This should be handled outside.
  49       mat[offset+(UpLo==Lower ? 0 : i)] = numext::real(mat[offset+(UpLo==Lower ? 0 : i)]);
  50       offset += UpLo==Lower ? size-i : (i+1);
  51     }
  52   }
  53 };
  54
  55 } // end namespace internal
  56
  57 #endif // EIGEN_RANK2UPDATE_H