Commit 9a5c2a5b authored by Pavel Kus's avatar Pavel Kus
Browse files

Added check of B-orthogonality of the generalized eigenvectors

previously the orthogonality check has been skipped for generalized EVP
parent 504027c8
...@@ -186,16 +186,28 @@ ...@@ -186,16 +186,28 @@
endif endif
! 2. Eigenvector orthogonality ! 2. Eigenvector orthogonality
!TODO for the generalized eigenvector problem, the orthogonality test has to be altered if(present(bs)) then
!TODO at the moment, it is skipped !for the generalized EVP, the eigenvectors should be B-orthogonal, not orthogonal
if(.not. present(bs)) then ! tmp2 = B * Z
! tmp1 = Z**T * Z tmp2(:,:) = 0.0_rck
#ifdef WITH_MPI
call scal_PRECISION_GEMM('N', 'N', na, nev, na, ONE, bs, 1, 1, &
sc_desc, z, 1, 1, sc_desc, ZERO, tmp2, 1, 1, sc_desc)
#else /* WITH_MPI */
call PRECISION_GEMM('N','N', na, nev, na, ONE, bs, na, z, na, ZERO, tmp2, na)
#endif /* WITH_MPI */
else
tmp2(:,:) = z(:,:)
endif
! tmp1 = Z**T * tmp2
! actually tmp1 = Z**T * Z for standard case and tmp1 = Z**T * B * Z for generalized
tmp1 = 0 tmp1 = 0
#ifdef WITH_MPI #ifdef WITH_MPI
call scal_PRECISION_GEMM(BLAS_TRANS_OR_CONJ, 'N', nev, nev, na, ONE, z, 1, 1, & call scal_PRECISION_GEMM(BLAS_TRANS_OR_CONJ, 'N', nev, nev, na, ONE, z, 1, 1, &
sc_desc, z, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc) sc_desc, tmp2, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
#else /* WITH_MPI */ #else /* WITH_MPI */
call PRECISION_GEMM(BLAS_TRANS_OR_CONJ,'N',nev,nev,na,ONE,z,na,z,na,ZERO,tmp1,na) call PRECISION_GEMM(BLAS_TRANS_OR_CONJ,'N',nev,nev,na,ONE,z,na,tmp2,na,ZERO,tmp1,na)
#endif /* WITH_MPI */ #endif /* WITH_MPI */
! First check, whether the elements on diagonal are 1 .. "normality" of the vectors ! First check, whether the elements on diagonal are 1 .. "normality" of the vectors
err = 0.0_rk err = 0.0_rk
...@@ -231,16 +243,16 @@ ...@@ -231,16 +243,16 @@
errmax = err errmax = err
#endif /* WITH_MPI */ #endif /* WITH_MPI */
if (myid==0) print *,'Error Orthogonality:',errmax if (myid==0) print *,'Error Orthogonality:',errmax
if (nev .ge. 2) then
if (errmax .gt. tol_orth .or. errmax .eq. 0.0_rk) then if (nev .ge. 2) then
status = 1 if (errmax .gt. tol_orth .or. errmax .eq. 0.0_rk) then
endif status = 1
else
if (errmax .gt. tol_orth) then
status = 1
endif
endif endif
endif ! skiping test of orthogonality for generalized eigenproblem else
if (errmax .gt. tol_orth) then
status = 1
endif
endif
end function end function
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment