In the following description, the loop for \(N_R\) is omitted for simplicity and instead of \(G_{i j}(z_k)\), the \(N_L\)-dimensional vector \({\bf x}_{k}\) is obtained by using the library.
The names of the routines is defined as follows.
komega_bicg_init, komega_cocg_init, komega_cg_c_init, komega_cg_r_init
Set the initial conditions such as the allocation of variables used in the library.
komega_bicg_update, komega_cocg_update, komega_cg_c_update, komega_cg_r_update
These routines are called in the iteration to update the solution vectors.
komega_bicg_finalize, komega_cocg_finalize, komega_cg_c_finalize, komega_cg_r_finalize
Release the allocated vectors in the library.
komega_bicg_getcoef, komega_cocg_getcoef, komega_cg_c_getcoef, komega_cg_r_getcoef
Get the \(\alpha\), \(\beta\), \(z_{\rm seed}\), \({\bf r}^{\rm L}\) conserved at each iteration.
komega_bicg_getvec, komega_cocg_getvec, komega_cg_c_getvec, komega_cg_r_getvec
Get the vectors \({\boldsymbol r}\), \({\boldsymbol r}^{\rm old}\), \({\tilde {\boldsymbol r}}\), \({\tilde {\boldsymbol r}}^{\rm old}\).
komega_bicg_restart, komega_cocg_restart, komega_cg_c_restart, CG_R_restart
Note
Give the vector size \(N_H\) corresponding to the size of the Hilbert space and the number of the frequency \(z\).
Allocate the two vectors (in the case of BiCG method, four vectors) with the size of \(N_H\).
Give the function for the Hamiltonian-vector production.
Allocate the solution vectors. It is noted that the length of each solution vector is not always equal to \(N_H\). In fact, the its length in the previous section is \(N_L\). In this case, the length of the (bi-)conjugate gradient vector \({\bf p}_k (k=1,\cdots N_z)\) also becomes \(N_L\). We have to prepare a code for projecting \(N_H\)-dimensional vector onto \(N_L\)dimensional space.
If the result converges (or a breakdown occurs), komega_*_update return the first element of status as a negative integer. Therefore, please exit loop when status(1) < 0 .
The 2-norm is used for the convergence check in the routine komega_*_update. Therefore, if 2-norms of residual vectors at all shift points becomes smaller than threshold, this routine assumes the result is converged.
We can obtain the history of \(\alpha, \beta, {\bf r}^{\rm L}\) for restarting calculation. In this case, itermax must not be 0.
Allocate \({\boldsymbol v}_{1 2}\), \({\boldsymbol v}_{1 3}\), \({\boldsymbol v}_2\), \({\boldsymbol v}_3\), \(\{{\bf x}_k\}, {\bf r}^{\rm L}\) \({\boldsymbol v}_2 = {\boldsymbol \varphi_j}\)
komega_bicg_init(N_H, N_L, N_z, x, z, itermax, threshold) start
Allocate \({\boldsymbol v}_3\), \({\boldsymbol v}_5\), \(\{\pi_k\}\) , \(\{\pi_k^{\rm old}\}\), \(\{{\bf p}_k\}\)
Copy \(\{z_k\}\)
If itermax \(\neq\) 0 , allocate arrays to store \(\alpha\), \(\beta\), and:math:{bf r}^{rm L} at each iteration.
\({\boldsymbol v}_4 = {\boldsymbol v}_2^*\) (an arbitrary vector), \({\boldsymbol v}_3 = {\boldsymbol v}_5 = {\bf 0}\),
\({\bf p}_{k} = {\bf x}_k = {\bf 0}(k=1 \cdots N_z),\; \pi_k=\pi_k^{\rm old} = 1(k=1 \cdots N_z)\)
\(\rho = \infty,\; \alpha = 1,\; z_{\rm seed}=0\)
( \({\boldsymbol v}_2 \equiv {\boldsymbol r}\), \({\boldsymbol v}_3 \equiv {\boldsymbol r}^{\rm old}\), \({\boldsymbol v}_4 \equiv {\tilde {\boldsymbol r}}\), \({\boldsymbol v}_5 \equiv {\tilde {\boldsymbol r}}^{\rm old}\). )
komega_bicg_init finish
do iteration
\({\bf r}^{\rm L} = {\hat P}^\dagger {\boldsymbol v}_2\)
\({\boldsymbol v}_{1 2} = {\hat H} {\boldsymbol v}_2\), \({\boldsymbol v}_{1 4} = {\hat H} {\boldsymbol v}_4\) [Or \(({\boldsymbol v}_{1 2}, {\boldsymbol v}_{1 4}) = {\hat H} ({\boldsymbol v}_2, {\boldsymbol v}_4)\) ]
komega_bicg_update(v_12, v_2, v_14, v_4, x, r_small, status) start
\(\circ\) Seed equation
\(\rho^{\rm old} = \rho,\; \rho = {\boldsymbol v}_4^* \cdot {\boldsymbol v}_2\)
\(\beta = \rho / \rho^{\rm old}\)
\({\boldsymbol v}_{1 2} = z_{\rm seed} {\boldsymbol v}_2 - {\boldsymbol v}_{1 2}\), \({\boldsymbol v}_{1 4} = z_{\rm seed}^* {\boldsymbol v}_4 - {\boldsymbol v}_{1 4}\)
\(\alpha^{\rm old} = \alpha,\; \alpha = \frac{\rho}{{\boldsymbol v}_3^* \cdot {\boldsymbol v}_{1 2} - \beta \rho / \alpha }\)
\(\circ\) Shifted equation
do \(k = 1 \cdots N_z\)
\(\pi_k^{\rm new} = [1+\alpha(z_k-z_{\rm seed})]\pi_k - \frac{\alpha \beta}{\alpha^{\rm old}}(\pi_k^{\rm old} - \pi_k)\)
\({\bf p}_{k} = \frac{1}{\pi_k} {\bf r}^{\rm L} + \frac{\pi^{\rm old}_k \pi^{\rm old}_k}{\pi_k \pi_k} \beta {\bf p}_{k}\)
\({\bf x}_{k} = {\bf x}_{k} + \frac{\pi_k}{\pi_k^{\rm new}} \alpha {\bf p}_{k}\)
\(\pi_k^{\rm old} = \pi_k\), \(\pi_k = \pi_k^{\rm new}\)
end do \(k\)
\({\boldsymbol v}_{1 2} = \left( 1 + \frac{\alpha \beta}{\alpha^{\rm old}} \right) {\boldsymbol v}_2 - \alpha {\boldsymbol v}_{1 2} - \frac{\alpha \beta}{\alpha^{\rm old}} {\boldsymbol v}_3\), \({\boldsymbol v}_3 = {\boldsymbol v}_2,\; {\boldsymbol v}_2 = {\boldsymbol v}_{1 2}\)
\({\boldsymbol v}_{1 4} = \left( 1 + \frac{\alpha^* \beta^*}{\alpha^{{\rm old}*}} \right) {\boldsymbol v}_4 - \alpha^* {\boldsymbol v}_{1 4} - \frac{\alpha^* \beta^*}{\alpha^{{\rm old} *}} {\boldsymbol v}_5\), \({\boldsymbol v}_5 = {\boldsymbol v}_4,\; {\boldsymbol v}_4 = {\boldsymbol v}_{1 4}\)
\(\circ\) Seed switch
Search \(k\) which gives the smallest \(|\pi_k|\) . \(\rightarrow z_{\rm seed},\; \pi_{\rm seed},\; \pi_{\rm seed}^{\rm old}\)
\({\boldsymbol v}_2 = {\boldsymbol v}_2 / \pi_{\rm seed}\), \({\boldsymbol v}_3 = {\boldsymbol v}_3 / \pi_{\rm seed}^{\rm old}\), \({\boldsymbol v}_4 = {\boldsymbol v}_4 / \pi_{\rm seed}^{*}\), \({\boldsymbol v}_5 = {\boldsymbol v}_5 / \pi_{\rm seed}^{\rm old *}\)
\(\alpha = (\pi_{\rm seed}^{\rm old} / \pi_{\rm seed}) \alpha\), \(\rho = \rho / (\pi_{\rm seed}^{\rm old} \pi_{\rm seed}^{\rm old})\)
\(\{\pi_k = \pi_k / \pi_{\rm seed},\; \pi_k^{\rm old} = \pi_k^{\rm old} / \pi_{\rm seed}^{\rm old}\}\)
komega_bicg_update finish
if(status(1) < 0 (This indicates \(|{\boldsymbol v}_2| <\) Threshold)) exit
end do iteration
komega_bicg_finalize start
Deallocate \({\boldsymbol v}_4\), \({\boldsymbol v}_5\), \(\{\pi_k\}\), \(\{\pi_k^{\rm old}\}\), \(\{{\bf p}_k\}\)
komega_bicg_finalize finish
Allocate \({\boldsymbol v}_1\), \({\boldsymbol v}_2\), \(\{{\bf x}_k\}, {\bf r}^{\rm L}\) \({\boldsymbol v}_2 = {\boldsymbol \varphi_j}\)
komega_cocg_init(N_H, N_L, N_z, x, z, itermax, threshold) start
Allocate \({\boldsymbol v}_3\), \(\{\pi_k\}\), \(\{\pi_k^{\rm old}\}\), \(\{{\bf p}_k\}\)
Copy \(\{z_k\}\)
If itermax \(\neq\) 0 , allocate arrays to store \(\alpha\), \(\beta\), and \({\bf r}^{\rm L}\) .
\({\boldsymbol v}_3 = {\bf 0}\),
\({\bf p}_{k} = {\bf x}_k = {\bf 0}(k=1 \cdots N_z),\; \pi_k=\pi_k^{\rm old} = 1(k=1 \cdots N_z)\)
\(\rho = \infty,\; \alpha = 1,\; \beta=0,\; z_{\rm seed}=0\)
( \({\boldsymbol v}_2 \equiv {\boldsymbol r}\), \({\boldsymbol v}_3 \equiv {\boldsymbol r}^{\rm old}\). )
komega_cocg_init finish
do iteration
\({\bf r}^{\rm L} = {\hat P}^\dagger {\boldsymbol v}_2\)
\({\boldsymbol v}_1 = {\hat H} {\boldsymbol v}_2\)
komega_cocg_update(v_1, v_2, x, r_small, status) start
\(\circ\) Seed equationw
\(\rho^{\rm old} = \rho,\; \rho = {\boldsymbol v}_2 \cdot {\boldsymbol v}_2\)
\(\beta = \rho / \rho^{\rm old}\)
\({\boldsymbol v}_1 = z_{\rm seed} {\boldsymbol v}_2 - {\boldsymbol v}_1\)
\(\alpha^{\rm old} = \alpha,\; \alpha = \frac{\rho}{{\boldsymbol v}_2 \cdot {\boldsymbol v}_1 - \beta \rho / \alpha }\)
\(\circ\) Shifted equations
do \(k = 1 \cdots N_z\)
\(\pi_k^{\rm new} = [1+\alpha(z_k-z_{\rm seed})]\pi_k - \frac{\alpha \beta}{\alpha^{\rm old}}(\pi_k^{\rm old} - \pi_k)\)
\({\bf p}_{k} = \frac{1}{\pi_k} {\bf r}^{\rm L} + \frac{\pi^{\rm old}_k \pi^{\rm old}_k}{\pi_k \pi_k} \beta {\bf p}_{k}\)
\({\bf x}_{k} = {\bf x}_{k} + \frac{\pi_k}{\pi_k^{\rm new}} \alpha {\bf p}_{k}\)
\(\pi_k^{\rm old} = \pi_k,\; \pi_k = \pi_k^{\rm new}\)
end do \(k\)
\({\boldsymbol v}_1 = \left( 1 + \frac{\alpha \beta}{\alpha^{\rm old}} \right) {\boldsymbol v}_2 - \alpha {\boldsymbol v}_1 - \frac{\alpha \beta}{\alpha^{\rm old}} {\boldsymbol v}_3\)
\({\boldsymbol v}_3 = {\boldsymbol v}_2\), \({\boldsymbol v}_2 = {\boldsymbol v}_1\)
\(\circ\) Seed switch
Search \(k\) which gives the smallest |pi_k| . \(\rightarrow z_{\rm seed},\; \pi_{\rm seed},\; \pi_{\rm seed}^{\rm old}\)
\({\boldsymbol v}_2 = {\boldsymbol v}_2 / \pi_{\rm seed}\), \({\boldsymbol v}_3 = {\boldsymbol v}_3 / \pi_{\rm seed}^{\rm old}\)
\(\alpha = (\pi_{\rm seed}^{\rm old} / \pi_{\rm seed}) \alpha\), \(\rho = \rho / (\pi_{\rm seed}^{\rm old} \pi_{\rm seed}^{\rm old})\)
\(\{\pi_k = \pi_k / \pi_{\rm seed},\; \pi_k^{\rm old} = \pi_k^{\rm old} / \pi_{\rm seed}^{\rm old}\}\)
komega_cocg_update finish
if(status(1) < 0 (This indicates \(|{\boldsymbol v}_2| <\) Threshold.)) exit
end do iteration
komega_cocg_finalize start
Deallocate \({\boldsymbol v}_3\), \(\{\pi_k\}\), \(\{\pi_k^{\rm old}\}\), \(\{{\bf p}_k\}\)
komega_cocg_finalize finish
The workflow is the same as that of the shifted COCG library.