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 <mpi.h>
48 #include <math.h>
49 #include <string.h>
50
51 #include <elpa/elpa.h>
52 #include <assert.h>
53
54 #define assert_elpa_ok(x) assert(x == ELPA_OK)
55
main(int argc,char ** argv)56 int main(int argc, char** argv) {
57 /* matrix dimensions */
58 const int na = 1000;
59 const int nev = 500;
60 const int nblk = 16;
61
62 /* mpi */
63 int myid, nprocs;
64 int na_cols, na_rows;
65 int np_cols, np_rows;
66 int my_prow, my_pcol;
67 int mpi_comm;
68
69 /* blacs */
70 int my_blacs_ctxt, sc_desc[9], info;
71
72 /* The Matrix */
73 double *a, *as, *z;
74 double *ev;
75
76 int error, status;
77
78 elpa_t handle;
79
80 int value;
81 MPI_Init(&argc, &argv);
82 MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
83 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
84
85
86 for (np_cols = (int) sqrt((double) nprocs); np_cols > 1; np_cols--) {
87 if (nprocs % np_cols == 0) {
88 break;
89 }
90 }
91
92 np_rows = nprocs/np_cols;
93
94 /* set up blacs */
95 /* convert communicators before */
96 mpi_comm = MPI_Comm_c2f(MPI_COMM_WORLD);
97 set_up_blacsgrid_f(mpi_comm, np_rows, np_cols, 'C', &my_blacs_ctxt, &my_prow, &my_pcol);
98 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);
99
100 /* allocate the matrices needed for elpa */
101 a = calloc(na_rows*na_cols, sizeof(double));
102 z = calloc(na_rows*na_cols, sizeof(double));
103 as = calloc(na_rows*na_cols, sizeof(double));
104 ev = calloc(na, sizeof(double));
105
106 // TODO: prepare properly
107 memset(a, 0, na_rows * na_cols * sizeof(double));
108 //prepare_matrix_real_double_f(na, myid, na_rows, na_cols, sc_desc, a, z, as);
109
110 if (elpa_init(20170403) != ELPA_OK) {
111 fprintf(stderr, "Error: ELPA API version not supported");
112 exit(1);
113 }
114
115 if(myid == 0) printf("init done\n");
116
117 handle = elpa_allocate(&error);
118 assert_elpa_ok(error);
119
120 /* Set parameters */
121 elpa_set(handle, "na", na, &error);
122 assert_elpa_ok(error);
123
124 elpa_set(handle, "nev", nev, &error);
125 assert_elpa_ok(error);
126
127 elpa_set(handle, "local_nrows", na_rows, &error);
128 assert_elpa_ok(error);
129
130 elpa_set(handle, "local_ncols", na_cols, &error);
131 assert_elpa_ok(error);
132
133 elpa_set(handle, "nblk", nblk, &error);
134 assert_elpa_ok(error);
135
136 elpa_set(handle, "mpi_comm_parent", MPI_Comm_c2f(MPI_COMM_WORLD), &error);
137 assert_elpa_ok(error);
138
139 elpa_set(handle, "process_row", my_prow, &error);
140 assert_elpa_ok(error);
141
142 elpa_set(handle, "process_col", my_pcol, &error);
143 assert_elpa_ok(error);
144 //
145 /* Setup */
146 assert_elpa_ok(elpa_setup(handle));
147
148 if(myid == 0) printf("setup done\n");
149 /* Set tunables */
150 elpa_set(handle, "solver", ELPA_SOLVER_1STAGE, &error);
151 // elpa_set(handle, "solver", ELPA_SOLVER_2STAGE, &error);
152 assert_elpa_ok(error);
153
154 // elpa_set(handle, "real_kernel", TEST_KERNEL, &error);
155 // assert_elpa_ok(error);
156
157 if(myid == 0) printf("solve..\n");
158
159 /* Solve EV problem */
160 elpa_eigenvectors(handle, a, ev, z, &error);
161 assert_elpa_ok(error);
162 if(myid == 0) printf("solve done \n");
163
164 // for(int i = 0; i < na; i++)
165 // printf("%lf, ", ev[i]);
166 // printf("\n");
167
168 elpa_deallocate(handle);
169 elpa_uninit();
170
171
172 /* check the results */
173 // status = check_correctness_real_double_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
174
175 // if (status !=0){
176 // printf("The computed EVs are not correct !\n");
177 // }
178 // if (status ==0){
179 // printf("All ok!\n");
180 // }
181
182 free(a);
183 free(z);
184 free(as);
185 free(ev);
186
187 MPI_Finalize();
188
189 return !!status;
190 }
191