elpa2_test_complex_c_version.c 6.78 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>
Andreas Marek's avatar
Andreas Marek committed
54
#include <test/shared_sources/generated.h>
Andreas Marek's avatar
Andreas Marek committed
55
#include <complex.h>
Andreas Marek's avatar
Andreas Marek committed
56

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

65
   int status;
Andreas Marek's avatar
Andreas Marek committed
66

67
   int np_cols, np_rows, np_colsStart;
Andreas Marek's avatar
Andreas Marek committed
68

69
   int my_blacs_ctxt, nprow, npcol, my_prow, my_pcol;
Andreas Marek's avatar
Andreas Marek committed
70

71
   int mpierr;
Andreas Marek's avatar
Andreas Marek committed
72

73 74
   int my_mpi_comm_world;
   int mpi_comm_rows, mpi_comm_cols;
Andreas Marek's avatar
Andreas Marek committed
75

76
   int info, *sc_desc;
Andreas Marek's avatar
Andreas Marek committed
77

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

Andreas Marek's avatar
Andreas Marek committed
81 82 83
   complex double *a, *z, *as, *tmp1, *tmp2;

   double *ev, *xr;
Andreas Marek's avatar
Andreas Marek committed
84

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
   int THIS_COMPLEX_ELPA_KERNEL_API;
90
#ifdef WITH_MPI
91 92 93
   MPI_Init(&argc, &argv);
   MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
94 95 96 97 98
#else
   nprocs = 1;
   myid =0;
   MPI_COMM_WORLD=1;
#endif
99 100 101
   na = 1000;
   nev = 500;
   nblk = 16;
Andreas Marek's avatar
Andreas Marek committed
102

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

114
     printf("\n");
Andreas Marek's avatar
Andreas Marek committed
115

116
   }
Andreas Marek's avatar
Andreas Marek committed
117

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

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

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

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

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

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

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

173
   /* allocate the matrices needed for elpa */
Andreas Marek's avatar
Andreas Marek committed
174 175 176 177 178 179
   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
180 181 182 183 184
   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
185 186 187 188


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

Andreas Marek's avatar
Andreas Marek committed
189 190
   tmp1  = malloc(na_rows*na_cols*sizeof(complex double));
   tmp2 = malloc(na_rows*na_cols*sizeof(complex double));
Andreas Marek's avatar
Andreas Marek committed
191 192 193

   iseed = malloc(4096*sizeof(int));

Andreas Marek's avatar
Andreas Marek committed
194 195 196
   prepare_matrix_complex_from_fortran(na, myid, na_rows, na_cols, sc_desc, iseed, xr, a, z, as);

   free(xr);
Andreas Marek's avatar
Andreas Marek committed
197 198 199

   if (myid == 0) {
     printf("\n");
Andreas Marek's avatar
Andreas Marek committed
200
     printf("Entering ELPA 2stage complex solver\n");
Andreas Marek's avatar
Andreas Marek committed
201 202
     printf("\n");
   }
203
#ifdef WITH_MPI
Andreas Marek's avatar
Andreas Marek committed
204
   mpierr = MPI_Barrier(MPI_COMM_WORLD);
205
#endif
Andreas Marek's avatar
Andreas Marek committed
206 207
   THIS_COMPLEX_ELPA_KERNEL_API = ELPA2_COMPLEX_KERNEL_GENERIC;
   success = elpa_solve_evp_complex_2stage(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);
Andreas Marek's avatar
Andreas Marek committed
208 209 210

   if (success != 1) {
     printf("error in ELPA solve \n");
211
#ifdef WITH_MPI
Andreas Marek's avatar
Andreas Marek committed
212
     mpierr = MPI_Abort(MPI_COMM_WORLD, 99);
213
#endif
Andreas Marek's avatar
Andreas Marek committed
214 215 216 217 218
   }


   if (myid == 0) {
     printf("\n");
Andreas Marek's avatar
Andreas Marek committed
219
     printf("2stage ELPA complex solver complete\n");
Andreas Marek's avatar
Andreas Marek committed
220 221 222
     printf("\n");
   }

223
   /* check the results */
Andreas Marek's avatar
Andreas Marek committed
224
   status = check_correctness_complex_from_fortran(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid, tmp1, tmp2);
Andreas Marek's avatar
Andreas Marek committed
225 226 227 228 229

   if (status !=0){
     printf("The computed EVs are not correct !\n");
   }
   if (status ==0){
Andreas Marek's avatar
Andreas Marek committed
230 231 232
     if (myid == 0) {
       printf("All ok!\n");
     }
Andreas Marek's avatar
Andreas Marek committed
233 234
   }

235 236 237 238 239 240 241
   free(sc_desc);
   free(a);
   free(z);
   free(as);

   free(tmp1);
   free(tmp2);
242
#ifdef WITH_MPI
243
   MPI_Finalize();
244
#endif
245
   return 0;
Andreas Marek's avatar
Andreas Marek committed
246
}