mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-18 14:34:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

Simplify Rank2Update.h

Browse Source
This commit is contained in:
Chen-Pang He 2012-09-08 22:20:05 +08:00
parent 17c746523e
commit dac5a8a37d
1 changed files with 14 additions and 54 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

68
blas/Rank2Update.h
View File

@ -16,34 +16,16 @@ namespace internal {
* 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>
struct rank2_update_selector
{
static void run(Index size, Scalar* mat, Index stride, const Scalar* _u, const Scalar* _v, Scalar alpha)
static void run(Index size, Scalar* mat, Index stride, const Scalar* u, const Scalar* v, Scalar alpha)
{
Map<const Matrix<Scalar,Dynamic,1> > u(_u, size), v(_v, size);
typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap;
for (Index i=0; i<size; ++i)
{
Map<Matrix<Scalar,Dynamic,1> >(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)
{
Map<const Matrix<Scalar,Dynamic,1> > u(_u, size), v(_v, size);
for (Index i=0; i<size; ++i)
{
Map<Matrix<Scalar,Dynamic,1> >(mat+(stride+1)*i, size-i) +=
conj(alpha) * conj(_u[i]) * v.tail(size-i)
+ alpha * conj(_v[i]) * u.tail(size-i);
Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1)) +=
conj(alpha) * conj(u[i]) * OtherMap(v+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1))
+ alpha * conj(v[i]) * OtherMap(u+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1));
}
}
};
@ -52,42 +34,20 @@ struct rank2_update_selector<Scalar,Index,Lower>
* The matrix is in packed form.
*/
template<typename Scalar, typename Index, int UpLo>
struct packed_rank2_update_selector;
template<typename Scalar, typename Index>
struct packed_rank2_update_selector<Scalar,Index,Upper>
struct packed_rank2_update_selector
{
static void run(Index size, Scalar* mat, const Scalar* _u, const Scalar* _v, Scalar alpha)
static void run(Index size, Scalar* mat, const Scalar* u, const Scalar* v, Scalar alpha)
{
Map<const Matrix<Scalar,Dynamic,1> > u(_u, size), v(_v, size);
typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap;
Index offset = 0;
for (Index i=0; i<size; ++i)
{
Map<Matrix<Scalar,Dynamic,1> >(mat+offset, i+1) +=
conj(alpha) * conj(_u[i]) * v.head(i+1)
+ alpha * conj(_v[i]) * u.head(i+1);
Map<Matrix<Scalar,Dynamic,1> >(mat+offset, UpLo==Lower ? size-i : (i+1)) +=
conj(alpha) * conj(u[i]) * OtherMap(v+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1))
+ alpha * conj(v[i]) * OtherMap(u+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1));
//FIXME This should be handled outside.
mat[offset+i] = real(mat[offset+i]);
offset += i+1;
}
}
};
template<typename Scalar, typename Index>
struct packed_rank2_update_selector<Scalar,Index,Lower>
{
static void run(Index size, Scalar* mat, const Scalar* _u, const Scalar* _v, Scalar alpha)
{
Map<const Matrix<Scalar,Dynamic,1> > u(_u, size), v(_v, size);
Index offset = 0;
for (Index i=0; i<size; ++i)
{
Map<Matrix<Scalar,Dynamic,1> >(mat+offset, size-i) +=
conj(alpha) * conj(_u[i]) * v.tail(size-i)
+ alpha * conj(_v[i]) * u.tail(size-i);
//FIXME This should be handled outside.
mat[offset] = real(mat[offset]);
offset += size-i;
mat[offset+(UpLo==Lower ? 0 : i)] = real(mat[offset+(UpLo==Lower ? 0 : i)]);
offset += UpLo==Lower ? size-i : (i+1);
}
}
};