...

Commits (5)
 ... ... @@ -91,7 +91,7 @@ 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-4) real(kind=C_DATATYPE_KIND) :: h_mat(4, nb+3) real(kind=C_DATATYPE_KIND) :: tau1, tau2, tau3, tau4 #endif /* REALCASE==1 */ ... ... @@ -115,14 +115,21 @@ #endif /* COMPLEXCASE==1 */ integer(kind=ik) :: i print *, "SIMPLE BLOCK4, nb, nq, ldq, ldh ", nb, nq, ldq, ldh print *, "Q:", q(1:ldq,1:nb+3) print *, "HH:", hh(1:ldh,1:6) ! Calculate dot product of the two Householder vectors #if REALCASE==1 h_mat(:,:) = 0.0 h_mat(1,4) = -1.0 h_mat(2,3) = -1.0 h_mat(3,2) = -1.0 h_mat(4,1) = -1.0 h_mat(1,5:nb+3) = -hh(2:nb, 1) h_mat(2,4:nb+2) = -hh(2:nb, 2) h_mat(3,3:nb+1) = -hh(2:nb, 3) h_mat(4,2:nb) = -hh(2:nb, 4) s_1_2 = hh(2,2) s_1_3 = hh(3,3) s_2_3 = hh(2,3) ... ... @@ -152,15 +159,7 @@ s_1_4 = s_1_4 + hh(i-3,1) * hh(i,4) enddo #endif #if COMPLEXCASE==1 stop !s = conjg(hh(2,2))*1.0 !do i=3,nb ! s = s+(conjg(hh(i,2))*hh(i-1,1)) !enddo #endif ! Do the Householder transformations a_1_1(1:nq) = q(1:nq,4) ... ... @@ -175,7 +174,6 @@ h_4_2 = hh(3,4) h_4_1 = hh(4,4) #if REALCASE == 1 w(1:nq) = a_3_1(1:nq) * h_4_3 + a_4_1(1:nq) w(1:nq) = a_2_1(1:nq) * h_4_2 + w(1:nq) w(1:nq) = a_1_1(1:nq) * h_4_1 + w(1:nq) ... ... @@ -186,25 +184,12 @@ y(1:nq) = a_1_1(1:nq) * h_2_1 + a_2_1(1:nq) x(1:nq) = a_1_1(1:nq) #endif #if COMPLEXCASE==1 stop !y(1:nq) = q(1:nq,1) + q(1:nq,2)*conjg(hh(2,2)) #endif do i=5,nb #if REALCASE == 1 h1 = hh(i-3,1) h2 = hh(i-2,2) h3 = hh(i-1,3) h4 = hh(i ,4) #endif #if COMPLEXCASE==1 stop ! h1 = conjg(hh(i-1,1)) ! h2 = conjg(hh(i,2)) #endif x(1:nq) = x(1:nq) + q(1:nq,i) * h1 y(1:nq) = y(1:nq) + q(1:nq,i) * h2 ... ... @@ -216,16 +201,9 @@ h2 = hh(nb-1,2) h3 = hh(nb ,3) #if REALCASE==1 x(1:nq) = x(1:nq) + q(1:nq,nb+1) * h1 y(1:nq) = y(1:nq) + q(1:nq,nb+1) * h2 z(1:nq) = z(1:nq) + q(1:nq,nb+1) * h3 #endif #if COMPLEXCASE==1 stop !x(1:nq) = x(1:nq) + q(1:nq,nb+1)*conjg(hh(nb,1)) #endif h1 = hh(nb-1,1) h2 = hh(nb ,2) ... ... @@ -237,7 +215,6 @@ x(1:nq) = x(1:nq) + q(1:nq,nb+3) * h1 ! Rank-1 update tau1 = hh(1,1) tau2 = hh(1,2) ... ... @@ -270,69 +247,13 @@ w(1:nq) = w(1:nq) * h1 - ( z(1:nq) * h4 + y(1:nq) * h3 + x(1:nq) * h2) q(1:nq,1) = q(1:nq,1) - w(1:nq) h4 = hh(2,4) q(1:nq,2) = q(1:nq,2) - (w(1:nq) * h4 + z(1:nq)) h3 = hh(2,3) h4 = hh(3,4) q(1:nq,3) = q(1:nq,3) - y(1:nq) q(1:nq,3) = -( z(1:nq) * h3) + q(1:nq,3) q(1:nq,3) = -( w(1:nq) * h4) + q(1:nq,3) h2 = hh(2,2) h3 = hh(3,3) h4 = hh(4,4) q(1:nq,4) = q(1:nq,4) - x(1:nq) q(1:nq,4) = -(y(1:nq) * h2) + q(1:nq,4) q(1:nq,4) = -(z(1:nq) * h3) + q(1:nq,4) q(1:nq,4) = -(w(1:nq) * h4) + q(1:nq,4) w_comb(:,1) = x w_comb(:,2) = y w_comb(:,3) = z w_comb(:,4) = w ! do i=5,nb ! h1 = hh(i-3,1) ! h2 = hh(i-2,2) ! h3 = hh(i-1,3) ! h4 = hh(i ,4) ! ! h_comb = (/-h1, -h2, -h3, -h4/) ! ! q(1:nq,i) = matmul(w_comb, h_comb) + q(1:nq,i) ! ! enddo !h_mat(:,:) = 0.0 h_mat(1,1:nb-4) = -hh(2:nb-3, 1) h_mat(2,1:nb-4) = -hh(3:nb-2, 2) h_mat(3,1:nb-4) = -hh(4:nb-1, 3) h_mat(4,1:nb-4) = -hh(5:nb, 4) q(1:nq, 5:nb) = matmul(w_comb, h_mat) + q(1:nq, 5:nb) h1 = hh(nb-2,1) h2 = hh(nb-1,2) h3 = hh(nb ,3) q(1:nq,nb+1) = -(x(1:nq) * h1) + q(1:nq,nb+1) q(1:nq,nb+1) = -(y(1:nq) * h2) + q(1:nq,nb+1) q(1:nq,nb+1) = -(z(1:nq) * h3) + q(1:nq,nb+1) h1 = hh(nb-1,1) h2 = hh(nb ,2) q(1:nq,nb+2) = - (x(1:nq) * h1) + q(1:nq,nb+2) q(1:nq,nb+2) = - (y(1:nq) * h2) + q(1:nq,nb+2) h1 = hh(nb,1) q(1:nq,nb+3) = - (x(1:nq) * h1) + q(1:nq,nb+3) q(1:nq, 1:nb+3) = matmul(w_comb, h_mat) + q(1:nq, 1:nb+3) end subroutine ... ...
 ... ... @@ -114,8 +114,8 @@ !TODO remove print *, "SIMPLE BLOCK4, nb, nq, ldq, ldh ", nb, nq, ldq, ldh print *, "Q:", q(1:ldq,1:nb+3) print *, "HH:", hh(1:ldh,1:6) !print *, "Q:", q(1:ldq,1:nb+3) !print *, "HH:", hh(1:ldh,1:6) ! call the blas kernel for future comparison ... ...