elpa1_test_real_c_version.c 7.26 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 52
/*  */
/*     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>

53 54 55
#include "test/shared_sources/generated.h"

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_REAL
77
   double *a, *z, *as, *ev, *tmp1, *tmp2;
78 79 80
#else
   float *a, *z, *as, *ev, *tmp1, *tmp2;
#endif
81
   int *iseed;
Andreas Marek's avatar
Andreas Marek committed
82

83
   int success;
Andreas Marek's avatar
Andreas Marek committed
84

85 86 87
   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
88

89 90 91
   na = 1000;
   nev = 500;
   nblk = 16;
Andreas Marek's avatar
Andreas Marek committed
92

93 94 95 96 97 98 99 100
   if (myid == 0) {
     printf("This is the c version of an ELPA test-programm\n");
     printf("\n");
     printf("It will call the 1stage ELPA real solver for an\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");
Andreas Marek's avatar
Andreas Marek committed
101 102
     printf("as it's Fortran counterpart. It's only purpose is to show how \n");
     printf("to evoke ELPA1 from a c programm\n");
103
     printf("\n");
104 105 106 107 108 109
#ifdef DOUBLE_PRECISION_REAL
     printf("The double precision version of ELPA1 is used\n");
#else
     printf("The single precision version of ELPA1 is used\n");
#endif
     printf("\n");
Andreas Marek's avatar
Andreas Marek committed
110

111
   }
Andreas Marek's avatar
Andreas Marek committed
112

113
   status = 0;
Andreas Marek's avatar
Andreas Marek committed
114

115 116 117 118 119 120 121
   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
122

123
   np_rows = nprocs/np_cols;
Andreas Marek's avatar
Andreas Marek committed
124

125 126 127 128
   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
129

130 131
   /* set up blacs */
   /* convert communicators before */
Andreas Marek's avatar
Andreas Marek committed
132 133 134 135 136 137 138 139 140
   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");
   }

141 142
   /* 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
143
   my_mpi_comm_world = MPI_Comm_c2f(MPI_COMM_WORLD);
144
   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
145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161

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

162
   /* allocate the matrices needed for elpa */
Andreas Marek's avatar
Andreas Marek committed
163 164 165 166 167
   if (myid == 0) {
     printf("\n");
     printf("Allocating matrices with na_rows=%d and na_cols=%d\n",na_rows, na_cols);
     printf("\n");
   }
168
#ifdef DOUBLE_PRECISION_REAL
Andreas Marek's avatar
Andreas Marek committed
169 170 171 172 173 174 175 176 177
   a  = malloc(na_rows*na_cols*sizeof(double));
   z  = malloc(na_rows*na_cols*sizeof(double));
   as = malloc(na_rows*na_cols*sizeof(double));


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

   tmp1  = malloc(na_rows*na_cols*sizeof(double));
   tmp2 = malloc(na_rows*na_cols*sizeof(double));
178 179 180 181
#else
   a  = malloc(na_rows*na_cols*sizeof(float));
   z  = malloc(na_rows*na_cols*sizeof(float));
   as = malloc(na_rows*na_cols*sizeof(float));
Andreas Marek's avatar
Andreas Marek committed
182 183


184
   ev = malloc(na*sizeof(float));
Andreas Marek's avatar
Andreas Marek committed
185

186 187 188 189 190 191 192 193 194 195
   tmp1  = malloc(na_rows*na_cols*sizeof(float));
   tmp2 = malloc(na_rows*na_cols*sizeof(float));
#endif

   iseed = malloc(4096*sizeof(int));
#ifdef DOUBLE_PRECISION_REAL
   prepare_matrix_real_from_fortran_double_precision(na, myid, na_rows, na_cols, sc_desc, iseed, a, z, as);
#else
   prepare_matrix_real_from_fortran_single_precision(na, myid, na_rows, na_cols, sc_desc, iseed, a, z, as);
#endif
Andreas Marek's avatar
Andreas Marek committed
196 197 198 199 200 201 202
   if (myid == 0) {
     printf("\n");
     printf("Entering ELPA 1stage real solver\n");
     printf("\n");
   }

   mpierr = MPI_Barrier(MPI_COMM_WORLD);
203 204 205 206 207
#ifdef DOUBLE_PRECISION_REAL
   success = elpa_solve_evp_real_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_real_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
208 209 210 211 212 213 214 215 216 217 218 219
   if (success != 1) {
     printf("error in ELPA solve \n");
     mpierr = MPI_Abort(MPI_COMM_WORLD, 99);
   }


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

220
#ifdef DOUBLE_PRECISION_REAL
221
   /* check the results */
222 223 224 225
   status = check_correctness_real_from_fortran_double_precision(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid, tmp1, tmp2);
#else
   status = check_correctness_real_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
226 227 228 229 230 231 232
   if (status !=0){
     printf("The computed EVs are not correct !\n");
   }
   if (status ==0){
     printf("All ok!\n");
   }

233 234 235 236 237 238 239
   free(sc_desc);
   free(a);
   free(z);
   free(as);

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

241
   MPI_Finalize();
Andreas Marek's avatar
Andreas Marek committed
242

243
   return 0;
Andreas Marek's avatar
Andreas Marek committed
244
}