GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-22 17:49:36 +08:00
Code Issues Packages Projects Releases Wiki Activity

Implement rank-2 update for packed matrices.

Browse Source
This commit is contained in:
Chen-Pang He 2012-09-08 17:29:44 +08:00
parent b5f9bec8ac
commit e4e7585a24
4 changed files with 146 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
blas/CMakeLists.txt
View File

@ -18,10 +18,10 @@ if(EIGEN_Fortran_COMPILER_WORKS)
set(EigenBlas_SRCS ${EigenBlas_SRCS}
complexdots.f
srotm.f srotmg.f drotm.f drotmg.f
lsame.f chpr2.f dspmv.f dtpsv.f ssbmv.f sspr.f stpmv.f
zhpr2.f chbmv.f chpr.f ctpmv.f dspr2.f sspmv.f stpsv.f
zhbmv.f zhpr.f ztpmv.f chpmv.f ctpsv.f dsbmv.f dspr.f dtpmv.f sspr2.f
zhpmv.f ztpsv.f
lsame.f dspmv.f dtpsv.f ssbmv.f sspr.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
zhpmv.f ztpsv.f
dtbmv.f stbmv.f ctbmv.f ztbmv.f
)
else()

54
blas/Rank2Update.h
View File

@ -28,9 +28,8 @@ struct rank2_update_selector<Scalar,Index,Upper>
for (Index i=0; i<size; ++i)
{
Map<PlainVector>(mat+stride*i, i+1) +=
conj(alpha) * conj(_u[i]) * v.head(i+1)
+ alpha * conj(_v[i]) * u.head(i+1);
Map<PlainVector>(mat+stride*i, i+1) += conj(alpha) * conj(_u[i]) * v.head(i+1)
+ alpha * conj(_v[i]) * u.head(i+1);
}
}
};
@ -45,9 +44,52 @@ struct rank2_update_selector<Scalar,Index,Lower>
for (Index i=0; i<size; ++i)
{
Map<PlainVector>(mat+(stride+1)*i, size-i) +=
conj(alpha) * conj(_u[i]) * v.tail(size-i)
+ alpha * conj(_v[i]) * u.tail(size-i);
Map<PlainVector>(mat+(stride+1)*i, size-i) += conj(alpha) * conj(_u[i]) * v.tail(size-i)
+ alpha * conj(_v[i]) * u.tail(size-i);
}
}
};
/* Optimized selfadjoint matrix += alpha * uv' + conj(alpha)*vu'
* 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>
{
static void run(Index size, Scalar* mat, const Scalar* _u, const Scalar* _v, Scalar alpha)
{
typedef Matrix<Scalar,Dynamic,1> PlainVector;
Map<const PlainVector> u(_u, size), v(_v, size);
Index offset = 0;
for (Index i=0; i<size; ++i)
{
offset += i;
Map<PlainVector>(mat+offset, i+1) += conj(alpha) * conj(_u[i]) * v.head(i+1)
+ alpha * conj(_v[i]) * u.head(i+1);
mat[offset+i] = real(mat[offset+i]);
}
}
};
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)
{
typedef Matrix<Scalar,Dynamic,1> PlainVector;
Map<const PlainVector> u(_u, size), v(_v, size);
Index offset = 0;
for (Index i=0; i<size; ++i)
{
Map<PlainVector>(mat+offset, size-i) += conj(alpha) * conj(_u[i]) * v.tail(size-i)
+ alpha * conj(_v[i]) * u.tail(size-i);
mat[offset] = real(mat[offset]);
offset += size-i;
}
}
};

51
blas/level2_cplx_impl.h
View File

@ -120,10 +120,53 @@ int EIGEN_BLAS_FUNC(hemv)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa
* where alpha is a scalar, x and y are n element vectors and A is an
* n by n hermitian matrix, supplied in packed form.
*/
// int EIGEN_BLAS_FUNC(hpr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *x, int *incx, RealScalar *y, int *incy, RealScalar *ap)
// {
// return 1;
// }
int EIGEN_BLAS_FUNC(hpr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pap)
{
typedef void (*functype)(int, Scalar*, const 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::packed_rank2_update_selector<Scalar,int,Upper>::run);
func[LO] = (internal::packed_rank2_update_selector<Scalar,int,Lower>::run);
init = true;
}
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* y = reinterpret_cast<Scalar*>(py);
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;
else if(*incy==0) info = 7;
if(info)
return xerbla_(SCALAR_SUFFIX_UP"HPR2 ",&info,6);
if(alpha==Scalar(0))
return 1;
Scalar* x_cpy = get_compact_vector(x, *n, *incx);
Scalar* y_cpy = get_compact_vector(y, *n, *incy);
int code = UPLO(*uplo);
if(code>=2 || func[code]==0)
return 0;
func[code](*n, ap, x_cpy, y_cpy, alpha);
if(x_cpy!=x) delete[] x_cpy;
if(y_cpy!=y) delete[] y_cpy;
return 1;
}
/** ZHER performs the hermitian rank 1 operation
*

51
blas/level2_real_impl.h
View File

@ -243,10 +243,53 @@ int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px
* where alpha is a scalar, x and y are n element vectors and A is an
* n by n symmetric matrix, supplied in packed form.
*/
// int EIGEN_BLAS_FUNC(spr2)(char *uplo, int *n, RealScalar *alpha, RealScalar *x, int *incx, RealScalar *y, int *incy, RealScalar *ap)
// {
// return 1;
// }
int EIGEN_BLAS_FUNC(spr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pap)
{
typedef void (*functype)(int, Scalar*, const 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::packed_rank2_update_selector<Scalar,int,Upper>::run);
func[LO] = (internal::packed_rank2_update_selector<Scalar,int,Lower>::run);
init = true;
}
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* y = reinterpret_cast<Scalar*>(py);
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;
else if(*incy==0) info = 7;
if(info)
return xerbla_(SCALAR_SUFFIX_UP"SPR2 ",&info,6);
if(alpha==Scalar(0))
return 1;
Scalar* x_cpy = get_compact_vector(x, *n, *incx);
Scalar* y_cpy = get_compact_vector(y, *n, *incy);
int code = UPLO(*uplo);
if(code>=2 || func[code]==0)
return 0;
func[code](*n, ap, x_cpy, y_cpy, alpha);
if(x_cpy!=x) delete[] x_cpy;
if(y_cpy!=y) delete[] y_cpy;
return 1;
}
/** DGER performs the rank 1 operation
*