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