elpa1_test_real_c_version.c 6.02 KB
Newer Older
Andreas Marek's avatar
Andreas Marek committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 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
/*     This file is part of ELPA. */
/*  */
/*     The ELPA library was originally created by the ELPA consortium, */
/*     consisting of the following organizations: */
/*  */
/*     - Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG), */
/*     - 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: */
/*     http://elpa.rzg.mpg.de/ */
/*  */
/*     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>

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

56
   int na, nev, nblk;
Andreas Marek's avatar
Andreas Marek committed
57

58
   int status;
Andreas Marek's avatar
Andreas Marek committed
59

60
   int np_cols, np_rows, np_colsStart;
Andreas Marek's avatar
Andreas Marek committed
61

62
   int my_blacs_ctxt, nprow, npcol, my_prow, my_pcol;
Andreas Marek's avatar
Andreas Marek committed
63

64
   int mpierr;
Andreas Marek's avatar
Andreas Marek committed
65

66 67
   int my_mpi_comm_world;
   int mpi_comm_rows, mpi_comm_cols;
Andreas Marek's avatar
Andreas Marek committed
68

69
   int info, *sc_desc;
Andreas Marek's avatar
Andreas Marek committed
70

71 72
   int na_rows, na_cols;
   double startVal;
Andreas Marek's avatar
Andreas Marek committed
73

74
   double *a, *z, *as, *ev, *tmp1, *tmp2;
Andreas Marek's avatar
Andreas Marek committed
75

76
   int *iseed;
Andreas Marek's avatar
Andreas Marek committed
77

78
   int success;
Andreas Marek's avatar
Andreas Marek committed
79

80 81 82
   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
83

84 85 86
   na = 1000;
   nev = 500;
   nblk = 16;
Andreas Marek's avatar
Andreas Marek committed
87

88 89 90 91 92 93 94 95
   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
96 97
     printf("as it's Fortran counterpart. It's only purpose is to show how \n");
     printf("to evoke ELPA1 from a c programm\n");
98
     printf("\n");
Andreas Marek's avatar
Andreas Marek committed
99

100
   }
Andreas Marek's avatar
Andreas Marek committed
101

102
   status = 0;
Andreas Marek's avatar
Andreas Marek committed
103

104 105 106 107 108 109 110
   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
111

112
   np_rows = nprocs/np_cols;
Andreas Marek's avatar
Andreas Marek committed
113

114 115 116 117
   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
118

119 120
   /* set up blacs */
   /* convert communicators before */
Andreas Marek's avatar
Andreas Marek committed
121 122 123 124 125 126 127 128 129
   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");
   }

130 131
   /* 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
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
   my_mpi_comm_world = MPI_Comm_c2f(MPI_COMM_WORLD);
   mpierr = elpa_get_communicators(my_mpi_comm_world, my_prow, my_pcol, &mpi_comm_rows, &mpi_comm_cols);

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

151
   /* allocate the matrices needed for elpa */
Andreas Marek's avatar
Andreas Marek committed
152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 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");
   }

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

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

   prepare_matrix_real_from_fortran(na, myid, na_rows, na_cols, sc_desc, iseed, a, z, as);

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

   mpierr = MPI_Barrier(MPI_COMM_WORLD);

180
   success = elpa_solve_evp_real_1stage(na, nev, a, na_rows, ev, z, na_rows, nblk, na_cols, mpi_comm_rows, mpi_comm_cols);
Andreas Marek's avatar
Andreas Marek committed
181 182 183 184 185 186 187 188 189 190 191 192 193

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

194
   /* check the results */
Andreas Marek's avatar
Andreas Marek committed
195 196 197 198 199 200 201 202 203
   status = check_correctness_real_from_fortran(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid, tmp1, tmp2);

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

204 205 206 207 208 209 210
   free(sc_desc);
   free(a);
   free(z);
   free(as);

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

212
   MPI_Finalize();
Andreas Marek's avatar
Andreas Marek committed
213

214
   return 0;
Andreas Marek's avatar
Andreas Marek committed
215
}