legacy_complex_2stage_c_version.c 7.56 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
#include <math.h>

53
#include <elpa/elpa_legacy.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, 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 81
   int na_rows, na_cols;
   double startVal;
82
#ifdef DOUBLE_PRECISION_COMPLEX
83 84
   complex double *a, *z, *as;
   double *ev;
85
#else
86 87
   complex *a, *z, *as;
   float *ev;
88
#endif
Andreas Marek's avatar
Andreas Marek committed
89

90
   int success;
Andreas Marek's avatar
Andreas Marek committed
91

92
   int THIS_COMPLEX_ELPA_KERNEL_API, useGPU;
93

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

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

118
     printf("\n");
Andreas Marek's avatar
Andreas Marek committed
119

120 121 122 123 124
#ifdef DOUBLE_PRECISION_COMPLEX
    printf(" Double precision version of ELPA2 is used. \n");
#else
    printf(" Single precision version of ELPA2 is used. \n");
#endif
125
   }
Andreas Marek's avatar
Andreas Marek committed
126

127
   status = 0;
Andreas Marek's avatar
Andreas Marek committed
128

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

137
   np_rows = nprocs/np_cols;
Andreas Marek's avatar
Andreas Marek committed
138

139 140 141 142
   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
143

144 145
   /* set up blacs */
   /* convert communicators before */
146
#ifdef WITH_MPI
Andreas Marek's avatar
Andreas Marek committed
147
   my_mpi_comm_world = MPI_Comm_c2f(MPI_COMM_WORLD);
148 149 150
#else
   my_mpi_comm_world = 1;
#endif
151
   set_up_blacsgrid_f(my_mpi_comm_world, np_rows, np_cols, 'C', &my_blacs_ctxt, &my_prow, &my_pcol);
Andreas Marek's avatar
Andreas Marek committed
152 153 154 155 156 157 158

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

159 160
   /* get the ELPA row and col communicators. */
   /* These are NOT usable in C without calling the MPI_Comm_f2c function on them !! */
161
#ifdef WITH_MPI
Andreas Marek's avatar
Andreas Marek committed
162
   my_mpi_comm_world = MPI_Comm_c2f(MPI_COMM_WORLD);
163
#endif
164
   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
165 166 167 168 169 170 171 172 173

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

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

174
   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);
Andreas Marek's avatar
Andreas Marek committed
175 176 177 178 179 180 181

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

182
   /* allocate the matrices needed for elpa */
Andreas Marek's avatar
Andreas Marek committed
183 184 185 186 187
   if (myid == 0) {
     printf("\n");
     printf("Allocating matrices with na_rows=%d and na_cols=%d\n",na_rows, na_cols);
     printf("\n");
   }
188
#ifdef DOUBLE_PRECISION_COMPLEX
Andreas Marek's avatar
Andreas Marek committed
189 190 191
   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));
Andreas Marek's avatar
Andreas Marek committed
192
   ev = malloc(na*sizeof(double));
193 194 195 196 197 198 199
#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));
   ev = malloc(na*sizeof(float));
#endif
#ifdef DOUBLE_PRECISION_COMPLEX
200
   prepare_matrix_random_complex_double_f(na, myid, na_rows, na_cols, sc_desc, a, z, as);
201
#else
202
   prepare_matrix_random_complex_single_f(na, myid, na_rows, na_cols, sc_desc, a, z, as);
203
#endif
Andreas Marek's avatar
Andreas Marek committed
204

Andreas Marek's avatar
Andreas Marek committed
205 206
   if (myid == 0) {
     printf("\n");
Andreas Marek's avatar
Andreas Marek committed
207
     printf("Entering ELPA 2stage complex solver\n");
Andreas Marek's avatar
Andreas Marek committed
208 209
     printf("\n");
   }
210
#ifdef WITH_MPI
Andreas Marek's avatar
Andreas Marek committed
211
   mpierr = MPI_Barrier(MPI_COMM_WORLD);
212
#endif
213
   useGPU = 0;
214
   THIS_COMPLEX_ELPA_KERNEL_API = ELPA_2STAGE_COMPLEX_DEFAULT;
215
#ifdef DOUBLE_PRECISION_COMPLEX
216
   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, useGPU);
217
#else
218
   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, useGPU);
219
#endif
Andreas Marek's avatar
Andreas Marek committed
220 221 222

   if (success != 1) {
     printf("error in ELPA solve \n");
223
#ifdef WITH_MPI
Andreas Marek's avatar
Andreas Marek committed
224
     mpierr = MPI_Abort(MPI_COMM_WORLD, 99);
225 226
#else
     exit(99);
227
#endif
Andreas Marek's avatar
Andreas Marek committed
228 229 230 231 232
   }


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

237
   /* check the results */
238
#ifdef DOUBLE_PRECISION_COMPLEX
239
   status = check_correctness_evp_numeric_residuals_complex_double_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
240
#else
241
   status = check_correctness_evp_numeric_residuals_complex_single_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
242
#endif
Andreas Marek's avatar
Andreas Marek committed
243 244 245 246
   if (status !=0){
     printf("The computed EVs are not correct !\n");
   }
   if (status ==0){
Andreas Marek's avatar
Andreas Marek committed
247 248 249
     if (myid == 0) {
       printf("All ok!\n");
     }
Andreas Marek's avatar
Andreas Marek committed
250 251
   }

252 253 254 255
   free(sc_desc);
   free(a);
   free(z);
   free(as);
256
   free(ev);
257

258
#ifdef WITH_MPI
259
   MPI_Finalize();
260
#endif
261
   return 0;
Andreas Marek's avatar
Andreas Marek committed
262
}