mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-11-27 06:30:28 +08:00
Code Issues Packages Projects Releases Wiki Activity

Implement rank-1 update for self-adjoint packed matrices.

Browse Source
This commit is contained in:
Chen-Pang He 2012-09-08 22:51:40 +08:00
parent dac5a8a37d
commit 1b8f416408
4 changed files with 93 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
blas/CMakeLists.txt
View File

@ -18,9 +18,9 @@ if(EIGEN_Fortran_COMPILER_WORKS)
set(EigenBlas_SRCS ${EigenBlas_SRCS}
complexdots.f
srotm.f srotmg.f drotm.f drotmg.f
lsame.f dspmv.f dtpsv.f ssbmv.f sspr.f stpmv.f
lsame.f dspmv.f dtpsv.f ssbmv.f stpmv.f
chbmv.f chpr.f ctpmv.f sspmv.f stpsv.f
zhbmv.f zhpr.f ztpmv.f chpmv.f ctpsv.f dsbmv.f dspr.f dtpmv.f
zhbmv.f zhpr.f ztpmv.f chpmv.f ctpsv.f dsbmv.f dtpmv.f
zhpmv.f ztpsv.f
dtbmv.f stbmv.f ctbmv.f ztbmv.f
)

47
blas/SelfadjointPackedProduct.h Normal file
View File

@ -0,0 +1,47 @@
// 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_SELFADJOINT_PACKED_PRODUCT_H
#define EIGEN_SELFADJOINT_PACKED_PRODUCT_H
namespace internal {
/* Optimized matrix += alpha * uv'
* The matrix is in packed form.
*
* FIXME I always fail tests for complex self-adjoint matrices.
*
******* FATAL ERROR - PARAMETER NUMBER 6 WAS CHANGED INCORRECTLY *******
******* xHPR FAILED ON CALL NUMBER:
2: xHPR ('U', 1, 0.0, X, 1, AP)
*/
template<typename Scalar, typename Index, int UpLo>
struct selfadjoint_packed_rank1_update
{
static void run(Index size, Scalar* mat, const Scalar* vec, Scalar alpha)
{
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, UpLo==Lower ? size-i : (i+1))
+= alpha * conj(vec[i]) * OtherMap(vec+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1));
//FIXME This should be handled outside.
mat[offset+(UpLo==Lower ? 0 : i)] = real(mat[offset+(UpLo==Lower ? 0 : i)]);
offset += UpLo==Lower ? size-i : (i+1);
}
}
};
//TODO struct selfadjoint_packed_product_selector
} // end namespace internal
#endif // EIGEN_SELFADJOINT_PACKED_PRODUCT_H

1
blas/common.h
View File

@ -76,6 +76,7 @@ namespace Eigen {
#include "BandTriangularSolver.h"
#include "GeneralRank1Update.h"
#include "Rank2Update.h"
#include "SelfadjointPackedProduct.h"
}
using namespace Eigen;

47
blas/level2_real_impl.h
View File

@ -231,10 +231,49 @@ int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px
* where alpha is a real scalar, x is an n element vector and A is an
* n by n symmetric matrix, supplied in packed form.
*/
// int EIGEN_BLAS_FUNC(spr)(char *uplo, int *n, Scalar *alpha, Scalar *x, int *incx, Scalar *ap)
// {
// return 1;
// }
int EIGEN_BLAS_FUNC(spr)(char *uplo, int *n, Scalar *palpha, Scalar *px, int *incx, Scalar *pap)
{
typedef void (*functype)(int, 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::selfadjoint_packed_rank1_update<Scalar,int,Upper>::run);
func[LO] = (internal::selfadjoint_packed_rank1_update<Scalar,int,Lower>::run);
init = true;
}
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* ap = reinterpret_cast<Scalar*>(pap);
Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
int info = 0;
if(UPLO(*uplo)==INVALID) info = 1;
else if(*n<0) info = 2;
else if(*incx==0) info = 5;
if(info)
return xerbla_(SCALAR_SUFFIX_UP"SPR ",&info,6);
if(alpha==Scalar(0))
return 1;
Scalar* x_cpy = get_compact_vector(x, *n, *incx);
int code = UPLO(*uplo);
if(code>=2 || func[code]==0)
return 0;
func[code](*n, ap, x_cpy, alpha);
if(x_cpy!=x) delete[] x_cpy;
return 1;
}
/** DSPR2 performs the symmetric rank 2 operation
*