Skip to content
Snippets Groups Projects
Select Git revision
  • a88b277b3112ee5c052a880585ac811aad1c717c
  • master default protected
  • OmpSs-FreeParametersWithFunctions
  • OmpSs-WithFunctions
  • OmpSs-FPGA
  • OmpSs-GPU
  • OmpSs
7 results

memRequest.cpp

Blame
  • clffft.cpp 15.19 KiB
    /* Copyright 2004 INAF - Osservatorio Astronomico di Cagliari
    
       Licensed under the Apache License, Version 2.0 (the "License");
       you may not use this file except in compliance with the License.
       You may obtain a copy of the License at
    
           http://www.apache.org/licenses/LICENSE-2.0
    
       Unless required by applicable law or agreed to in writing, software
       distributed under the License is distributed on an "AS IS" BASIS,
       WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
       See the License for the specific language governing permissions and
       limitations under the License.
     */
    
    /*! \file clffft.cpp
     *
     * \brief C++ implementation of LFFFT functions.
     */
    #include <cstdio>
    #include <exception>
    #include <fstream>
    #include <regex>
    #include <string>
    
    #ifndef INCLUDE_TYPES_H_
    #include "../include/types.h"
    #endif
    
    #ifndef INCLUDE_PARSERS_H_
    #include "../include/Parsers.h"
    #endif
    
    #ifndef INCLUDE_LOGGING_H_
    #include "../include/logging.h"
    #endif
    
    #ifndef INCLUDE_CONFIGURATION_H_
    #include "../include/Configuration.h"
    #endif
    
    #ifndef INCLUDE_COMMONS_H_
    #include "../include/Commons.h"
    #endif
    
    #ifndef INCLUDE_SPH_SUBS_H_
    #include "../include/sph_subs.h"
    #endif
    
    #ifndef INCLUDE_TFRFME_H_
    #include "../include/tfrfme.h"
    #endif
    
    #ifndef INCLUDE_TRA_SUBS_H_
    #include "../include/tra_subs.h"
    #endif
    
    using namespace std;
    
    /*! \brief C++ implementation of LFFFT
     *
     *  \param data_file: `string` Name of the input data file.
     *  \param output_path: `string` Directory to write the output files in.
     */
    void lffft(string data_file, string output_path) {
      const dcomplex uim = 0.0 + 1.0 * I;
      const double sq2i = 1.0 / sqrt(2.0);
      const dcomplex sq2iti = sq2i * uim;
    
      fstream tlfff, tlfft;
      double ****zpv = NULL;
      dcomplex *ac = NULL, *ws = NULL, *wsl = NULL;
      dcomplex **am0m = NULL;
      dcomplex **amd = NULL;
      int **indam = NULL;
      dcomplex *tmsm = NULL, *tmse = NULL, **tms = NULL;
      dcomplex **_wsum = NULL;
      int jft, jss, jtw;
      int is, le, nvam = 0;
      double vks, exris;
      CIL *cil = new CIL();
      CCR *ccr = new CCR();
    
      int num_lines = 0;
      string *file_lines = load_file(data_file, &num_lines);
      regex re = regex("-?[0-9]+");
      smatch m;
      string str_target = file_lines[0];
      for (int mi = 0; mi < 3; mi++) {
        regex_search(str_target, m, re);
        if (mi == 0) jft = stoi(m.str());
        else if (mi == 1) jss = stoi(m.str());
        else if (mi == 2) jtw = stoi(m.str());
        str_target = m.suffix().str();
      } // mi loop
      string ttms_name = output_path + "/c_TTMS";
      fstream ttms;
      ttms.open(ttms_name, ios::in | ios::binary);
      if (ttms.is_open()) {
        ttms.read(reinterpret_cast<char *>(&is), sizeof(int));
        ttms.read(reinterpret_cast<char *>(&le), sizeof(int));
        ttms.read(reinterpret_cast<char *>(&vks), sizeof(double));
        ttms.read(reinterpret_cast<char *>(&exris), sizeof(double));
        cil->le = le;
        cil->nlem = le * (le + 2);
        cil->nlemt = cil->nlem + cil->nlem;
        if (is >= 2222) { // label 120
          tms = new dcomplex*[le];
          for (int ti = 0; ti < le; ti++) tms[ti] = new dcomplex[3]();
          // QUESTION|WARNING: original code uses LM without defining it. Where does it come from?
          int lm = le;
          for (int i = 0; i < lm; i++) {
    	double vreal, vimag;
    	ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
    	ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
    	tms[i][0] = vreal + vimag * I;
    	ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
    	ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
    	tms[i][1] = vreal + vimag * I;
    	ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
    	ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
    	tms[i][2] = vreal + vimag * I;
          } // i loop
        } else if (is >= 1111) { // label 125
          tmsm = new dcomplex[le]();
          tmse = new dcomplex[le]();
          for (int i = 0; i < le; i++) {
    	double vreal, vimag;
    	ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
    	ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
    	tmsm[i] = vreal + vimag * I;
    	ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
    	ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
    	tmse[i] = vreal + vimag * I;
          } // i loop
        } else if (is >= 0) { // label 135
          am0m = new dcomplex*[cil->nlemt];
          for (int ai = 0; ai < cil->nlemt; ai++) am0m[ai] = new dcomplex[cil->nlemt]();
          for (int i = 0; i < cil->nlemt; i++) {
    	for (int j = 0; j < cil->nlemt; j++) {
    	  double vreal, vimag;
    	  ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
    	  ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
    	  am0m[i][j] = vreal + vimag * I;
    	} // j loop
          } // i loop
        } else if (is < 0) {
          nvam = le * le + (le * (le + 1) * (le * 2 + 1)) / 3;
          amd = new dcomplex*[nvam];
          for (int ai = 0; ai < nvam; ai++) amd[ai] = new dcomplex[4]();
          for (int i = 0; i < nvam; i++) {
    	for (int j = 0; j < 4; j++) {
    	  double vreal, vimag;
    	  ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
    	  ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
    	  amd[i][j] = vreal + vimag * I;
    	} // j loop
          } // i loop
          indam = new int*[le];
          int vint;
          for (int ii = 0; ii < le; ii++) indam[ii] = new int[le]();
          for (int i = 0; i < le; i++) {
    	for (int j = 0; j < le; j++) {
    	  ttms.read(reinterpret_cast<char *>(&vint), sizeof(int));
    	  indam[i][j] = vint;
    	} // j loop
          } // i loop
          ttms.read(reinterpret_cast<char *>(&vint), sizeof(int));
          cil->mxmpo = vint;
          cil->mxim = vint * 2 - 1;
        }
        // label 150
        ttms.close();
        TFRFME *tfrfme = NULL;
        string binary_name;
        if (jss != 1) binary_name = output_path + "/c_TFRFME.hd5";
        else binary_name = output_path + "/c_TWS";
        tfrfme = TFRFME::from_binary(binary_name, "HDF5");
        if (tfrfme != NULL) {
          int lmode, lm, nkv, nxv, nyv, nzv;
          double vk, exri, an, ff, tra;
          double spd, frsh, exril;
          lmode = (int)tfrfme->get_param("lmode");
          lm = (int)tfrfme->get_param("lm");
          nkv = (int)tfrfme->get_param("nkv");
          nxv = (int)tfrfme->get_param("nxv");
          nyv = (int)tfrfme->get_param("nyv");
          nzv = (int)tfrfme->get_param("nzv");
          _wsum = tfrfme->get_matrix();
          double *_xv = tfrfme->get_x();
          double *_yv = tfrfme->get_y();
          double *_zv = tfrfme->get_z();
          if (lm >= le) {
    	vk = tfrfme->get_param("vk");
    	exri = tfrfme->get_param("exri");
    	an = tfrfme->get_param("an");
    	ff = tfrfme->get_param("ff");
    	tra = tfrfme->get_param("tra");
    	if (vk == vks && exri == exris) {
    	  spd = tfrfme->get_param("spd");
    	  frsh = tfrfme->get_param("frsh");
    	  exril = tfrfme->get_param("exril");
    	  bool goto160 = false;
    	  if (jft <= 0) {
    	    zpv = new double***[le];
    	    for (int zi = 0; zi < le; zi++) {
    	      zpv[zi] = new double**[3];
    	      for (int zj = 0; zj < 3; zj++) {
    		zpv[zi][zj] = new double*[2];
    		for (int zk = 0; zk < 2; zk++) zpv[zi][zj][zk] = new double[2]();
    	      } // zj loop
    	    } // zi loop
    	    thdps(le, zpv);
    	    double exdc = exri * exri;
    	    double sqk = vk * vk * exdc;
    	    ccr->cof = 1.0 / sqk;
    	    ccr->cimu = ccr->cof / sqrt(2.0);
    	    if (jss != 1) {
    	      string tlfff_name = output_path + "/c_TLFFF";
    	      tlfff.open(tlfff_name.c_str(), ios::out | ios::binary);
    	      if (tlfff.is_open()) {
    		tlfff.write(reinterpret_cast<char *>(&lmode), sizeof(int));
    		tlfff.write(reinterpret_cast<char *>(&le), sizeof(int));
    		tlfff.write(reinterpret_cast<char *>(&nkv), sizeof(int));
    		tlfff.write(reinterpret_cast<char *>(&nxv), sizeof(int));
    		tlfff.write(reinterpret_cast<char *>(&nyv), sizeof(int));
    		tlfff.write(reinterpret_cast<char *>(&nzv), sizeof(int));
    		tlfff.write(reinterpret_cast<char *>(&vk), sizeof(double));
    		tlfff.write(reinterpret_cast<char *>(&exri), sizeof(double));
    		tlfff.write(reinterpret_cast<char *>(&an), sizeof(double));
    		tlfff.write(reinterpret_cast<char *>(&ff), sizeof(double));
    		tlfff.write(reinterpret_cast<char *>(&tra), sizeof(double));
    		tlfff.write(reinterpret_cast<char *>(&spd), sizeof(double));
    		tlfff.write(reinterpret_cast<char *>(&frsh), sizeof(double));
    		tlfff.write(reinterpret_cast<char *>(&exril), sizeof(double));
    		for (int i = 0; i < nxv; i++) {
    		  double x = _xv[i];
    		  tlfff.write(reinterpret_cast<char *>(&x), sizeof(double));
    		}
    		for (int i = 0; i < nyv; i++) {
    		  double y = _yv[i];
    		  tlfff.write(reinterpret_cast<char *>(&y), sizeof(double));
    		}
    		for (int i = 0; i < nzv; i++) {
    		  double z = _zv[i];
    		  tlfff.write(reinterpret_cast<char *>(&z), sizeof(double));
    		}
    		if (jft < 0) goto160 = true;
    	      } else { // Should never happen.
    		printf("ERROR: could not open TLFFF file.\n");
    	      }
    	    }
    	  }
    	  // label 155
    	  if (!goto160) {
    	    if (jss != 1) {
    	      string tlfft_name = output_path + "/c_TLFFT";
    	      tlfft.open(tlfft_name.c_str(), ios::out | ios::binary);
    	      if (tlfft.is_open()) {
    		tlfft.write(reinterpret_cast<char *>(&lmode), sizeof(int));
    		tlfft.write(reinterpret_cast<char *>(&le), sizeof(int));
    		tlfft.write(reinterpret_cast<char *>(&nkv), sizeof(int));
    		tlfft.write(reinterpret_cast<char *>(&nxv), sizeof(int));
    		tlfft.write(reinterpret_cast<char *>(&nyv), sizeof(int));
    		tlfft.write(reinterpret_cast<char *>(&nzv), sizeof(int));
    		tlfft.write(reinterpret_cast<char *>(&vk), sizeof(double));
    		tlfft.write(reinterpret_cast<char *>(&exri), sizeof(double));
    		tlfft.write(reinterpret_cast<char *>(&an), sizeof(double));
    		tlfft.write(reinterpret_cast<char *>(&ff), sizeof(double));
    		tlfft.write(reinterpret_cast<char *>(&tra), sizeof(double));
    		tlfft.write(reinterpret_cast<char *>(&spd), sizeof(double));
    		tlfft.write(reinterpret_cast<char *>(&frsh), sizeof(double));
    		tlfft.write(reinterpret_cast<char *>(&exril), sizeof(double));
    		for (int i = 0; i < nxv; i++) {
    		  double x = _xv[i];
    		  tlfft.write(reinterpret_cast<char *>(&x), sizeof(double));
    		}
    		for (int i = 0; i < nyv; i++) {
    		  double y = _yv[i];
    		  tlfft.write(reinterpret_cast<char *>(&y), sizeof(double));
    		}
    		for (int i = 0; i < nzv; i++) {
    		  double z = _zv[i];
    		  tlfft.write(reinterpret_cast<char *>(&z), sizeof(double));
    		}
    	      } else { // Should never happen.
    		printf("ERROR: could not open TLFFT file.\n");
    	      }
    	    }
    	  }
    	  // label 160
    	  const int nlmm = lm * (lm + 2);
    	  const int nlmmt = nlmm + nlmm;
    	  ws = new dcomplex[nlmmt]();
    	  if (lm > le) wsl = new dcomplex[nlmmt]();
    	  // FORTRAN writes two output formatted files without opening them
    	  // explicitly. It is assumed thay can be opened here.
    	  string out66_name = output_path + "/c_out66.txt";
    	  string out67_name = output_path + "/c_out67.txt";
    	  FILE *output66 = fopen(out66_name.c_str(), "w");
    	  FILE *output67 = fopen(out67_name.c_str(), "w");
    	  for (int iz475 = 0; iz475 < nzv; iz475++) {
    	    for (int iy475 = 0; iy475 < nyv; iy475++) {
    	      for (int ix475 = 0; ix475 < nxv; ix475++) {
    		for (int i = 0; i < nlmmt; i++) {
    		  //double vreal, vimag;
    		  //binary_input.read(reinterpret_cast<char *>(&vreal), sizeof(double));
    		  //binary_input.read(reinterpret_cast<char *>(&vimag), sizeof(double));
    		  int row = i;
    		  int col = (nyv * nxv * iz475) + (nxv * iy475) + ix475;
    		  dcomplex value = _wsum[row][col];
    		  if (lm <= le) {
    		    ws[i] = value;
    		  } else { // label 170
    		    wsl[i] = value;
    		    for (int i175 = 0; i175 < cil->nlem; i175++) {
    		      int ie = i175 + cil->nlem;
    		      int iel = i175 + nlmm;
    		      ws[i175] = wsl[i175];
    		      ws[ie] = wsl[iel];
    		    } // i175 loop
    		  }
    		} // i loop
    		// label 180
    		if (is != 2222) {
    		  if (is != 1111) {
    		    if (is > 0) { // Goes to 305
    		      ac = new dcomplex[cil->nlemt]();
    		      camp(ac, am0m, ws, cil);
    		      // Goes to 445
    		    } else if (is < 0) { // Goes to 405
    		      ac = new dcomplex[cil->nlemt]();
    		      czamp(ac, amd, indam, ws, cil);
    		      // Goes to 445
    		    }
    		  } else {
    		    ac = new dcomplex[cil->nlemt]();
    		    samp(ac, tmsm, tmse, ws, cil);
    		    // Goes to 445
    		  }
    		} else {
    		  ac = new dcomplex[cil->nlemt]();
    		  sampoa(ac, tms, ws, cil);
    		  // Goes to 445
    		}
    		bool goto475 = false;
    		// label 445
    		if (jft <= 0) {
    		  double *fffe = new double[3]();
    		  double *fffs = new double[3]();
    		  ffrf(zpv, ac, ws, fffe, fffs, cil, ccr);
    		  if (jss == 1) {
    		    // Writes to 66
    		    fprintf(
    			    output66, " %18.16lE%18.16lE%18.16lE\n",
    			    fffe[0], fffs[0], fffe[0] - fffs[0]
    			    );
    		    fprintf(
    			    output66, " %18.16lE%18.16lE%18.16lE\n",
    			    fffe[1], fffs[1], fffe[1] - fffs[1]
    			    );
    		    fprintf(
    			    output66, " %18.16lE%18.16lE%18.16lE\n",
    			    fffe[2], fffs[2], fffe[2] - fffs[2]
    			    );
    		  } else { // label 450
    		    for (int i = 0; i < 3; i++) {
    		      double value = fffe[i] - fffs[i];
    		      tlfff.write(reinterpret_cast<char *>(&value), sizeof(double));
    		    }
    		    if (jtw == 1) {
    		      // Writes to 66
    		      fprintf(
    			      output66, " %5d%4d%4d%15.4lE%15.4lE%15.4lE\n",
    			      ix475 + 1, iy475 + 1, iz475 + 1,
    			      fffe[0] - fffs[0], fffe[1] - fffs[1], fffe[2] - fffs[2]
    			      );
    		    }
    		  }
    		  if (jft < 0) goto475 = true;
    		  delete[] fffe;
    		  delete[] fffs;
    		}
    		// label 460
    		if (!goto475) {
    		  double *ffte = new double[3]();
    		  double *ffts = new double[3]();
    		  ffrt(ac, ws, ffte, ffts, cil);
    		  if (jss == 1) {
    		    // Writes to 67
    		    fprintf(
    			    output67, " %18.16lE%18.16lE%18.16lE\n",
    			    ffte[0], ffts[0], ffte[0] - ffts[0]
    			    );
    		    fprintf(
    			    output67, " %18.16lE%18.16lE%18.16lE\n",
    			    ffte[1], ffts[1], ffte[1] - ffts[1]
    			    );
    		    fprintf(
    			    output67, " %18.16lE%18.16lE%18.16lE\n",
    			    ffte[2], ffts[2], ffte[2] - ffts[2]
    			    );
    		  } else { // label 470
    		    for (int i = 0; i < 3; i++) {
    		      double value = ffte[i] - ffts[i];
    		      tlfft.write(reinterpret_cast<char *>(&value), sizeof(double));
    		    }
    		    if (jtw == 1) {
    		      // Writes to 67
    		      fprintf(
    			      output67, " %5d%4d%4d%15.4lE%15.4lE%15.4lE\n",
    			      ix475 + 1, iy475 + 1, iz475 + 1,
    			      ffte[0] - ffts[0], ffte[1] - ffts[1], ffte[2] - ffts[2]
    			      );
    		    }
    		  }
    		  delete[] ffte;
    		  delete[] ffts;
    		}
    	      } // ix475 loop
    	    } // iy475 loop
    	  } // iz475 loop
    	  if (jss != 1) {
    	    if (jft <= 0) tlfff.close();
    	    if (jft >= 0) tlfft.close();
    	  }
    	  fclose(output66);
    	  fclose(output67);
    	}
          }
          delete tfrfme;
        } else {
          printf("ERROR: could not open binary input file %s.\n", binary_name.c_str());
        }
      } else {
        printf("ERROR: could not open TTMS file.\n");
      }
      
      // Clean up memory
      if (ac != NULL) delete[] ac;
      if (ws != NULL) delete[] ws;
      if (wsl != NULL) delete[] wsl;
      if (tmsm != NULL) delete[] tmsm;
      if (tmse != NULL) delete[] tmse;
      if (tms != NULL) {
        for (int ti = le - 1; ti > -1; ti--) delete[] tms[ti];
        delete[] tms;
      }
      if (am0m != NULL) {
        for (int ai = cil->nlemt - 1; ai > -1; ai--) delete[] am0m[ai];
        delete[] am0m;
      }
      if (amd != NULL) {
        for (int ai = nvam - 1; ai > -1; ai--) delete[] amd[ai];
        delete[] amd;
      }
      if (indam != NULL) {
        for (int ii = le - 1; ii > -1; ii--) delete[] indam[ii];
        delete[] indam;
      }
      if (zpv != NULL) {
        for (int zi = le - 1; zi > -1; zi--) {
          for (int zj = 2; zj > -1; zj--) {
    	for (int zk = 1; zk > -1; zk--) delete[] zpv[zi][zj][zk];
    	delete[] zpv[zi][zj];
          } // zj loop
          delete[] zpv[zi];
        } // zi loop
        delete[] zpv;
      }
      delete cil;
      delete ccr;
      delete[] file_lines;
      printf("LFFT: Done.\n");
    }