1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_SELFADJOINTRANK2UPTADE_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_SELFADJOINTRANK2UPTADE_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/* Optimized selfadjoint matrix += alpha * uv' + conj(alpha)*vu' 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * It corresponds to the Level2 syr2 BLAS routine 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar, typename Index, typename UType, typename VType, int UpLo> 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct selfadjoint_rank2_update_selector; 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar, typename Index, typename UType, typename VType> 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct selfadjoint_rank2_update_selector<Scalar,Index,UType,VType,Lower> 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static void run(Scalar* mat, Index stride, const UType& u, const VType& v, const Scalar& alpha) 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = u.size(); 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i=0; i<size; ++i) 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+i, size-i) += 337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez (numext::conj(alpha) * numext::conj(u.coeff(i))) * v.tail(size-i) 347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez + (alpha * numext::conj(v.coeff(i))) * u.tail(size-i); 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar, typename Index, typename UType, typename VType> 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct selfadjoint_rank2_update_selector<Scalar,Index,UType,VType,Upper> 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static void run(Scalar* mat, Index stride, const UType& u, const VType& v, const Scalar& alpha) 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = u.size(); 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i=0; i<size; ++i) 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i, i+1) += 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez (numext::conj(alpha) * numext::conj(u.coeff(i))) * v.head(i+1) 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez + (alpha * numext::conj(v.coeff(i))) * u.head(i+1); 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<bool Cond, typename T> struct conj_expr_if 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : conditional<!Cond, const T&, 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CwiseUnaryOp<scalar_conjugate_op<typename traits<T>::Scalar>,T> > {}; 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, unsigned int UpLo> 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename DerivedU, typename DerivedV> 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathSelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo> 617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez::rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, const Scalar& alpha) 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef internal::blas_traits<DerivedU> UBlasTraits; 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename UBlasTraits::DirectLinearAccessType ActualUType; 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::remove_all<ActualUType>::type _ActualUType; 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::add_const_on_value_type<ActualUType>::type actualU = UBlasTraits::extract(u.derived()); 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef internal::blas_traits<DerivedV> VBlasTraits; 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename VBlasTraits::DirectLinearAccessType ActualVType; 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::remove_all<ActualVType>::type _ActualVType; 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::add_const_on_value_type<ActualVType>::type actualV = VBlasTraits::extract(v.derived()); 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // If MatrixType is row major, then we use the routine for lower triangular in the upper triangular case and 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // vice versa, and take the complex conjugate of all coefficients and vector entries. 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 }; 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar actualAlpha = alpha * UBlasTraits::extractScalarFactor(u.derived()) 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * numext::conj(VBlasTraits::extractScalarFactor(v.derived())); 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (IsRowMajor) 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez actualAlpha = numext::conj(actualAlpha); 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::selfadjoint_rank2_update_selector<Scalar, Index, 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::remove_all<typename internal::conj_expr_if<IsRowMajor ^ UBlasTraits::NeedToConjugate,_ActualUType>::type>::type, 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::remove_all<typename internal::conj_expr_if<IsRowMajor ^ VBlasTraits::NeedToConjugate,_ActualVType>::type>::type, 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (IsRowMajor ? int(UpLo==Upper ? Lower : Upper) : UpLo)> 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::run(_expression().const_cast_derived().data(),_expression().outerStride(),actualU,actualV,actualAlpha); 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_SELFADJOINTRANK2UPTADE_H 94