From 47b243ef08f4101de3d936f2373343e67eaa4dd5 Mon Sep 17 00:00:00 2001 From: "Field G. Van Zee" Date: Wed, 23 Jul 2014 13:41:13 -0500 Subject: [PATCH] Call setid for early return from herk/her2k. Details: - Added setid call (to zero imaginary parts of diagonal elements) to early return branches of herk_front() and her2k_front() for cases where alpha is zero. Thanks to Murtaza Ali for suggesting this fix. - Comment update. --- frame/3/her2k/bli_her2k_front.c | 11 +++++++++-- frame/3/herk/bli_herk_front.c | 11 +++++++++-- 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/frame/3/her2k/bli_her2k_front.c b/frame/3/her2k/bli_her2k_front.c index 141c5a17d..dd0ed5321 100644 --- a/frame/3/her2k/bli_her2k_front.c +++ b/frame/3/her2k/bli_her2k_front.c @@ -52,10 +52,12 @@ void bli_her2k_front( obj_t* alpha, if ( bli_error_checking_is_enabled() ) bli_her2k_check( alpha, a, b, beta, c ); - // If alpha is zero, scale by beta and return. + // If alpha is zero, scale by beta, zero the imaginary components of + // the diagonal elements, and return. if ( bli_obj_equals( alpha, &BLIS_ZERO ) ) { bli_scalm( beta, c ); + bli_setid( &BLIS_ZERO, c ); return; } @@ -139,7 +141,12 @@ void bli_her2k_front( obj_t* alpha, #endif - // Explicitly set diagonal elements' imaginary components to zero. + // The Hermitian rank-2k product was computed as A*B'+B*A', even for + // the diagonal elements. Mathematically, the imaginary components of + // diagonal elements of a Hermitian rank-2k product should always be + // zero. However, in practice, they sometimes accumulate meaningless + // non-zero values. To prevent this, we explicitly set those values + // to zero before returning. bli_setid( &BLIS_ZERO, &c_local ); } diff --git a/frame/3/herk/bli_herk_front.c b/frame/3/herk/bli_herk_front.c index 6733ba310..0b0b23dcd 100644 --- a/frame/3/herk/bli_herk_front.c +++ b/frame/3/herk/bli_herk_front.c @@ -48,10 +48,12 @@ void bli_herk_front( obj_t* alpha, if ( bli_error_checking_is_enabled() ) bli_herk_check( alpha, a, beta, c ); - // If alpha is zero, scale by beta and return. + // If alpha is zero, scale by beta, zero the imaginary components of + // the diagonal elements, and return. if ( bli_obj_equals( alpha, &BLIS_ZERO ) ) { bli_scalm( beta, c ); + bli_setid( &BLIS_ZERO, c ); return; } @@ -93,7 +95,12 @@ void bli_herk_front( obj_t* alpha, bli_herk_thrinfo_free_paths( infos, n_threads ); - // Explicitly set diagonal elements' imaginary components to zero. + // The Hermitian rank-k product was computed as A*A', even for the + // diagonal elements. Mathematically, the imaginary components of + // diagonal elements of a Hermitian rank-k product should always be + // zero. However, in practice, they sometimes accumulate meaningless + // non-zero values. To prevent this, we explicitly set those values + // to zero before returning. bli_setid( &BLIS_ZERO, &c_local ); }