[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] [getfem-commits] branch master updated: Make gmm qr_fac
From: |
Konstantinos Poulios |
Subject: |
[Getfem-commits] [getfem-commits] branch master updated: Make gmm qr_factor more consistent with lapack for real square matrices |
Date: |
Wed, 20 Dec 2023 18:29:15 -0500 |
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 <logari81@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));
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] [getfem-commits] branch master updated: Make gmm qr_factor more consistent with lapack for real square matrices,
Konstantinos Poulios <=