From 7f57876f8b0d130dc0dcc5c60b76f8ad1ef8b2b3 Mon Sep 17 00:00:00 2001 From: Pavel Kus Date: Wed, 12 Jun 2019 12:18:50 +0200 Subject: [PATCH] replacing the HH transformation by matrix multiplication --- src/elpa2/kernels/blas_block4_template.F90 | 54 +--------------------- 1 file changed, 1 insertion(+), 53 deletions(-) diff --git a/src/elpa2/kernels/blas_block4_template.F90 b/src/elpa2/kernels/blas_block4_template.F90 index e94d3806..7f16a0f9 100644 --- a/src/elpa2/kernels/blas_block4_template.F90 +++ b/src/elpa2/kernels/blas_block4_template.F90 @@ -118,59 +118,7 @@ ! TODO lets see how much does it matter s_mat = matmul(h_mat, transpose(h_mat)) - ! Do the Householder transformations - a_1_1(1:nq) = q(1:nq,4) - a_2_1(1:nq) = q(1:nq,3) - a_3_1(1:nq) = q(1:nq,2) - a_4_1(1:nq) = q(1:nq,1) - - h_2_1 = hh(2,2) - h_3_2 = hh(2,3) - h_3_1 = hh(3,3) - h_4_3 = hh(2,4) - h_4_2 = hh(3,4) - h_4_1 = hh(4,4) - - w_comb(1:nq,4) = a_3_1(1:nq) * h_4_3 + a_4_1(1:nq) - w_comb(1:nq,4) = a_2_1(1:nq) * h_4_2 + w_comb(1:nq,4) - w_comb(1:nq,4) = a_1_1(1:nq) * h_4_1 + w_comb(1:nq,4) - - w_comb(1:nq,3) = a_2_1(1:nq) * h_3_2 + a_3_1(1:nq) - w_comb(1:nq,3) = a_1_1(1:nq) * h_3_1 + w_comb(1:nq,3) - - w_comb(1:nq,2) = a_1_1(1:nq) * h_2_1 + a_2_1(1:nq) - - w_comb(1:nq,1) = a_1_1(1:nq) - - do i=5,nb - h1 = hh(i-3,1) - h2 = hh(i-2,2) - h3 = hh(i-1,3) - h4 = hh(i ,4) - - w_comb(1:nq,1) = w_comb(1:nq,1) + q(1:nq,i) * h1 - w_comb(1:nq,2) = w_comb(1:nq,2) + q(1:nq,i) * h2 - w_comb(1:nq,3) = w_comb(1:nq,3) + q(1:nq,i) * h3 - w_comb(1:nq,4) = w_comb(1:nq,4) + q(1:nq,i) * h4 - enddo - - h1 = hh(nb-2,1) - h2 = hh(nb-1,2) - h3 = hh(nb ,3) - - w_comb(1:nq,1) = w_comb(1:nq,1) + q(1:nq,nb+1) * h1 - w_comb(1:nq,2) = w_comb(1:nq,2) + q(1:nq,nb+1) * h2 - w_comb(1:nq,3) = w_comb(1:nq,3) + q(1:nq,nb+1) * h3 - - h1 = hh(nb-1,1) - h2 = hh(nb ,2) - - w_comb(1:nq,1) = w_comb(1:nq,1) + q(1:nq,nb+2) * h1 - w_comb(1:nq,2) = w_comb(1:nq,2) + q(1:nq,nb+2) * h2 - - h1 = hh(nb,1) - - w_comb(1:nq,1) = w_comb(1:nq,1) + q(1:nq,nb+3) * h1 + w_comb = matmul(q(1:ldq, 1:nb+3), -transpose(h_mat)) ! Rank-1 update tau1 = hh(1,1) -- GitLab