/** \file helpers.c
 * Malloc-, Realloc-, Free-wrappers and other small routines.
 * 
 * This file contains small helpful routines that do such omnipresent tasks
 * as allocation memory Malloc(), Malloci(), Mallocii(), reallocating Realloc(),
 * Realloci(), Reallocii() or disallocating Free() with checks error messages
 * in the case of failure. Speed Measuring is initiated InitSpeedMeasure() by setting
 * SpeedStruct entries to zero,
 * completed for output to file CompleteSpeedMeasure() and actually performed
 * SpeedMeasure() herein.\n
 * RemoveEverything() frees all memory that was alloceted during the course of the
 * programme.\n
 * StartParallel() initializes process groups and ranks, while WaitForOtherProcs()
 * waits for the other running processes between the end of init phase and the begin
 * of the calculation.\n
 * Finally, GetRGB() does a colour translation of an RGB value from integer but it's
 * not used anymore.
 * 
  Project: ParallelCarParrinello
 \author Jan Hamaekers
 \date 2000

  File: helpers.c
  $Id: helpers.c,v 1.39 2007-02-12 09:44:55 foo Exp $
*/

#include <stdlib.h>
#include <stdio.h>
#include <stddef.h>
#include <math.h>
#include <string.h>
#include "data.h"
#include "errors.h"
#include "helpers.h"
#include "grad.h"
#include "pseudo.h"
#include "ions.h"
#include "output.h"
#include "run.h"
#include "gramsch.h"
#include "mymath.h"
#include "myfft.h"
#include "perturbed.h"

/** Output of a debug message to stderr.
 * \param *P Problem at hand, points to ParallelSimulationData#me
 * \param output output string
 */
#ifdef HAVE_DEBUG
void debug_in(struct Problem *P, const char *output, const char *file, const int line) {
  if (output) fprintf(stderr,"(%i) DEBUG: in %s at line %i: %s\n",P->Par.me, file, line, output);
}
#else
void debug_in(struct Problem *P, const char *output, const char *file, const int line) {}  // print nothing
#endif

/** Malloc with error check.
 * malloc() wrapper, printing error on null pointer
 * \param size size of memory to be allocated
 * \param output error message to be given on failure
 */
void* Malloc(size_t size, const char* output)
{
  void* dummy = malloc(size);
  if (!dummy)
    Error(MallocError, output);
  return dummy;
}

/** Malloc string array and set its length to the allocated length.
 * \param *output message if malloc fails
 * \param size number of chars to alloc for \a *buffer
 * \return pointer to string array
 */
char* MallocString(size_t size, const char* output) {
  int i;
  char *buffer;
  buffer = Malloc(sizeof(char) * (size+1), output); // alloc
  for (i=0;i<size;i++)  // reset
    buffer[i] = i % 2 == 0 ? 'p': 'c';
  buffer[size] = '\0'; // and set length marker on its end
  return(buffer);
}


/** Malloc with error check and 1 int.
 * malloc() wrapper, printing error on null pointer
 * \param size size of memory to be allocated
 * \param output error message to be given on failure (containing one %d or similar)
 * \param i integer value to be put into error message (%d in output)
 */
void* Malloci(size_t size, const char* output, int i)
{
  char dummyoutput[MAXDUMMYSTRING];
  void* dummy = (void *) malloc(size);
  if (!dummy) {
    sprintf(dummyoutput,output,i);
    Error(MallocError, dummyoutput);
  }
  return dummy;
}

/** Malloc with error check and 2 ints.
 * malloc() wrapper, printing error on null pointer
 * \param size size of memory to be allocated
 * \param output error message to be given on failure (containing two %d or similar)
 * \param i first integer value to be put into error message (%d in output)
 * \param j second integer value to be put into error message (%d in output)
 */
void* Mallocii(size_t size, const char* output, int i, int j)
{
  char dummyoutput[MAXDUMMYSTRING];
  void* dummy = (void *) malloc(size);
  if (!dummy) {
    sprintf(dummyoutput,output,i,j);
    Error(MallocError, dummyoutput);
  }
  return dummy;
}

/** Change size of allocated memory range.
 * Wrapper for realloc(), giving error message on failure
 * \param pointer points to allocated memory
 * \param size new desired size
 * \param output error message on failure
 */
void* Realloc(void* pointer, size_t size, const char* output)
{
  void *dummy = (void *) realloc(pointer, size);
  if (!dummy)
    Error(MallocError, output);
  return dummy;
}

/** Change size of allocated memory range, 1 int in error msg.
 * Wrapper for realloc(), giving error message on failure
 * \param pointer points to allocated memory
 * \param size new desired size
 * \param output error message on failure (containing one %d or similar)
 * \param i integer in error message
 */
void* Realloci(void* pointer, size_t size, const char* output, int i)
{
  char dummyoutput[MAXDUMMYSTRING];
  void *dummy = (void *) realloc(pointer, size);
  if (!dummy) {
    sprintf(dummyoutput, output, i);
    Error(MallocError, dummyoutput);
  }
  return dummy;
}

/** Change size of allocated memory range, 2 int in error msg.
 * Wrapper for realloc(), giving error message on failure
 * \param pointer points to allocated memory
 * \param size new desired size
 * \param output error message on failure (containing two %d or similar)
 * \param i first integer in error message
 * \param j second integer in error message
 */
void* Reallocii(void* pointer, size_t size, const char* output, int i, int j)
{
  char dummyoutput[MAXDUMMYSTRING];
  void *dummy = (void *) realloc(pointer, size);
  if (!dummy) {
    sprintf(dummyoutput,output,i,j);
    Error(MallocError, dummyoutput);
  }
  return dummy;
}

/** Free memory with error check.
 * Wrapper for free(), free is only called on given non-null pointer
 * \param *ptr pointer to allocated memory region
 * \param *output message to print on failure
 */
void Free (void *ptr, const char* output) {
  if (ptr != NULL)
    free(ptr);
  else 
    fprintf(stderr,"Free not necessary: %s\n", output);
}

/** Translates a given integer into color code
 * The value is interpreted by powers of the range into an rgb set of
 * values.
 * \param i specifies value within the range
 * \param max	specifies color range
 * \param rgb pointer to double[3] 
 */
void GetRGB(int i, unsigned int max, double* rgb) {
  int basis,basis2,basis3;
  int maxRGB,RGBint;
  int rgbint[3];
  i++;
  if (max <= 1) max = 2; /* keine Division durch 0 */
  basis = ceil(pow((double)(max+1),1./3.));
  basis2 = basis*basis;
  basis3 = basis2*basis;
  maxRGB = basis3-1;
  /* Koordinatentrans von 1..max -> 1..maxRGB */
  RGBint = floor(1. + (double)(maxRGB - 1)/(double)(max-1)*(double)(i-1));
  /* Transferiere RGBint jetzt in Basisdarstellung */
  rgbint[2] = RGBint/basis2;
  RGBint %= basis2;
  rgbint[1] = RGBint/basis;
  rgbint[0] = RGBint%basis;
  /* Jetzt Basisdarstellung in rgb darstellung */
  rgb[2] = rgbint[2]*(1./(basis-1));
  rgb[1] = rgbint[1]*(1./(basis-1));
  rgb[0] = rgbint[0]*(1./(basis-1));
}

/** Waits for other processes to finish.
 * Simply waits until all processes reach this routine for the MPI_Gather to work 
 * \param *P Problem at hand
 * \param out	1 - output ready note, 0 - no output
 */
void WaitForOtherProcs(struct Problem *P, int out) {
  /* Warte, bis alle fertig sind! */
  int buffint = 0;
  MPI_Gather(&buffint, 0, MPI_INT, &buffint, 0, MPI_INT, 0, P->Par.comm);
  if(out && P->Call.out[MinOut]) fprintf(stderr,"\nAll procs ready: %i\n",P->Par.me);
}

/** Initialization for parallel computing.
 * determines process ids and defines mpi groups (aka Communicators), being
 * the ones working on the same Psi (when doing fft) and those working on the
 * same section (transposed GramSchmidt)
 * \param *P Problem at hand
 * \param **argv is actually unnecessary
 */
void StartParallel (struct Problem *P, char **argv) {
  int info;
  char processor_name[MPI_MAX_PROCESSOR_NAME];
  char dummy[MAXDUMMYSTRING];
  P->Par.procs = P->Par.proc[PEPsi]*P->Par.proc[PEGamma];
  P->Par.me = -1; /* eigene Position noch unklar */
  /* MPI Initialisierung */
  argv = argv; /* wird nicht benoetigt */
  /* mpi Gruppen definieren */
  P->Par.world = MPI_COMM_WORLD;
  MPI_Comm_dup(P->Par.world, &P->Par.comm);
  MPI_Comm_size(P->Par.comm, &info);
  if (info != P->Par.procs) {
    sprintf(dummy,"StartParallel. Wrong number of processes started %i != %i",info, P->Par.procs);
    Error(SomeError, dummy);
  }
  MPI_Comm_rank(P->Par.comm, &P->Par.mytid);
  MPI_Get_processor_name(processor_name, &info);
  P->Par.me = P->Par.mytid;
  P->Par.mypos[PEGamma] = P->Par.me % P->Par.proc[PEGamma];
  P->Par.mypos[PEPsi] = P->Par.me / P->Par.proc[PEGamma];
  if(P->Call.out[NormalOut]) fprintf(stderr,"Process %d(%i,%i) on %s \n", P->Par.mytid, P->Par.mypos[PEPsi], P->Par.mypos[PEGamma], processor_name);
  if (P->Par.me == 0) { /* Master */
    if (P->Par.procs > 9999)
      fprintf(stderr, "(%d)Warning: procs>9999 will cause trouble with output files\n", P->Par.mytid);
  }
  if(P->Call.out[MinOut]) fprintf(stderr,"Info: %i processes started\n", P->Par.procs); /* Prozesse gestartet */


  /* Erstellung der Communicatoren */
  switch (P->Lat.Psi.Use) {
  case UseSpinDouble:
    P->Par.my_color_comm_ST = (int)SpinDouble;
    P->Par.me_comm_ST = P->Par.me;
    P->Par.Max_my_color_comm_ST = 1;
    P->Par.Max_me_comm_ST = P->Par.procs;
    break;
  case UseSpinUpDown:
    P->Par.my_color_comm_ST = P->Par.me / (P->Par.procs/2);
    P->Par.me_comm_ST = P->Par.me % (P->Par.procs/2);
    P->Par.Max_my_color_comm_ST = 2;
    P->Par.Max_me_comm_ST = P->Par.procs/2;
    break;
  }
  MPI_Comm_split (P->Par.comm, P->Par.my_color_comm_ST, P->Par.me_comm_ST, &P->Par.comm_ST);
  if (P->Lat.Psi.Use == UseSpinUpDown) {
    MPI_Intercomm_create ( P->Par.comm_ST, 0, P->Par.comm, 
			   (P->Par.my_color_comm_ST ? 0 : P->Par.procs/2)
			   , InterTag1, &P->Par.comm_STInter);
  }

  /* Alle Procs mit gleichen Psis (Sind alle an einer fft beteiligt)*/
  P->Par.my_color_comm_ST_Psi = P->Par.me_comm_ST / P->Par.proc[PEGamma];
  P->Par.me_comm_ST_Psi = P->Par.me_comm_ST % P->Par.proc[PEGamma];
  P->Par.Max_my_color_comm_ST_Psi =  P->Par.Max_me_comm_ST/P->Par.proc[PEGamma];
  P->Par.Max_me_comm_ST_Psi = P->Par.proc[PEGamma];
  MPI_Comm_split (P->Par.comm_ST, P->Par.my_color_comm_ST_Psi, P->Par.me_comm_ST_Psi, &P->Par.comm_ST_Psi);

  /* Alle Procs mit gleichen PsisAbschnitt (Transposed Psi - GramSch)*/
  P->Par.my_color_comm_ST_PsiT = P->Par.me_comm_ST % P->Par.proc[PEGamma];
  P->Par.me_comm_ST_PsiT = P->Par.me_comm_ST / P->Par.proc[PEGamma];
  P->Par.Max_my_color_comm_ST_PsiT = P->Par.proc[PEGamma]; 
  P->Par.Max_me_comm_ST_PsiT = P->Par.Max_me_comm_ST/P->Par.proc[PEGamma];
  MPI_Comm_split (P->Par.comm_ST, P->Par.my_color_comm_ST_PsiT, P->Par.me_comm_ST_PsiT, &P->Par.comm_ST_PsiT);
}


/** Removes everything from memory.
 * Frees memory and exits program
 * \param *P Problem at hand
 */
void RemoveEverything(struct Problem *P)
{
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  int i,d, type;
  if (P->Call.out[NormalOut]) fprintf(stderr,"(%i)RemoveEverything !\n",P->Par.me);
  RemoveGradients(P, &P->Grad);
  RemovePseudoRead(P);
  RemoveIonsRead(&P->Ion);
  fft_3d_destroy_plan(P,P->Lat.plan);
  if (P->Files.MeOutVis) {
    Free(P->Files.PosTemp, "RemoveEverything: P->Files.PosTemp");
    Free(P->Files.OutputPosType, "RemoveEverything: P->Files.OutputPosType");
    //Free(P->Files.OutVisStep, "RemoveEverything: P->Files.OutVisStep");
  }
  Free(P->Files.mainname, "RemoveEverything: P->Files.mainname");
  Free(P->Files.mainpath, "RemoveEverything: P->Files.mainpath");
  Free(P->Call.MainParameterFile, "RemoveEverything: P->Call.MainParameterFile");
	if (P->Call.ForcesFile != NULL)
		Free(P->Call.ForcesFile, "RemoveEverything: P->Call.ForcesFile");
  Free(Lat->MaxNoOfnFields, "RemoveEverything: Lat->MaxNoOfnFields");
  for (i=0; i < P->Lat.MaxLevel; i++) {
    Free(Lat->Lev[i].AllMaxG, "RemoveEverything: Lat->Lev[i].AllMaxG");
    Free(Lat->NFields[i], "RemoveEverything: Lat->NFields");
    Free(Lat->Lev[i].GArray, "RemoveEverything: Lat->Lev[i].GArray");
    Free(Lat->Lev[i].HashG, "RemoveEverything: Lat->Lev[i].HashG");
    if (Lat->Lev[i].MaxDoubleG)
      Free(Lat->Lev[i].DoubleG, "RemoveEverything: Lat->Lev[i].DoubleG");
    if (i != 0) {
      Free(Lat->Lev[i].PosFactorUp, "RemoveEverything: Lat->Lev[i].PosFactorUp");
    } else {
      for (d=0; d < (P->R.DoPerturbation == 1 ? MaxDensityTypes : MaxInitDensityTypes); d++) {
      	if (Lat->Lev[i].Dens->DensityArray[d] != NULL) Free(Lat->Lev[i].Dens->DensityArray[d], "RemoveEverything: Lat->Lev[i].Dens->DensityArray[d]");
      	if (Lat->Lev[i].Dens->DensityCArray[d] != NULL) Free(Lat->Lev[i].Dens->DensityCArray[d], "RemoveEverything: Lat->Lev[i].Dens->DensityCArray[d]");
      }
    }
    /* if (i != Lat->MaxLevel-1) */
      Free(Lat->Lev[i].Dens, "RemoveEverything: Lat->Lev[i].Dens");
  }
  for (i=1; i < Lat->MaxLevel; i++) {
    Free(Lat->Lev[i].LPsi->LocalPsi, "RemoveEverything: Lat->Lev[i].LPsi->LocalPsi");
    Free(Lat->Lev[i].LPsi->OldLocalPsi, "RemoveEverything: Lat->Lev[i].LPsi->OldLocalPsi");
    if (i == 1) {
      Free(Lat->Lev[i].LPsi->PsiDat, "RemoveEverything: Lat->Lev[i].LPsi->PsiDat");
      Free(Lat->Lev[i].LPsi->OldPsiDat, "RemoveEverything: Lat->Lev[i].LPsi->OldPsiDat");
    }
    Free(Lat->Lev[i].LPsi, "RemoveEverything: Lat->Lev[i].LPsi");
  }
  for (i=0;i<Lat->Psi.NoOfTotalPsis;i++) {
    Free(Psi->lambda[i], "RemoveEverything: Psi->lambda[i]");
  }  
  for (i=0;i<Lat->Psi.NoOfPsis;i++) {
    Free(Psi->Overlap[i], "RemoveEverything: Psi->Overlap[i]");
  }
  Free(Psi->lambda, "RemoveEverything: Psi->lambda");
  Free(Psi->Overlap, "RemoveEverything: Psi->Overlap");
  Free(Psi->AllPsiStatus, "RemoveEverything: Psi->AllPsiStatus");
  Free(Psi->AllPsiStatusForSort, "RemoveEverything: Psi->AllPsiStatusForSort");
  Free(Psi->LocalPsiStatus, "RemoveEverything: Psi->LocalPsiStatus");
  Free(Psi->AllLocalNo, "RemoveEverything: Psi->AllLocalNo");
  Free(Psi->RealAllLocalNo, "RemoveEverything: Psi->RealAllLocalNo");
  Free(Psi->TempSendA, "RemoveEverything: Psi->TempSendA");
  Free(Psi->AddData, "RemoveEverything: Psi->AddData");
  Free(Psi->AllActualLocalPsiNo, "RemoveEverything: Psi->AllActualLocalPsiNo");
  Free(Psi->AllOldActualLocalPsiNo, "RemoveEverything: Psi->AllOldActualLocalPsiNo");
  Free(Lat->Lev, "RemoveEverything: Lat->Lev");
  Free(Lat->LevelSizes, "RemoveEverything: Lat->LevelSizes");
  Free(Lat->NFields, "RemoveEverything: Lat->NFields");
  if (Lat->RT.Use == UseRT) {
    Free(Lat->RT.Coeff, "RemoveEverything: Lat->RT.Coeff");
    Free(Lat->RT.RTLaplaceS, "RemoveEverything: Lat->RT.RTLaplaceS");
    Free(Lat->RT.RTLaplace0, "RemoveEverything: Lat->RT.RTLaplace0");
    for (i=0; i< MAXRTPOSFAC; i++) Free(Lat->RT.PosFactor[i], "RemoveEverything: Lat->RT.PosFactor[i]");
    for (i=0; i< MAXRTARRAYS; i++) {
      Free(Lat->RT.DensityC[i], "RemoveEverything: Lat->RT.DensityC[i]");
    }
    Free(Lat->RT.TempC, "RemoveEverything: Lat->RT.TempC");
  }
  for (type=Occupied;type<Extra;type++)
    for (i=0; i<MAXALLPSIENERGY; i++) {
      Free(Lat->Energy[type].PsiEnergy[i], "RemoveEverything: Lat->Energy[type].PsiEnergy[i]");
    }
  for (i=Occupied;i<Extra;i++)
    Free(P->R.MinimisationName[i], "RemoveEverything: P->R.MinimisationName[i]");
  Free(P->R.MinimisationName, "RemoveEverything: P->R.MinimisationName");
  MPI_Comm_free(&P->Par.comm_ST_PsiT);
  //debug(P,"comm_ST_PsiT");  
  MPI_Comm_free(&P->Par.comm_ST_Psi);
  //debug(P,"comm_ST_Psi");  
  if (Psi->Use == UseSpinUpDown) MPI_Comm_free(&P->Par.comm_STInter);
  MPI_Comm_free(&P->Par.comm_ST);
  //debug(P,"comm_ST");  
  MPI_Comm_free(&P->Par.comm);
  //debug(P,"comm");  
  Free(P, "RemoveEverything: P");
  FreeMPI_OnePsiElement();
  //debug(P,"OnePsiElement");  
  
  // Free string names from Files structure
  Free(P->Files.filename, "RemoveEverything: P->Files.filename");
  Free(P->Files.default_path, "RemoveEverything: P->Files.default_path");
  Free(P->Files.pseudopot_path, "RemoveEverything: P->Files.pseudopot_path");
}

static const char suffixspeed[] = ".speed";		//!< suffix of the FileData#SpeedFile where the timings are written to


/** Initializes the timer array.
 * Sets every time, average, min/max and deviation to zero.
 * SpeedStep is set to 1.
 * Opens Filedata::SpeedFile for appended writing
 * \param *P Problem at hand
 */
void InitSpeedMeasure(struct Problem *P)
{
  struct FileData *F = &P->Files;
  struct SpeedStruct *S = &P->Speed;
  int i;
  F->SpeedFile = NULL;
  if (P->Par.me == 0) 
    OpenFile(P, &F->SpeedFile, suffixspeed, "a",P->Call.out[ReadOut]);
  for (i=0; i < MAXTIMETYPES; i++)
    S->SpeedStep[i] = 1;
  SetArrayToDouble0(S->time1,MAXTIMETYPES);
  SetArrayToDouble0(S->time2,MAXTIMETYPES);
  SetArrayToDouble0(S->time,MAXTIMETYPES);
  SetArrayToDouble0(S->average,MAXTIMETYPES);
  SetArrayToDouble0(S->min,MAXTIMETYPES); 
  SetArrayToDouble0(S->max,MAXTIMETYPES);
  SetArrayToDouble0(S->stddev,MAXTIMETYPES); 
  S->InitSteps = 0;
  S->LevUpSteps = 0;
  S->Steps = 0;
}

/** Measure speed of a routine.
 * Points to SpeedStruct of the Problem, gets current time and depending on
 * \a TOT (whether we start or stop the watch) and the timing group \a TT
 * puts time in SpeedStruct::time1 (start) or in SpeedStruct::time2 (stop)
 * and in SpeedStruct::time the difference between the two is added.
 * \param *P Problem at hand containing SpeedStruct
 * \param TT TimeTypes
 * \param TOT TimeDotypes
 */
void SpeedMeasure(struct Problem *P, enum TimeTypes TT, enum TimeDoTypes TOT)
{ 
  struct SpeedStruct *S = &P->Speed;
  double timeA = GetTime();
  switch (TOT) {
  case StartTimeDo:
    S->time1[TT] = timeA;
    break;
  case StopTimeDo:
    S->time2[TT] = timeA;
    S->time[TT] += (S->time2[TT]-S->time1[TT]);
    break;
  } 
}  

/** Final Measuring and overall time.
 * Does the last measurement calculation and combines results from all processes
 * to calculate min, max, average and writes it to FileData::SpeedFile
 * \param *P Problem at hand
 */
void CompleteSpeedMeasure(struct Problem *P)
{ /* gibt letzte Messung aus und benoetigte Gesamtzeit */
  struct SpeedStruct *S = &P->Speed;
  struct FileData *F = &P->Files;
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  struct RunStruct *R = &P->R;
  struct Ions *I = &P->Ion;
  double average[MAXTIMETYPES];
  double stddev[MAXTIMETYPES];
  double time[MAXTIMETYPES];
  int i;
  S->InitSteps += S->LevUpSteps;
  for (i=0; i<MAXTIMETYPES; i++) 
    switch ((enum TimeTypes)i) {
    case InitSimTime:
    case InitGramSchTime:
    case InitLocTime:
    case InitNonLocTime:
    case InitDensityTime:
    case LocFTime:
    case NonLocFTime:
    case EwaldTime:
    case GapTime:
      S->SpeedStep[i] = S->InitSteps;
      break;
    case SimTime:
    case GramSchTime:
    case LocTime:
    case NonLocTime:
    case DensityTime:
    case WannierTime:
    case ReadnWriteTime:
      S->SpeedStep[i] = S->Steps*P->Par.proc[0]/(double)P->Lat.Psi.GlobalNo[PsiMaxNo];
      break;
    default: 
      S->SpeedStep[i] = 1;
    }
  S->time[LevSMaxG] = R->LevS->MaxG;
  MPI_Allreduce( S->time, time, MAXTIMETYPES, MPI_DOUBLE, MPI_SUM, P->Par.comm);
  for (i=0; i<MAXTIMETYPES; i++) 
    time[i] /= (double)P->Par.procs;
  time[LevSMaxG] = R->LevS->TotalAllMaxG;
  for (i=0; i<MAXTIMETYPES; i++) 
    average[i] = S->time[i]/S->SpeedStep[i];
  MPI_Allreduce( average, S->average, MAXTIMETYPES, MPI_DOUBLE, MPI_SUM, P->Par.comm);
  for (i=0; i<MAXTIMETYPES; i++) 
    S->average[i] /= (double)P->Par.procs;
  S->average[LevSMaxG] = R->LevS->TotalAllMaxG/(double)P->Par.proc[1];
  for (i=0; i<MAXTIMETYPES; i++) 
    stddev[i] = (average[i]-S->average[i])*(average[i]-S->average[i]);
  MPI_Allreduce( stddev, S->stddev, MAXTIMETYPES, MPI_DOUBLE, MPI_SUM, P->Par.comm);
  MPI_Allreduce( &stddev[LevSMaxG], &S->stddev[LevSMaxG], 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
  for (i=0; i<MAXTIMETYPES-1; i++) 
    S->stddev[i] /= (double)P->Par.procs;
  S->stddev[LevSMaxG] /= (double)P->Par.proc[1];
  for (i=0; i<MAXTIMETYPES; i++) 
    S->stddev[i] = sqrt(S->stddev[i]);
  
  MPI_Allreduce( average, S->max, MAXTIMETYPES, MPI_DOUBLE, MPI_MAX, P->Par.comm);
  MPI_Allreduce( average, S->min, MAXTIMETYPES, MPI_DOUBLE, MPI_MIN, P->Par.comm);
  if (P->Par.me != 0) return;
  
  fprintf(F->SpeedFile, "LevNo\tMaxG\tRMaxG\tMaxN\tN[0]\tN[1]\tN[2]\tRLevSStep\n");
  for (i=0; i < Lat->MaxLevel; i++)
    fprintf(F->SpeedFile, "%i\t%i\t%i\t%i\t%i\t%i\t%i\t%i\n", Lat->Lev[i].LevelNo, Lat->Lev[i].TotalAllMaxG, Lat->Lev[i].TotalRealAllMaxG, Lat->Lev[i].MaxN, Lat->Lev[i].N[0], Lat->Lev[i].N[1], Lat->Lev[i].N[2], Lat->Lev[i].Step);
  fprintf(F->SpeedFile, "procs\tproc[0]\tproc[1]\tPsiMax\tPsiDoub\tPsiUp\tPsiDown\tNRStep\tTotalIons\n");
  fprintf(F->SpeedFile, "%i\t%i\t%i\t%i\t%i+%i\t%i+%i\t%i+%i\t%i\t%i\n", P->Par.procs, P->Par.proc[0], P->Par.proc[1], Psi->GlobalNo[PsiMaxNo], Psi->GlobalNo[PsiMaxNoDouble], Psi->GlobalNo[PsiMaxAdd], Psi->GlobalNo[PsiMaxNoUp], Psi->GlobalNo[PsiMaxAdd], Psi->GlobalNo[PsiMaxNoDown], Psi->GlobalNo[PsiMaxAdd], R->NewRStep, I->Max_TotalIons);
  
  fprintf(F->SpeedFile, "Type\tSimTime\t\tInitSimTime\tInitTime\tInitGramSchTime\tInitLocTime\tInitNonLocTime\tInitDensTime\tGramSchTime\tLocTime\t\tNonLocTime\tDensTime\tLocFTime\tNonLocFTime\tEwaldTime\tGapTime\tCurrDensTime\tWannierTime\tReadnWriteTime\tLevSMaxG\n");
  fprintf(F->SpeedFile, "Steps");
  for(i=0;i<MAXTIMETYPES;i++)
    fprintf(F->SpeedFile, "\t%e", S->SpeedStep[i]);
  fprintf(F->SpeedFile, "\n");

  fprintf(F->SpeedFile, "total");
  for(i=0;i<MAXTIMETYPES;i++)
    fprintf(F->SpeedFile, "\t%e", time[i]);
  fprintf(F->SpeedFile, "\n");
  
  fprintf(F->SpeedFile, "average");
  for(i=0;i<MAXTIMETYPES;i++)
    fprintf(F->SpeedFile, "\t%e", S->average[i]);
  fprintf(F->SpeedFile, "\n");

  fprintf(F->SpeedFile, "stddev");
  for(i=0;i<MAXTIMETYPES;i++)
    fprintf(F->SpeedFile, "\t%e", S->stddev[i]);
  fprintf(F->SpeedFile, "\n");

  fprintf(F->SpeedFile, "min");
  for(i=0;i<MAXTIMETYPES;i++)
    fprintf(F->SpeedFile, "\t%e", S->min[i]);
  fprintf(F->SpeedFile, "\n");

  fprintf(F->SpeedFile, "max");
  for(i=0;i<MAXTIMETYPES;i++)
    fprintf(F->SpeedFile, "\t%e", S->max[i]);
  fprintf(F->SpeedFile, "\n");

  fclose(F->SpeedFile);
}