Commit 8489f58b authored by Yingzhou Li's avatar Yingzhou Li Committed by Victor Yu

Added omm estimated upper bound in function and modified documentation accordingly

parent d80714fa
No preview for this file type
......@@ -741,6 +741,10 @@ holds information for conducting operations in ELSI-RCI.
& \verbb+total_energy+ & real double & null
& Total energy of last iteration \\
\toprule
\multirow{1}{*}{\shortstack{Solver (OMM)}}
& \verbb+omm_est_ub+ & real double & null
& Estimated upper bound for the spectrum \\
\midrule
\multirow{3}{*}{\shortstack{Solver\\(PPCG)}}
& \verbb+ppcg_sbsize+ & integer & $\min(8,\verbb+n_state+)$
& PPCG block size \\
......@@ -1145,7 +1149,7 @@ $A$ and store the results in $B$ as,
\subsubsection{Remark:} If the overlapping matrix is Hermitian,
then we only need to implement $B \leftarrow SA$.
\subsubsection{Required by solver:} All
\subsubsection{Required by solver:} All (if not \verbb+ovlp_is_unit+)
\subsection{RCI\_P\_MULTI}
......@@ -1498,27 +1502,30 @@ tasks.
\verbb+ovlp_is_unit, tol_iter, verbose+
\subsubsection{Required task:} \RCINULL, \RCISTOP, \RCICONVERGE,
\RCIALLOCATE, \RCIDEALLOCATE, \RCIHMULTI, \RCISMULTI,
\RCIALLOCATE, \RCIDEALLOCATE, \RCIHMULTI,
\RCIPMULTI, \RCICOPY, \RCISUBCOPY, \RCISUBCOL, \RCICOLSCALE, \RCIAXPY,
\RCICOLNORM, \RCIGEMM, \RCIHEEV / \RCIHEGV, \RCIPOTRF, \RCITRSM
\subsubsection{Other task:} None
\subsubsection{Other task:}
\RCISMULTI (if not \verbb+ovlp_is_unit+)
\subsection{OMM}
\subsubsection{Required configuration:} \verbb+solver+ ==
\verbb+RCI_SOLVER_OMM+
\verbb+RCI_SOLVER_OMM+, \verbb+omm_est_ub+
\subsubsection{Optional configuration:} \verbb+n_basis, n_state,+
\verbb+ovlp_is_unit, tol_iter, verbose+
\subsubsection{Required task:} \RCINULL, \RCISTOP, \RCICONVERGE,
\RCIALLOCATE, \RCIDEALLOCATE, \RCIHMULTI, \RCISMULTI,
\RCIALLOCATE, \RCIDEALLOCATE, \RCIHMULTI,
\RCIPMULTI, \RCICOPY, \RCITRACE, \RCIDOT, \RCIAXPY,
\RCIGEMM
\subsubsection{Other task:} \RCIHEEV / \RCIHEGV(if eigenvectors is needed)
\subsubsection{Other task:}
\RCISMULTI (if not \verbb+ovlp_is_unit+),
\RCIHEEV / \RCIHEGV(if eigenvectors is needed)
\subsection{PPCG}
......@@ -1531,11 +1538,12 @@ tasks.
\verbb+ppcg_tol_lock+
\subsubsection{Required task:} \RCINULL, \RCISTOP, \RCICONVERGE,
\RCIALLOCATE, \RCIDEALLOCATE, \RCIHMULTI, \RCISMULTI, \RCIPMULTI,
\RCIALLOCATE, \RCIDEALLOCATE, \RCIHMULTI, \RCIPMULTI,
\RCICOPY, \RCISUBCOPY, \RCISUBCOL, \RCISUBROW, \RCITRACE, \RCICOLSCALE,
\RCIAXPY, \RCICOLNORM, \RCIGEMM, \RCIHEEV / \RCIHEGV, \RCIPOTRF, \RCITRSM
\subsubsection{Other task:} None
\subsubsection{Other task:}
\RCISMULTI (if not \verbb+ovlp_is_unit+)
\subsection{ChebFilter}
......@@ -1547,11 +1555,12 @@ tasks.
\verbb+ovlp_is_unit, tol_iter, verbose, cheb_max_inneriter+
\subsubsection{Required task:} \RCINULL, \RCISTOP, \RCICONVERGE,
\RCIALLOCATE, \RCIDEALLOCATE, \RCIHMULTI, \RCISMULTI, \RCIPMULTI,
\RCIALLOCATE, \RCIDEALLOCATE, \RCIHMULTI, \RCIPMULTI,
\RCICOPY, \RCISUBCOPY, \RCISUBCOL, \RCISCALE,
\RCIAXPY, \RCIGEMM
\subsubsection{Other task:} None
\subsubsection{Other task:}
\RCISMULTI (if not \verbb+ovlp_is_unit+)
......
......@@ -195,6 +195,7 @@ contains
integer, save :: icg ! CG step num.
real(r8), save :: est_ub ! estimated upper bound of spectrum
real(r8), save :: lambda ! CG step size
real(r8), save :: lambda_d ! lambda denominator
real(r8), save :: lambda_n ! lambda numerator
......@@ -239,6 +240,7 @@ contains
n = r_h%n_state
max_n = r_h%max_n
ovlp_is_unit = r_h%ovlp_is_unit
est_ub = r_h%omm_est_ub
icg = 0
lambda = 0.0_r8
call rci_op_null(task)
......@@ -250,27 +252,39 @@ contains
ijob = ijob + 1
return
end if
! -- HW = C'*HC
! -calculate the overlap matrix in WF basis: SW=C^T*S*C
! -- SC = S*C
if (ijob == SID_INIT + 2) then
call rci_op_gemm(iS, task, 'C', 'N', n, n, m, 1.0_r8, &
MID_C, m, MID_HC, m, 0.0_r8, MID_HW, max_n)
if (ovlp_is_unit) then
call rci_op_null(task)
else
call rci_op_s_multi(iS, task, 'N', m, n, MID_C, MID_SC)
end if
ijob = ijob + 1
return
end if
! -calculate the overlap matrix in WF basis: SW=C^T*S*C
! -- SC = S*C
if (ijob == SID_INIT + 3) then
if (ovlp_is_unit) then
call rci_op_null(task)
call rci_op_axpy(iS, task, m, n, -est_ub, MID_C, m, MID_HC, m)
else
call rci_op_s_multi(iS, task, 'N', m, n, MID_C, MID_SC)
call rci_op_axpy(iS, task, m, n, -est_ub, MID_SC, m, MID_HC, m)
end if
ijob = ijob + 1
return
end if
! -- SW = C'*SC
! -- HW = C'*HC
if (ijob == SID_INIT + 4) then
call rci_op_gemm(iS, task, 'C', 'N', n, n, m, 1.0_r8, &
MID_C, m, MID_HC, m, 0.0_r8, MID_HW, max_n)
ijob = ijob + 1
return
end if
! -- SW = C'*SC
if (ijob == SID_INIT + 5) then
if (ovlp_is_unit) then
call rci_op_gemm(iS, task, 'C', 'N', n, n, m, 1.0_r8, &
MID_C, m, MID_C, m, 0.0_r8, MID_SW, max_n)
......@@ -285,20 +299,20 @@ contains
! -calculate the gradient: G=2*(2*HC-SC*HW-HC*SW)
! Note: It is required that G = 0 as the initial
! -- G = 4*HC
if (ijob == SID_INIT + 5) then
if (ijob == SID_INIT + 6) then
call rci_op_axpy(iS, task, m, n, 4.0_r8, MID_HC, m, MID_G, m)
ijob = ijob + 1
return
end if
! -- G = -2 HC*SW + G
if (ijob == SID_INIT + 6) then
if (ijob == SID_INIT + 7) then
call rci_op_gemm(iS, task, 'N', 'N', m, n, n, -2.0_r8, &
MID_HC, m, MID_SW, max_n, 1.0_r8, MID_G, m)
ijob = ijob + 1
return
end if
! -- G = -2 SC*HW + G
if (ijob == SID_INIT + 7) then
if (ijob == SID_INIT + 8) then
if (ovlp_is_unit) then
call rci_op_gemm(iS, task, 'N', 'N', m, n, n, -2.0_r8, &
MID_C, m, MID_HW, max_n, 1.0_r8, MID_G, m)
......@@ -312,14 +326,14 @@ contains
! -calculate the preconditioned gradient by premultiplying G by P
! -- PG = P*G
if (ijob == SID_INIT + 8) then
if (ijob == SID_INIT + 9) then
call rci_op_p_multi(iS, task, 'N', m, n, MID_G, MID_PG)
ijob = ijob + 1
return
end if
! - calculate HWd = G'*H*C = PG'*HC
if (ijob == SID_INIT + 9) then
if (ijob == SID_INIT + 10) then
call rci_op_gemm(iS, task, 'C', 'N', n, n, m, 1.0_r8, &
MID_PG, m, MID_HC, m, 0.0_r8, MID_HWd, max_n)
ijob = ijob + 1
......@@ -327,7 +341,7 @@ contains
end if
! - calculate SWd = G'*S*C = PG'*SC
if (ijob == SID_INIT + 10) then
if (ijob == SID_INIT + 11) then
if (ovlp_is_unit) then
call rci_op_gemm(iS, task, 'C', 'N', n, n, m, 1.0_r8, &
MID_PG, m, MID_C, m, 0.0_r8, MID_SWd, max_n)
......@@ -341,22 +355,11 @@ contains
! - calculate HWdd = G'*H*G = PG'*(H*PG)
! -- HD = H*PG
if (ijob == SID_INIT + 11) then
call rci_op_h_multi(iS, task, 'N', m, n, MID_PG, MID_HD)
ijob = ijob + 1
return
end if
! -- HWdd = PG'*HD
if (ijob == SID_INIT + 12) then
call rci_op_gemm(iS, task, 'C', 'N', n, n, m, 1.0_r8, &
MID_PG, m, MID_HD, m, 0.0_r8, &
MID_HWdd, max_n)
call rci_op_h_multi(iS, task, 'N', m, n, MID_PG, MID_HD)
ijob = ijob + 1
return
end if
! SWdd=G^T*S*G
! - calculate SWdd = G'*S*G = PG'*(S*PG)
! -- SD = S*PG
if (ijob == SID_INIT + 13) then
if (ovlp_is_unit) then
......@@ -367,8 +370,30 @@ contains
ijob = ijob + 1
return
end if
! -- SWdd = PG'*SD
if (ijob == SID_INIT + 14) then
if (ovlp_is_unit) then
call rci_op_axpy(iS, task, m, n, -est_ub, MID_PG, m, MID_HD, m)
else
call rci_op_axpy(iS, task, m, n, -est_ub, MID_SD, m, MID_HD, m)
end if
ijob = ijob + 1
return
end if
! -- HWdd = PG'*HD
if (ijob == SID_INIT + 15) then
call rci_op_gemm(iS, task, 'C', 'N', n, n, m, 1.0_r8, &
MID_PG, m, MID_HD, m, 0.0_r8, &
MID_HWdd, max_n)
ijob = ijob + 1
return
end if
! SWdd=G^T*S*G
! - calculate SWdd = G'*S*G = PG'*(S*PG)
! -- SWdd = PG'*SD
if (ijob == SID_INIT + 16) then
if (ovlp_is_unit) then
call rci_op_gemm(iS, task, 'C', 'N', n, n, m, 1.0_r8, &
MID_PG, m, MID_PG, m, 0.0_r8, &
......@@ -579,29 +604,41 @@ contains
ijob = ijob + 1
return
end if
! -- HWdd = D'*HD
! -- SD = S*D
if (ijob == SID_ITER + 8) then
call rci_op_gemm(iS, task, 'C', 'N', n, n, m, 1.0_r8, &
MID_D, m, MID_HD, m, 0.0_r8, &
MID_HWdd, max_n)
if (ovlp_is_unit) then
call rci_op_null(task)
else
call rci_op_s_multi(iS, task, 'N', m, n, MID_D, MID_SD)
end if
ijob = ijob + 1
return
end if
! SWdd=D'*S*D
! - calculate SWdd = D'*S*D
! -- SD = S*D
if (ijob == SID_ITER + 9) then
if (ovlp_is_unit) then
call rci_op_null(task)
call rci_op_axpy(iS, task, m, n, -est_ub, MID_D, m, MID_HD, m)
else
call rci_op_s_multi(iS, task, 'N', m, n, MID_D, MID_SD)
call rci_op_axpy(iS, task, m, n, -est_ub, MID_SD, m, MID_HD, m)
end if
ijob = ijob + 1
return
end if
! -- SWdd = D'*SD
! -- HWdd = D'*HD
if (ijob == SID_ITER + 10) then
call rci_op_gemm(iS, task, 'C', 'N', n, n, m, 1.0_r8, &
MID_D, m, MID_HD, m, 0.0_r8, &
MID_HWdd, max_n)
ijob = ijob + 1
return
end if
! SWdd=D'*S*D
! - calculate SWdd = D'*S*D
! -- SWdd = D'*SD
if (ijob == SID_ITER + 11) then
if (ovlp_is_unit) then
call rci_op_gemm(iS, task, 'C', 'N', n, n, m, 1.0_r8, &
MID_D, m, MID_D, m, 0.0_r8, &
......@@ -861,6 +898,7 @@ contains
! - Psi = Psi*G
if (ijob == SID_FINISH + 4) then
resvec(1:n) = resvec(1:n) + est_ub
r_h%total_energy = sum(resvec(1:n))
if (r_h%verbose > 0) then
write(*,*) 'Eigenvalues are:'
......
......@@ -32,6 +32,9 @@ module ELSI_RCI_DATATYPE
integer(i4) :: total_iter
real(r8) :: total_energy
! OMM
real(r8) :: omm_est_ub
! PPCG
integer(i4) :: ppcg_sbsize
integer(i4) :: ppcg_rrstep
......@@ -42,6 +45,7 @@ module ELSI_RCI_DATATYPE
real(r8) :: cheb_est_ub
integer(i4) :: cheb_max_inneriter
end type
! RCI Instruction type
......
......@@ -42,6 +42,9 @@ contains
r_h%verbose = verbose
! OMM default
r_h%omm_est_ub = 20
! PPCG default
r_h%ppcg_tol_lock = tol_iter/10
r_h%ppcg_sbsize = min(8,n_state)
......
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