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