Line data Source code
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 :
44 : #include "config-f90.h"
45 :
46 : #include <stdio.h>
47 : #include <stdlib.h>
48 : #ifdef WITH_MPI
49 : #include <mpi.h>
50 : #endif
51 : #include <math.h>
52 :
53 : #include <elpa/elpa_legacy.h>
54 : #include <test/shared/generated.h>
55 :
56 96 : int main(int argc, char** argv) {
57 : int myid;
58 : int nprocs;
59 : #ifndef WITH_MPI
60 : int MPI_COMM_WORLD;
61 : #endif
62 : int na, nev, nblk;
63 :
64 : int status;
65 :
66 : int np_cols, np_rows, np_colsStart;
67 :
68 : int my_blacs_ctxt, my_prow, my_pcol;
69 :
70 : int mpierr;
71 :
72 : int my_mpi_comm_world;
73 : int mpi_comm_rows, mpi_comm_cols;
74 :
75 : int info, *sc_desc;
76 :
77 : int na_rows, na_cols;
78 : float startVal;
79 :
80 : float *a, *z, *as, *ev;
81 :
82 : int i;
83 : int useGPU;
84 : int success;
85 :
86 : int useQr, THIS_REAL_ELPA_KERNEL_API;
87 : #ifdef WITH_MPI
88 64 : MPI_Init(&argc, &argv);
89 64 : MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
90 64 : MPI_Comm_rank(MPI_COMM_WORLD, &myid);
91 : #else
92 32 : nprocs = 1;
93 32 : myid=0;
94 32 : MPI_COMM_WORLD=1;
95 : #endif
96 96 : na = 1000;
97 96 : nev = 500;
98 96 : nblk = 16;
99 :
100 96 : if (myid == 0) {
101 64 : printf("This is the c version of an ELPA test-programm\n");
102 64 : printf("\n");
103 64 : printf("It will call the ELPA real solver for an\n");
104 64 : printf("of matrix size %d. It will compute %d eigenvalues\n",na,nev);
105 64 : printf("and uses a blocksize of %d\n",nblk);
106 64 : printf("\n");
107 64 : printf("This is an example program with much less functionality\n");
108 64 : printf("as it's Fortran counterpart. It's only purpose is to show how \n");
109 64 : printf("to evoke ELPA1 from a c programm\n");
110 64 : printf("\n");
111 :
112 : }
113 :
114 96 : status = 0;
115 :
116 96 : startVal = sqrt((float) nprocs);
117 96 : np_colsStart = (int) round(startVal);
118 96 : for (np_cols=np_colsStart;np_cols>1;np_cols--){
119 0 : if (nprocs %np_cols ==0){
120 0 : break;
121 : }
122 : }
123 :
124 96 : np_rows = nprocs/np_cols;
125 :
126 96 : if (myid == 0) {
127 64 : printf("\n");
128 64 : printf("Number of processor rows %d, cols %d, total %d \n",np_rows,np_cols,nprocs);
129 : }
130 :
131 : /* set up blacs */
132 : /* convert communicators before */
133 : #ifdef WITH_MPI
134 64 : my_mpi_comm_world = MPI_Comm_c2f(MPI_COMM_WORLD);
135 : #else
136 32 : my_mpi_comm_world = 1;
137 : #endif
138 96 : set_up_blacsgrid_f(my_mpi_comm_world, np_rows, np_cols, 'C', &my_blacs_ctxt, &my_prow, &my_pcol);
139 :
140 96 : if (myid == 0) {
141 64 : printf("\n");
142 64 : printf("Past BLACS_Gridinfo...\n");
143 64 : printf("\n");
144 : }
145 :
146 : /* get the ELPA row and col communicators. */
147 : /* These are NOT usable in C without calling the MPI_Comm_f2c function on them !! */
148 : #ifdef WITH_MPI
149 64 : my_mpi_comm_world = MPI_Comm_c2f(MPI_COMM_WORLD);
150 : #endif
151 96 : mpierr = elpa_get_communicators(my_mpi_comm_world, my_prow, my_pcol, &mpi_comm_rows, &mpi_comm_cols);
152 :
153 96 : if (myid == 0) {
154 64 : printf("\n");
155 64 : printf("Past split communicator setup for rows and columns...\n");
156 64 : printf("\n");
157 : }
158 :
159 96 : sc_desc = malloc(9*sizeof(int));
160 :
161 96 : 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);
162 :
163 96 : if (myid == 0) {
164 64 : printf("\n");
165 64 : printf("Past scalapack descriptor setup...\n");
166 64 : printf("\n");
167 : }
168 :
169 : /* allocate the matrices needed for elpa */
170 96 : if (myid == 0) {
171 64 : printf("\n");
172 64 : printf("Allocating matrices with na_rows=%d and na_cols=%d\n",na_rows, na_cols);
173 64 : printf("\n");
174 : }
175 :
176 96 : a = malloc(na_rows*na_cols*sizeof(float));
177 96 : z = malloc(na_rows*na_cols*sizeof(float));
178 96 : as = malloc(na_rows*na_cols*sizeof(float));
179 96 : ev = malloc(na*sizeof(float));
180 :
181 96 : prepare_matrix_random_real_single_f(na, myid, na_rows, na_cols, sc_desc, a, z, as);
182 :
183 96 : if (myid == 0) {
184 64 : printf("\n");
185 64 : printf("Entering ELPA 1stage real solver\n");
186 64 : printf("\n");
187 : }
188 : #ifdef WITH_MPI
189 64 : mpierr = MPI_Barrier(MPI_COMM_WORLD);
190 : #endif
191 96 : useQr = 0;
192 96 : useGPU = 0;
193 96 : THIS_REAL_ELPA_KERNEL_API = ELPA_2STAGE_REAL_DEFAULT;
194 :
195 96 : success = elpa_solve_evp_real_single(na, nev, a, na_rows, ev, z, na_rows, nblk, na_cols, mpi_comm_rows, mpi_comm_cols, my_mpi_comm_world, THIS_REAL_ELPA_KERNEL_API, useQr, useGPU, "1stage");
196 :
197 96 : if (success != 1) {
198 0 : printf("error in ELPA solve \n");
199 : #ifdef WITH_MPI
200 0 : mpierr = MPI_Abort(MPI_COMM_WORLD, 99);
201 : #else
202 0 : exit(99);
203 : #endif
204 : }
205 :
206 :
207 96 : if (myid == 0) {
208 64 : printf("\n");
209 64 : printf("1stage ELPA real solver complete\n");
210 64 : printf("\n");
211 : }
212 :
213 64000096 : for (i=0;i<na_rows*na_cols;i++){
214 64000000 : a[i] = as[i];
215 64000000 : z[i] = as[i];
216 : }
217 :
218 96 : if (myid == 0) {
219 64 : printf("\n");
220 64 : printf("Entering ELPA 2stage real solver\n");
221 64 : printf("\n");
222 : }
223 : #ifdef WITH_MPI
224 64 : mpierr = MPI_Barrier(MPI_COMM_WORLD);
225 : #endif
226 96 : useGPU = 0;
227 96 : useQr = 0;
228 96 : THIS_REAL_ELPA_KERNEL_API = ELPA_2STAGE_REAL_DEFAULT;
229 :
230 96 : success = elpa_solve_evp_real_single(na, nev, a, na_rows, ev, z, na_rows, nblk, na_cols, mpi_comm_rows, mpi_comm_cols, my_mpi_comm_world, THIS_REAL_ELPA_KERNEL_API, useQr, useGPU, "2stage");
231 :
232 96 : if (success != 1) {
233 0 : printf("error in ELPA solve \n");
234 : #ifdef WITH_MPI
235 0 : mpierr = MPI_Abort(MPI_COMM_WORLD, 99);
236 : #else
237 0 : exit(99);
238 : #endif
239 : }
240 96 : if (myid == 0) {
241 64 : printf("\n");
242 64 : printf("2stage ELPA real solver complete\n");
243 64 : printf("\n");
244 : }
245 :
246 64000096 : for (i=0;i<na_rows*na_cols;i++){
247 64000000 : a[i] = as[i];
248 64000000 : z[i] = as[i];
249 : }
250 :
251 96 : if (myid == 0) {
252 64 : printf("\n");
253 64 : printf("Entering auto-chosen ELPA real solver\n");
254 64 : printf("\n");
255 : }
256 : #ifdef WITH_MPI
257 64 : mpierr = MPI_Barrier(MPI_COMM_WORLD);
258 : #endif
259 96 : useQr = 0;
260 96 : useGPU = 0;
261 96 : THIS_REAL_ELPA_KERNEL_API = ELPA_2STAGE_REAL_DEFAULT;
262 :
263 96 : success = elpa_solve_evp_real_single(na, nev, a, na_rows, ev, z, na_rows, nblk, na_cols, mpi_comm_rows, mpi_comm_cols, my_mpi_comm_world, THIS_REAL_ELPA_KERNEL_API, useQr, useGPU, "auto");
264 :
265 96 : if (success != 1) {
266 0 : printf("error in ELPA solve \n");
267 : #ifdef WITH_MPI
268 0 : mpierr = MPI_Abort(MPI_COMM_WORLD, 99);
269 : #else
270 0 : exit(99);
271 : #endif
272 : }
273 96 : if (myid == 0) {
274 64 : printf("\n");
275 64 : printf("Auto-chosen ELPA real solver complete\n");
276 64 : printf("\n");
277 : }
278 :
279 : /* check the results */
280 96 : status = check_correctness_evp_numeric_residuals_real_single_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
281 :
282 96 : if (status !=0){
283 0 : printf("The computed EVs are not correct !\n");
284 : }
285 96 : if (status ==0){
286 96 : if (myid ==0) {
287 64 : printf("All ok!\n");
288 : }
289 : }
290 :
291 96 : free(sc_desc);
292 96 : free(a);
293 96 : free(z);
294 96 : free(as);
295 96 : free(ev);
296 :
297 : #ifdef WITH_MPI
298 64 : MPI_Finalize();
299 : #endif
300 96 : return 0;
301 : }
|
