forked from etmc/tmLQCD
/
eigenvalues.c
138 lines (126 loc) · 4.21 KB
/
eigenvalues.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
/*******************************************************************************
* $Id$
*
* Here we compute the nr_of_eigenvalues lowest eigenvalues
* of (gamma5*D)^2. Therefore we use the arnoldi routines.
*
* The computed eigenvalues are stored in g_eigenvalues
* and the computed eigenvectors in g_ev
*
* inout:
* nr_of_eigenvalues: input: Number of eigenvalues to compute
* output: Number of computed eigenvalues
* input:
* crylov_space_dimension: Dimension of crylov space dimension
* to be used in the arnoldi routines
* operator_flag: Choose if we want to use D_Wilson
* or D_Overlap
*
* Autor: Carsten Urbach <urbach@ifh.de>
*
*******************************************************************************/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "global.h"
#include "su3.h"
#include "linalg_eo.h"
#include "start.h"
#include "tm_operators.h"
#include "solver/solver.h"
#include "solver/jdher.h"
#ifdef MPI
#include "solver/pjdher.h"
#endif
#include "eigenvalues.h"
/*********************************************************
*
* We need here another function Qsqr_psi, representing
* (gamma5*D)^2, because this is used in the CG solver
*
* It is not identical to Q_sqr_psi and not externally
* accessible.
*
*********************************************************/
spinor *eigenvectors = NULL;
double * eigenvls = NULL;
int eigenvalues_for_cg_computed = 0;
void eigenvalues(int * nr_of_eigenvalues, const int operator_flag,
const int max_iterations, const double precision) {
static spinor * eigenvectors_ = NULL;
static int allocated = 0;
spinor *temp_field, *temp_field_ = NULL, *aux_ = NULL;
/**********************
* For Jacobi-Davidson
**********************/
int verbosity = 5, converged = 0, blocksize = 1, blockwise = 0;
int solver_it_max = 50, j_max, j_min;
/*int it_max = 10000;*/
complex *eigv_ = NULL, *eigv;
double decay_min = 1.7, decay_max = 1.5, prec,
threshold_min = 1.e-3, threshold_max = 5.e-2;
static int v0dim = 0;
/**********************
* General variables
**********************/
int returncode=0;
if(g_proc_id == g_stdio_proc) printf("\nNumber of lowest eigenvalues to compute = %d\n\n",(*nr_of_eigenvalues));
eigenvalues_for_cg_computed = 1;
if(g_proc_id == g_stdio_proc) printf("Using Jacobi-Davidson method! \n");
if((*nr_of_eigenvalues) < 8){
j_max = 15;
j_min = 8;
}
else{
j_max = 2*(*nr_of_eigenvalues);
j_min = (*nr_of_eigenvalues);
}
if(precision < 1.e-14){
prec = 1.e-14;
}
else{
prec = precision;
}
g_mu = 0.00343;
/* prec = 1.e-10; */
if(allocated == 0) {
allocated = 1;
#if (defined SSE || defined SSE2 || defined SSE3)
eigenvectors_ = calloc((VOLUMEPLUSRAND)/2*(*nr_of_eigenvalues)+1, sizeof(spinor));
eigenvectors = (spinor *)(((unsigned long int)(eigenvectors_)+ALIGN_BASE)&~ALIGN_BASE);
temp_field_ = calloc((VOLUMEPLUSRAND)/2+1, sizeof(spinor));
temp_field = (spinor *)(((unsigned long int)(temp_field_)+ALIGN_BASE)&~ALIGN_BASE);
#else
eigenvectors_= calloc((VOLUMEPLUSRAND)/2*(*nr_of_eigenvalues), sizeof(spinor));
temp_field_ = calloc((VOLUMEPLUSRAND)/2, sizeof(spinor));
eigenvectors = eigenvectors_;
temp_field = temp_field_;
#endif
eigenvls = (double*)malloc((*nr_of_eigenvalues)*sizeof(double));
}
/* compute minimal eigenvalues */
#ifdef MPI
pjdher((VOLUME)/2*sizeof(spinor)/sizeof(complex), (VOLUMEPLUSRAND)/2*sizeof(spinor)/sizeof(complex),
0., prec,
(*nr_of_eigenvalues), j_max, j_min,
max_iterations, blocksize, blockwise, v0dim, (complex*) eigenvectors,
CG, solver_it_max,
threshold_min, decay_min, verbosity,
&converged, (complex*) eigenvectors, eigenvls,
&returncode, JD_MINIMAL, 1,
&Qtm_pm_psi);
#else
jdher((VOLUME)/2*sizeof(spinor)/sizeof(complex),
0., prec,
(*nr_of_eigenvalues), j_max, j_min,
max_iterations, blocksize, blockwise, v0dim, (complex*) eigenvectors,
BICGSTAB, solver_it_max,
threshold_min, decay_min, verbosity,
&converged, (complex*) eigenvectors, eigenvls,
&returncode, JD_MINIMAL, 1,
&Qtm_pm_psi);
#endif
(*nr_of_eigenvalues) = converged;
v0dim = converged;
free(eigenvls);
}