text:
/*****************************************************************\
__
/ /
/ /                     __  __
/ /______    _______    / / / / ________   __       __
/ ______  \  /_____  \  / / / / / _____  | / /      / /
/ /      | / _______| / / / / / / /____/ / / /      / /
/ /      / / / _____  / / / / / / _______/ / /      / /
/ /      / / / /____/ / / / / / / |______  / |______/ /
/_/      /_/ |________/ / / / /  \_______/  \_______  /
/_/ /_/                     / /
/ /
High Level Game Framework            /_/

---------------------------------------------------------------

Copyright (c) 2007-2011 - Rodrigo Braz Monteiro.
This file is subject to the terms of halley_license.txt.

\*****************************************************************/

#pragma once

#include "halley/text/halleystring.h"
#include "halley/file/path.h"
#include <memory>
#include <functional>
#include <halley/concurrency/future.h>
#include <gsl/gsl>
#include "metadata.h"

namespace Halley {
	enum class AssetType;

	class ResourceDataReader {
	public:
		virtual ~ResourceDataReader() {}
		virtual size_t size() const = 0;
		virtual int read(gsl::span<gsl::byte> dst) = 0;
		virtual void seek(int64_t pos, int whence) = 0;
		virtual size_t tell() const = 0;
		virtual void close() = 0;

		Bytes readAll();
	};

	class ResourceData {
	public:
		ResourceData(String path);
		virtual ~ResourceData() {}

		String getPath() const { return path; }

	private:
		String path;
	};

	class ResourceDataStatic : public ResourceData {
	public:
		ResourceDataStatic(String path);
		ResourceDataStatic(const void* data, size_t size, String path, bool owning = true);

		void set(const void* data, size_t size, bool owning = true);
		bool isLoaded() const;

		const void* getData() const;
		gsl::span<const gsl::byte> getSpan() const;
		size_t getSize() const;
		String getString() const;
		void inflate();

		static std::unique_ptr<ResourceDataStatic> loadFromFileSystem(Path path);
		void writeToFileSystem(String path) const;

	private:
		std::shared_ptr<const char> data;
		size_t size = 0;
		bool loaded;
	};

	typedef std::function<std::unique_ptr<ResourceDataReader>()> ResourceDataMakeReader;
	class ResourceDataStream : public ResourceData {
	public:
		ResourceDataStream(String path, ResourceDataMakeReader makeReader);
		std::unique_ptr<ResourceDataReader> getReader() const { return make(); }

	private:
		ResourceDataMakeReader make;
	};

	class IResourceLocator
	{
	public:
		virtual ~IResourceLocator() {}
		virtual const Metadata& getMetaData(const String& resource, AssetType type) const = 0;
		virtual std::unique_ptr<ResourceDataStatic> getStatic(const String& asset, AssetType type) = 0;
		virtual std::unique_ptr<ResourceDataStream> getStream(const String& asset, AssetType type) = 0;
	};

	enum class ResourceLoadPriority {
		Low = 0,
		Normal = 1,
		High = 2
	};

	class HalleyAPI;
	class Metadata;

	class ResourceLoader
	{
		friend class ResourceCollectionBase;

	public:
		const String& getName() const { return name; }
		ResourceLoadPriority getPriority() const { return priority; }
		const HalleyAPI& getAPI() const { return *api; }
		const Metadata& getMeta() const { return *metadata; }

		std::unique_ptr<ResourceDataStatic> getStatic();
		std::unique_ptr<ResourceDataStream> getStream();
		Future<std::unique_ptr<ResourceDataStatic>> getAsync() const;

	private:
		ResourceLoader(ResourceLoader&& loader) noexcept;
		ResourceLoader(IResourceLocator& locator, const String& name, AssetType type, ResourceLoadPriority priority, const HalleyAPI* api);
		~ResourceLoader();

		IResourceLocator& locator;
		String name;
		AssetType type;
		ResourceLoadPriority priority;
		const HalleyAPI* api;
		const Metadata* metadata;
		bool loaded = false;
	};

}

max_stars_repo_name: lye/halley

max_issues_repo_name: lye/halley

max_forks_repo_name: lye/halley


text:
#ifndef DATA_PROCESS_H
#define DATA_PROCESS_H

#include <stdio.h>
#include <string.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_statistics.h>
#include <math.h>

void support_degree_generator(gsl_matrix*, double*);
void eigenvec_calc(gsl_vector*, gsl_matrix*, gsl_matrix*);
void principal_comp_calc(gsl_matrix*, gsl_matrix*, gsl_matrix*);
void contri_rate_calc_kth(gsl_vector*, double*);
void major_contri_calc(double*, double*);
void integ_supp_score_calc(double*, gsl_matrix*, gsl_vector*);
void elliminate_incorrect_data(gsl_vector*, int*);
void weight_coeff_calc(gsl_vector*, int*, double*);
double fused_output(double*, double*);

#endif

ext: h

max_stars_repo_name: karanbirsandhu/nu-sense

max_issues_repo_name: karanbirsandhu/nu-sense

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include "airsar.h"
#include "asf_meta.h"
#include <ctype.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>

#define SQR(x) (x*x)

static char *get_airsar(char *buf, char *str)
{
  char *p, *q;

  static char value[51];
  memset(value,0,51);

  q = (char *) CALLOC(51,sizeof(char));
  p = strstr(buf, str);
  if (p) {
    strncpy(q, p, 50);
    strcpy(value, q+strlen(str));
    //printf("%s: %s\n", str, value);
  }
  else
    strcpy(value, "");

  FREE(q);

  return value;
}

airsar_header *read_airsar_header(const char *dataFile)
{
  airsar_header *header = NULL;
  FILE *fp;
  char buf[4400], *value;
  double version;
printf("\n\ndataFile: %s\n\n", dataFile);
  // Allocate memory and file handling
  value = (char *) MALLOC(sizeof(char)*25);
  header = (airsar_header *) CALLOC(1, sizeof(airsar_header));
  fp = FOPEN(dataFile, "r");
  if (fgets(buf, 4400, fp) == NULL)
    asfPrintError("Could not read general header\n");
  FCLOSE(fp);

  // Check for the processor version
  // We take this as an indicator that we actually deal with AirSAR data
  version = atof(get_airsar(buf, "JPL AIRCRAFT SAR PROCESSOR VERSION"));

  // Read general header
  if (version > 0.0) {
    header->record_length = atoi(get_airsar(buf, "RECORD LENGTH IN BYTES ="));
    header->number_records = atoi(get_airsar(buf, "NUMBER OF HEADER RECORDS ="));
    header->sample_count =
      atoi(get_airsar(buf, "NUMBER OF SAMPLES PER RECORD ="));
    header->line_count =
      atoi(get_airsar(buf, "NUMBER OF LINES IN IMAGE ="));
    sprintf(header->processor, "%s",
      trim_spaces(get_airsar(buf, "JPL AIRCRAFT SAR PROCESSOR VERSION")));
    value = trim_spaces(get_airsar(buf, "DATA TYPE ="));
    if (strcmp(value, "INTEGER*2") == 0)
      header->data_type = INTEGER16;
    sprintf(header->range_projection, "%s",
	    trim_spaces(get_airsar(buf, "RANGE PROJECTION =")));
    header->x_pixel_size =
      atof(get_airsar(buf, "RANGE PIXEL SPACING (METERS) ="));
    header->y_pixel_size =
      atof(get_airsar(buf, "AZIMUTH PIXEL SPACING (METERS) ="));
    header->first_data_offset =
      atoi(get_airsar(buf, "BYTE OFFSET OF FIRST DATA RECORD ="));
    if (header->first_data_offset == 0)
      header->first_data_offset = header->record_length*header->number_records;
    header->parameter_header_offset =
      atoi(get_airsar(buf, "BYTE OFFSET OF PARAMETER HEADER ="));
    header->calibration_header_offset =
      atoi(get_airsar(buf, "BYTE OFFSET OF CALIBRATION HEADER ="));
    header->dem_header_offset =
      atoi(get_airsar(buf, "BYTE OFFSET OF DEM HEADER ="));
  }

  return header;
}

airsar_param_header *read_airsar_params(const char *dataFile)
{
  airsar_param_header *params=NULL;
  airsar_header *header = read_airsar_header(dataFile);
  FILE *fp;
  char *buf;
  int size;

  // Allocate memory and file handling
  params = (airsar_param_header *) CALLOC(1, sizeof(airsar_param_header));
  if (header->calibration_header_offset > 0)
    size = header->calibration_header_offset - header->parameter_header_offset;
  else if (header->dem_header_offset > 0)
    size = header->dem_header_offset - header->parameter_header_offset;
  else
    size = header->first_data_offset - header->parameter_header_offset;
  buf = (char *) MALLOC(sizeof(char)*size);
  fp = FOPEN(dataFile, "r");
  FSEEK(fp, header->parameter_header_offset, 1);
  if (fgets(buf, size, fp) == NULL)
    asfPrintError("Could not read parameter header\n");
  FCLOSE(fp);

  // Read parameter header
  sprintf(params->site_name, "%s",
	  trim_spaces(get_airsar(buf, "SITE NAME")));
  sprintf(params->cct_type, "%s",
	  trim_spaces(get_airsar(buf, "CCT TYPE")));
  params->cct_id = atoi(get_airsar(buf, "CCT ID"));
  params->start_lat =
    atof(get_airsar(buf,
		    "LATITUDE AT START OF SCENE (DEGREES)"));
  params->start_lon =
    atof(get_airsar(buf,
		    "LONGITUDE AT START OF SCENE (DEGREES)"));
  params->end_lat =
    atof(get_airsar(buf,
		    "LATITUDE AT END OF SCENE (DEGREES)"));
  params->end_lon =
    atof(get_airsar(buf,
		    "LONGITUDE AT END OF SCENE (DEGREES)"));
  sprintf(params->acquisition_date, "%s",
	  trim_spaces(get_airsar(buf,
				 "DATE OF ACQUISITION (GMT)")));
  params->acquisition_seconds =
    atof(get_airsar(buf,
		    "TIME OF ACQUISITION: SECONDS IN DAY"));
  sprintf(params->frequencies, "%s",
	  trim_spaces(get_airsar(buf,
				 "FREQUENCIES COLLECTED")));
  params->prf =
    atof(get_airsar(buf, "PRF AT START OF TRANSFER (HZ)"));
  params->range_sampling_rate =
    atof(get_airsar(buf, "SAMPLING RATE (MHZ)"));
  params->chirp_bandwidth =
    atof(get_airsar(buf, "CHIRP BANDWIDTH (MHZ)"));
  params->pulse_length =
    atof(get_airsar(buf, "PULSE LENGTH (MICROSECONDS)"));
  params->wavelength =
    atof(get_airsar(buf, "PROCESSOR WAVELENGTH (METERS)"));
  params->near_slant_range =
    atof(get_airsar(buf, "NEAR SLANT RANGE (METERS)"));
  params->far_slant_range =
    atof(get_airsar(buf, "FAR SLANT RANGE (METERS)"));
  params->near_look_angle =
    atof(get_airsar(buf, "NEAR LOOK ANGLE (DEGREES)"));
  params->far_look_angle =
    atof(get_airsar(buf, "FAR LOOK ANGLE (DEGEES)"));
  params->azimuth_look_count =
    atoi(get_airsar(buf,
		    "NUMBER OF LOOKS PROCESSED IN AZIMUTH"));
  params->range_look_count =
    atoi(get_airsar(buf,
		    "NUMBER OF LOOKS PROCESSED IN RANGE"));
  params->deskewed =
    atoi(get_airsar(buf,
		    "DESKEW FLAG (1=DESKEWED, 2=NOT DESKEWED)"));
  params->sr_sample_spacing =
    atof(get_airsar(buf,
		    "SLANT RANGE SAMPLE SPACING (METERS)"));
  params->slant_range_resolution =
    atof(get_airsar(buf,
		    "NOMINAL SLANT RANGE RESOLUTION (METERS)"));
  params->azimuth_sample_spacing =
    atof(get_airsar(buf, "AZIMUTH SAMPLE SPACING (METERS)"));
  params->azimuth_resolution =
    atof(get_airsar(buf,
		    "NOMINAL AZIMUTH RESOLUTION (METERS)"));
  params->center_lat =    atof(get_airsar(buf, "IMAGE CENTER LATITUDE (DEGREES)"));
  params->center_lon =
    atof(get_airsar(buf,
		    "IMAGE CENTER LONGITUDE (DEGREES)"));
  params->scale_factor =
    atof(get_airsar(buf, "GENERAL SCALE FACTOR"));
  params->cal_factor_hh =
    atof(get_airsar(buf, "CALIBRATION FACTOR APPLIED, DB, HH"));
  params->cal_factor_hv =
    atof(get_airsar(buf, "CALIBRATION FACTOR APPLIED, DB, HV"));
  params->cal_factor_vh =
    atof(get_airsar(buf, "CALIBRATION FACTOR APPLIED, DB, VH"));
  params->cal_factor_vv =
    atof(get_airsar(buf, "CALIBRATION FACTOR APPLIED, DB, VV"));
  params->gps_altitude =
    atof(get_airsar(buf, "GPS ALTITUDE, M"));
  params->lat_peg_point =
    atof(get_airsar(buf, "LATITUDE OF PEG POINT"));
  params->lon_peg_point =
    atof(get_airsar(buf, "LONGITUDE OF PEG POINT"));
  params->head_peg_point =
    atof(get_airsar(buf, "HEADING AT PEG POINT"));
  params->along_track_offset =
    atof(get_airsar(buf, "ALONG-TRACK OFFSET S0  (M)"));
  params->cross_track_offset =
    atof(get_airsar(buf, "CROSS-TRACK OFFSET C0  (M)"));

  return params;
}

airsar_dem_header *read_airsar_dem(const char *dataFile)
{
  airsar_dem_header *dem = NULL;
  airsar_header *header = read_airsar_header(dataFile);
  FILE *fp;
  char *buf;
  int size;

  // Allocate memory and file handling
  dem = (airsar_dem_header *) CALLOC(1, sizeof(airsar_dem_header));
  size = header->first_data_offset - header->dem_header_offset;
  buf = (char *) MALLOC(sizeof(char)*size);
  fp = FOPEN(dataFile, "r");
  FSEEK(fp, header->dem_header_offset, 1);
  if (fgets(buf, size, fp) == NULL)
    asfPrintError("Could not read parameter header\n");
  FCLOSE(fp);

  // Read DEM header
  dem->elevation_offset = atof(get_airsar(buf, "ELEVATION OFFSET (M) ="));
  dem->elevation_increment =
    atof(get_airsar(buf, "ELEVATION INCREMENT (M) ="));
  dem->corner1_lat = atof(get_airsar(buf, "LATITUDE OF CORNER 1 ="));
  dem->corner1_lon = atof(get_airsar(buf, "LONGITUDE OF CORNER 1 ="));
  dem->corner2_lat = atof(get_airsar(buf, "LATITUDE OF CORNER 2 ="));
  dem->corner2_lon = atof(get_airsar(buf, "LONGITUDE OF CORNER 2 ="));
  dem->corner3_lat = atof(get_airsar(buf, "LATITUDE OF CORNER 3 ="));
  dem->corner3_lon = atof(get_airsar(buf, "LONGITUDE OF CORNER 3 ="));
  dem->corner4_lat = atof(get_airsar(buf, "LATITUDE OF CORNER 4 ="));
  dem->corner4_lon = atof(get_airsar(buf, "LONGITUDE OF CORNER 4 ="));
  dem->lat_peg_point =
    atof(get_airsar(buf, "LATITUDE OF PEG POINT ="));
  dem->lon_peg_point =
    atof(get_airsar(buf, "LONGITUDE OF PEG POINT ="));
  dem->head_peg_point =
    atof(get_airsar(buf, "HEADING AT PEG POINT (DEGREES) ="));
  dem->along_track_offset =
    atof(get_airsar(buf, "ALONG-TRACK OFFSET S0 (M) ="));
  dem->cross_track_offset =
    atof(get_airsar(buf, "CROSS-TRACK OFFSET C0 (M) ="));

  return dem;
}

airsar_cal_header *read_airsar_cal(const char *dataFile)
{
  airsar_cal_header *cal = NULL;
  airsar_header *header = read_airsar_header(dataFile);
  FILE *fp;
  char *buf;
  int size;

  // Allocate memory and file handling
  cal = (airsar_cal_header *) CALLOC(1, sizeof(airsar_cal_header));
  size = header->first_data_offset - header->calibration_header_offset;
  buf = (char *) MALLOC(sizeof(char)*size);
  fp = FOPEN(dataFile, "r");
  FSEEK(fp, header->calibration_header_offset, 1);
  if (fgets(buf, size, fp) == NULL)
    asfPrintError("Could not read parameter header\n");
  FCLOSE(fp);

  // Read DEM header
  cal->scale_factor =
    atof(get_airsar(buf, "GENERAL SCALE FACTOR (dB)"));
  cal->hh_amp_cal_factor =
    atof(get_airsar(buf, "HH AMPLITUDE CALIBRATION FACTOR (dB)"));
  cal->hv_amp_cal_factor =
    atof(get_airsar(buf, "HV AMPLITUDE CALIBRATION FACTOR (dB)"));
  cal->vh_amp_cal_factor =
    atof(get_airsar(buf, "VH AMPLITUDE CALIBRATION FACTOR (dB)"));
  cal->vv_amp_cal_factor =
    atof(get_airsar(buf, "VV AMPLITUDE CALIBRATION FACTOR (dB)"));
  cal->hh_phase_cal_factor =
    atof(get_airsar(buf, "HH PHASE CALIBRATION FACTOR (DEGREES)"));
  cal->hv_phase_cal_factor =
    atof(get_airsar(buf, "HV PHASE CALIBRATION FACTOR (DEGREES)"));
  cal->vh_phase_cal_factor =
    atof(get_airsar(buf, "VH PHASE CALIBRATION FACTOR (DEGREES)"));
  cal->vv_phase_cal_factor =
    atof(get_airsar(buf, "VV PHASE CALIBRATION FACTOR (DEGREES)"));
  cal->hh_noise_sigma0 =
    atof(get_airsar(buf, "HH NOISE EQUIVALENT SIGMA ZERO (dB)"));
  cal->hv_noise_sigma0 =
    atof(get_airsar(buf, "HV NOISE EQUIVALENT SIGMA ZERO (dB)"));
  cal->vv_noise_sigma0 =
    atof(get_airsar(buf, "VV NOISE EQUIVALENT SIGMA ZERO (dB)"));
  cal->byte_offset_hh_corr =
    atof(get_airsar(buf, "BYTE OFFSET TO HH CORRECTION VECTOR"));
  cal->byte_offset_hv_corr =
    atof(get_airsar(buf, "BYTE OFFSET TO HV CORRECTION VECTOR"));
  cal->byte_offset_vv_corr =
    atof(get_airsar(buf, "BYTE OFFSET TO VV CORRECTION VECTOR"));
  cal->num_bytes_corr =
    atof(get_airsar(buf, "NUMBER OF BYTES IN CORRECTION VECTORS"));

  return cal;
}

static airsar_general *read_airsar_general(const char *inBaseName)
{
  // Read general metadata file
  airsar_general *general=NULL;
  char line[256]="", *value, *p, *q;

  general = (airsar_general *) CALLOC(1, sizeof(airsar_general));

  char *metaFile = MALLOC(sizeof(char)*(strlen(inBaseName)+20));
  if (!fileExists(inBaseName))
    sprintf(metaFile, "%s_meta.airsar", inBaseName);
  else {
    strcpy(metaFile, inBaseName);
    q = strstr(inBaseName, "_meta.airsar");
    if (q)
      *q = '\0';
  }

  FILE *fpIn = FOPEN(metaFile, "r");

  while (NULL != fgets(line, 255, fpIn)) {
    p = strchr(line, '=');
    if (p) {
      value = p+1;
      if (strncmp(line, "SCENE ID", 8) == 0)
        sprintf(general->scene_id, "%s", value);
      else if (strncmp(line, "Name of Flight-line", 19) == 0)
        sprintf(general->flight_line, "%s", value);
      else if (strncmp(line, "Date of Acquistion", 18) == 0)
        sprintf(general->date_acquisition, "%s", value);
      else if (strncmp(line, "Date of Processing", 18) == 0)
        sprintf(general->date_processing, "%s", value);
      else if (strncmp(line, "Radar Projection (for highest freq product)",
               43) == 0)
        sprintf(general->radar_projection, "%s", value);
      else if (strncmp(line, "Width in Km(for highest freq product)", 37) == 0)
        general->width = atof(value);
      else if (strncmp(line, "Length in Km(for highest freq product)", 38) == 0)
        general->length = atof(value);
      else if (strncmp(line, "Range pixel spacing", 19) == 0)
        general->range_pixel_spacing = atof(value);
      else if (strncmp(line, "Azimuth pixel spacing", 21) == 0)
        general->azimuth_pixel_spacing = atof(value);
      else if (strncmp(line, "Corner 1 latitude in degrees", 28) == 0)
        general->corner1_lat = atof(value);
      else if (strncmp(line, "Corner 1 longitude in degrees", 29) == 0)
        general->corner1_lon = atof(value);
      else if (strncmp(line, "Corner 2 latitude in degrees", 28) == 0)
        general->corner2_lat = atof(value);
      else if (strncmp(line, "Corner 2 longitude in degrees", 29) == 0)
        general->corner2_lon = atof(value);
      else if (strncmp(line, "Corner 3 latitude in degrees", 28) == 0)
        general->corner3_lat = atof(value);
      else if (strncmp(line, "Corner 3 longitude in degrees", 29) == 0)
        general->corner3_lon = atof(value);
      else if (strncmp(line, "Corner 4 latitude in degrees", 28) == 0)
        general->corner4_lat = atof(value);
      else if (strncmp(line, "Corner 4 longitude in degrees", 29) == 0)
        general->corner4_lon = atof(value);
      else if (strncmp(line, "C-band Radar Bandwidth in MHZ", 29) == 0)
        general->c_bandwidth = atof(value);
      else if (strncmp(line, "L-band Radar Bandwidth in MHZ", 29) == 0)
        general->l_bandwidth = atof(value);
      else if (strncmp(line, "P-band Radar Bandwidth in MHZ", 29) == 0)
        general->p_bandwidth = atof(value);
      else if (strncmp(line, "AIRSAR Mode", 11) == 0)
        sprintf(general->mode, "%s", value);
      else if (strncmp(line, "C-band Polarimetric data", 24) == 0) {
        if (strncmp(value, "Yes", 3) == 0)
            general->c_pol_data = TRUE;
        else if (strncmp(value, "None", 4) == 0)
            general->c_pol_data = FALSE;
      }
      else if (strncmp(line, "L-band Polarimetric data", 24) == 0) {
        if (strncmp(value, "Yes", 3) == 0)
            general->l_pol_data = TRUE;
        else if (strncmp(value, "None", 4) == 0)
            general->l_pol_data = FALSE;
      }
      else if (strncmp(line, "P-band Polarimetric data", 24) == 0) {
        if (strncmp(value, "Yes", 3) == 0)
            general->p_pol_data = TRUE;
        else if (strncmp(value, "None", 4) == 0)
            general->p_pol_data = FALSE;
      }
      else if (strncmp(line, "C-band Cross Track Interferometric data", 39) == 0) {
        if (strncmp(value, "Yes", 3) == 0)
            general->c_cross_data = TRUE;
        else if (strncmp(value, "None", 4) == 0)
            general->c_cross_data = FALSE;
      }
      else if (strncmp(line, "L-band Cross Track Interferometric data", 39) == 0) {
        if (strncmp(value, "Yes", 3) == 0)
            general->l_cross_data = TRUE;
        else if (strncmp(value, "None", 4) == 0)
            general->l_cross_data = FALSE;
      }
      else if (strncmp(line, "C-band Along Track Interferometric data", 39) == 0) {
        if (strncmp(value, "Yes", 3) == 0)
            general->c_along_data = TRUE;
        else if (strncmp(value, "None", 4) == 0)
        general->c_along_data = FALSE;
      }
      else if (strncmp(line, "L-band Along Track Interferometric data", 39) == 0) {
        if (strncmp(value, "Yes", 3) == 0)
            general->l_along_data = TRUE;
        else if (strncmp(value, "None", 4) == 0)
            general->l_along_data = FALSE;
      }
      else if (strncmp(line, "Interferometry Baseline Length", 30) == 0)
        sprintf(general->baseline, "%s", value);
      else if (strncmp(line, "Single Frequency/Channel band 1", 31) == 0)
        sprintf(general->frequency_band1, "%s", value);
      else if (strncmp(line, "Single Frequency/Channel band 2", 31) == 0)
        sprintf(general->frequency_band2, "%s", value);
      else if (strncmp(line, "Single Frequency/Channel band 3", 31) == 0)
        sprintf(general->frequency_band3, "%s", value);
      else if (strncmp(line, "Data Format for band 1", 22) == 0)
        sprintf(general->format_band1, "%s", value);
      else if (strncmp(line, "Data Format for band 2", 22) == 0)
        sprintf(general->format_band2, "%s", value);
      else if (strncmp(line, "Data Format for band 3", 22) == 0)
        sprintf(general->format_band3, "%s", value);
    }
  }
  FCLOSE(fpIn);

  free(metaFile);
  return general;
}

meta_parameters *import_airsar_meta(const char *dataName,
				    const char *inBaseName, int force)
{
  airsar_general *general = NULL;
  airsar_header *header = read_airsar_header(dataName);
  airsar_param_header *params = read_airsar_params(dataName);
  airsar_dem_header *dem = NULL;
  airsar_cal_header *cal = NULL;

  if (!force || inBaseName)
    general = read_airsar_general(inBaseName);

  if (!force && (general->c_cross_data || general->l_cross_data)) {
    // Figure out the whether we have a DEM header
    char *demFile;
    demFile = (char *) MALLOC(sizeof(char)*1024);
    int found_c_dem = FALSE, found_l_dem = FALSE, found_p_dem = FALSE;

    // The C-Band DEM is the preferred file for extracting the metadata.
    // Look for that first. If no DEM is around then error out.
    sprintf(demFile, "%s_c.demi2", inBaseName);
    if (fileExists(demFile))
      found_c_dem = TRUE;
    sprintf(demFile, "%s_l.demi2", inBaseName);
    if (fileExists(demFile))
      found_l_dem = TRUE;
    sprintf(demFile, "%s_p.demi2", inBaseName);
    if (fileExists(demFile))
      found_p_dem = TRUE;
    if (!found_c_dem && !found_l_dem && !found_p_dem)
      return NULL;

    // Assign the correct DEM for generating the AirSAR metadata
    if (found_c_dem)
      sprintf(demFile, "%s_c.demi2", inBaseName);
    else if (found_l_dem)
      sprintf(demFile, "%s_l.demi2", inBaseName);
    else if (found_p_dem)
      sprintf(demFile, "%s_p.demi2", inBaseName);

    dem = read_airsar_dem(demFile);
  }
  if (!dem && header->dem_header_offset > 0)
    dem = read_airsar_dem(dataName);
  if (header->calibration_header_offset > 0)
    cal = read_airsar_cal(dataName);

  meta_parameters *ret = airsar2meta(header, params, dem);

  // For old AirSAR data, we need to extract the corner coordinates out of the
  // general header file. If that is not available than we can't do any
  // geocoding later.
  if (strcmp_case(ret->general->sensor, "AIRSAR") == 0 &&
      strstr(ret->general->processor, "3.56")) {
    if (inBaseName) {
      FREE(ret->airsar);
      ret->airsar = NULL;
      ret->location->lat_start_near_range = general->corner1_lat;
      ret->location->lon_start_near_range = general->corner1_lon;
      ret->location->lat_start_far_range = general->corner2_lat;
      ret->location->lon_start_far_range = general->corner2_lon;
      ret->location->lat_end_near_range = general->corner4_lat;
      ret->location->lon_end_near_range = general->corner4_lon;
      ret->location->lat_end_far_range = general->corner3_lat;
      ret->location->lon_end_far_range = general->corner3_lon;
      double lat, lon;
      double yLine, xSamp;
      meta_get_lineSamp(ret, 37.7697, -122.5218, 0.0, &yLine, &xSamp);
      printf("center: line: %.3lf, sample: %.3lf\n", yLine, xSamp);
      meta_get_latLon(ret, 0, 0, 0.0, &lat, &lon);
      printf("start near - lat: %.4lf, lon: %.4lf\n", lat, lon);
      meta_get_lineSamp(ret, ret->location->lat_start_near_range,
			ret->location->lon_start_near_range, 0.0,
			&yLine, &xSamp);
      printf("start near - line: %.3lf, sample: %.3lf, lat: %.4lf, lon: %.4lf\n"
	     ,yLine, xSamp, ret->location->lat_start_near_range,
	     ret->location->lon_start_near_range);
      meta_get_latLon(ret, 0, ret->general->sample_count, 0.0, &lat, &lon);
      printf("start far - lat: %.4lf, lon: %.4lf\n", lat, lon);
      meta_get_lineSamp(ret, ret->location->lat_start_far_range,
			ret->location->lon_start_far_range, 0.0,
			&yLine, &xSamp);
      printf("start far - line: %.3lf, sample: %.3lf, lat: %.4lf, lon: %.4lf\n",
	     yLine, xSamp, ret->location->lat_start_far_range,
	     ret->location->lon_start_far_range);
      meta_get_latLon(ret, ret->general->line_count, 0, 0.0, &lat, &lon);
      printf("end near - lat: %.4lf, lon: %.4lf\n", lat, lon);
      meta_get_lineSamp(ret, ret->location->lat_end_near_range,
			ret->location->lon_end_near_range, 0.0,
			&yLine, &xSamp);
      printf("end near - line: %.3lf, sample: %.3lf, lat: %.4lf, lon: %.4lf\n",
	     yLine, xSamp, ret->location->lat_end_near_range,
	     ret->location->lon_end_near_range);
      meta_get_latLon(ret, ret->general->line_count,
		      ret->general->sample_count, 0.0, &lat, &lon);
      printf("end far - lat: %.4lf, lon: %.4lf\n", lat, lon);
      meta_get_lineSamp(ret, ret->location->lat_end_far_range,
			ret->location->lon_end_far_range, 0.0,
			&yLine, &xSamp);
      printf("end far - line: %.3lf, sample: %.3lf, lat: %.4lf, lon: %.4lf\n",
	     yLine, xSamp, ret->location->lat_end_far_range,
	     ret->location->lon_end_far_range);

    }
    else
      asfPrintWarning("No main header file available to extract location "
		      "information. Can't geocode this data later\n");
  }

  if (general)
    FREE(general);
  FREE(header);
  FREE(params);
  if (dem)
    FREE(dem);

  return ret;
}

//static void fudge_airsar_params(meta_parameters *meta);

int ingest_insar_data(const char *inBaseName, const char *outBaseName,
		      char band)
{
  airsar_header *header;
  meta_parameters *metaIn, *metaOut;
  FILE *fpIn, *fpOut;
  char *inFile=NULL, *outFile=NULL;
  int ii, kk, line_offset, ret=FALSE;
  float *floatBuf=NULL;

  // Generate metadata file
  inFile = (char *) MALLOC(sizeof(char)*255);
  outFile = (char *) MALLOC(sizeof(char)*255);

  // Ingest the DEM
  sprintf(inFile, "%s_%c.demi2", inBaseName, band);
  if (!fileExists(inFile))
    asfPrintStatus("   Could not find DEM file (%s) ...\n", inFile);
  else {
    asfPrintStatus("   Ingesting DEM ...\n");
    sprintf(outFile, "%s_%c_dem.img", outBaseName, band);
    header = read_airsar_header(inFile);
    metaIn = import_airsar_meta(inFile, inBaseName, FALSE);
    line_offset = header->first_data_offset/metaIn->general->sample_count/2;
    metaIn->general->line_count += line_offset;
    metaOut = import_airsar_meta(inFile, inBaseName, FALSE);
    metaIn->general->data_type = INTEGER16;
    metaOut->general->data_type = REAL32;
    floatBuf = (float *) MALLOC(sizeof(float)*metaIn->general->sample_count);
    fpIn = FOPEN(inFile, "rb");
    fpOut = FOPEN(outFile, "wb");
    for (ii=0; ii<metaOut->general->line_count; ii++) {
      get_float_line(fpIn, metaIn, ii+line_offset, floatBuf);
      for (kk=0; kk<metaIn->general->sample_count; kk++)
	floatBuf[kk] = floatBuf[kk]*metaIn->airsar->elevation_increment +
	  metaIn->airsar->elevation_offset;
      put_float_line(fpOut, metaOut, ii, floatBuf);
      asfLineMeter(ii, metaOut->general->line_count);
    }
    FCLOSE(fpIn);
    FCLOSE(fpOut);
    metaOut->general->image_data_type = DEM;
    strcpy(metaOut->general->bands, "DEM");
    //fudge_airsar_params(metaOut);
    meta_write(metaOut, outFile);
    FREE(header);
    ret = TRUE;
  }

  // Ingest amplitude image
  sprintf(inFile, "%s_%c.vvi2", inBaseName, band);
  if (!fileExists(inFile))
    asfPrintStatus("   Could not find amplitude file (%s) ...\n", inFile);
  else {
    asfPrintStatus("   Ingesting amplitude image ...\n");
    sprintf(outFile, "%s_%c_vv.img", outBaseName, band);
    header = read_airsar_header(inFile);
    line_offset = header->first_data_offset/metaIn->general->sample_count/2;
    metaIn->general->line_count = metaOut->general->line_count + line_offset;
    metaIn->general->data_type = INTEGER16;
    metaOut->general->data_type = REAL32;
    floatBuf = (float *) MALLOC(sizeof(float)*metaIn->general->sample_count);
    fpIn = FOPEN(inFile, "rb");
    fpOut = FOPEN(outFile, "wb");
    for (ii=0; ii<metaOut->general->line_count; ii++) {
      get_float_line(fpIn, metaIn, ii+line_offset, floatBuf);
      put_float_line(fpOut, metaOut, ii, floatBuf);
      asfLineMeter(ii, metaOut->general->line_count);
    }
    FCLOSE(fpIn);
    FCLOSE(fpOut);
    metaOut->general->image_data_type = AMPLITUDE_IMAGE;
    strcpy(metaOut->general->bands, "AMP");
    //fudge_airsar_params(metaOut);
    meta_write(metaOut, outFile);
    FREE(header);
    ret = TRUE;
  }

  // Ingest coherence image
  sprintf(inFile, "%s_%c.corgr", inBaseName, band);
  if (!fileExists(inFile))
    asfPrintStatus("   Could not find coherence image (%s) ...\n", inFile);
  else {
    asfPrintStatus("   Ingesting coherence image ...\n");
    sprintf(outFile, "%s_%c_coh.img", outBaseName, band);
    header = read_airsar_header(inFile);
    line_offset = header->first_data_offset / metaIn->general->sample_count;
    metaIn->general->line_count = metaOut->general->line_count + line_offset;
    metaIn->general->data_type = BYTE;
    metaOut->general->data_type = REAL32;
    floatBuf = (float *) MALLOC(sizeof(float)*metaIn->general->sample_count);
    fpIn = FOPEN(inFile, "rb");
    fpOut = FOPEN(outFile, "wb");
    for (ii=0; ii<metaOut->general->line_count; ii++) {
      get_float_line(fpIn, metaIn, ii+line_offset, floatBuf);
      put_float_line(fpOut, metaOut, ii, floatBuf);
      asfLineMeter(ii, metaOut->general->line_count);
    }
    FCLOSE(fpIn);
    FCLOSE(fpOut);
    metaOut->general->image_data_type = COHERENCE_IMAGE;
    strcpy(metaOut->general->bands, "COH");
    //fudge_airsar_params(metaOut);
    meta_write(metaOut, outFile);
    FREE(header);
    ret = TRUE;
  }

  // Clean up
  if (floatBuf)
    FREE(floatBuf);
  if (inFile)
    FREE(inFile);
  if (outFile)
    FREE(outFile);

  return ret;
}

static int sign(char byteBuf)
{
  if (byteBuf < 0)
    return -1;
  else
    return 1;
}

int ingest_polsar_data(const char *inBaseName, const char *outBaseName,
		       radiometry_t radiometry, char band)
{
  FILE *fpIn, *fpOut;
  char *inFile, *outFile;
  int ii, kk, ret;
  char *byteBuf;
  float *power, *shh_amp, *shh_phase, *shv_amp, *shv_phase, *svh_amp;
  float *svh_phase, *svv_amp, *svv_phase;
  float total_power, ysca, amp, phase;
  float m11, m12, m13, m14, m22, m23, m24, m33, m34, m44;
  complexFloat cpx;

  // Allocate memory
  inFile = (char *) MALLOC(sizeof(char)*255);
  outFile = (char *) MALLOC(sizeof(char)*255);

  // Ingest polarimetric data
  sprintf(inFile, "%s_%c.datgr", inBaseName, band);
  if (!fileExists(inFile))
    sprintf(inFile, "%s_%c.dat", inBaseName, band);
  if (!fileExists(inFile)) {
    asfPrintStatus("   Cound not find polarimetric data set (%s_%c) ...\n",
		   inBaseName, band);
    return FALSE;
  }
  else {
    meta_parameters *meta = import_airsar_meta(inFile, inBaseName, FALSE);
    meta->general->data_type = REAL32;
    meta->general->image_data_type = POLARIMETRIC_IMAGE;
    meta->general->band_count = 9;
    if (radiometry == r_AMP)
      strcpy(meta->general->bands,
	     "AMP,AMP_HH,PHASE_HH,AMP_HV,PHASE_HV,AMP_VH,PHASE_VH,"\
	     "AMP_VV,PHASE_VV");
    else if (radiometry == r_SIGMA)
      strcpy(meta->general->bands,
	     "AMP,SIGMA-AMP-HH,SIGMA-PHASE-HH,SIGMA-AMP-HV,SIGMA-PHASE-HV,"\
	     "SIGMA-AMP-VH,SIGMA-PHASE-VH,SIGMA-AMP-VV,SIGMA-PHASE-VV");
    else if (radiometry == r_SIGMA_DB)
      strcpy(meta->general->bands,
	     "AMP,SIGMA_DB-AMP-HH,SIGMA_DB-PHASE-HH,SIGMA_DB-AMP-HV,"\
	     "SIGMA_DB-PHASE-HV,SIGMA_DB-AMP-VH,SIGMA_DB-PHASE-VH,"\
	     "SIGMA_DB-AMP-VV,SIGMA_DB-PHASE-VV");
    power = (float *) MALLOC(sizeof(float)*meta->general->sample_count);
    shh_amp = (float *) MALLOC(sizeof(float)*meta->general->sample_count);
    shh_phase = (float *) MALLOC(sizeof(float)*meta->general->sample_count);
    shv_amp = (float *) MALLOC(sizeof(float)*meta->general->sample_count);
    shv_phase = (float *) MALLOC(sizeof(float)*meta->general->sample_count);
    svh_amp = (float *) MALLOC(sizeof(float)*meta->general->sample_count);
    svh_phase = (float *) MALLOC(sizeof(float)*meta->general->sample_count);
    svv_amp = (float *) MALLOC(sizeof(float)*meta->general->sample_count);
    svv_phase = (float *) MALLOC(sizeof(float)*meta->general->sample_count);
    byteBuf = (char *) MALLOC(sizeof(char)*10);
    airsar_header *header = read_airsar_header(inFile);
//    airsar_param_header *params = read_airsar_params(inFile);
    long offset = header->first_data_offset;
    sprintf(outFile, "%s_%c.img", outBaseName, band);
    fpIn = FOPEN(inFile, "rb");
    fpOut = FOPEN(outFile, "wb");
    FSEEK(fpIn, offset, SEEK_SET);
    for (ii=0; ii<meta->general->line_count; ii++) {
      for (kk=0; kk<meta->general->sample_count; kk++) {
	FREAD(byteBuf, sizeof(char), 10, fpIn);
	// Scale is always 1.0 according to Bruce Chapman
	m11 = ((float)byteBuf[1]/254.0 + 1.5) * pow(2, byteBuf[0]);
	m12 = (float)byteBuf[2] * m11 / 127.0;
	m13 = sign(byteBuf[3]) * SQR((float)byteBuf[3] / 127.0) * m11;
	m14 = sign(byteBuf[4]) * SQR((float)byteBuf[4] / 127.0) * m11;
	m23 = sign(byteBuf[5]) * SQR((float)byteBuf[5] / 127.0) * m11;
	m24 = sign(byteBuf[6]) * SQR((float)byteBuf[6] / 127.0) * m11;
	m33 = (float)byteBuf[7] * m11 / 127.0;
	m34 = (float)byteBuf[8] * m11 / 127.0;
	m44 = (float)byteBuf[9] * m11 / 127.0;
	m22 = 1 - m33 -m44;
	total_power =
	  ((float)byteBuf[1]/254.0 + 1.5) * pow(2, byteBuf[0]);
	ysca = 2.0 * sqrt(total_power);
	power[kk] = sqrt(total_power);
	cpx.real = (float)byteBuf[2] * ysca / 127.0;
	cpx.imag = (float)byteBuf[3] * ysca / 127.0;
	amp = sqrt(cpx.real*cpx.real + cpx.imag*cpx.imag);
	phase = atan2(cpx.imag, cpx.real);
	if (radiometry == r_AMP) {
	  shh_amp[kk] = amp;
	  shh_phase[kk] = phase;
	}
	else if (radiometry == r_SIGMA) {
	  shh_amp[kk] = amp*amp;
	  shh_phase[kk] = phase;
	}
	else if (radiometry == r_SIGMA_DB) {
	  shh_amp[kk] = amp;
	  shh_phase[kk] = phase;
	}
	cpx.real = (float)byteBuf[4] * ysca / 127.0;
	cpx.imag = (float)byteBuf[5] * ysca / 127.0;
	amp = sqrt(cpx.real*cpx.real + cpx.imag*cpx.imag);
	phase = atan2(cpx.imag, cpx.real);
	if (radiometry == r_AMP) {
	  shv_amp[kk] = amp;
	  shv_phase[kk] = phase;
	}
	else if (radiometry == r_SIGMA) {
	  shv_amp[kk] = amp*amp;
	  shv_phase[kk] = phase;
	}
	else if (radiometry == r_SIGMA_DB) {
	  shv_amp[kk] = amp;
	  shv_phase[kk] = phase;
	}
	cpx.real = (float)byteBuf[6] * ysca / 127.0;
	cpx.imag = (float)byteBuf[7] * ysca / 127.0;
	amp = sqrt(cpx.real*cpx.real + cpx.imag*cpx.imag);
	phase = atan2(cpx.imag, cpx.real);
	if (radiometry == r_AMP) {
	  svh_amp[kk] = amp;
	  svh_phase[kk] = phase;
	}
	else if (radiometry == r_SIGMA) {
	  svh_amp[kk] = amp*amp;
	  svh_phase[kk] = phase;
	}
	else if (radiometry == r_SIGMA_DB) {
	  svh_amp[kk] = amp;
	  svh_phase[kk] = phase;
	}
	cpx.real = (float)byteBuf[8] * ysca / 127.0;
	cpx.imag = (float)byteBuf[9] * ysca / 127.0;
	amp = sqrt(cpx.real*cpx.real + cpx.imag*cpx.imag);
	phase = atan2(cpx.imag, cpx.real);
	if (radiometry == r_AMP) {
	  svv_amp[kk] = amp;
	  svv_phase[kk] = phase;
	}
	else if (radiometry == r_SIGMA) {
	  svv_amp[kk] = amp*amp;
	  svv_phase[kk] = phase;
	}
	else if (radiometry == r_SIGMA_DB) {
	  svv_amp[kk] = amp;
	  svv_phase[kk] = phase;
	}
      }
      put_band_float_line(fpOut, meta, 0, ii, power);
      put_band_float_line(fpOut, meta, 1, ii, shh_amp);
      put_band_float_line(fpOut, meta, 2, ii, shh_phase);
      put_band_float_line(fpOut, meta, 3, ii, shv_amp);
      put_band_float_line(fpOut, meta, 4, ii, shv_phase);
      put_band_float_line(fpOut, meta, 5, ii, svh_amp);
      put_band_float_line(fpOut, meta, 6, ii, svh_phase);
      put_band_float_line(fpOut, meta, 7, ii, svv_amp);
      put_band_float_line(fpOut, meta, 8, ii, svv_phase);
      asfLineMeter(ii, meta->general->line_count);
    }
    FCLOSE(fpIn);
    FCLOSE(fpOut);
    meta_write(meta, outFile);

    // Clean up
    FREE(power);
    FREE(shh_amp);
    FREE(shh_phase);
    FREE(shv_amp);
    FREE(shv_phase);
    FREE(svh_amp);
    FREE(svh_phase);
    FREE(svv_amp);
    FREE(svv_phase);
    FREE(inFile);
    FREE(outFile);
    FREE(byteBuf);

    ret = TRUE;
  }

  return ret;
}

void import_airsar(const char *inBaseName, radiometry_t radiometry,
		   const char *outBaseName)
{
  char inFile[1024];
  int found_c_dem = FALSE, found_l_dem = FALSE;

  // Check radiometry
  if (radiometry != r_AMP &&
      radiometry != r_SIGMA && radiometry != r_SIGMA_DB) {
    asfPrintWarning("Radiometry other than AMPLITUDE and SIGMA is not "
		    "supported for AirSAR data.\nDefaulting back to "
		    "AMPLITUDE.\n");
    radiometry = r_AMP;
  }

  airsar_general *general = read_airsar_general(inBaseName);

  // Check for existence of DEMs.
  // The C-Band DEM is the preferred file for extracting the metadata.
  // Look for that first. If no DEM is around then error out.
  if (general->c_cross_data || general->l_cross_data) {
    sprintf(inFile, "%s_c.demi2", inBaseName);
    if (fileExists(inFile))
      found_c_dem = TRUE;
    sprintf(inFile, "%s_l.demi2", inBaseName);
    if (fileExists(inFile))
      found_l_dem = TRUE;
    if (!found_c_dem && found_l_dem)
      asfPrintWarning("Could not find C-band DEM.\nRequired for most reliable "
		      "metadata extraction\n");
    if (!found_c_dem && !found_l_dem)
      asfPrintError("Could not find any DEM.\nCan't reliably import this "
		    "Airsar data.\n");
  }

  // Check for interferometric data
  if (general->c_cross_data) {
    asfPrintStatus("\n   Ingesting C-band cross track interferometric data ..."
           "\n\n");
    ingest_insar_data(inBaseName, outBaseName, 'c');
  }
  if (general->l_cross_data) {
    asfPrintStatus("\n   Ingesting L-band cross track interferometric data ..."
           "\n\n");
    ingest_insar_data(inBaseName, outBaseName, 'l');
  }

  // Kept out the along-track interferometric data for the moment.
  // Only a few data sets were acquired that way and we have no real way
  // to verify the results.

  // Check for polarimetric data
  if (general->c_pol_data) {
    asfPrintStatus("\n   Ingesting C-band polarimetric data ...\n\n");
    ingest_polsar_data(inBaseName, outBaseName, radiometry, 'c');
  }
  if (general->l_pol_data) {
    asfPrintStatus("\n   Ingesting L-band polarimetric data ...\n\n");
    ingest_polsar_data(inBaseName, outBaseName, radiometry, 'l');
  }
  if (general->p_pol_data) {
    asfPrintStatus("\n   Ingesting P-band polarimetric data ...\n\n");
    ingest_polsar_data(inBaseName, outBaseName, radiometry, 'p');
  }

  FREE(general);
}

// The purpose of this code is to refine the values of the along
// and cross track offsets, so that the meta_get_latLon() value
// returned at the corners of the image matches what the correct
// corner lat/lon values are, from the metadata.  Why do we believe
// the corners locations more than the given along/cross track
// offsets?  From the data sets we've been looking at, it seems like
// the corner coords are good, but the offsets aren't...

// We just do a 2-d minimization of the total error (in line/sample
// space -- not lat/lon (minimizing in lat/lon coords would favor
// minimizing the lon difference near the poles)) between the
// corner coordinates.
struct fudge_airsar_params {
        meta_parameters *meta;
};

static int
getObjective(const gsl_vector *x, void *params, gsl_vector *f)
{
  double c0 = gsl_vector_get(x,0);
  double s0 = gsl_vector_get(x,1);

  if (!meta_is_valid_double(c0) || !meta_is_valid_double(s0)) {
    // This does happen sometimes, when we've already found the root
    return GSL_FAILURE;
  }

  struct fudge_airsar_params *p = (struct fudge_airsar_params *)params;
  meta_parameters *meta = p->meta;

  int nl = meta->general->line_count;
  int ns = meta->general->sample_count;

  double old_c0 = meta->airsar->cross_track_offset;
  double old_s0 = meta->airsar->along_track_offset;
  meta->airsar->cross_track_offset = c0;
  meta->airsar->along_track_offset = s0;

  double line, samp, err = 0.;

  meta_get_lineSamp(meta, meta->location->lat_start_near_range,
                          meta->location->lon_start_near_range,
                    0., &line, &samp);
  err += hypot(line,samp);

  meta_get_lineSamp(meta, meta->location->lat_start_far_range,
                          meta->location->lon_start_far_range,
                    0., &line, &samp);
  err += hypot(line,samp-(double)ns);

  meta_get_lineSamp(meta, meta->location->lat_end_near_range,
                          meta->location->lon_end_near_range,
                    0., &line, &samp);
  err += hypot(line-(double)nl,samp);

  meta_get_lineSamp(meta, meta->location->lat_end_far_range,
                          meta->location->lon_end_far_range,
                    0., &line, &samp);
  err += hypot(line-(double)nl,samp-(double)ns);

  //printf("getObjective> [%f,%f] -> %f\n", c0, s0, err);

  gsl_vector_set(f,0,err);
  gsl_vector_set(f,1,err);

  meta->airsar->cross_track_offset = old_c0;
  meta->airsar->along_track_offset = old_s0;

  return GSL_SUCCESS;
}
/*
static void coarse_search(double c0_extent_min, double c0_extent_max,
                          double s0_extent_min, double s0_extent_max,
                          double *c0_min, double *s0_min,
                          meta_parameters *meta)
{
    struct fudge_airsar_params params;
    params.meta = meta;

    double the_min = 9999999;
    double min_c0=99, min_s0=99;
    int i,j,k=6;
    double c0_extent = c0_extent_max - c0_extent_min;
    double s0_extent = s0_extent_max - s0_extent_min;
    gsl_vector *v = gsl_vector_alloc(2);
    gsl_vector *u = gsl_vector_alloc(2);
    //printf("           ");
    //for (j = 0; j <= k; ++j) {
    //    double s0 = s0_extent_min + ((double)j)/k*s0_extent;
    //    printf("%9.3f ", s0);
    //}
    //printf("\n           ");
    //for (j = 0; j <= k; ++j)
    //    printf("--------- ");
    //printf("\n");
    for (i = 0; i <= k; ++i) {
        double c0 = c0_extent_min + ((double)i)/k*c0_extent;
        //printf("%9.3f | ", c0);

        for (j = 0; j <= k; ++j) {
            double s0 = s0_extent_min + ((double)j)/k*s0_extent;

            gsl_vector_set(v, 0, c0);
            gsl_vector_set(v, 1, s0);
            getObjective(v,(void*)(&params), u);
            double n = gsl_vector_get(u,0);
            //printf("%9.3f ", n);
            if (n<the_min) {
                the_min=n;
                min_c0=gsl_vector_get(v,0);
                min_s0=gsl_vector_get(v,1);
            }
        }
        //printf("\n");
    }

    *c0_min = min_c0;
    *s0_min = min_s0;

    gsl_vector_free(v);
    gsl_vector_free(u);
}
*/
/*
static void
generate_start(meta_parameters *meta, double c0, double s0,
               double *start_c0, double *start_s0)
{
    int i;

    double extent_c0_min = -100000. + c0;
    double extent_c0_max = 100000. + c0;

    double extent_s0_min = -100000. + s0;
    double extent_s0_max = 100000. + s0;

    double c0_range = extent_c0_max - extent_c0_min;
    double s0_range = extent_s0_max - extent_s0_min;

    for (i=0; i<12; ++i)
    //for (i=0; i<4; ++i)
    {
        coarse_search(extent_c0_min, extent_c0_max,
                      extent_s0_min, extent_s0_max,
                      start_c0, start_s0, meta);

        c0_range /= 3.;
        s0_range /= 3.;

        extent_c0_min = *start_c0 - c0_range/2.;
        extent_c0_max = *start_c0 + c0_range/2.;

        extent_s0_min = *start_s0 - s0_range/2.;
        extent_s0_max = *start_s0 + s0_range/2.;

        //printf("refining search to region: cross: (%9.3f,%9.3f)\n"
        //       "                           along: (%9.3f,%9.3f)\n",
        //       extent_c0_min, extent_c0_max,
        //       extent_s0_min, extent_s0_max);
    }
}
*/
/*
static void show_error(meta_parameters *meta, const char *descrip)
{
  int nl = meta->general->line_count;
  int ns = meta->general->sample_count;
  double line, samp, err=0;

  asfPrintStatus("Computing known corner coordinates vs. calculated "
                 "metadata differences:\n");

  meta_get_lineSamp(meta, meta->location->lat_start_near_range,
                          meta->location->lon_start_near_range,
                    0., &line, &samp);
  asfPrintStatus("  Start Near: %f (%d) %f (%d)\n", line, 0, samp, 0);
  err += hypot(line, samp);

  meta_get_lineSamp(meta, meta->location->lat_start_far_range,
                          meta->location->lon_start_far_range,
                    0., &line, &samp);
  asfPrintStatus("  Start Far: %f (%d) %f (%d)\n", line, 0, samp, ns);
  err += hypot(line, samp-(double)ns);

  meta_get_lineSamp(meta, meta->location->lat_end_near_range,
                          meta->location->lon_end_near_range,
                    0., &line, &samp);
  asfPrintStatus("  End Near: %f (%d) %f (%d)\n", line, nl, samp, 0);
  err += hypot(line-(double)nl, samp);

  meta_get_lineSamp(meta, meta->location->lat_end_far_range,
                          meta->location->lon_end_far_range,
                    0., &line, &samp);
  asfPrintStatus("  End Far: %f (%d) %f (%d)\n", line, nl, samp, ns);
  err += hypot(line-(double)nl, samp-(double)ns);

  asfPrintStatus("%s, average corner error: %.2f pixels\n", descrip, err/4.);
}
*/
/*
static void fudge_airsar_params(meta_parameters *meta)
{
  int status, iter = 0, max_iter = 1000;
  const gsl_multiroot_fsolver_type *T;
  gsl_multiroot_fsolver *s;
  gsl_error_handler_t *prev;
  struct fudge_airsar_params params;
  const size_t n = 2;
  double c0_initial, s0_initial;
  double out_c0, out_s0;

  params.meta = meta;

  asfPrintStatus("Refining airsar cross-track and along-track offsets to "
         "match corner locations.\n");

  show_error(meta, "Prior to airsar geolocation refinement");

  asfPrintStatus("Prior to coarse refinement (original metadata values):\n");
  asfPrintStatus("  cross track offset: %.1fm\n",
                 meta->airsar->cross_track_offset);
  asfPrintStatus("  along track offset: %.1fm\n",
                 meta->airsar->along_track_offset);
  generate_start(meta, meta->airsar->cross_track_offset,
                 meta->airsar->along_track_offset, &c0_initial, &s0_initial);
  asfPrintStatus("Starting iterative search with:\n");
  asfPrintStatus("  cross track offset: %.1fm\n", c0_initial);
  asfPrintStatus("  along track offset: %.1fm\n", s0_initial);

  gsl_multiroot_function F = {&getObjective, n, &params};
  gsl_vector *x = gsl_vector_alloc(n);

  gsl_vector_set (x, 0, c0_initial);
  gsl_vector_set (x, 1, s0_initial);

  T = gsl_multiroot_fsolver_hybrid;
  s = gsl_multiroot_fsolver_alloc(T, n);
  gsl_multiroot_fsolver_set(s, &F, x);

  prev = gsl_set_error_handler_off();

  do {
    ++iter;
    status = gsl_multiroot_fsolver_iterate(s);

    // abort if stuck
    if (status) break;

    status = gsl_multiroot_test_residual (s->f, 1e-8);
  } while (status == GSL_CONTINUE && iter < max_iter);

  // we allow GSL_ENOPROG (not making progress), since often the coarse
  // search ends up at or very close to the minimum
  if (status == GSL_SUCCESS || status == GSL_ENOPROG) {
    asfPrintStatus("Converged after %d iteration%s.\n",
                   iter, iter==1 ? "" : "s");
    out_c0 = gsl_vector_get(s->x, 0);
    out_s0 = gsl_vector_get(s->x, 1);

    asfPrintStatus("Final values:\n");
    asfPrintStatus("  cross track offset: %.1fm [adjusted by %.1fm]\n",
                   out_c0, fabs(out_c0-meta->airsar->cross_track_offset));
    asfPrintStatus("  along track offset: %.1fm [adjusted by %.1fm]\n",
                   out_s0, fabs(out_s0-meta->airsar->along_track_offset));

    meta->airsar->cross_track_offset = out_c0;
    meta->airsar->along_track_offset = out_s0;

    show_error(meta, "After airsar geolocation refinement");
  }
  else {
    asfPrintStatus("After %d iterations, failed to converge:\n  %s\n"
                   "Metadata parameters left unchanged.\n",
                   iter, gsl_strerror(status));
    out_c0 = out_s0 = 0;
  }

  gsl_multiroot_fsolver_free(s);
  gsl_vector_free(x);
  gsl_set_error_handler(prev);

}
*/

void read_meta_airsar(char *inBaseName, char *outBaseName)
{
  airsar_general *general = read_airsar_general(inBaseName);
  char *inFile = (char *) MALLOC(sizeof(char)*255);
  char *outFile = (char *) MALLOC(sizeof(char)*255);
  airsar_header *header;
  meta_parameters *meta;
  int line_offset;

  if (general->c_cross_data) {
    sprintf(inFile, "%s_c.demi2", inBaseName);
    if (fileExists(inFile)) {
      asfPrintStatus("C-band DEM ...\n");
      sprintf(outFile, "%s_c_dem.xml", outBaseName);
      meta = import_airsar_meta(inFile, inBaseName, FALSE);
      meta->general->data_type = REAL32;
      meta->general->image_data_type = DEM;
      strcpy(meta->general->bands, "DEM");
      meta_write_xml(meta, outFile);
      meta_free(meta);
    }
    sprintf(inFile, "%s_c.vvi2", inBaseName);
    if (fileExists(inFile)) {
      asfPrintStatus("C-band amplitude image ...\n");
      sprintf(outFile, "%s_c_vv.xml", outBaseName);
      meta = import_airsar_meta(inFile, inBaseName, FALSE);
      meta->general->data_type = REAL32;
      meta->general->image_data_type = AMPLITUDE_IMAGE;
      strcpy(meta->general->bands, "AMP");
      meta_write_xml(meta, outFile);
      meta_free(meta);
    }
    sprintf(inFile, "%s_c.corgr", inBaseName);
    if (fileExists(inFile)) {
      asfPrintStatus("C-band coherence image ...\n");
      sprintf(outFile, "%s_c_coh.xml", outBaseName);
      meta = import_airsar_meta(inFile, inBaseName, FALSE);
      meta->general->data_type = REAL32;
      meta->general->image_data_type = COHERENCE_IMAGE;
      strcpy(meta->general->bands, "COH");
      meta_write_xml(meta, outFile);
      meta_free(meta);
    }
  }
  if (general->l_cross_data) {
    sprintf(inFile, "%s_l.demi2", inBaseName);
    if (fileExists(inFile)) {
      asfPrintStatus("L-band DEM ...\n");
      sprintf(outFile, "%s_l_dem.xml", outBaseName);
      meta = import_airsar_meta(inFile, inBaseName, FALSE);
      meta->general->data_type = REAL32;
      meta->general->image_data_type = DEM;
      strcpy(meta->general->bands, "DEM");
      meta_write_xml(meta, outFile);
      meta_free(meta);
    }
    sprintf(inFile, "%s_l.vvi2", inBaseName);
    if (fileExists(inFile)) {
      asfPrintStatus("L-band amplitude image ...\n");
      sprintf(outFile, "%s_l_vv.xml", outBaseName);
      meta = import_airsar_meta(inFile, inBaseName, FALSE);
      meta->general->data_type = REAL32;
      meta->general->image_data_type = AMPLITUDE_IMAGE;
      strcpy(meta->general->bands, "AMP");
      meta_write_xml(meta, outFile);
      meta_free(meta);
    }
    sprintf(inFile, "%s_l.corgr", inBaseName);
    if (fileExists(inFile)) {
      asfPrintStatus("L-band coherence image ...\n");
      sprintf(outFile, "%s_l_coh.xml", outBaseName);
      meta = import_airsar_meta(inFile, inBaseName, FALSE);
      meta->general->data_type = REAL32;
      meta->general->image_data_type = COHERENCE_IMAGE;
      strcpy(meta->general->bands, "COH");
      meta_write_xml(meta, outFile);
      meta_free(meta);
    }
  }
  if (general->c_pol_data) {
    sprintf(inFile, "%s_c.datgr", inBaseName);
    if (!fileExists(inFile))
      sprintf(inFile, "%s_c.dat", inBaseName);
    if (fileExists(inFile)) {
      asfPrintStatus("C-band polarimetric data set ...\n");
      meta = import_airsar_meta(inFile, inBaseName, FALSE);
      meta->general->data_type = REAL32;
      meta->general->image_data_type = POLARIMETRIC_IMAGE;
      meta->general->band_count = 9;
      strcpy(meta->general->bands,
	     "AMP,AMP_HH,PHASE_HH,AMP_HV,PHASE_HV,AMP_VH,PHASE_VH,"	\
	     "AMP_VV,PHASE_VV");
      sprintf(outFile, "%s_c.xml", outBaseName);
      meta_write_xml(meta, outFile);
      meta_free(meta);
    }
  }
  if (general->l_pol_data) {
    sprintf(inFile, "%s_l.datgr", inBaseName);
    if (!fileExists(inFile))
      sprintf(inFile, "%s_l.dat", inBaseName);
    if (fileExists(inFile)) {
      asfPrintStatus("L-band polarimetric data set ...\n");
      meta = import_airsar_meta(inFile, inBaseName, FALSE);
      meta->general->data_type = REAL32;
      meta->general->image_data_type = POLARIMETRIC_IMAGE;
      meta->general->band_count = 9;
      strcpy(meta->general->bands,
	     "AMP,AMP_HH,PHASE_HH,AMP_HV,PHASE_HV,AMP_VH,PHASE_VH,"	\
	     "AMP_VV,PHASE_VV");
      sprintf(outFile, "%s_l.xml", outBaseName);
      meta_write_xml(meta, outFile);
      meta_free(meta);
    }
  }
  if (general->p_pol_data) {
    sprintf(inFile, "%s_p.datgr", inBaseName);
    if (!fileExists(inFile))
      sprintf(inFile, "%s_p.dat", inBaseName);
    if (fileExists(inFile)) {
      asfPrintStatus("P-band polarimetric data set ...\n");
      meta = import_airsar_meta(inFile, inBaseName, FALSE);
      meta->general->data_type = REAL32;
      meta->general->image_data_type = POLARIMETRIC_IMAGE;
      meta->general->band_count = 9;
      strcpy(meta->general->bands,
	     "AMP,AMP_HH,PHASE_HH,AMP_HV,PHASE_HV,AMP_VH,PHASE_VH,"	\
	     "AMP_VV,PHASE_VV");
      sprintf(outFile, "%s_p.xml", outBaseName);
      meta_write_xml(meta, outFile);
      meta_free(meta);
    }
  }

  FREE(general);
}

ext: c

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

void
setupInitialSeed(
		 int id,
		 gsl_rng *rgen );

void createInitialSample(
			 double *x,
			 long N,
			 gsl_rng *rgen,
			 int state_dim );

max_stars_repo_name: heinekmp/AIRPF

max_issues_repo_name: heinekmp/AIRPF

max_forks_repo_name: heinekmp/AIRPF


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <time.h>
#define PI 3.1416
#define USAGE "./mcmc_solar.x n_steps n_burn"

void likelihood(float *chi, float *y_obs, float *y_model, int n_puntos);
void my_model(float *mine, float *anio ,float a, float b, float c, float d, int n_puntos);
void mcmc(int n_steps, float *a, float *b, float *c, float *d, float *l, float *x_obs, float *y_obs,int n_puntos, float *mine);
void min(float *mine, float *x_obs, float *a, float *b, float *c, float *d, float *l, int n_steps, int n_puntos);

void leer(float *anio, float *manc, int n_puntos);
float *reserva(int n_puntos);
void print_array(float *a, float *b, float *c, float *d , float *l, int n_steps, int n_burn);

int main(int argc, char **argv){

	srand48(time(NULL));

	//Valores que entran
	int n_steps=atof(argv[1]);
	int n_burn=atof(argv[2]);
	//Inicializacion de los arrays que guardan los datos
	int n_puntos= 4141;
	float *manc;
	float *anio;
	anio=reserva(n_puntos);
	manc=reserva(n_puntos);

	leer(anio, manc, n_puntos);

	//Array con los resultados del ajuste
	float *mine;
	mine=reserva(n_puntos);

	//Arrays que guardan los datos de las variables
	float *a;
	float *b;
	float *c;
	float *d;
	float *l;
	a=reserva(n_steps);
	b=reserva(n_steps);
	c=reserva(n_steps);
	d=reserva(n_steps);
	l=reserva(n_steps);
	//Primer paso en la cadena
	a[0]=drand48();
	b[0]=drand48();
	c[0]=drand48();
	d[0]=drand48();
	float chi;
	my_model(mine,anio,a[0],b[0],c[0],d[0],n_puntos);
	likelihood(&chi,manc,mine,n_puntos);
	l[0]=chi;

	//Iteraciones
	mcmc(n_steps,a,b,c,d,l,anio,manc,n_puntos,mine);
	min(mine,anio,a,b,c,d,l,n_steps,n_puntos);

	print_array(a,b,c,d,l,n_steps,n_burn);

	free(mine);
	free(anio);
	free(manc);
	free(a);
	free(b);
	free(c);
	free(d);

	return 0;
}

//OBTENCION DEL MEJOR AJUSTE
void min(float *mine, float *x_obs, float *a, float *b, float *c, float *d, float *l, int n_steps, int n_puntos){
	int i;
	float temp, a2, b2, c2, d2;
	temp=l[0];
	a2=a[0];
	b2=b[0];
	c2=c[0];
	d2=d[0];
	for(i=1;i<n_steps;i++){
		if(temp > l[i] ){
			temp=l[i];
			a2=a[i];
			b2=b[i];
			c2=c[i];
			d2=d[i];
		}
	}
	my_model(mine,x_obs,a2,b2,c2,d2,n_puntos);
}

//ITERACIONES
void mcmc(int n_steps, float *a, float *b, float *c, float *d, float *l, float *x_obs, float *y_obs,int n_puntos, float *mine){
	int i;
	float a2,b2,c2,d2, chi, chi2, alpha, beta;
	const gsl_rng_type * T;
	gsl_rng * r;
	gsl_rng_env_setup();
	T = gsl_rng_default;
	r = gsl_rng_alloc (T);

    for(i=0;i<(n_steps-1);i++){
		a2 = a[i] + gsl_ran_gaussian(r, 0.1);
		b2 = b[i] + gsl_ran_gaussian(r, 0.1);
		c2 = c[i] + gsl_ran_gaussian(r, 0.1);
		d2 = c[i] + gsl_ran_gaussian(r, 0.1);

		float *y_init;
		float *y_prime;
		y_init=reserva(n_puntos);
		y_prime=reserva(n_puntos);

		my_model(y_init,x_obs,a[i],b[i],c[i],d[i],n_puntos);
		my_model(y_prime,x_obs,a2,b2,c2,d2,n_puntos);

		likelihood(&chi,y_obs,y_init,n_puntos);
		likelihood(&chi2,y_obs,y_prime,n_puntos);

		alpha = -chi2 + chi;
		if(alpha >= 0.0){
			a[i+1]=a2;
			b[i+1]=b2;
			c[i+1]=c2;
			d[i+1]=d2;
			l[i+1]=chi2;
		}
		else{
			beta = drand48();
			if(alpha < -6000.0 ){
				a[i+1]=a[i];
				b[i+1]=b[i];
				c[i+1]=c[i];
				d[i+1]=d[i];
				l[i+1]=chi;
			}
			else if(beta <= exp(alpha)){
				a[i+1]=a2;
				b[i+1]=b2;
				c[i+1]=c2;
				d[i+1]=d2;
				l[i+1]=chi2;
			}
			else{
				a[i+1]=a[i];
				b[i+1]=b[i];
				c[i+1]=c[i];
				d[i+1]=d[i];
				l[i+1]=chi;
			}
		}

		free(y_init);
		free(y_prime);
	}

  gsl_rng_free (r);
}

//LECTURA DEL ARCHIVO
void leer(float *anio, float *manc, int n_puntos){

	float *mes;
	mes=reserva(n_puntos);

	FILE *in;
	int i;
	in=fopen("monthrg.dat","r");
	if(!in){
    printf("problems opening the file %s\n","monthrg.dat");
    exit(1);
	}
	float temp1;
	float temp2;
	for(i=0;i<n_puntos;i++){
	  do{
		  fscanf(in, "%f %f %f %f %f\n", &anio[i], &mes[i], &temp1, &manc[i] , &temp2);
	  }while(manc[i] == -99.0);
	}
	fclose(in);

	//Calculo del anio como decimal
	for(i=0;i<n_puntos;i++){
		anio[i] = anio[i] + (mes[i]-1.0)/12.0;
	}
	free(mes);
}

//IMPRESION DE LOS ARRAYS
void print_array(float *a, float *b, float *c, float *d , float *l, int n_steps, int n_burn){
  int i;
  for(i=n_burn;i<n_steps;i++){
    printf("%f %f %f %f %f\n", a[i], b[i], c[i] , d[i], l[i]);
  }
}

//EVALUAR LOS PUNTOS CON EL AJUSTE
void my_model(float *mine, float *anio ,float a, float b, float c, float d, int n_puntos){
	int i;
	for(i=0;i<n_puntos;i++){
		mine[i] = a * cos((2*PI/d)*anio[i] + b) + c;
	}
}

//RESERVA DE MEMORIA
float *reserva(int n_puntos){
  float *array;
  int i;
  if(!(array = malloc(n_puntos * sizeof(float)))){
    printf("Problema en reserva\n");
    exit(1);
  }
  for(i=0;i<n_puntos;i++){
    array[i] = 0.0;
  }
  return array;
}

//CALCULO DE LA VARIACION DEL AJUSTE
void likelihood(float *chi, float *y_obs, float *y_model, int n_puntos){
	int i;
	float temp=0;
	for(i=0;i<n_puntos;i++){
		temp=temp+pow(y_obs[i]-y_model[i],2.0);
	}
	*chi=(1.0/2.0)*temp;
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <gsl/gsl_sf_expint.h>
#include "SE_trans.h"

#define ALPHA 0.99
#define DOMD 3.14

double f(double t)
{
  return 0.25 * ((1+t) * log1p(t) + (1-t) * log1p(-t) - 2 * M_LN2) / M_LN2;
}

double G(double x)
{
  return 0.25 * x * SE_trans_div(-1, 1, x) / M_LN2;
}

double J(int j, double h, double x)
{
  return h * (0.5 + M_1_PI * gsl_sf_Si(M_PI * (x / h - j)));
}

double Fapp(double t, int m)
{
  double h = sqrt(M_PI*DOMD / (ALPHA*m));

  int j;
  double sum1 = 0;
  double sum2 = 0;

  for (j = -m; j <  0; j++) {
    sum1 += G(j*h) * J(j, h, SE_trans_inv(-1, 1, t));
  }
  for (j =  m; j >= 0; j--) {
    sum2 += G(j*h) * J(j, h, SE_trans_inv(-1, 1, t));
  }

  return (sum1 + sum2);
}

int main()
{
  int i;
  int n;
  double t;
  double err, maxerr;
  clock_t start, end;
  double time;

  int STEP = 1000;

  for (n = 3; n <= 150; n += 6) {
    start = clock();

    maxerr = 0;
    for (i = - STEP + 1; i <= STEP - 1; i++) {
      t = (double)i / (double)(STEP);

      err = fabs(f(t) - Fapp(t, n));
      maxerr = fmax(err, maxerr);
    }
    end = clock();

    time = (double)(end - start) / CLOCKS_PER_SEC;
    printf("%d\t%e\t%e\n", n, err, time);
  }

  return EXIT_SUCCESS;
}

max_stars_repo_name: okayamat/sinc-indef

max_issues_repo_name: okayamat/sinc-indef

max_forks_repo_name: okayamat/sinc-indef


text:
// Diagonalization Routines using LAPACK(E)
#include <stdio.h>
#include <stdlib.h>

#include <emscripten.h>
#include <lapacke.h>

// void ssyevr_(char *jobz, char *range, char *uplo, int *n, float *a, int *lda,
//              float *vl, float *vu, int *il, int *iu, float *abstol, int *m,
//              float *w, float *z, int *ldz, int *isuppz, float *work, int *lwork,
//              int *iwork, int *liwork, int *info);

void EMSCRIPTEN_KEEPALIVE diagonalize(float *coords, int natoms, int neig, float *evals, float *evecs);

void EMSCRIPTEN_KEEPALIVE diagonalize(float *mtx, int natoms, int neig, float *evals, float *evecs)
{

    // int i;  // debug

    // Parameters
    lapack_int info;

    lapack_int n = natoms;
    lapack_int lda = n;
    lapack_int ldz = neig;

    // printf("natoms = %d | neig = %d \n", natoms, neig);
    lapack_int il = 1;     // 1st eigval idx
    lapack_int iu = neig; // last idx.

    float abstol = -1.0; // default machine value.

    float vl = 0.0;
    float vu = 0.0;
    lapack_int m;
    lapack_int isuppz[natoms];

    // Solve for eigenvalues/eigenvectors
    info = LAPACKE_ssyevr(
        LAPACK_ROW_MAJOR, 'V', 'I', 'U', n, mtx, lda,
        vl, vu, il, iu, abstol, &m, evals, evecs, ldz, isuppz);

    if (info > 0)
    {
        printf("The algorithm failed to compute eigenvalues.\n");
        exit(1);
    }

    // printf("\n The total number of eigenvalues found:%2i\n", m);
    // for (i = 0; i < neig; i++)
    // {
    //     printf("eval %d = %f\n", i, evals[i]);
    // }
}

max_stars_repo_name: JoaoRodrigues/nmaJS

max_issues_repo_name: JoaoRodrigues/nmaJS

max_forks_repo_name: JoaoRodrigues/nmaJS


text:
static const char help[] =
"Build and view a tiny three-triangle mesh using DMPlex, and integrate a\n"
"scalar function over it.  Option prefixes tny_ and plex_view_.\n\n";

/* either build DMPlex via call to DMPlexCreateFromCellList()
    ./tiny
or by directly setting cones "by hand":
    ./tiny -tny_by_hand

compute an integral by symmetric Gauss quadrature:
    ./tiny -tny_integrate_f
    ./tiny -tny_integrate_f -tny_quaddegree 1|[2]|3

compare these views:
    ./tiny -dm_view
    ./tiny -section_view
    ./tiny -v_vec_view

and -plex_view_xxx options (see plexview.h)

parallel refinement works:
    mpiexec -n 2 ./tiny -dm_refine 1 -plex_view_ranges -plex_view_coords

FIXME interpolation of f(x,y) by P2

FIXME interpolate f(x,y) by P2 and integrate it

FIXME add option -tny_element P1|P2|P3
*/

#include <petsc.h>
#include "../../quadrature.h"
#include "plexview.h"

// Describe the mesh "triangle style" with separate numbering for cells and vertices.
static const int    dim = 2,
                    ncell = 3,
                    nvert = 5,
                    cells[9] = {0, 3, 2,  // 9 = ncell * (dim+1)
                                0, 2, 1,
                                2, 3, 4};
static const double coordverts[10] = {0.0, 0.0,  // 10 = nvert * dim
                                      0.0, 1.0,
                                      0.5, 1.0,
                                      1.0, 0.0,
                                      1.0, 1.0};

// Describe same mesh, but directly as DMPlex, i.e. by giving cell and
// edge cones in DAG.  These values are generated by DMPlexCreateFromCellList()
// internally, so they are redundant if we use that create method.
static const int npoint = 15,
                 ccone[3][3] = {{8,9,10},
                                {10,11,12},
                                {9,13,14}},
                 econe[7][2] = {{3,6},
                                {5,6},
                                {3,5},
                                {4,5},
                                {3,4},
                                {6,7},
                                {5,7}};

static double f(double x, double y) {
    return exp(- x - 2 * y);
}

extern PetscErrorCode CreateMeshByHand(DM*);
extern PetscErrorCode CreateCoordinateSectionByHand(DM*);
extern PetscErrorCode CreateSectionP2(DM,PetscSection*);
extern PetscErrorCode IntegrateF(DM,double(double,double),int,double*);
extern PetscErrorCode EvaluateFP2(DM,PetscSection,double(double,double),Vec*);

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    DM            dmplex;
    PetscSection  p2section;
    Vec           v;
    PetscBool     by_hand = PETSC_FALSE,
                  integrate_f = PETSC_FALSE;
    int           quaddegree = 2;

    PetscInitialize(&argc,&argv,NULL,help);
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "tny_", "options for tiny", "");CHKERRQ(ierr);
    ierr = PetscOptionsBool("-by_hand",
        "use by-hand construction",
        "tiny.c", by_hand, &by_hand, NULL);CHKERRQ(ierr);
    ierr = PetscOptionsBool("-integrate_f",
        "integrate f(x,y) by summing over cells and using quadrature",
        "tiny.c", integrate_f, &integrate_f, NULL);CHKERRQ(ierr);
    ierr = PetscOptionsInt("-quaddegree",
        "use this quadrature degree for numerical integrations",
        "tiny.c", quaddegree, &quaddegree, NULL);CHKERRQ(ierr);
    ierr = PetscOptionsEnd();

    // create the DMPlex mesh
    if (by_hand) {
        ierr = CreateMeshByHand(&dmplex); CHKERRQ(ierr);
        ierr = CreateCoordinateSectionByHand(&dmplex); CHKERRQ(ierr);
    } else {
        PetscMPIInt rank;
        ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
        if (rank == 0) { // create mesh on rank 0
            ierr = DMPlexCreateFromCellList(PETSC_COMM_WORLD,
                dim,ncell,nvert,dim+1,
                PETSC_TRUE,  // "interpolate" flag; TRUE means "topologically-interpolate"
                             // i.e. create edges (1D) from vertices (0D) and cells (2D)
                cells,dim,coordverts,
                &dmplex); CHKERRQ(ierr);
        } else { // empty mesh on rank > 0
            ierr = DMPlexCreateFromCellList(PETSC_COMM_WORLD,
                dim,0,0,dim+1,PETSC_TRUE,NULL,dim,NULL,&dmplex); CHKERRQ(ierr);
        }
    }

    // distribute mesh over processes using default partitioner
    {
        DM  distributedMesh = NULL;
        // overlap of 0 is appropriate to P2 etc. FEM:
        ierr = DMPlexDistribute(dmplex, 0, NULL, &distributedMesh);CHKERRQ(ierr);
        if (distributedMesh) {
          ierr = DMDestroy(&dmplex);CHKERRQ(ierr);
          dmplex  = distributedMesh;
        }
    }

    // reset names before viewing
    ierr = PetscObjectSetName((PetscObject)dmplex, "tiny mesh"); CHKERRQ(ierr);
    {
        DM           cdm;
        PetscSection csection;
        ierr = DMGetCoordinateDM(dmplex, &cdm); CHKERRQ(ierr);
        ierr = DMGetDefaultSection(cdm, &csection); CHKERRQ(ierr);
        ierr = PetscObjectSetName((PetscObject)csection, "vertex coordinate section"); CHKERRQ(ierr);
    }
    ierr = DMSetFromOptions(dmplex); CHKERRQ(ierr);
    ierr = DMViewFromOptions(dmplex, NULL, "-dm_view"); CHKERRQ(ierr);
    ierr = PlexViewFromOptions(dmplex); CHKERRQ(ierr);

    // integrate f(x,y) by summing over cells and using quadrature
    if (integrate_f) {
        double       integral;
        const double intfexact = (1.0 - exp(-1.0)) * 0.5 * (1.0 - exp (-2.0));
        ierr = IntegrateF(dmplex,f,quaddegree,&integral); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_WORLD,
            "integral of f(x,y) is %.12f with error %.2e\n",
            integral,fabs(integral - intfexact)); CHKERRQ(ierr);
    }

    // create nodes (degrees of freedom) for P2 elements using PetscSection
    ierr = CreateSectionP2(dmplex,&p2section); CHKERRQ(ierr);
    ierr = DMSetDefaultSection(dmplex, p2section); CHKERRQ(ierr);
    ierr = PetscObjectViewFromOptions((PetscObject)p2section,NULL,"-tny_view_section"); CHKERRQ(ierr);

    // put function f(x,y) into v by local calculations using coordinates and p2section
    ierr = DMCreateGlobalVector(dmplex, &v); CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject)v, "v"); CHKERRQ(ierr);
    ierr = EvaluateFP2(dmplex,p2section,f,&v); CHKERRQ(ierr);
    ierr = PetscObjectViewFromOptions((PetscObject)v,NULL,"-v_vec_view"); CHKERRQ(ierr);

    VecDestroy(&v);  PetscSectionDestroy(&p2section);  DMDestroy(&dmplex);
    return PetscFinalize();
}

/* This function is essentially equivalent to using DMPlexCreateFromCellList().
Note that rank 0 gets the actual mesh and other ranks get an empty mesh.
See the implementations of
    DMPlexBuildFromCellList_Private()
    DMPlexCreateFromCellListParallel()
    DMPlexInterpolate()
    DMPlexBuildCoordinates_Private()      */
PetscErrorCode CreateMeshByHand(DM *dmplex) {
    PetscErrorCode ierr;
    int           j;
    PetscMPIInt   rank;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
    ierr = DMPlexCreate(PETSC_COMM_WORLD,dmplex); CHKERRQ(ierr);
    ierr = DMSetDimension(*dmplex,dim); CHKERRQ(ierr);
    if (rank == 0) {
        // set the total number of points (npoint = ncell + nvert + nedges)
        ierr = DMPlexSetChart(*dmplex, 0, npoint); CHKERRQ(ierr);
        // the points are cells, vertices, edges in that order
        // we only set cones for cells and edges
        for (j = 0; j < ncell; j++) {
            ierr = DMPlexSetConeSize(*dmplex, j, dim+1); CHKERRQ(ierr);
        }
        for (j = ncell + nvert; j < npoint; j++) {
            ierr = DMPlexSetConeSize(*dmplex, j, dim); CHKERRQ(ierr);
        }
        ierr = DMSetUp(*dmplex);
        for (j = 0; j < ncell; j++) {
            ierr = DMPlexSetCone(*dmplex, j, ccone[j]); CHKERRQ(ierr);
        }
        for (j = ncell + nvert; j < npoint; j++) {
            ierr = DMPlexSetCone(*dmplex, j, econe[j-ncell-nvert]); CHKERRQ(ierr);
        }
    } else {
        ierr = DMPlexSetChart(*dmplex, 0, 0); CHKERRQ(ierr);
    }
    // with cones we have only upward directions and no labels for the strata
    // (note: both Symmetrize & Stratify are required, and they must be in this order
    ierr = DMPlexSymmetrize(*dmplex); CHKERRQ(ierr);
    ierr = DMPlexStratify(*dmplex); CHKERRQ(ierr);
    return 0;
}

// Set up a PetscSection which holds vertex coordinates.
PetscErrorCode CreateCoordinateSectionByHand(DM *dmplex) {
    PetscErrorCode ierr;
    PetscSection  coordSection;
    DM            cdm;
    Vec           coordinates;
    double        *acoord;
    int           j, d, dim, vertexstart, vertexend;
    // you have to setup the PetscSection returned by DMGetCoordinateSection() first,
    // or else the Vec returned by DMCreateLocalVector() has zero size
    // (and thus seg faults)
    ierr = DMGetDimension(*dmplex, &dim); CHKERRQ(ierr);
    ierr = DMGetCoordinateSection(*dmplex, &coordSection); CHKERRQ(ierr);
    ierr = DMPlexGetDepthStratum(*dmplex, 0, &vertexstart, &vertexend); CHKERRQ(ierr);
    ierr = PetscSectionSetNumFields(coordSection, 1); CHKERRQ(ierr);
    ierr = PetscSectionSetFieldComponents(coordSection, 0, dim); CHKERRQ(ierr);
    ierr = PetscSectionSetChart(coordSection, vertexstart, vertexend); CHKERRQ(ierr);
    for (j = vertexstart; j < vertexend; j++) {
        ierr = PetscSectionSetDof(coordSection, j, dim); CHKERRQ(ierr);
        ierr = PetscSectionSetFieldDof(coordSection, j, 0, dim); CHKERRQ(ierr);
    }
    ierr = PetscSectionSetUp(coordSection); CHKERRQ(ierr);
    // now we can actually set up the coordinate Vec
    ierr = DMGetCoordinateDM(*dmplex, &cdm); CHKERRQ(ierr);
    ierr = DMCreateLocalVector(cdm, &coordinates); CHKERRQ(ierr);
    ierr = VecSetBlockSize(coordinates,dim); CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates"); CHKERRQ(ierr);
    ierr = VecGetArray(coordinates, &acoord); CHKERRQ(ierr);
    for (j = 0; j < vertexend-vertexstart; j++) {
        for (d = 0; d < dim; ++d) {
            acoord[j*dim+d] = coordverts[j*dim+d];
        }
    }
    ierr = VecRestoreArray(coordinates, &acoord); CHKERRQ(ierr);
    // finally we tell the DM that it has coordinates
    ierr = DMSetCoordinatesLocal(*dmplex, coordinates); CHKERRQ(ierr);
    VecDestroy(&coordinates);
    return 0;
}

// use quadrature to integrate a function by summing over cells
PetscErrorCode IntegrateF(DM dmplex, double f(double x,double y),
                          int quaddegree, double *fint) {
    PetscErrorCode ierr;
    DM              cdm;
    Vec             coords;
    const Quad2DTri q = symmgauss[quaddegree-1];
    const double    *acoords;
    double          x[3], y[3], xr, yr, dx1, dx2, dy1, dy2, absdetJ,
                    csum, fintloc;
    int             numpts, *pts = NULL, voff, j, p, l, r,
                    vertstart, vertend, cellstart, cellend;
    MPI_Comm        comm;

    ierr = DMGetCoordinateDM(dmplex, &cdm); CHKERRQ(ierr);
    ierr = DMGetCoordinatesLocal(dmplex, &coords); CHKERRQ(ierr);
    ierr = DMPlexGetHeightStratum(dmplex, 0, &cellstart, &cellend); CHKERRQ(ierr);
    ierr = DMPlexGetDepthStratum(dmplex, 0, &vertstart, &vertend); CHKERRQ(ierr);
    ierr = VecGetArrayRead(coords, &acoords); CHKERRQ(ierr);
    // integral by sum over cells
    fintloc = 0.0;
    for (j = cellstart; j < cellend; j++) {
        ierr = DMPlexGetTransitiveClosure(dmplex, j, PETSC_TRUE, &numpts, &pts);
        if (numpts != 7) {
            SETERRQ(PETSC_COMM_WORLD,1,"wrong: assume closure of triangle has 7 points\n");
        }
        // record vertex coordinates
        for (l = 0; l < 3; l++) { // loop through vertex points
            p = 8 + 2 * l; // p=0,1 are cell info; p=2,...,7 are edge info; omit orientations
            if ((pts[p] < vertstart) || (pts[p] >= vertend))  {
                SETERRQ(PETSC_COMM_WORLD,2,"wrong: pts[p] should be a vertex\n");
            }
            voff = pts[p] - vertstart;
            x[l] = acoords[2*voff+0];
            y[l] = acoords[2*voff+1];
        }
        // geometry of element (cell)
        dx1 = x[1] - x[0];
        dx2 = x[2] - x[0];
        dy1 = y[1] - y[0];
        dy2 = y[2] - y[0];
        absdetJ = fabs(dx1 * dy2 - dx2 * dy1);
        // sum over quadrature points on cell
        csum = 0.0;
        for (r = 0; r < q.n; r++) {
            xr = x[0] + dx1 * q.xi[r] + dx2 * q.eta[r];
            yr = y[0] + dy1 * q.xi[r] + dy2 * q.eta[r];
            csum += q.w[r] * f(xr,yr);
        }
        // add cell contribution
        fintloc += absdetJ * csum;
        ierr = DMPlexRestoreTransitiveClosure(dmplex, j, PETSC_TRUE, &numpts, &pts); CHKERRQ(ierr);
    }
    ierr = VecRestoreArrayRead(coords, &acoords); CHKERRQ(ierr);

    ierr = PetscObjectGetComm((PetscObject)dmplex,&comm); CHKERRQ(ierr);
    ierr = MPI_Allreduce(&fintloc,fint,1,MPI_DOUBLE,MPI_SUM,comm); CHKERRQ(ierr);
    return 0;
}

PetscErrorCode CreateSectionP2(DM dmplex, PetscSection *section) {
    PetscErrorCode ierr;
    int  j, pstart, pend, vertexstart, edgeend;
    ierr = PetscSectionCreate(PETSC_COMM_WORLD,section); CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject)*section, "P2 scalar section"); CHKERRQ(ierr);
    ierr = PetscSectionSetNumFields(*section, 1); CHKERRQ(ierr);
    ierr = DMPlexGetChart(dmplex, &pstart, &pend); CHKERRQ(ierr);
    ierr = PetscSectionSetChart(*section, pstart, pend); CHKERRQ(ierr);
    ierr = DMPlexGetDepthStratum(dmplex, 0, &vertexstart, NULL); CHKERRQ(ierr);
    ierr = DMPlexGetDepthStratum(dmplex, 1, NULL, &edgeend); CHKERRQ(ierr);
    for (j = pstart; j < pend; ++j) {
        if (j < vertexstart) {
            ierr = PetscSectionSetDof(*section, j, 0); CHKERRQ(ierr);
            ierr = PetscSectionSetFieldDof(*section, j, 0, 0); CHKERRQ(ierr);
        } else {
            ierr = PetscSectionSetDof(*section, j, 1); CHKERRQ(ierr);
            ierr = PetscSectionSetFieldDof(*section, j, 0, 1); CHKERRQ(ierr);
        }
    }
    ierr = PetscSectionSetUp(*section); CHKERRQ(ierr);
    return 0;
}

PetscErrorCode EvaluateFP2(DM dmplex, PetscSection section,
                                double f(double x,double y), Vec *v) {
    PetscErrorCode ierr;
    DM           cdm;
    Vec          vloc, coords;
    double       *avloc, x, y;
    const double *acoords;
    int          numpts, *pts = NULL, dof, off, j, p,
                 vertstart, vertend, edgestart, edgeend, cellstart, cellend;

    // we put values in a local vector, thus redundantly on overlap
    ierr = DMGetLocalVector(dmplex, &vloc); CHKERRQ(ierr);
    ierr = VecGetArray(vloc, &avloc); CHKERRQ(ierr);
    // need coordinates of P2 nodes (i.e. both vertices and edges)
    ierr = DMGetCoordinateDM(dmplex, &cdm); CHKERRQ(ierr);
    ierr = DMGetCoordinatesLocal(dmplex,&coords); CHKERRQ(ierr);
    ierr = DMPlexGetHeightStratum(dmplex, 0, &cellstart, &cellend); CHKERRQ(ierr);
    ierr = DMPlexGetDepthStratum(dmplex, 0, &vertstart, &vertend); CHKERRQ(ierr);
    ierr = DMPlexGetDepthStratum(dmplex, 1, &edgestart, &edgeend); CHKERRQ(ierr);
    ierr = VecGetArrayRead(coords, &acoords); CHKERRQ(ierr);
    for (j = cellstart; j < cellend; j++) {
        ierr = DMPlexGetTransitiveClosure(dmplex, j, PETSC_TRUE, &numpts, &pts);
        for (p = 0; p < numpts*2; p += 2) {   // omit orientations
            PetscSectionGetDof(section, pts[p], &dof);
            if (dof > 0) {
                // compute (x,y) for vertex or edge center from coords
                if (pts[p] < vertstart) {
                    SETERRQ(PETSC_COMM_WORLD,1,"cell center computation not implemented\n");
                } else if (pts[p] < edgestart) {
                    int voff;
                    voff = pts[p] - vertstart;
                    x = acoords[2*voff+0];
                    y = acoords[2*voff+1];
                } else { // pts[p] is an edge ...
                    const int *vpts;
                    int       voff[2];
                    ierr = DMPlexGetCone(dmplex, pts[p], &vpts); CHKERRQ(ierr);
                    voff[0] = vpts[0] - vertstart;
                    voff[1] = vpts[1] - vertstart;
                    x = 0.5 * (acoords[2*voff[0]+0] + acoords[2*voff[1]+0]);
                    y = 0.5 * (acoords[2*voff[0]+1] + acoords[2*voff[1]+1]);
                }
                // get index "off" into Vec vloc based on P2 scalar section
                PetscSectionGetOffset(section, pts[p], &off);
                avloc[off] = f(x,y);
            }
        }
        ierr = DMPlexRestoreTransitiveClosure(dmplex, j, PETSC_TRUE, &numpts, &pts); CHKERRQ(ierr);
    }
    ierr = VecRestoreArray(vloc, &avloc); CHKERRQ(ierr);
    ierr = VecRestoreArrayRead(coords, &acoords); CHKERRQ(ierr);

    // now we want v global, i.e. only values not duplicate ghosts
    ierr = DMLocalToGlobalBegin(dmplex,vloc,INSERT_VALUES,*v); CHKERRQ(ierr);
    ierr = DMLocalToGlobalEnd(dmplex,vloc,INSERT_VALUES,*v); CHKERRQ(ierr);
    ierr = DMRestoreLocalVector(dmplex, &vloc); CHKERRQ(ierr);
    return 0;
}

// FIXME InterpolateFP2()


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
#include <timerlib.h>
#ifndef NOBLAS
#ifdef MKL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#endif
#include "gdrdgemm.h"

#define NMAT 2048
void gdr_check_and_restart(double a[][NMAT],
			   double b[][NMAT],
			   double c[][NMAT])
{
    int try =0;
    static int initialized = 0;
    if (initialized) return;
    while(1){
	int i,j;
	for(i=0;i<NMAT;i++){
	    for(j=0;j<NMAT;j++){
		a[i][j]=0;
		b[i][j]=i*NMAT+j;
		c[i][j]=0;
	    }
	}
	for(i=0;i<NMAT;i++)a[i][i]=1;
	cblas_dgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans,
		     NMAT,NMAT, NMAT, 1.0, a, NMAT, b, NMAT, 0.0, c, NMAT);
	cblas_dgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans,
		     NMAT,NMAT, NMAT, 1.0, a, NMAT, b, NMAT, 0.0, c, NMAT);
	mygdrdgemm(NMAT, NMAT, NMAT, 1.0, (double*)a, NMAT,
		   (double*)b, NMAT, 0.0, (double*) c, NMAT);
	int err = 0;
	for(i=0;i<NMAT;i++){
	    for(j=0;j<NMAT;j++){
		if (b[i][j] != c[i][j]){
		    err ++;
		}
	    }
	}
	if (err == 0){
	    fprintf(stderr,"gdr_check_and_restart passed %d\n", try);
	    initialized=1;
	    return;
	}
	try++;
	fprintf(stderr, "gdr_check_and_restart, err=%d try=%d\n", err, try);
	gdr_free();
	gdr_init();
    }
}

void gdrblas_dgemm
(
#ifndef MKL
   const enum CBLAS_ORDER             ORDER,
   const enum CBLAS_TRANSPOSE             TRANSA,
   const enum CBLAS_TRANSPOSE             TRANSB,
#else
   const CBLAS_ORDER             ORDER,
   const CBLAS_TRANSPOSE             TRANSA,
   const CBLAS_TRANSPOSE             TRANSB,
#endif
   const int                        M,
   const int                        N,
   const int                        K,
   const double                     ALPHA,
   const double *                   A,
   const int                        LDA,
   const double *                   B,
   const int                        LDB,
   const double                     BETA,
   double *                         C,
   const int                        LDC
)
{
    int NOTA, NOTB,gdrdoneflag;
    double alpha = ALPHA, beta = BETA;
    int F77M=M, F77N=N, F77K=K, F77lda=LDA, F77ldb=LDB, F77ldc=LDC;

    if(      TRANSA == CblasNoTrans ){
	NOTA = 1;
    }else{
	NOTA = 0;
    }
    if(      TRANSB == CblasNoTrans ){
	NOTB = 1;
    }else{
	NOTB = 0;
    }

    if( ORDER == CblasColMajor ){
	gdrdoneflag=0;
	gdr_dgemm_(&NOTA,&NOTB,&F77M,&F77N,&F77K, &alpha, &beta,
		   &F77lda, &F77ldb, &F77ldc,A,B,C, &gdrdoneflag);
    }  else  {
	gdrdoneflag=0;
	gdr_dgemm_(&NOTA,&NOTB,&F77N,&F77M,&F77K, &alpha, &beta,
		   &F77ldb, &F77lda, &F77ldc,B,A,C, &gdrdoneflag);
    }

    if(gdrdoneflag!=1){
	cblas_dgemm(ORDER, TRANSA, TRANSB,
		    M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC);
    }
    /*
     * End of HPL_dgemm
     */
}
void dumpcmat(int m, int n, int nc, double c[][nc])
{
    static int callcount = 0;
    static FILE* fid;
    if (callcount == 0){
	fid = fopen("/tmp/matdata", "w");
    }
    callcount ++;
    if(callcount < 8){
	fprintf(fid,"\nPrint CMAT callcount=%d\n",callcount);
	int i, j;
	for(i=0;i<m;i++){
	    fprintf(fid,"\ni=%d\n", i);
	    for(j=0;j<n;j++){
		if ((j%8)==0 )fprintf(fid,"\n%5d:", j);
		fprintf(fid," %20.12e",c[i][j]);
	    }
	}
    }
}
double touchcmat(int m, int n, int nc, double c[][nc])
{
    double sum=0;
    int i, j;
    for(i=0;i<m;i++){
	for(j=0;j<n;j++){
	    sum += c[i][j]*c[i][j];
	}
    }
    return sum;
}

double ssum = 0.0;

void mygdrdgemm(int m,
		int n,
		int k,
		double alpha,
		double * a,
		int na,
		double * b,
		int nb,
		double beta,
		double * c,
		int nc)
{
    int nota=1, notb=1;
    int gdrdoneflag = 0;
    static int first_call = 1;
    if (first_call){
	gdrdgemm_set_procname(MP_myprocid());
	first_call=0;
	init_current_time();
    }
    //    char str[128];
    //    sprintf(str,"before sums = %25.20e %25.20e %25.20e",
    //	    touchcmat(m,k,na,a),  touchcmat(k,n,nb,b), touchcmat(m,n,nc,c));
    //    MP_message(str);
    //    dprintf(9,"mygdrdgemm omp_max_threads=%d procs=%d\n",
    //	    omp_get_max_threads(),omp_get_num_procs());
    double zero=0.0;
    //    int tmp=0;
    //    gdr_dgemm_(&nota, &notb, &n, &m, &k, &zero, &beta, &nb, &na, &nc,
    //	       b, a, c, &tmp);
    gdr_dgemm_(&nota, &notb, &n, &m, &k, &alpha, &beta, &nb, &na, &nc,
	       b, a, c, &gdrdoneflag);
    //	cblas_dgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans,
    //		     m,n, k, alpha, a, na, b, nb, beta, c, nc);
    //	gdrdoneflag=1;
    //    fprintf(stderr,"gdrflag=%d\n", gdrdoneflag);
    if(gdrdoneflag!=1){
	cblas_dgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans,
		     m,n, k, alpha, a, na, b, nb, beta, c, nc);
    }
    //    dumpcmat(m,n,nc,c);
    //    sprintf(str,"after sums = %25.20e %25.20e %25.20e",
    //	    touchcmat(m,k,na,a),  touchcmat(k,n,nb,b), touchcmat(m,n,nc,c));
    //    MP_message(str);
}


text:
#ifndef STDAFX__H
#define STDAFX__H

#if defined(_WIN32)

#include <SDKDDKVer.h>

#if !defined(_STL_EXTRA_DISABLED_WARNINGS)
#define _STL_EXTRA_DISABLED_WARNINGS 4061 4324 4365 4514 4571 4582 4583 4623 4625 4626 4710 4774 4820 4987 5026 5027 5039
#endif

#if !defined(_SCL_SECURE_NO_WARNINGS)
#define _SCL_SECURE_NO_WARNINGS 1
#endif

#if !defined(_CRT_SECURE_NO_WARNINGS)
#define _CRT_SECURE_NO_WARNINGS 1
#endif

#if !defined(_SILENCE_CXX17_ITERATOR_BASE_CLASS_DEPRECATION_WARNING)
#define _SILENCE_CXX17_ITERATOR_BASE_CLASS_DEPRECATION_WARNING 1
#endif

#if !defined(_SILENCE_CXX17_OLD_ALLOCATOR_MEMBERS_DEPRECATION_WARNING)
#define _SILENCE_CXX17_OLD_ALLOCATOR_MEMBERS_DEPRECATION_WARNING 1
#endif

#if !defined(_SILENCE_CXX17_CODECVT_HEADER_DEPRECATION_WARNING)
#define _SILENCE_CXX17_CODECVT_HEADER_DEPRECATION_WARNING 1
#endif

#define STRICT
#define NOMINMAX

#pragma warning(disable: 4571) // warning C4571: Informational: catch(...) semantics changed since Visual C++ 7.1; structured exceptions (SEH) are no longer caught
#pragma warning(disable: 4668) // warning C4668: '%s' is not defined as a preprocessor macro, replacing with '0' for '#if/#elif'
#pragma warning(disable: 4710) // warning C4710: '%s': function not inlined
#pragma warning(disable: 4711) // warning C4711: function '%s' selected for automatic inline expansion
#pragma warning(disable: 4820) // warning C4820: '%s': '%d' bytes padding added after data member '%s'
#pragma warning(disable: 5045) // warning C5045: Compiler will insert Spectre mitigation for memory load if /Qspectre switch specified

#include <Windows.h>

#pragma warning(push)
#pragma warning(disable: 4005) // warning C4005: '%s': macro redefinition
#include <Winternl.h>
#include <ntstatus.h>
#pragma warning(pop)

#else // defined(_WIN32)

#include <cpuid.h>

#endif

#if defined(_MSC_VER)

// currently broken (triggers on iterators and other objects that are usually unnamed)
#pragma warning(disable: 26444) // warning C26444: Avoid unnamed objects with custom construction and destruction (es.84: http://go.microsoft.com/fwlink/?linkid=862923).

// disable for everything
#pragma warning(disable: 4061) // warning C4061: enumerator '%s' in switch of enum '%s' is not explicitly handled by a case label
#pragma warning(disable: 4324) // warning C4234: structure was padded due to alignment specifier
#pragma warning(disable: 4514) // warning C4514: '%s': unreferenced inline function has been removed
#pragma warning(disable: 4623) // warning C4623: '%s': default constructor was implicitly defined as deleted
#pragma warning(disable: 4625) // warning C4625: '%s': copy constructor was implicitly defined as deleted
#pragma warning(disable: 4626) // warning C4626: '%s': assignment operator was implicitly defined as deleted
#pragma warning(disable: 4710) // warning C4710: '%s': function not inlined
#pragma warning(disable: 4820) // warning C4820: '%s': '%d' bytes padding added after data member '%s'
#pragma warning(disable: 5026) // warning C5026: '%s': move constructor was implicitly defined as deleted
#pragma warning(disable: 5027) // warning C5027: '%s': move assignment operator was implicitly defined as deleted

#pragma warning(disable: 26412) // warning C26412: Do not dereference an invalid pointer (lifetimes rule 1). 'return of %s' was invalidated at line %d by 'no initialization'.
#pragma warning(disable: 26426) // warning C26426: Global initializer calls a non-constexpr function '%s' (i.22: http://go.microsoft.com/fwlink/?linkid=853919).
#pragma warning(disable: 26481) // warning C26481: Don't use pointer arithmetic. Use span instead. (bounds.1: http://go.microsoft.com/fwlink/p/?LinkID=620413)
#pragma warning(disable: 26482) // warning C26482: Only index into arrays using constant expressions (bounds.2: http://go.microsoft.com/fwlink/p/?LinkID=620414).
#pragma warning(disable: 26485) // warning C26485: Expression '%s::`vbtable'': No array to pointer decay. (bounds.3: http://go.microsoft.com/fwlink/p/?LinkID=620415)
#pragma warning(disable: 26490) // warning C26490: Don't use reinterpret_cast. (type.1: http://go.microsoft.com/fwlink/p/?LinkID=620417)
#pragma warning(disable: 26499) // warning C26499: Could not find any lifetime tracking information for '%s'

// disable for standard headers
#pragma warning(push)
#pragma warning(disable: 26400) // warning C26400: Do not assign the result of an allocation or a function call with an owner<T> return value to a raw pointer, use owner<T> instead. (i.11 http://go.microsoft.com/fwlink/?linkid=845474)
#pragma warning(disable: 26401) // warning C26401: Do not delete a raw pointer that is not an owner<T>. (i.11: http://go.microsoft.com/fwlink/?linkid=845474)
#pragma warning(disable: 26408) // warning C26408: Avoid malloc() and free(), prefer the nothrow version of new with delete. (r.10 http://go.microsoft.com/fwlink/?linkid=845483)
#pragma warning(disable: 26409) // warning C26409: Avoid calling new and delete explicitly, use std::make_unique<T> instead. (r.11 http://go.microsoft.com/fwlink/?linkid=845485)
#pragma warning(disable: 26411) // warning C26411: The parameter '%s' is a reference to unique pointer and it is never reassigned or reset, use T* or T& instead. (r.33 http://go.microsoft.com/fwlink/?linkid=845479)
#pragma warning(disable: 26412) // warning C26412: Do not dereference an invalid pointer (lifetimes rule 1). 'return of %s' was invalidated at line %d by 'end of function scope (local lifetimes end)'.
#pragma warning(disable: 26413) // warning C26413: Do not dereference nullptr (lifetimes rule 2). 'nullptr' was pointed to nullptr at line %d.
#pragma warning(disable: 26423) // warning C26423: The allocation was not directly assigned to an owner.
#pragma warning(disable: 26424) // warning C26424: Failing to delete or assign ownership of allocation at line %d.
#pragma warning(disable: 26425) // warning C26425: Assigning '%s' to a static variable.
#pragma warning(disable: 26444) // warning C26444: Avoid unnamed objects with custom construction and destruction (es.84: http://go.microsoft.com/fwlink/?linkid=862923).
#pragma warning(disable: 26461) // warning C26461: The reference argument '%s' for function %s can be marked as const. (con.3: https://go.microsoft.com/fwlink/p/?LinkID=786684)
#pragma warning(disable: 26471) // warning C26471: Don't use reinterpret_cast. A cast from void* can use static_cast. (type.1: http://go.microsoft.com/fwlink/p/?LinkID=620417).
#pragma warning(disable: 26481) // warning C26481: Don't use pointer arithmetic. Use span instead. (bounds.1: http://go.microsoft.com/fwlink/p/?LinkID=620413)
#pragma warning(disable: 26482) // warning C26482: Only index into arrays using constant expressions. (bounds.2: http://go.microsoft.com/fwlink/p/?LinkID=620414)
#pragma warning(disable: 26490) // warning C26490: Don't use reinterpret_cast. (type.1: http://go.microsoft.com/fwlink/p/?LinkID=620417)
#pragma warning(disable: 26493) // warning C26493: Don't use C-style casts that would perform a static_cast downcast, const_cast, or reinterpret_cast. (type.4: http://go.microsoft.com/fwlink/p/?LinkID=620420)
#pragma warning(disable: 26494) // warning C26494: Variable '%s' is uninitialized. Always initialize an object. (type.5: http://go.microsoft.com/fwlink/p/?LinkID=620421)
#pragma warning(disable: 26495) // warning C26495: Variable '%s' is uninitialized. Always initialize a member variable. (type.6: http://go.microsoft.com/fwlink/p/?LinkID=620422)
#pragma warning(disable: 26496) // warning C26496: Variable '%s' is assigned only once, mark it as const. (con.4: https://go.microsoft.com/fwlink/p/?LinkID=784969)
#pragma warning(disable: 26497) // warning C26497: This function %s could be marked constexpr if compile-time evaluation is desired. (f.4: https://go.microsoft.com/fwlink/p/?LinkID=784970)

#else

// linux warnings go here

#endif

#include <cstddef>
#include <map>
#include <set>
#include <vector>
#include <string>
#include <algorithm>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <type_traits>
#include <utility>
#include <memory>
#include <tuple>
#include <thread>
#include <cstdlib>
#include <codecvt>

#if defined(_MSC_VER)

// disable additionally for third-party libraries
#pragma warning(push)
#pragma warning(disable: 4456) // warning C4456: declaration of '%s' hides previous local declaration
#pragma warning(disable: 4458) // warning C4458: declaration of '%s' hides class member
#pragma warning(disable: 4459) // warning C4459: declaration of '%s' hides global declaration

#pragma warning(disable: 26429) // warning C26429: Symbol '%s' is never tested for nullness, it can be marked as not_null (f.23: http://go.microsoft.com/fwlink/?linkid=853921).
#pragma warning(disable: 26432) // warning C26432: If you define or delete any default operation in the type '%s', define or delete them all (c.21: http://go.microsoft.com/fwlink/?linkid=853922).
#pragma warning(disable: 26433) // warning C26433: Function '%s' should be marked with 'override' (c.128: http://go.microsoft.com/fwlink/?linkid=853923).
#pragma warning(disable: 26434) // warning C26434: Function '%s' hides a non-virtual function '%s' (c.128: http://go.microsoft.com/fwlink/?linkid=853923).
#pragma warning(disable: 26436) // warning C26436: The type '%s' with a virtual function needs either public, virtual or protected non-virtual destructor (c.35: http://go.microsoft.com/fwlink/?linkid=853924).
#pragma warning(disable: 26439) // warning C26439: This kind of function may not throw. Declare it 'noexcept' (f.6: http://go.microsoft.com/fwlink/?linkid=853927).
#pragma warning(disable: 26440) // warning C26440: Function '%s' can be declared 'noexcept' (f.6: http://go.microsoft.com/fwlink/?linkid=853927).
#pragma warning(disable: 26443) // warning C26443: Overriding destructor should not use explicit 'override' or 'virtual' specifiers (c.128: http://go.microsoft.com/fwlink/?linkid=853923).
#pragma warning(disable: 26462) // warning C26462: The value pointed to by '%s' is assigned only once, mark it as a pointer to const (con.4: https://go.microsoft.com/fwlink/p/?LinkID=784969).
#pragma warning(disable: 26472) // warning C26472: Don't use a static_cast for arithmetic conversions. Use brace initialization, gsl::narrow_cast or gsl::narow (type.1: http://go.microsoft.com/fwlink/p/?LinkID=620417).
#pragma warning(disable: 26474) // warning C26474: Don't cast between pointer types when the conversion could be implicit (type.1: http://go.microsoft.com/fwlink/p/?LinkID=620417).
#pragma warning(disable: 26491) // warning C26491: Don't use static_cast downcasts (type.2: http://go.microsoft.com/fwlink/p/?LinkID=620418).
#pragma warning(disable: 26498) // warning C26498: The function '%s' is constexpr, mark variable '%s' constexpr if compile-time evaluation is desired (con.5: https://go.microsoft.com/fwlink/p/?LinkID=784974).

#endif

#ifndef BOOST_CONFIG_SUPPRESS_OUTDATED_MESSAGE
#define BOOST_CONFIG_SUPPRESS_OUTDATED_MESSAGE
#endif

#include <boost/algorithm/string.hpp>
#include <boost/xpressive/xpressive.hpp>

#include <gsl/gsl>

#include <fmt/format.h>

#if defined(_MSC_VER)

#pragma warning(pop)
#pragma warning(pop)

#else

// linux warning restoration goes here

#endif

#endif

max_stars_repo_name: DrPizza/jsm

max_issues_repo_name: DrPizza/jsm

max_forks_repo_name: DrPizza/jsm


text:
/* ode-initval/odeiv.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program 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
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>

gsl_odeiv2_step *
gsl_odeiv2_step_alloc (const gsl_odeiv2_step_type * T, size_t dim)
{
  gsl_odeiv2_step *s = (gsl_odeiv2_step *) malloc (sizeof (gsl_odeiv2_step));

  if (s == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for ode struct", GSL_ENOMEM);
    };

  s->type = T;
  s->dimension = dim;

  s->state = s->type->alloc (dim);

  if (s->state == 0)
    {
      free (s);                 /* exception in constructor, avoid memory leak */

      GSL_ERROR_NULL ("failed to allocate space for ode state", GSL_ENOMEM);
    };

  return s;
}

const char *
gsl_odeiv2_step_name (const gsl_odeiv2_step * s)
{
  return s->type->name;
}

unsigned int
gsl_odeiv2_step_order (const gsl_odeiv2_step * s)
{
  return s->type->order (s->state);
}

int
gsl_odeiv2_step_apply (gsl_odeiv2_step * s,
                       double t,
                       double h,
                       double y[],
                       double yerr[],
                       const double dydt_in[],
                       double dydt_out[], const gsl_odeiv2_system * dydt)
{
  return s->type->apply (s->state, s->dimension, t, h, y, yerr, dydt_in,
                         dydt_out, dydt);
}

int
gsl_odeiv2_step_reset (gsl_odeiv2_step * s)
{
  return s->type->reset (s->state, s->dimension);
}

void
gsl_odeiv2_step_free (gsl_odeiv2_step * s)
{
  RETURN_IF_NULL (s);
  s->type->free (s->state);
  free (s);
}

int
gsl_odeiv2_step_set_driver (gsl_odeiv2_step * s, const gsl_odeiv2_driver * d)
{
  if (d != NULL)
    {
      s->type->set_driver (s->state, d);
    }
  else
    {
      GSL_ERROR ("driver pointer is null", GSL_EFAULT);
    }

  return GSL_SUCCESS;
}

max_stars_repo_name: igormcoelho/optstats

max_issues_repo_name: igormcoelho/optstats

max_forks_repo_name: igormcoelho/optstats