legacy_single_complex_driver_c_version.c 8.07 KB
Newer Older
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
/*     This file is part of ELPA. */
/*  */
/*     The ELPA library was originally created by the ELPA consortium, */
/*     consisting of the following organizations: */
/*  */
/*     - Max Planck Computing and Data Facility (MPCDF), formerly known as */
/*       Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG), */
/*     - Bergische Universität Wuppertal, Lehrstuhl für angewandte */
/*       Informatik, */
/*     - Technische Universität München, Lehrstuhl für Informatik mit */
/*       Schwerpunkt Wissenschaftliches Rechnen , */
/*     - Fritz-Haber-Institut, Berlin, Abt. Theorie, */
/*     - Max-Plack-Institut für Mathematik in den Naturwissenschaften, */
/*       Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition, */
/*       and */
/*     - IBM Deutschland GmbH */
/*  */
/*  */
/*     More information can be found here: */
/*     http://elpa.mpcdf.mpg.de/ */
/*  */
/*     ELPA is free software: you can redistribute it and/or modify */
/*     it under the terms of the version 3 of the license of the */
/*     GNU Lesser General Public License as published by the Free */
/*     Software Foundation. */
/*  */
/*     ELPA is distributed in the hope that it will be useful, */
/*     but WITHOUT ANY WARRANTY; without even the implied warranty of */
/*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the */
/*     GNU Lesser General Public License for more details. */
/*  */
/*     You should have received a copy of the GNU Lesser General Public License */
/*     along with ELPA.  If not, see <http://www.gnu.org/licenses/> */
/*  */
/*     ELPA reflects a substantial effort on the part of the original */
/*     ELPA consortium, and we ask you to respect the spirit of the */
/*     license that we chose: i.e., please contribute any changes you */
/*     may have back to the original ELPA library distribution, and keep */
/*     any derivatives of ELPA under the same license that we chose for */
/*     the original distribution, the GNU Lesser General Public License. */
/*  */
/*  */

#include "config-f90.h"

#include <stdio.h>
#include <stdlib.h>
#ifdef WITH_MPI
#include <mpi.h>
#endif
#include <math.h>

53
#include <elpa/elpa_legacy.h>
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
#include <test/shared/generated.h>
#include <complex.h>

int main(int argc, char** argv) {
   int myid;
   int nprocs;
#ifndef WITH_MPI
   int MPI_COMM_WORLD;
#endif
   int na, nev, nblk;

   int status;

   int np_cols, np_rows, np_colsStart;

69
   int my_blacs_ctxt, my_prow, my_pcol;
70
71
72
73
74
75
76
77
78

   int mpierr;

   int my_mpi_comm_world;
   int mpi_comm_rows, mpi_comm_cols;

   int info, *sc_desc;

   int na_rows, na_cols;
Andreas Marek's avatar
Andreas Marek committed
79
   float startVal;
80

Andreas Marek's avatar
Andreas Marek committed
81
   complex float *a, *z, *as;
82

Andreas Marek's avatar
Andreas Marek committed
83
   float *ev;
84
85
86
87

   int success;
   int i;

88
   int useGPU, THIS_COMPLEX_ELPA_KERNEL_API;
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
#ifdef WITH_MPI
   MPI_Init(&argc, &argv);
   MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
#else
   nprocs = 1;
   myid =0;
   MPI_COMM_WORLD=1;
#endif
   na = 1000;
   nev = 500;
   nblk = 16;

   if (myid == 0) {
     printf("This is the c version of an ELPA test-programm\n");
     printf("\n");
     printf("It will call the ELPA complex solver for a matrix\n");
     printf("of matrix size %d. It will compute %d eigenvalues\n",na,nev);
     printf("and uses a blocksize of %d\n",nblk);
     printf("\n");
     printf("This is an example program with much less functionality\n");
     printf("as it's Fortran counterpart. It's only purpose is to show how \n");
     printf("to evoke ELPA1 from a c programm\n");

     printf("\n");

   }

   status = 0;

Andreas Marek's avatar
Andreas Marek committed
119
   startVal = sqrt((float) nprocs);
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
   np_colsStart = (int) round(startVal);
   for (np_cols=np_colsStart;np_cols>1;np_cols--){
     if (nprocs %np_cols ==0){
     break;
     }
   }

   np_rows = nprocs/np_cols;

   if (myid == 0) {
     printf("\n");
     printf("Number of processor rows %d, cols %d, total %d \n",np_rows,np_cols,nprocs);
   }

   /* set up blacs */
   /* convert communicators before */
#ifdef WITH_MPI
   my_mpi_comm_world = MPI_Comm_c2f(MPI_COMM_WORLD);
#else
   my_mpi_comm_world = 1;
#endif
141
   set_up_blacsgrid_f(my_mpi_comm_world, np_rows, np_cols, 'C', &my_blacs_ctxt, &my_prow, &my_pcol);
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163

   if (myid == 0) {
     printf("\n");
     printf("Past BLACS_Gridinfo...\n");
     printf("\n");
   }

   /* get the ELPA row and col communicators. */
   /* These are NOT usable in C without calling the MPI_Comm_f2c function on them !! */
#ifdef WITH_MPI
   my_mpi_comm_world = MPI_Comm_c2f(MPI_COMM_WORLD);
#endif
   mpierr = elpa_get_communicators(my_mpi_comm_world, my_prow, my_pcol, &mpi_comm_rows, &mpi_comm_cols);

   if (myid == 0) {
     printf("\n");
     printf("Past split communicator setup for rows and columns...\n");
     printf("\n");
   }

   sc_desc = malloc(9*sizeof(int));

164
   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
167
168
169
170
171
172
173
174
175
176
177
178

   if (myid == 0) {
     printf("\n");
     printf("Past scalapack descriptor setup...\n");
     printf("\n");
   }

   /* allocate the matrices needed for elpa */
   if (myid == 0) {
     printf("\n");
     printf("Allocating matrices with na_rows=%d and na_cols=%d\n",na_rows, na_cols);
     printf("\n");
   }

Andreas Marek's avatar
Andreas Marek committed
179
180
181
182
   a  = malloc(na_rows*na_cols*sizeof(complex float));
   z  = malloc(na_rows*na_cols*sizeof(complex float));
   as = malloc(na_rows*na_cols*sizeof(complex float));
   ev = malloc(na*sizeof(float));
183

Andreas Marek's avatar
Andreas Marek committed
184
   prepare_matrix_complex_single_f(na, myid, na_rows, na_cols, sc_desc, a, z, as);
185
186
187
188
189
190
191
192
193

   if (myid == 0) {
     printf("\n");
     printf("Entering ELPA 1stage complex solver\n");
     printf("\n");
   }
#ifdef WITH_MPI
   mpierr = MPI_Barrier(MPI_COMM_WORLD);
#endif
194
   useGPU = 0;
195
   THIS_COMPLEX_ELPA_KERNEL_API = ELPA_2STAGE_COMPLEX_DEFAULT;
Andreas Marek's avatar
Andreas Marek committed
196
   success = elpa_solve_evp_complex_single(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
199
200
201

   if (success != 1) {
     printf("error in ELPA solve \n");
#ifdef WITH_MPI
     mpierr = MPI_Abort(MPI_COMM_WORLD, 99);
202
203
#else
     exit(99);
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
#endif
   }


   if (myid == 0) {
     printf("\n");
     printf("1stage ELPA complex solver complete\n");
     printf("\n");
   }

   for (i=0;i<na_rows*na_cols;i++){
      a[i] = as[i];
      z[i] = as[i];
   }
   if (myid == 0) {
     printf("\n");
     printf("Entering ELPA 2stage complex solver\n");
     printf("\n");
   }
#ifdef WITH_MPI
   mpierr = MPI_Barrier(MPI_COMM_WORLD);
#endif
226
   useGPU =0;
227
   THIS_COMPLEX_ELPA_KERNEL_API = ELPA_2STAGE_COMPLEX_DEFAULT;
Andreas Marek's avatar
Andreas Marek committed
228
   success = elpa_solve_evp_complex_single(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
231
232
233

   if (success != 1) {
     printf("error in ELPA solve \n");
#ifdef WITH_MPI
     mpierr = MPI_Abort(MPI_COMM_WORLD, 99);
234
235
#else
     exit(99);
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
#endif
   }

   if (myid == 0) {
     printf("\n");
     printf("2stage ELPA complex solver complete\n");
     printf("\n");
   }

   for (i=0;i<na_rows*na_cols;i++){
      a[i] = as[i];
      z[i] = as[i];
   }
   if (myid == 0) {
     printf("\n");
     printf("Entering auto-chosen ELPA complex solver\n");
     printf("\n");
   }
#ifdef WITH_MPI
   mpierr = MPI_Barrier(MPI_COMM_WORLD);
#endif
257
   useGPU = 0;
258
   THIS_COMPLEX_ELPA_KERNEL_API = ELPA_2STAGE_COMPLEX_DEFAULT;
Andreas Marek's avatar
Andreas Marek committed
259
   success = elpa_solve_evp_complex_single(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
262
263
264

   if (success != 1) {
     printf("error in ELPA solve \n");
#ifdef WITH_MPI
     mpierr = MPI_Abort(MPI_COMM_WORLD, 99);
265
266
#else
     exit(99);
267
268
269
270
271
272
273
274
275
276
#endif
   }

   if (myid == 0) {
     printf("\n");
     printf("Auto-chosen ELPA complex solver complete\n");
     printf("\n");
   }

   /* check the results */
Andreas Marek's avatar
Andreas Marek committed
277
   status = check_correctness_complex_single_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
278
279
280
281
282
283
284
285
286
287
288
289
290
291

   if (status !=0){
     printf("The computed EVs are not correct !\n");
   }
   if (status ==0){
     if (myid == 0) {
       printf("All ok!\n");
     }
   }

   free(sc_desc);
   free(a);
   free(z);
   free(as);
292
   free(ev);
293

294
295
296
297
298
#ifdef WITH_MPI
   MPI_Finalize();
#endif
   return 0;
}