CHROMA
ritz_array.h
Go to the documentation of this file.
1 #ifndef __ritz_array_h__
2 #define __ritz_array_h__
3 
4 #include "chromabase.h"
5 #include "linearop.h"
6 
7 namespace Chroma {
8 
9 //! Minimizes the Ritz functional with a CG based algorithm
10 /*!
11  * \ingroup eig
12  *
13  * This subroutine minimizes the Ritz functional with a CG based
14  * algorithm to find the n-th lowest eigenvalue of a hermitian A
15 
16  * lambda_n = min_z <z|Az>/<z|z>
17 
18  * In the subspace orthogonal to the (n-1) lower eigenvectors of A
19 
20  * This routine can handle the original version or be a part of a
21  * Kalkreuter Simma algorithm.
22  *
23  * The algorithm has been modified that there is now only one target
24  * residue -- the relative one.
25  *
26  * so the original stopping criterion is now:
27  *
28  * || g || < Rsd_r * lambda
29  *
30  * The two key features for Kalkreuter - Simma over the normal procedure
31  * are:
32  * i) Adaptive termination of the procedure (controlled from caller)
33  *
34  * The iteration will stop if the gradient has decreased by at least
35  * a given factor of gamma^{-1} which according to the paper is O(10)
36  *
37  * ii) The initial error estimate is now the "delta_cycle" criterion.
38  * basically this scales the norm of || g^2 ||. So the iteration will
39  * also stop if delta_cycle || gamma^2 || < Rsd_r mu.
40  *
41  * To run in non K-S mode, set gamma = 1, delta_cycle = 1.
42  *
43  * For details of the Kalkreuter Simma algorithm see
44  * hep-lat/9507023
45  *
46  * Algorithm:
47 
48  * Apsi[0] := A . Psi[0] ;
49  * mu[0] := < Psi[0] | A . Psi[0] > ; Initial Ritz value
50  * Note: we assume < Psi[0] |Psi[0] > = 1
51  * p[1] = g[0] := ( A - mu[0] ) Psi[0] ; Initial direction
52  * if converged then return
53  * FOR k FROM 1 TO MaxCG DO CG iterations
54  * s1 = (mu[k-1] + < p[k] | A p[k] >/p[k]|^2) / 2;
55  * s2 = (mu[k-1] - < p[k] | A p[k] >/p[k]|^2) / 2;
56  * s3 = |g[k-1]|^2 / |p[k]|;
57  * Compute a and cos(theta), sin(theta)
58  * (see DESY internal report, September 1994, by Bunk, Jansen,
59  * Luescher and Simma)
60  * lambda = mu[k] = mu[k-1] - 2 a sin^2(theta);
61  * Psi[k] = cos(theta) Psi[k-1] + sin(theta) p[k]/|p[k]|;
62  * Apsi[k] = cos(theta) Apsi[k-1] + sin(theta) A p[k]/|p[k]|;
63  * g[k] = Apsi[k] - mu[k] Psi[k]
64  *
65  * if now converged, then exit
66  *
67  * b[k+1] := cos(theta) |g[k]|^2 / |g[k-1]|^2;
68  * p[k+1] := g[k] + b[k+1] (p[k] - Psi[k] < Psi[k] | p[k] >); New direction
69 
70  * Arguments:
71 
72  * A The hermitian Matrix as a lin op (Read)
73  * lambda Eigenvalue to find (Write)
74  * Psi_all Eigenvectors (Modify)
75  * N_eig Current Eigenvalue number (Read)
76  * RsdR_r (relative) residue (Read)
77  * n_renorm Renormalize every n_renorm iter.(Read)
78  * n_min Minimum no of CG iters to do (Read)
79  * n_max Maximum no of CG iters to do (Read)
80  * n_count Number of CG iteration (done) (Write)
81  * Kalk_Sim Use Kalkreuter-Simma criterion (Read)
82  * delta_cycle The initial error estimate (Read)
83  * gamma "Convergence factor" gamma required (Read)
84  * ProjApsiP flag for projecting A.psi (Read)
85 
86  * Local Variables:
87 
88  * psi New eigenstd::vector
89  * p Direction std::vector
90  * Apsi Temporary for A.psi
91  * Ap Temporary for A.p, and other
92  * mu Ritz functional value
93  * g2 | g[k] |**2
94  * p2 | p[k] |**2
95  * k CG iteration counter
96  * b beta[k+1]
97  * and some others...
98  */
99 
100 template < typename T >
101 void Ritz_t(const LinearOperatorArray<T>& A, // Herm Pos Def
102  Real& lambda, // Current E-value
103  multi2d<T>& psi_all, // E-std::vector array
104  int N_eig, // Current evec index
105  const Real& Rsd_r, // Target relative residue
106  const Real& Rsd_a, // Target absolute residue
107  const Real& zero_cutoff, // If ev slips below this we consider
108  // if to be zero
109  int n_renorm, // Renormalise frequency
110  int n_min, int n_max, // minimum / maximum no of iters to do
111  int MaxCG, // Max iters after which we bomb
112  bool ProjApsiP, // Project A (?) -- user option
113  int& n_count, // No of iters actually taken
114  Real& final_grad, // Final gradient
115  bool Kalk_Sim, // Are we in Kalk Simma mode?
116  const Real& delta_cycle, // Initial error estimate (KS mode)
117  const Real& gamma_factor) ; // Convergence factor Gamma
118 
119 
120 void Ritz(const LinearOperatorArray<LatticeFermion>& A, // Herm Pos Def
121  Real& lambda, // Current E-value
122  multi2d<LatticeFermion>& psi_all, // E-std::vector array
123  int N_eig, // Current evec index
124  const Real& Rsd_r, // Target relative residue
125  const Real& Rsd_a, // Target absolute residue
126  const Real& zero_cutoff, // if ev slips below this we consider
127  // if to be zero.
128  int n_renorm, // Renormalise frequency
129  int n_min, int n_max, // minimum / maximum no of iters to do
130  int MaxCG, // Maximum no if iters after which we bomb
131  bool ProjApsiP, // Project A (?) -- user option
132  int& n_count, // No of iters actually taken
133  Real& final_grad, // Final Gradient
134  bool Kalk_Sim, // Are we in Kalk Simma mode?
135  const Real& delta_cycle, // Initial error estimate (KS mode)
136  const Real& gamma_factor) ; // Convergence factor Gamma
137 
138 } // end namespace Chroma
139 
140 #endif
Primary include file for CHROMA library code.
Linear Operator to arrays.
Definition: linearop.h:61
void Ritz_t(const LinearOperator< T > &A, Real &lambda, multi1d< T > &psi_all, int N_eig, const Real &Rsd_r, const Real &Rsd_a, const Real &zero_cutoff, int n_renorm, int n_min, int n_max, int MaxCG, bool ProjApsiP, int &n_count, Real &final_grad, bool Kalk_Sim, const Real &delta_cycle, const Real &gamma_factor)
Minimizes the Ritz functional with a CG based algorithm.
Definition: ritz.cc:103
Linear Operators.
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
int n_count
Definition: invbicg.cc:78
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > enum InvType invType const multi1d< Real > int MaxCG
Definition: pbg5p_w.cc:32
void Ritz(const LinearOperator< LatticeFermion > &A, Real &lambda, multi1d< LatticeFermion > &psi_all, int N_eig, const Real &Rsd_r, const Real &Rsd_a, const Real &zero_cutoff, int n_renorm, int n_min, int n_max, int MaxCG, bool ProjApsiP, int &n_count, Real &final_grad, bool Kalk_Sim, const Real &delta_cycle, const Real &gamma_factor)
Definition: ritz.cc:559
A(A, psi, r, Ncb, PLUS)