This is an automated email from the git hooks/post-receive script. logari81 pushed a commit to branch master in repository getfem.
The following commit(s) were added to refs/heads/master by this push: new 27cc7174 Make gmm qr_factor more consistent with lapack for real square matrices 27cc7174 is described below commit 27cc7174e4ab7187452e331fa320e46158f93d54 Author: Konstantinos Poulios <logar...@gmail.com> AuthorDate: Thu Dec 21 00:28:59 2023 +0100 Make gmm qr_factor more consistent with lapack for real square matrices - lapack skips the last reflection in dlargf.f for real square matrices - fixes make_gmm_test fails when compiled with the lapack interface --- src/gmm/gmm_dense_qr.h | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/src/gmm/gmm_dense_qr.h b/src/gmm/gmm_dense_qr.h index 3b266868..4bdde208 100644 --- a/src/gmm/gmm_dense_qr.h +++ b/src/gmm/gmm_dense_qr.h @@ -54,8 +54,11 @@ namespace gmm { GMM_ASSERT2(m >= n, "dimensions mismatch"); std::vector<T> W(m), V(m); - - for (size_type j = 0; j < n; ++j) { + // skip the last reflection for *real* square matrices (m-1 limit) + // lapack does the same in dlarfg.f but not in zlarfg.f + const size_type jmax = (m==n && !is_complex(T())) ? m-1 + : n; + for (size_type j = 0; j < jmax; ++j) { sub_interval SUBI(j, m-j), SUBJ(j, n-j); V.resize(m-j); W.resize(n-j); @@ -80,7 +83,10 @@ namespace gmm { if (m == 0) return; std::vector<T> V(m), W(mat_nrows(A)); V[0] = T(1); - for (size_type j = 0; j < n; ++j) { + // assume that the last reflection was skipped for *real* square matrices + const size_type jmax = (m==n && !is_complex(T())) ? m-1 + : n; + for (size_type j = 0; j < jmax; ++j) { V.resize(m-j); for (size_type i = j+1; i < m; ++i) V[i-j] = QR(i, j); @@ -102,7 +108,10 @@ namespace gmm { if (m == 0) return; std::vector<T> V(m), W(mat_ncols(A)); V[0] = T(1); - for (size_type j = 0; j < n; ++j) { + // assume that the last reflection was skipped for *real* square matrices + const size_type jmax = (m==n && !is_complex(T())) ? m-1 + : n; + for (size_type j = 0; j < jmax; ++j) { V.resize(m-j); for (size_type i = j+1; i < m; ++i) V[i-j] = QR(i, j); row_house_update(sub_matrix(A, sub_interval(j, m-j), @@ -123,8 +132,11 @@ namespace gmm { std::vector<T> W(m); dense_matrix<T> VV(m, n); - - for (size_type j = 0; j < n; ++j) { + // skip the last reflection for *real* square matrices (m-1 limit) + // lapack does the same in dlarfg.f but not in zlarfg.f + const size_type jmax = (m==n && !is_complex(T())) ? m-1 + : n; + for (size_type j = 0; j < jmax; ++j) { sub_interval SUBI(j, m-j), SUBJ(j, n-j); for (size_type i = j; i < m; ++i) VV(i,j) = Q(i, j); @@ -136,8 +148,10 @@ namespace gmm { gmm::copy(sub_matrix(Q, sub_interval(0, n), sub_interval(0, n)), R); gmm::copy(identity_matrix(), Q); - - for (size_type j = n-1; j != size_type(-1); --j) { + // assume that the last reflection was skipped for *real* square matrices + const size_type jstart = (m==n && !is_complex(T())) ? m-2 + : n-1; + for (size_type j = jstart; j != size_type(-1); --j) { sub_interval SUBI(j, m-j), SUBJ(j, n-j); row_house_update(sub_matrix(Q, SUBI, SUBJ), sub_vector(mat_col(VV,j), SUBI), sub_vector(W, SUBJ));