1 /* This file is part of ELPA.
2
3 The ELPA library was originally created by the ELPA consortium,
4 consisting of the following organizations:
5
6 - Max Planck Computing and Data Facility (MPCDF), formerly known as
7 Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
8 - Bergische Universität Wuppertal, Lehrstuhl für angewandte
9 Informatik,
10 - Technische Universität München, Lehrstuhl für Informatik mit
11 Schwerpunkt Wissenschaftliches Rechnen ,
12 - Fritz-Haber-Institut, Berlin, Abt. Theorie,
13 - Max-Plack-Institut für Mathematik in den Naturwissenschaften,
14 Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
15 and
16 - IBM Deutschland GmbH
17
18
19 More information can be found here:
20 http://elpa.mpcdf.mpg.de/
21
22 ELPA is free software: you can redistribute it and/or modify
23 it under the terms of the version 3 of the license of the
24 GNU Lesser General Public License as published by the Free
25 Software Foundation.
26
27 ELPA is distributed in the hope that it will be useful,
28 but WITHOUT ANY WARRANTY; without even the implied warranty of
29 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 GNU Lesser General Public License for more details.
31
32 You should have received a copy of the GNU Lesser General Public License
33 along with ELPA. If not, see <http://www.gnu.org/licenses/>
34
35 ELPA reflects a substantial effort on the part of the original
36 ELPA consortium, and we ask you to respect the spirit of the
37 license that we chose: i.e., please contribute any changes you
38 may have back to the original ELPA library distribution, and keep
39 any derivatives of ELPA under the same license that we chose for
40 the original distribution, the GNU Lesser General Public License.
41 */
42
43 #include "config.h"
44
45 #include <stdio.h>
46 #include <stdlib.h>
47 #include <string.h>
48 #ifdef WITH_MPI
49 #include <mpi.h>
50 #endif
51 #include <math.h>
52
53 #include <elpa/elpa.h>
54 #include <assert.h>
55
56 #include "test/shared/generated.h"
57
58 #if !(defined(TEST_REAL) ^ defined(TEST_COMPLEX))
59 #error "define exactly one of TEST_REAL or TEST_COMPLEX"
60 #endif
61
62 #if !(defined(TEST_SINGLE) ^ defined(TEST_DOUBLE))
63 #error "define exactly one of TEST_SINGLE or TEST_DOUBLE"
64 #endif
65
66 #if !(defined(TEST_SOLVER_1STAGE) ^ defined(TEST_SOLVER_2STAGE))
67 #error "define exactly one of TEST_SOLVER_1STAGE or TEST_SOLVER_2STAGE"
68 #endif
69
70 #ifdef TEST_GENERALIZED_DECOMP_EIGENPROBLEM
71 #define TEST_GENERALIZED_EIGENPROBLEM
72 #endif
73
74 #ifdef TEST_SINGLE
75 # define EV_TYPE float
76 # ifdef TEST_REAL
77 # define MATRIX_TYPE float
78 # define PREPARE_MATRIX_RANDOM prepare_matrix_random_real_single_f
79 # define PREPARE_MATRIX_RANDOM_SPD prepare_matrix_random_spd_real_single_f
80 # define CHECK_CORRECTNESS_EVP_NUMERIC_RESIDUALS check_correctness_evp_numeric_residuals_real_single_f
81 # define CHECK_CORRECTNESS_EVP_GEN_NUMERIC_RESIDUALS check_correctness_evp_gen_numeric_residuals_real_single_f
82 # else
83 # define MATRIX_TYPE complex float
84 # define PREPARE_MATRIX_RANDOM prepare_matrix_random_complex_single_f
85 # define PREPARE_MATRIX_RANDOM_SPD prepare_matrix_random_spd_complex_single_f
86 # define CHECK_CORRECTNESS_EVP_NUMERIC_RESIDUALS check_correctness_evp_numeric_residuals_complex_single_f
87 # define CHECK_CORRECTNESS_EVP_GEN_NUMERIC_RESIDUALS check_correctness_evp_gen_numeric_residuals_complex_single_f
88 # endif
89 #else
90 # define EV_TYPE double
91 # ifdef TEST_REAL
92 # define MATRIX_TYPE double
93 # define PREPARE_MATRIX_RANDOM prepare_matrix_random_real_double_f
94 # define PREPARE_MATRIX_RANDOM_SPD prepare_matrix_random_spd_real_double_f
95 # define CHECK_CORRECTNESS_EVP_NUMERIC_RESIDUALS check_correctness_evp_numeric_residuals_real_double_f
96 # define CHECK_CORRECTNESS_EVP_GEN_NUMERIC_RESIDUALS check_correctness_evp_gen_numeric_residuals_real_double_f
97 # else
98 # define MATRIX_TYPE complex double
99 # define PREPARE_MATRIX_RANDOM prepare_matrix_random_complex_double_f
100 # define PREPARE_MATRIX_RANDOM_SPD prepare_matrix_random_spd_complex_double_f
101 # define CHECK_CORRECTNESS_EVP_NUMERIC_RESIDUALS check_correctness_evp_numeric_residuals_complex_double_f
102 # define CHECK_CORRECTNESS_EVP_GEN_NUMERIC_RESIDUALS check_correctness_evp_gen_numeric_residuals_complex_double_f
103 # endif
104 #endif
105
106 #define assert_elpa_ok(x) assert(x == ELPA_OK)
107
108
main(int argc,char ** argv)109 int main(int argc, char** argv) {
110 /* matrix dimensions */
111 int na, nev, nblk;
112
113 /* mpi */
114 int myid, nprocs;
115 int na_cols, na_rows;
116 int np_cols, np_rows;
117 int my_prow, my_pcol;
118 int mpi_comm;
119 int provided_mpi_thread_level;
120
121 /* blacs */
122 int my_blacs_ctxt, sc_desc[9], info;
123
124 /* The Matrix */
125 MATRIX_TYPE *a, *as, *z, *b, *bs;
126 EV_TYPE *ev;
127
128 int error, status;
129
130 elpa_t handle;
131
132 int value;
133 #ifdef WITH_MPI
134 #ifndef WITH_OPENMP
135 MPI_Init(&argc, &argv);
136 #else
137 MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &provided_mpi_thread_level);
138
139 if (provided_mpi_thread_level != MPI_THREAD_MULTIPLE) {
140 fprintf(stderr, "MPI ERROR: MPI_THREAD_MULTIPLE is not provided on this system\n");
141 MPI_Finalize();
142 exit(77);
143 }
144 #endif
145
146 MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
147 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
148
149 #else
150 nprocs = 1;
151 myid = 0;
152 #endif
153
154 if (argc == 4) {
155 na = atoi(argv[1]);
156 nev = atoi(argv[2]);
157 nblk = atoi(argv[3]);
158 } else {
159 na = 500;
160 nev = 250;
161 nblk = 16;
162 }
163
164 for (np_cols = (int) sqrt((double) nprocs); np_cols > 1; np_cols--) {
165 if (nprocs % np_cols == 0) {
166 break;
167 }
168 }
169
170 np_rows = nprocs/np_cols;
171
172 /* set up blacs */
173 /* convert communicators before */
174 #ifdef WITH_MPI
175 mpi_comm = MPI_Comm_c2f(MPI_COMM_WORLD);
176 #else
177 mpi_comm = 0;
178 #endif
179 set_up_blacsgrid_f(mpi_comm, np_rows, np_cols, 'C', &my_blacs_ctxt, &my_prow, &my_pcol);
180 set_up_blacs_descriptor_f(na, nblk, my_prow, my_pcol, np_rows, np_cols, &na_rows, &na_cols, sc_desc, my_blacs_ctxt, &info);
181
182 /* allocate the matrices needed for elpa */
183 a = calloc(na_rows*na_cols, sizeof(MATRIX_TYPE));
184 z = calloc(na_rows*na_cols, sizeof(MATRIX_TYPE));
185 as = calloc(na_rows*na_cols, sizeof(MATRIX_TYPE));
186 ev = calloc(na, sizeof(EV_TYPE));
187
188 PREPARE_MATRIX_RANDOM(na, myid, na_rows, na_cols, sc_desc, a, z, as);
189
190 #if defined(TEST_GENERALIZED_EIGENPROBLEM)
191 b = calloc(na_rows*na_cols, sizeof(MATRIX_TYPE));
192 bs = calloc(na_rows*na_cols, sizeof(MATRIX_TYPE));
193 PREPARE_MATRIX_RANDOM_SPD(na, myid, na_rows, na_cols, sc_desc, b, z, bs, nblk, np_rows, np_cols, my_prow, my_pcol);
194 #endif
195
196 if (elpa_init(CURRENT_API_VERSION) != ELPA_OK) {
197 fprintf(stderr, "Error: ELPA API version not supported");
198 exit(1);
199 }
200
201 handle = elpa_allocate(&error);
202 //assert_elpa_ok(error);
203
204 /* Set parameters */
205 elpa_set(handle, "na", na, &error);
206 assert_elpa_ok(error);
207
208 elpa_set(handle, "nev", nev, &error);
209 assert_elpa_ok(error);
210
211 if (myid == 0) {
212 printf("Setting the matrix parameters na=%d, nev=%d \n",na,nev);
213 }
214 elpa_set(handle, "local_nrows", na_rows, &error);
215 assert_elpa_ok(error);
216
217 elpa_set(handle, "local_ncols", na_cols, &error);
218 assert_elpa_ok(error);
219
220 elpa_set(handle, "nblk", nblk, &error);
221 assert_elpa_ok(error);
222
223 #ifdef WITH_MPI
224 elpa_set(handle, "mpi_comm_parent", MPI_Comm_c2f(MPI_COMM_WORLD), &error);
225 assert_elpa_ok(error);
226
227 elpa_set(handle, "process_row", my_prow, &error);
228 assert_elpa_ok(error);
229
230 elpa_set(handle, "process_col", my_pcol, &error);
231 assert_elpa_ok(error);
232 #endif
233 #ifdef TEST_GENERALIZED_EIGENPROBLEM
234 elpa_set(handle, "blacs_context", my_blacs_ctxt, &error);
235 assert_elpa_ok(error);
236 #endif
237
238 /* Setup */
239 assert_elpa_ok(elpa_setup(handle));
240
241 /* Set tunables */
242 #ifdef TEST_SOLVER_1STAGE
243 elpa_set(handle, "solver", ELPA_SOLVER_1STAGE, &error);
244 #else
245 elpa_set(handle, "solver", ELPA_SOLVER_2STAGE, &error);
246 #endif
247 assert_elpa_ok(error);
248
249 elpa_set(handle, "gpu", TEST_GPU, &error);
250 assert_elpa_ok(error);
251
252 #if defined(TEST_SOLVE_2STAGE) && defined(TEST_KERNEL)
253 # ifdef TEST_COMPLEX
254 elpa_set(handle, "complex_kernel", TEST_KERNEL, &error);
255 # else
256 elpa_set(handle, "real_kernel", TEST_KERNEL, &error);
257 # endif
258 assert_elpa_ok(error);
259 #endif
260
261 elpa_get(handle, "solver", &value, &error);
262 if (myid == 0) {
263 printf("Solver is set to %d \n", value);
264 }
265
266 #if defined(TEST_GENERALIZED_EIGENPROBLEM)
267 elpa_generalized_eigenvectors(handle, a, b, ev, z, 0, &error);
268 #if defined(TEST_GENERALIZED_DECOMP_EIGENPROBLEM)
269 //a = as, so that the problem can be solved again
270 memcpy(a, as, na_rows * na_cols * sizeof(MATRIX_TYPE));
271 elpa_generalized_eigenvectors(handle, a, b, ev, z, 1, &error);
272 #endif
273 #else
274 /* Solve EV problem */
275 elpa_eigenvectors(handle, a, ev, z, &error);
276 #endif
277 assert_elpa_ok(error);
278
279 elpa_deallocate(handle, &error);
280 elpa_uninit(&error);
281
282 /* check the results */
283 #if defined(TEST_GENERALIZED_EIGENPROBLEM)
284 status = CHECK_CORRECTNESS_EVP_GEN_NUMERIC_RESIDUALS(na, nev, na_rows, na_cols, as, z, ev,
285 sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol, bs);
286 #else
287 status = CHECK_CORRECTNESS_EVP_NUMERIC_RESIDUALS(na, nev, na_rows, na_cols, as, z, ev,
288 sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol);
289 #endif
290
291 if (status !=0){
292 printf("The computed EVs are not correct !\n");
293 }
294 if (status ==0){
295 printf("All ok!\n");
296 }
297
298 free(a);
299 free(z);
300 free(as);
301 free(ev);
302 #if defined(TEST_GENERALIZED_EIGENPROBLEM)
303 free(b);
304 free(bs);
305 #endif
306
307 #ifdef WITH_MPI
308 MPI_Finalize();
309 #endif
310
311 return !!status;
312 }
313