...
 
Commits (4)
...@@ -81,21 +81,11 @@ ...@@ -81,21 +81,11 @@
real(kind=C_DATATYPE_KIND), intent(in) :: hh(1:ldh,1:6) real(kind=C_DATATYPE_KIND), intent(in) :: hh(1:ldh,1:6)
#endif #endif
real(kind=C_DATATYPE_KIND) :: h_2_1, h_3_2, h_3_1, h_4_3, h_4_2, h_4_1
real(kind=C_DATATYPE_KIND) :: a_1_1(nq), a_2_1(nq), a_3_1(nq), a_4_1(nq)
real(kind=C_DATATYPE_KIND) :: h1, h2, h3, h4
real(kind=C_DATATYPE_KIND) :: w(nq), z(nq), x(nq), y(nq)
real(kind=C_DATATYPE_KIND) :: w_orig(nq), z_orig(nq), x_orig(nq), y_orig(nq)
real(kind=C_DATATYPE_KIND) :: w_comb(nq, 4) real(kind=C_DATATYPE_KIND) :: w_comb(nq, 4)
real(kind=C_DATATYPE_KIND) :: h_comb(4)
real(kind=C_DATATYPE_KIND) :: h_mat(4, nb+3) real(kind=C_DATATYPE_KIND) :: h_mat(4, nb+3)
real(kind=C_DATATYPE_KIND) :: h4m(4, 4)
real(kind=C_DATATYPE_KIND) :: s_mat(4, 4) real(kind=C_DATATYPE_KIND) :: s_mat(4, 4)
real(kind=C_DATATYPE_KIND) :: tau1, tau2, tau3, tau4 integer(kind=ik) :: i, j, k
integer(kind=ik) :: i
! Calculate dot product of the two Householder vectors ! Calculate dot product of the two Householder vectors
...@@ -112,48 +102,33 @@ ...@@ -112,48 +102,33 @@
h_mat(3,3:nb+1) = -hh(2:nb, 3) h_mat(3,3:nb+1) = -hh(2:nb, 3)
h_mat(4,2:nb) = -hh(2:nb, 4) h_mat(4,2:nb) = -hh(2:nb, 4)
! TODO we actually need just the strictly upper triangle of s_mat ! TODO we do not need the diagonal, but how to do it with BLAS?
! TODO take care when changing to BLAS !s_mat = - matmul(h_mat, transpose(h_mat))
! TODO we do not even need diagonal, which might not be achievable by blas. call DSYRK('L', 'N', 4, nb+3, &
! TODO lets see how much does it matter -1.0_rk8, h_mat, 4, &
s_mat = matmul(h_mat, transpose(h_mat)) 0.0_rk8, s_mat, 4)
w_comb = matmul(q(1:ldq, 1:nb+3), -transpose(h_mat)) s_mat(1,1) = 1
s_mat(2,2) = 1
! Rank-1 update s_mat(3,3) = 1
tau1 = hh(1,1) s_mat(4,4) = 1
tau2 = hh(1,2)
tau3 = hh(1,3) ! w_comb = - matmul(q(1:nq, 1:nb+3), transpose(h_mat))
tau4 = hh(1,4) call DGEMM('N', 'T', nq, 4, nb+3, &
-1.0_rk8, q, ldq, &
! x_orig = x h_mat, 4, &
! y_orig = y 0.0_rk8, w_comb, nq)
! z_orig = z
! w_orig = w ! Rank-1 update
do i = 1, 4
h4m = 0.0 w_comb(1:nq,i) = matmul(w_comb(1:nq,1:i), hh(1,i) * s_mat(i,1:i))
enddo
h4m(1,1) = tau1
w_comb(1:nq,1) = w_comb(1:nq,1) * h4m(1,1) !q(1:nq, 1:nb+3) = matmul(w_comb, h_mat) + q(1:nq, 1:nb+3)
call DGEMM('N', 'N', nq, nb+3, 4, &
h4m(2,1) = - tau2 * s_mat(1,2) 1.0_rk8, w_comb, nq, &
h4m(2,2) = tau2 h_mat, 4, &
w_comb(1:nq,2) = w_comb(1:nq,1) * h4m(2,1) + w_comb(1:nq,2) * h4m(2,2) 1.0_rk8, q, ldq)
h4m(3,1) = - tau3 * s_mat(1,3)
h4m(3,2) = - tau3 * s_mat(2,3)
h4m(3,3) = tau3
w_comb(1:nq,3) = w_comb(1:nq,1) * h4m(3,1) + w_comb(1:nq,2) * h4m(3,2) + w_comb(1:nq,3) * h4m(3,3)
h4m(4,1) = - tau4 * s_mat(1,4)
h4m(4,2) = - tau4 * s_mat(2,4)
h4m(4,3) = - tau4 * s_mat(3,4)
h4m(4,4) = tau4
w_comb(1:nq,4) = w_comb(1:nq,1) * h4m(4,1) + w_comb(1:nq,2) * h4m(4,2) + w_comb(1:nq,3) * h4m(4,3) + w_comb(1:nq,4) * h4m(4,4)
q(1:nq, 1:nb+3) = matmul(w_comb, h_mat) + q(1:nq, 1:nb+3)
end subroutine end subroutine
......