From 3b1e8b7d64c563ef8ea785da8b45da4617c1cc86 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Thu, 25 Mar 2021 21:15:55 +0000 Subject: [PATCH 1/4] Fix default vector return conventions in rank-local SVD --- src/TiledArray/math/linalg/rank-local.cpp | 11 ++++++----- tests/CMakeLists.txt | 12 ++++++------ 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/src/TiledArray/math/linalg/rank-local.cpp b/src/TiledArray/math/linalg/rank-local.cpp index e26ec3da59..c74ee8ed52 100644 --- a/src/TiledArray/math/linalg/rank-local.cpp +++ b/src/TiledArray/math/linalg/rank-local.cpp @@ -143,18 +143,19 @@ template void svd(Matrix& A, std::vector& S, Matrix* U, Matrix* VT) { integer m = A.rows(); integer n = A.cols(); + integer k = std::min(m, n); T* a = A.data(); integer lda = A.rows(); - S.resize(std::min(m, n)); + S.resize(k); T* s = S.data(); auto jobu = lapack::Job::NoVec; T* u = nullptr; integer ldu = m; if (U) { - jobu = lapack::Job::AllVec; - U->resize(m, n); + jobu = lapack::Job::SomeVec; + U->resize(m, k); u = U->data(); ldu = U->rows(); } @@ -163,8 +164,8 @@ void svd(Matrix& A, std::vector& S, Matrix* U, Matrix* VT) { T* vt = nullptr; integer ldvt = n; if (VT) { - jobvt = lapack::Job::AllVec; - VT->resize(n, m); + jobvt = lapack::Job::SomeVec; + VT->resize(k, n); vt = VT->data(); ldvt = VT->rows(); } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index d75c0f8e5f..bab83c0b5c 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -88,11 +88,11 @@ set(ta_test_src_files ta_test.cpp reduce_task.cpp proc_grid.cpp dist_eval_contraction_eval.cpp - expressions.cpp - expressions_sparse.cpp - expressions_complex.cpp - expressions_btas.cpp - expressions_mixed.cpp + #expressions.cpp + #expressions_sparse.cpp + #expressions_complex.cpp + #expressions_btas.cpp + #expressions_mixed.cpp foreach.cpp solvers.cpp initializer_list.cpp @@ -102,7 +102,7 @@ set(ta_test_src_files ta_test.cpp tot_dist_array_part2.cpp random.cpp trace.cpp - tot_expressions.cpp + #tot_expressions.cpp annotation.cpp diagonal_array.cpp contraction_helpers.cpp From 343cc88a17b51d4773a20303047967ff6f69ebce Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Thu, 25 Mar 2021 21:44:45 +0000 Subject: [PATCH 2/4] Extended the flexibility of the rank-local SVD wrapper to allow for all allowed combinations of vector Job options provided by LAPACK --- src/TiledArray/math/linalg/rank-local.cpp | 55 ++++++++++++++++++----- src/TiledArray/math/linalg/rank-local.h | 11 ++++- 2 files changed, 55 insertions(+), 11 deletions(-) diff --git a/src/TiledArray/math/linalg/rank-local.cpp b/src/TiledArray/math/linalg/rank-local.cpp index c74ee8ed52..3cc10a99c8 100644 --- a/src/TiledArray/math/linalg/rank-local.cpp +++ b/src/TiledArray/math/linalg/rank-local.cpp @@ -140,7 +140,7 @@ void heig(Matrix& A, Matrix& B, std::vector& W) { } template -void svd(Matrix& A, std::vector& S, Matrix* U, Matrix* VT) { +void svd(Job jobu, Job jobvt, Matrix& A, std::vector& S, Matrix* U, Matrix* VT) { integer m = A.rows(); integer n = A.cols(); integer k = std::min(m, n); @@ -150,6 +150,7 @@ void svd(Matrix& A, std::vector& S, Matrix* U, Matrix* VT) { S.resize(k); T* s = S.data(); +#if 0 auto jobu = lapack::Job::NoVec; T* u = nullptr; integer ldu = m; @@ -169,6 +170,40 @@ void svd(Matrix& A, std::vector& S, Matrix* U, Matrix* VT) { vt = VT->data(); ldvt = VT->rows(); } +#else + T* u = nullptr; + T* vt = nullptr; + integer ldu = 1, ldvt = 1; + if( (jobu == Job::SomeVec or jobu == Job::AllVec) and (not U) ) + TA_LAPACK_ERROR("Requested out-of-place right singular vectors with null U input"); + if( (jobvt == Job::SomeVec or jobvt == Job::AllVec) and (not VT) ) + TA_LAPACK_ERROR("Requested out-of-place left singular vectors with null VT input"); + + if( jobu == Job::SomeVec ) { + U->resize(m, k); + u = U->data(); + ldu = m; + } + + if( jobu == Job::AllVec ) { + U->resize(m, m); + u = U->data(); + ldu = m; + } + + if( jobvt == Job::SomeVec ) { + VT->resize(k, n); + vt = VT->data(); + ldvt = k; + } + + if( jobvt == Job::AllVec ) { + VT->resize(n, n); + vt = VT->data(); + ldvt = n; + } + +#endif TA_LAPACK(gesvd, jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt); } @@ -195,15 +230,15 @@ void lu_inv(Matrix& A) { TA_LAPACK(getri, n, a, lda, ipiv.data()); } -#define TA_LAPACK_EXPLICIT(MATRIX, VECTOR) \ - template void cholesky(MATRIX&); \ - template void cholesky_linv(MATRIX&); \ - template void cholesky_solve(MATRIX&, MATRIX&); \ - template void cholesky_lsolve(Op, MATRIX&, MATRIX&); \ - template void heig(MATRIX&, VECTOR&); \ - template void heig(MATRIX&, MATRIX&, VECTOR&); \ - template void svd(MATRIX&, VECTOR&, MATRIX*, MATRIX*); \ - template void lu_solve(MATRIX&, MATRIX&); \ +#define TA_LAPACK_EXPLICIT(MATRIX, VECTOR) \ + template void cholesky(MATRIX&); \ + template void cholesky_linv(MATRIX&); \ + template void cholesky_solve(MATRIX&, MATRIX&); \ + template void cholesky_lsolve(Op, MATRIX&, MATRIX&); \ + template void heig(MATRIX&, VECTOR&); \ + template void heig(MATRIX&, MATRIX&, VECTOR&); \ + template void svd(Job,Job,MATRIX&, VECTOR&, MATRIX*, MATRIX*); \ + template void lu_solve(MATRIX&, MATRIX&); \ template void lu_inv(MATRIX&); TA_LAPACK_EXPLICIT(Matrix, std::vector); diff --git a/src/TiledArray/math/linalg/rank-local.h b/src/TiledArray/math/linalg/rank-local.h index 144c844e3c..87d4b44a56 100644 --- a/src/TiledArray/math/linalg/rank-local.h +++ b/src/TiledArray/math/linalg/rank-local.h @@ -10,6 +10,8 @@ namespace TiledArray::math::linalg::rank_local { +using Job = ::lapack::Job; + template using Matrix = ::Eigen::Matrix; @@ -35,7 +37,14 @@ template void heig(Matrix &A, Matrix &B, std::vector &W); template -void svd(Matrix &A, std::vector &S, Matrix *U, Matrix *VT); +void svd(Job jobu, Job jobvt, Matrix &A, std::vector &S, Matrix *U, Matrix *VT); + +template +void svd(Matrix &A, std::vector &S, Matrix *U, Matrix *VT) { + svd( U ? Job::SomeVec : Job::NoVec, + VT ? Job::SomeVec : Job::NoVec, + A, S, U, VT ); +} template void lu_solve(Matrix &A, Matrix &B); From f3e04149388a36f5e25278be39e028366f535d95 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Thu, 25 Mar 2021 21:49:24 +0000 Subject: [PATCH 3/4] Cleanup of old SVD code --- src/TiledArray/math/linalg/rank-local.cpp | 23 ----------------------- 1 file changed, 23 deletions(-) diff --git a/src/TiledArray/math/linalg/rank-local.cpp b/src/TiledArray/math/linalg/rank-local.cpp index 3cc10a99c8..5413ad0e5f 100644 --- a/src/TiledArray/math/linalg/rank-local.cpp +++ b/src/TiledArray/math/linalg/rank-local.cpp @@ -150,27 +150,6 @@ void svd(Job jobu, Job jobvt, Matrix& A, std::vector& S, Matrix* U, Mat S.resize(k); T* s = S.data(); -#if 0 - auto jobu = lapack::Job::NoVec; - T* u = nullptr; - integer ldu = m; - if (U) { - jobu = lapack::Job::SomeVec; - U->resize(m, k); - u = U->data(); - ldu = U->rows(); - } - - auto jobvt = lapack::Job::NoVec; - T* vt = nullptr; - integer ldvt = n; - if (VT) { - jobvt = lapack::Job::SomeVec; - VT->resize(k, n); - vt = VT->data(); - ldvt = VT->rows(); - } -#else T* u = nullptr; T* vt = nullptr; integer ldu = 1, ldvt = 1; @@ -203,8 +182,6 @@ void svd(Job jobu, Job jobvt, Matrix& A, std::vector& S, Matrix* U, Mat ldvt = n; } -#endif - TA_LAPACK(gesvd, jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt); } From cd5f534f67aca4a2dea54877a98c4aaaeba1605e Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Thu, 25 Mar 2021 15:12:22 -0700 Subject: [PATCH 4/4] Re-enable expressions tests --- tests/CMakeLists.txt | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index bab83c0b5c..d75c0f8e5f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -88,11 +88,11 @@ set(ta_test_src_files ta_test.cpp reduce_task.cpp proc_grid.cpp dist_eval_contraction_eval.cpp - #expressions.cpp - #expressions_sparse.cpp - #expressions_complex.cpp - #expressions_btas.cpp - #expressions_mixed.cpp + expressions.cpp + expressions_sparse.cpp + expressions_complex.cpp + expressions_btas.cpp + expressions_mixed.cpp foreach.cpp solvers.cpp initializer_list.cpp @@ -102,7 +102,7 @@ set(ta_test_src_files ta_test.cpp tot_dist_array_part2.cpp random.cpp trace.cpp - #tot_expressions.cpp + tot_expressions.cpp annotation.cpp diagonal_array.cpp contraction_helpers.cpp