devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
Chen-Pang He
2012-09-10 18:25:30 +08:00
parent 65caa40a3d
commit 04f315d692
7 changed files with 53 additions and 863 deletions
Download Patch File Download Diff File Expand all files Collapse all files

47
blas/level2_cplx_impl.h
View File

@@ -108,10 +108,49 @@ int EIGEN_BLAS_FUNC(hemv)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa
* where alpha is a real scalar, x is an n element vector and A is an
* n by n hermitian matrix, supplied in packed form.
*/
// int EIGEN_BLAS_FUNC(hpr)(char *uplo, int *n, RealScalar *alpha, RealScalar *x, int *incx, RealScalar *ap)
// {
// return 1;
// }
int EIGEN_BLAS_FUNC(hpr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pap)
{
typedef void (*functype)(int, Scalar*, const Scalar*, RealScalar);
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,ColMajor,Upper,false,Conj>::run);
func[LO] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run);
init = true;
}
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* ap = reinterpret_cast<Scalar*>(pap);
RealScalar alpha = *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"HPR ",&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;
}
/** ZHPR2 performs the hermitian rank 2 operation
*