elpa2_test_complex_c_version.c 7.74 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
48
49
50
51
/*  */
/*     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>
#include <mpi.h>
#include <math.h>

#include <elpa/elpa.h>
Andreas Marek's avatar
Andreas Marek committed
52
#include <test/shared_sources/generated.h>
Andreas Marek's avatar
Andreas Marek committed
53
#include <complex.h>
Andreas Marek's avatar
Andreas Marek committed
54

55
int main(int argc, char** argv) {
56
57
   int myid;
   int nprocs;
Andreas Marek's avatar
Andreas Marek committed
58

59
   int na, nev, nblk;
Andreas Marek's avatar
Andreas Marek committed
60

61
   int status;
Andreas Marek's avatar
Andreas Marek committed
62

63
   int np_cols, np_rows, np_colsStart;
Andreas Marek's avatar
Andreas Marek committed
64

65
   int my_blacs_ctxt, nprow, npcol, my_prow, my_pcol;
Andreas Marek's avatar
Andreas Marek committed
66

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

69
70
   int my_mpi_comm_world;
   int mpi_comm_rows, mpi_comm_cols;
Andreas Marek's avatar
Andreas Marek committed
71

72
   int info, *sc_desc;
Andreas Marek's avatar
Andreas Marek committed
73

74
75
   int na_rows, na_cols;
   double startVal;
76
#ifdef DOUBLE_PRECISION_COMPLEX
Andreas Marek's avatar
Andreas Marek committed
77
78
79
   complex double *a, *z, *as, *tmp1, *tmp2;

   double *ev, *xr;
80
81
#else
   complex *a, *z, *as, *tmp1, *tmp2;
Andreas Marek's avatar
Andreas Marek committed
82

83
84
   float *ev, *xr;
#endif
85
   int *iseed;
Andreas Marek's avatar
Andreas Marek committed
86

87
   int success;
Andreas Marek's avatar
Andreas Marek committed
88

Andreas Marek's avatar
Andreas Marek committed
89
90
   int THIS_COMPLEX_ELPA_KERNEL_API;

91
92
93
   MPI_Init(&argc, &argv);
   MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
Andreas Marek's avatar
Andreas Marek committed
94

95
96
97
   na = 1000;
   nev = 500;
   nblk = 16;
Andreas Marek's avatar
Andreas Marek committed
98

99
100
101
   if (myid == 0) {
     printf("This is the c version of an ELPA test-programm\n");
     printf("\n");
Andreas Marek's avatar
Andreas Marek committed
102
     printf("It will call the 1stage ELPA complex solver for a matrix\n");
103
104
105
106
     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
107
108
109
     printf("as it's Fortran counterpart. It's only purpose is to show how \n");
     printf("to evoke ELPA1 from a c programm\n");

110
     printf("\n");
Andreas Marek's avatar
Andreas Marek committed
111

112
113
114
115
116
#ifdef DOUBLE_PRECISION_COMPLEX
    printf(" Double precision version of ELPA2 is used. \n");
#else
    printf(" Single precision version of ELPA2 is used. \n");
#endif
117
   }
Andreas Marek's avatar
Andreas Marek committed
118

119
   status = 0;
Andreas Marek's avatar
Andreas Marek committed
120

121
122
123
124
125
126
127
   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
128

129
   np_rows = nprocs/np_cols;
Andreas Marek's avatar
Andreas Marek committed
130

131
132
133
134
   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
135

136
137
   /* set up blacs */
   /* convert communicators before */
Andreas Marek's avatar
Andreas Marek committed
138
139
140
141
142
143
144
145
146
   my_mpi_comm_world = MPI_Comm_c2f(MPI_COMM_WORLD);
   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");
   }

147
148
   /* get the ELPA row and col communicators. */
   /* These are NOT usable in C without calling the MPI_Comm_f2c function on them !! */
Andreas Marek's avatar
Andreas Marek committed
149
   my_mpi_comm_world = MPI_Comm_c2f(MPI_COMM_WORLD);
150
   mpierr = get_elpa_communicators(my_mpi_comm_world, my_prow, my_pcol, &mpi_comm_rows, &mpi_comm_cols);
Andreas Marek's avatar
Andreas Marek committed
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167

   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");
   }

168
   /* allocate the matrices needed for elpa */
Andreas Marek's avatar
Andreas Marek committed
169
170
171
172
173
   if (myid == 0) {
     printf("\n");
     printf("Allocating matrices with na_rows=%d and na_cols=%d\n",na_rows, na_cols);
     printf("\n");
   }
174
#ifdef DOUBLE_PRECISION_COMPLEX
Andreas Marek's avatar
Andreas Marek committed
175
176
177
178
179
   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
180
181
182
183


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

Andreas Marek's avatar
Andreas Marek committed
184
185
   tmp1  = malloc(na_rows*na_cols*sizeof(complex double));
   tmp2 = malloc(na_rows*na_cols*sizeof(complex double));
186
187
188
189
190
191
#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));

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


194
195
196
197
198
199
200
201
202
203
204
   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
205
206

   free(xr);
Andreas Marek's avatar
Andreas Marek committed
207
208
209

   if (myid == 0) {
     printf("\n");
Andreas Marek's avatar
Andreas Marek committed
210
     printf("Entering ELPA 2stage complex solver\n");
Andreas Marek's avatar
Andreas Marek committed
211
212
213
214
     printf("\n");
   }

   mpierr = MPI_Barrier(MPI_COMM_WORLD);
Andreas Marek's avatar
Andreas Marek committed
215
   THIS_COMPLEX_ELPA_KERNEL_API = ELPA2_COMPLEX_KERNEL_GENERIC;
216
217
218
219
220
#ifdef DOUBLE_PRECISION_COMPLEX
   success = elpa_solve_evp_complex_2stage_double_precision(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);
#else
   success = elpa_solve_evp_complex_2stage_single_precision(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);
#endif
Andreas Marek's avatar
Andreas Marek committed
221
222
223
224
225
226
227
228
229

   if (success != 1) {
     printf("error in ELPA solve \n");
     mpierr = MPI_Abort(MPI_COMM_WORLD, 99);
   }


   if (myid == 0) {
     printf("\n");
Andreas Marek's avatar
Andreas Marek committed
230
     printf("2stage ELPA complex solver complete\n");
Andreas Marek's avatar
Andreas Marek committed
231
232
233
     printf("\n");
   }

234
   /* check the results */
235
236
237
238
239
#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
240
241
242
243
   if (status !=0){
     printf("The computed EVs are not correct !\n");
   }
   if (status ==0){
Andreas Marek's avatar
Andreas Marek committed
244
245
246
     if (myid == 0) {
       printf("All ok!\n");
     }
Andreas Marek's avatar
Andreas Marek committed
247
248
   }

249
250
251
252
253
254
255
   free(sc_desc);
   free(a);
   free(z);
   free(as);

   free(tmp1);
   free(tmp2);
Andreas Marek's avatar
Andreas Marek committed
256

257
   MPI_Finalize();
Andreas Marek's avatar
Andreas Marek committed
258

259
   return 0;
Andreas Marek's avatar
Andreas Marek committed
260
}