elpa1_test_complex_c_version.c 7.85 KB
Newer Older
Andreas Marek's avatar
Andreas Marek committed
1
2
3
4
5
/*     This file is part of ELPA. */
/*  */
/*     The ELPA library was originally created by the ELPA consortium, */
/*     consisting of the following organizations: */
/*  */
6
7
/*     - Max Planck Computing and Data Facility (MPCDF), formerly known as */
/*       Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG), */
Andreas Marek's avatar
Andreas Marek committed
8
9
10
11
12
13
14
15
16
17
18
19
/*     - 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: */
20
/*     http://elpa.mpcdf.mpg.de/ */
Andreas Marek's avatar
Andreas Marek committed
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
/*  */
/*     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>
48
#ifdef WITH_MPI
Andreas Marek's avatar
Andreas Marek committed
49
#include <mpi.h>
50
#endif
Andreas Marek's avatar
Andreas Marek committed
51
52
53
#include <math.h>

#include <elpa/elpa.h>
54
#include <test/shared/generated.h>
Andreas Marek's avatar
Andreas Marek committed
55
#include <complex.h>
Andreas Marek's avatar
Andreas Marek committed
56

57
58
#define DOUBLE_PRECISION_COMPLEX 1

59
int main(int argc, char** argv) {
60
61
   int myid;
   int nprocs;
62
63
64
#ifndef WITH_MPI
   int MPI_COMM_WORLD;
#endif
65
   int na, nev, nblk;
Andreas Marek's avatar
Andreas Marek committed
66

67
   int status;
Andreas Marek's avatar
Andreas Marek committed
68

69
   int np_cols, np_rows, np_colsStart;
Andreas Marek's avatar
Andreas Marek committed
70

71
   int my_blacs_ctxt, nprow, npcol, my_prow, my_pcol;
Andreas Marek's avatar
Andreas Marek committed
72

73
   int mpierr;
Andreas Marek's avatar
Andreas Marek committed
74

75
76
   int my_mpi_comm_world;
   int mpi_comm_rows, mpi_comm_cols;
Andreas Marek's avatar
Andreas Marek committed
77

78
   int info, *sc_desc;
Andreas Marek's avatar
Andreas Marek committed
79

80
   int na_rows, na_cols;
81

82
   double startVal;
83
#ifdef DOUBLE_PRECISION_COMPLEX
Andreas Marek's avatar
Andreas Marek committed
84

Andreas Marek's avatar
Andreas Marek committed
85
86
87
   complex double *a, *z, *as, *tmp1, *tmp2;

   double *ev, *xr;
88
89
90
91
92
93
#else

   complex *a, *z, *as, *tmp1, *tmp2;

   float *ev, *xr;
#endif
Andreas Marek's avatar
Andreas Marek committed
94

95
   int *iseed;
Andreas Marek's avatar
Andreas Marek committed
96

97
   int success;
98
#ifdef WITH_MPI
99
100
101
   MPI_Init(&argc, &argv);
   MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
102
103
104
105
106
#else
   nprocs=1;
   myid=0;
   MPI_COMM_WORLD=1;
#endif
107
108
109
   na = 1000;
   nev = 500;
   nblk = 16;
Andreas Marek's avatar
Andreas Marek committed
110

111
112
113
   if (myid == 0) {
     printf("This is the c version of an ELPA test-programm\n");
     printf("\n");
Andreas Marek's avatar
Andreas Marek committed
114
     printf("It will call the 1stage ELPA complex solver for a matrix\n");
115
116
117
118
     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");
Andreas Marek's avatar
Andreas Marek committed
119
120
121
     printf("as it's Fortran counterpart. It's only purpose is to show how \n");
     printf("to evoke ELPA1 from a c programm\n");

122
     printf("\n");
Andreas Marek's avatar
Andreas Marek committed
123

124
125
126
127
128
129
#ifdef DOUBLE_PRECISION_COMPLEX
     printf("The double precision version of ELPA1 is used\n");
#else
     printf("The single precision version of ELPA1 is used\n");
#endif
     printf("\n");
130
   }
Andreas Marek's avatar
Andreas Marek committed
131

132
   status = 0;
Andreas Marek's avatar
Andreas Marek committed
133

134
135
136
137
138
139
140
   startVal = sqrt((double) nprocs);
   np_colsStart = (int) round(startVal);
   for (np_cols=np_colsStart;np_cols>1;np_cols--){
     if (nprocs %np_cols ==0){
     break;
     }
   }
Andreas Marek's avatar
Andreas Marek committed
141

142
   np_rows = nprocs/np_cols;
Andreas Marek's avatar
Andreas Marek committed
143

144
145
146
147
   if (myid == 0) {
     printf("\n");
     printf("Number of processor rows %d, cols %d, total %d \n",np_rows,np_cols,nprocs);
   }
Andreas Marek's avatar
Andreas Marek committed
148

149
150
   /* set up blacs */
   /* convert communicators before */
151
#ifdef WITH_MPI
Andreas Marek's avatar
Andreas Marek committed
152
   my_mpi_comm_world = MPI_Comm_c2f(MPI_COMM_WORLD);
153
154
155
#else
   my_mpi_comm_world = 1;
#endif
Andreas Marek's avatar
Andreas Marek committed
156
157
158
159
160
161
162
163
   set_up_blacsgrid_from_fortran(my_mpi_comm_world, &my_blacs_ctxt, &np_rows, &np_cols, &nprow, &npcol, &my_prow, &my_pcol);

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

164
165
   /* get the ELPA row and col communicators. */
   /* These are NOT usable in C without calling the MPI_Comm_f2c function on them !! */
166
#ifdef WITH_MPI
Andreas Marek's avatar
Andreas Marek committed
167
   my_mpi_comm_world = MPI_Comm_c2f(MPI_COMM_WORLD);
168
#endif
169
   mpierr = elpa_get_communicators(my_mpi_comm_world, my_prow, my_pcol, &mpi_comm_rows, &mpi_comm_cols);
Andreas Marek's avatar
Andreas Marek committed
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186

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

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

   set_up_blacs_descriptor_from_fortran(na, nblk, my_prow, my_pcol, np_rows, np_cols, &na_rows, &na_cols, sc_desc, my_blacs_ctxt, &info);

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

187
   /* allocate the matrices needed for elpa */
Andreas Marek's avatar
Andreas Marek committed
188
189
190
191
192
193
   if (myid == 0) {
     printf("\n");
     printf("Allocating matrices with na_rows=%d and na_cols=%d\n",na_rows, na_cols);
     printf("\n");
   }

194
#ifdef DOUBLE_PRECISION_COMPLEX
Andreas Marek's avatar
Andreas Marek committed
195
196
197
198
199
   a  = malloc(na_rows*na_cols*sizeof(complex double));
   z  = malloc(na_rows*na_cols*sizeof(complex double));
   as = malloc(na_rows*na_cols*sizeof(complex double));

   xr = malloc(na_rows*na_cols*sizeof(double));
Andreas Marek's avatar
Andreas Marek committed
200
201
202
203


   ev = malloc(na*sizeof(double));

Andreas Marek's avatar
Andreas Marek committed
204
205
   tmp1  = malloc(na_rows*na_cols*sizeof(complex double));
   tmp2 = malloc(na_rows*na_cols*sizeof(complex double));
206
207
208
209
#else
   a  = malloc(na_rows*na_cols*sizeof(complex));
   z  = malloc(na_rows*na_cols*sizeof(complex));
   as = malloc(na_rows*na_cols*sizeof(complex));
Andreas Marek's avatar
Andreas Marek committed
210

211
   xr = malloc(na_rows*na_cols*sizeof(float));
Andreas Marek's avatar
Andreas Marek committed
212

213
214
215
216
217
218
219
220
221
222
223
224
225

   ev = malloc(na*sizeof(float));

   tmp1  = malloc(na_rows*na_cols*sizeof(complex));
   tmp2 = malloc(na_rows*na_cols*sizeof(complex));
#endif

   iseed = malloc(4096*sizeof(int));
#ifdef DOUBLE_PRECISION_COMPLEX
   prepare_matrix_complex_from_fortran_double_precision(na, myid, na_rows, na_cols, sc_desc, iseed, xr, a, z, as);
#else
   prepare_matrix_complex_from_fortran_single_precision(na, myid, na_rows, na_cols, sc_desc, iseed, xr, a, z, as);
#endif
Andreas Marek's avatar
Andreas Marek committed
226
227

   free(xr);
Andreas Marek's avatar
Andreas Marek committed
228
229
230

   if (myid == 0) {
     printf("\n");
Andreas Marek's avatar
Andreas Marek committed
231
     printf("Entering ELPA 1stage complex solver\n");
Andreas Marek's avatar
Andreas Marek committed
232
233
     printf("\n");
   }
234
#ifdef WITH_MPI
Andreas Marek's avatar
Andreas Marek committed
235
   mpierr = MPI_Barrier(MPI_COMM_WORLD);
236
#endif
Andreas Marek's avatar
Andreas Marek committed
237

238
239
240
241
242
#ifdef DOUBLE_PRECISION_COMPLEX
   success = elpa_solve_evp_complex_1stage_double_precision(na, nev, a, na_rows, ev, z, na_rows, nblk, na_cols, mpi_comm_rows, mpi_comm_cols);
#else
   success = elpa_solve_evp_complex_1stage_single_precision(na, nev, a, na_rows, ev, z, na_rows, nblk, na_cols, mpi_comm_rows, mpi_comm_cols);
#endif
Andreas Marek's avatar
Andreas Marek committed
243
244
245

   if (success != 1) {
     printf("error in ELPA solve \n");
246
#ifdef WITH_MPI
Andreas Marek's avatar
Andreas Marek committed
247
     mpierr = MPI_Abort(MPI_COMM_WORLD, 99);
248
#endif
Andreas Marek's avatar
Andreas Marek committed
249
250
251
252
253
   }


   if (myid == 0) {
     printf("\n");
Andreas Marek's avatar
Andreas Marek committed
254
     printf("1stage ELPA complex solver complete\n");
Andreas Marek's avatar
Andreas Marek committed
255
256
257
     printf("\n");
   }

258
   /* check the results */
259
260
261
262
263
#ifdef DOUBLE_PRECISION_COMPLEX
   status = check_correctness_complex_from_fortran_double_precision(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid, tmp1, tmp2);
#else
   status = check_correctness_complex_from_fortran_single_precision(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid, tmp1, tmp2);
#endif
Andreas Marek's avatar
Andreas Marek committed
264
265
266
267
268
269
270
271

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

272
273
274
275
276
277
278
   free(sc_desc);
   free(a);
   free(z);
   free(as);

   free(tmp1);
   free(tmp2);
279
#ifdef WITH_MPI
280
   MPI_Finalize();
281
#endif
282
   return 0;
Andreas Marek's avatar
Andreas Marek committed
283
}