text:
/**
 * File: nicback_utils.c
 * Subroutine for the NICMOS background subtraction
 *
 */
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include "aXe_grism.h"
#include "aXe_utils.h"
#include "spc_FITScards.h"
#include "spce_sect.h"
#include "spce_fitting.h"
#include "spce_is_in.h"
#include "spc_back.h"
#include "spce_pathlength.h"
#include "trfit_utils.h"
#include "nicback_utils.h"

/**
 * Function: make_nicmos_back
 * The subroutine computes and stores the background image for a specific
 * NICMOS grism image. The non-masked pixels in the grism image are
 * correlated with their corresponding values in the master background image.
 * A linear relation is fitted to the value pairs. Deviating points are
 * excluded using kappa-sigma clipping.
 * Using the slope of the final fit, the background image is computed
 * and stored as a fits image.
 *
 * Parameters:
 * @param obs         - the observation for the grism image
 * @param msk_name    - the full name of the mask file
 * @param master_bck  - the full name of the master sky background
 * @param bck_name    - the full name of the background image
 *
 * Returns:
 * @return -
 */
void
make_nicmos_back(const observation * const obs, const char *msk_name,
		 const char *master_bck, char bck_name[], char plist_name[],
		 const char corr_bck[])
{
  px_point npixels;

  double skypix_frac;

  gsl_matrix *msk_img;
  gsl_matrix *mbck_img;
  gsl_matrix *corr_img=NULL;
  gsl_matrix *grism_bck;

  gsl_vector *lfit;

  fitbck_data *fbck_data;
  //int f_status=0;

  FITScards *cards;

  char  ID[MAXCHAR];

  // fix the extension name
  sprintf (ID, "BCK");

  // load the master background image
  fprintf(stdout,"Loading DATA from: %s...", master_bck);
  mbck_img = FITSimage_to_gsl(master_bck, 1, 1);
  fprintf(stdout,". Done.\n");

  if (strlen(corr_bck) > 0)
    {
      // load the correction background image
      fprintf(stdout,"Loading DATA from: %s...", corr_bck);
      corr_img = FITSimage_to_gsl(corr_bck, 1, 1);
      fprintf(stdout,". Done.\n");
    }

  // load the mask image
  fprintf(stdout,"Loading DATA from: %s...", msk_name);
  msk_img = FITSimage_to_gsl(msk_name, 2, 1);
  fprintf(stdout,". Done.\n");

  // get the dimension of
  // the grism image
  npixels = get_npixel(obs);

  // allocate memory for the fit data
  fbck_data = alloc_fitbck_data(npixels.x * npixels.y);

  // fill the fit data into the structure
  fill_fbck_data(obs, msk_img, mbck_img, fbck_data, plist_name, corr_img);

  // make the fit, including
  // kappa-sigma iterations
  lfit = make_ksig_linfit(fbck_data);

  // report the result of the fit onto the screen
  fprintf(stdout, "\nFitting result: c0 = %f +- %f , c1 = %f +- %f , chi^2 = %f\n\n",
	  gsl_vector_get(lfit, 1), gsl_vector_get(lfit, 3),
	  gsl_vector_get(lfit, 2), gsl_vector_get(lfit, 4), gsl_vector_get(lfit, 6));

  // compute the background for the grism frame
  grism_bck = compute_nicmos_back(mbck_img, gsl_vector_get(lfit, 2),
				  gsl_vector_get(lfit, 1), corr_img);

  // write the background image to a file
  fprintf(stdout,"Writing data to: %s...", bck_name);
  gsl_to_FITSimage (grism_bck, bck_name, 1, ID);
  fprintf(stdout,". Done.\n");

  // determine the fraction of pixels
  // used in the background determination
  skypix_frac = get_skypix_frac(fbck_data);

  // create the fits keywords
  cards = nicbck_info_to_FITScards(skypix_frac, gsl_vector_get(lfit, 2),
				   gsl_vector_get(lfit, 1));

  // transfer the fits keywords
  // to the background image
  put_FITS_cards (bck_name, 2, cards);

  // release memory
  free_fitbck_data(fbck_data);
  gsl_matrix_free(mbck_img);
  if (corr_img != NULL)
    gsl_matrix_free(corr_img);
  gsl_matrix_free(msk_img);
  gsl_matrix_free(grism_bck);
  gsl_vector_free(lfit);
  free_FITScards(cards);
}

/**
 * Function: alloc_fitbck_data
 */
fitbck_data *
alloc_fitbck_data(const int n_data)
{
  fitbck_data *fbck_data;

  fbck_data = (fitbck_data *)malloc(sizeof(fitbck_data));

  fbck_data->x_values = (double *)malloc(n_data * sizeof(double));
  fbck_data->y_values = (double *)malloc(n_data * sizeof(double));
  fbck_data->e_values = (double *)malloc(n_data * sizeof(double));
  fbck_data->x_pos = (int *)malloc(n_data * sizeof(int));
  fbck_data->y_pos = (int *)malloc(n_data * sizeof(int));

  fbck_data->n_data = n_data;

  return fbck_data;
}

/**
 * Function: get_skypix_frac
 * The subroutine computes the fraction of image pixels which were used for
 * determining the scaling factor of the master background and the computation
 * of the bacground image. Pixels not used were masked out, since they are part
 * of an object beam or rejected in the kappa-sigma clipping.
 *
 * Parameters:
 * @param fbck_data - the data used for the fit
 *
 * Returns:
 * @return skypix_frac - fraction of pixels used for background determination
 */
double
get_skypix_frac(const fitbck_data *fbck_data)
{
  int index;
  int nvalues=0;

  double skypix_frac=0.0;

  // go over all data values
  for (index=0; index < fbck_data->n_data; index++)
    {
      // check whether the weight is not NULL
      if (fbck_data->e_values[index])
	// enhance the number of used pixels
	nvalues++;
    }

  // compute the fration of used pixels
  skypix_frac = (double)nvalues / (double)fbck_data->n_data;

  // return the fraction
  return skypix_frac;
}

/**
 * Function: compute_nicmos_back
 * The subroutine applies a scaling factor to a master background
 * in order to get the background image for a specific grism image.
 *
 * Parameters:
 * @param mbck_img - the master background pixel
 * @param factor   - the scaling factor
 *
 * Returns:
 * @return grism_bck - the background image as gsl-matrix
 */
gsl_matrix *
compute_nicmos_back(const gsl_matrix *mbck_img, const double factor,
		    const double offset, const gsl_matrix *corr_img)
{
  gsl_matrix *grism_bck;

  int ii, jj;

  // allocate the space for the background image
  grism_bck = gsl_matrix_alloc(mbck_img->size1, mbck_img->size2);

  // go over each row
  for (ii=0; ii < (int)mbck_img->size1; ii++)
    {
      // go over each column
      for (jj=0; jj < (int)mbck_img->size2; jj++)
	{
	  // compute and set the value in the background image
	  if (corr_img != NULL)
	    // using the pedestal image
	    gsl_matrix_set(grism_bck, ii, jj,
	        (gsl_matrix_get(mbck_img, ii, jj) + gsl_matrix_get(corr_img, ii, jj)) * factor);
	  else
	    // using the fit values only
	    gsl_matrix_set(grism_bck, ii, jj, gsl_matrix_get(mbck_img, ii, jj) * factor + offset);
	}
    }

  // return the background image
  return grism_bck;
}

/**
 * Function: compute_nicmos_back
 * The subroutine iteratively fits a linear relation to the pairs of
 * grism image - master background pixel values. After each fit,
 * pairs with large deviations are determined and excluded.
 * The result of the final fit is returned.
 *
 * Parameters:
 * @param fbck_data - the data used for the fits
 *
 * Returns:
 * @return lfit - the results of the final fit
 */
gsl_vector *
make_ksig_linfit(fitbck_data *fbck_data)
{
  gsl_vector *lfit=NULL;

  int index=0;

  int clipped=1;

  // check whether iterations still
  // must be done and whether
  // the data was changed from clipping
  while (index < N_KSIGBCK_ITER && clipped)
    {
      if (check_fbck_data(fbck_data))
	{
	  // make a non-weighted linear fit
	  lfit = det_vector_linear(fbck_data->x_values, fbck_data->y_values, fbck_data->e_values,
				   fbck_data->n_data, 0);

	  // make a clipping iteration
	  clipped = clipp_fbck_data(fbck_data, lfit, N_KSIGBCK_KAPPA);

	  // inhance the clipping counter
	  index++;
	}
      else
	{
	  // release the space
	  // for the old vector
	  if (lfit)
	    gsl_vector_free(lfit);

	  // get a new dummy vector
	  lfit = get_fbck_defaults();

	  break;
	}
    }

  // check whether the scale is within the boundaries
  if (gsl_vector_get(lfit, 2) < BCK_SCALE_MIN || gsl_vector_get(lfit, 2) > BCK_SCALE_MAX)
    {

      fprintf(stdout, "aXe_NICBACK: Background scale is out of bounds. It is re-adjusted.\n");
      // release the space
      // for the old vector
      gsl_vector_free(lfit);

      // get a new dummy vector
      lfit = get_fbck_defaults();
    }

  // return the fit result
  return lfit;
}

int
check_fbck_data(fitbck_data *fbck_data)
{
  int check=1;
  int ix_min=60000;
  int ix_max=0;
  int iy_min=60000;
  int iy_max=0;
  int index;
  double frac;
  double x_ext, y_ext;

  // get the fraction of background pixel
  frac =  get_skypix_frac(fbck_data);

  if (frac < FRAC_BPIX_MIN)
    {
      check = 0;
      fprintf(stdout, "aXe_NICBACK: Not enough background pixels. Fraction: %e\n", frac);
    }
  else
    {
      // go over all data values
      for (index=0; index < fbck_data->n_data; index++)
	{
	  // check whether the weight is not NULL
	  if (fbck_data->e_values[index])
	    {
	      if (fbck_data->x_pos[index] > ix_max)
		ix_max = fbck_data->x_pos[index];
	      if (fbck_data->x_pos[index] < ix_min)
		ix_min = fbck_data->x_pos[index];
	      if (fbck_data->y_pos[index] > iy_max)
		iy_max = fbck_data->y_pos[index];
	      if (fbck_data->y_pos[index] < iy_min)
		iy_min = fbck_data->y_pos[index];
	    }
	}
      x_ext = (float)(ix_max-ix_min) / NPIX_DIM1;
      y_ext = (float)(iy_max-iy_min) / NPIX_DIM2;

      if (x_ext < AREA_EXTEND_MIN || y_ext < AREA_EXTEND_MIN)
	{
	  check = 0;
	  fprintf(stdout, "aXe_NICBACK: Not enough coverage. xcover,ycover: %e,%e\n", x_ext, y_ext);
	}
    }

  return check;
}

gsl_vector*
get_fbck_defaults()
{
  gsl_vector *ret;

  ret = gsl_vector_alloc(10);

  /*
  gsl_vector_set(ret, 0, xean);
  gsl_vector_set(ret, 1, c0);
  gsl_vector_set(ret, 2, c1);
  gsl_vector_set(ret, 3, cov00);
  gsl_vector_set(ret, 4, cov01);
  gsl_vector_set(ret, 5, cov11);
  gsl_vector_set(ret, 6, chisq);
  */

  gsl_vector_set(ret, 0, 0.0);
  gsl_vector_set(ret, 1, 0.0);
  gsl_vector_set(ret, 2, 1.0);
  gsl_vector_set(ret, 3, 1.0);
  gsl_vector_set(ret, 4, 0.0);
  gsl_vector_set(ret, 5, 1.0);
  gsl_vector_set(ret, 6, 0.0);

  return ret;
}

/**
 * Function: clipp_fbck_data
 * The subroutine applies kappa-sigma clipping to a set of data points.
 * The difference between the measured values and the expected values
 * according to linear regression are determined.
 * Data which is outside of the accepted interval is excluded by
 * setting the weight to '0.0'. If at least one value was clipped,
 * '1' is returned, otherwise '0' (if all values are within the
 * accepted interval.
 *
 * Parameters:
 * @param fbck_data - the data used for the fits
 * @param lfit      - the result of the fit
 * @param kappa     - the kappa value
 *
 * Returns:
 * @return clipped - boolean to indicate clipped values
 */
int
clipp_fbck_data(fitbck_data *fbck_data, const gsl_vector *lfit, const float kappa)
{
  int clipped=0;
  int nclipps=0;

  int ndata=0;
  int index;

  double dy_mean=0.0;
  double stdev=0.0;
  double absdev;
  double dy_act;

  double *y_diff;

  // allocate memory for the tmp-vector
  y_diff = (double *) malloc(fbck_data->n_data * sizeof(double));

  // go over all data values
  for (index=0; index < fbck_data->n_data; index++)
    {
      // check whether the weight is not NULL
      if (fbck_data->e_values[index])
	{
	  // compute the difference from the data point
	  // to the prediction from the fit
	  dy_act = (fbck_data->x_values[index] - gsl_vector_get(lfit, 0)) * gsl_vector_get(lfit, 2)
	    + gsl_vector_get(lfit, 1)-fbck_data->y_values[index];

	  // fill the tmp vector
	  y_diff[ndata] = dy_act;

	  // pre-comute the mean
	  dy_mean += dy_act;

	  // enhance the counter
	  // of non-NULL values
	  ndata++;
	}
    }

  // compute the mean offset
  dy_mean = dy_mean / (double)ndata;

  // compute the standard deviation
  stdev = comp_stdev_from_array(y_diff, ndata, dy_mean);

  // get the maximum allowed
  // deviation
  absdev = kappa*stdev;

  // go over the data
  for (index=0; index < fbck_data->n_data; index++)
    {
      // check whether the data has weigth
      if (fbck_data->e_values[index])
	{
	  // compute the absolute difference from the data point
	  // to the prediction from the fit
	  dy_act = fabs((fbck_data->x_values[index] - gsl_vector_get(lfit, 0)) * gsl_vector_get(lfit, 2)
			+ gsl_vector_get(lfit, 1)-fbck_data->y_values[index]);

	  // check whether the difference
	  // is larger than allowed
	  if (dy_act > absdev)
	    {
	      // make the weight NULL
	      fbck_data->e_values[index] = 0.0;

	      // set the clipped flagg
	      clipped=1;
	      nclipps++;
	    }
	}
    }

  fprintf(stdout, "Number of clipped pixels: %i\n", nclipps);

  // release memory
  free(y_diff);

  // return the integer
  // indicating clipped values
  return clipped;
}

/**
 * Function: comp_stdev_from_array
 * The subroutine computes the standard deviation
 * of a data vector, using the known mean value.
 *
 * Parameters:
 * @param data  - the data vector
 * @param ndata - number of data points to use
 * @param mean  - mean value of the data points
 *
 * Returns:
 * @return clipped - boolean to indicate clipped values
 */
double
comp_stdev_from_array(const double *data, const int ndata, const double mean)
{
  int index=0;

  double stdev=0.0;

  // go over the data array
  for (index=0; index < ndata; index++)
    {
      // compute and add the square difference to the mean
      stdev += (data[index] - mean) * (data[index] - mean);
    }

  // compute the average square difference
  if (ndata > 1)
    stdev /= ((double)ndata-1);

  // compute the final
  // standard deviation
  stdev = sqrt(stdev);

  // return the stdev
  return stdev;
}

/**
 * Function: make_nicmos_back
 * The subroutine creates and fills a structure with pairs of correpsonding
 * from a grism image and a master sky background. Dq-masked pixels values
 * and pixels which are part of an object beam are used, but marked with
 * a zero weight.
 *
 * Parameters:
 * @param obs       - the observation for the grism image
 * @param msk_image - the maks image values
 * @param mbck_img  - the master background values
 * @param fbck_data    - structure with the data pairs
 *
 * Returns:
 * @return -
 */
void
fill_fbck_data(const observation * const obs, const gsl_matrix *msk_img,
	       const gsl_matrix *mbck_img, fitbck_data *fbck_data, char plist_name[],
	       const gsl_matrix *corr_img)
{
  int ix=0;
  int iy=0;
  int index=0;
  char Buffer[10240];
  FILE *fout;

  int px_min=20;
  int px_max=180;
  int py_min=50;
  int py_max=245;
  //int px_min=0;
  //int px_max=255;
  //int py_min=0;
  //int py_max=255;
  //10:240,15:245[10:180,20:245][20:180,50:245]

  // intitialize the index for
  // the fit_data structure
  index = 0;

  double diffval;

  fout = fopen(plist_name, "w");

  // go over all pixels
  // in the grism image
  for (ix=0; ix < (int)obs->grism->size1; ix++)
    {
      for (iy=0; iy < (int)obs->grism->size2; iy++)
	{
	  // set the master background pixel as independent variable
	  fbck_data->x_values[index] = gsl_matrix_get(mbck_img, ix, iy);

	  if (corr_img != NULL)
	    // subtract the correction value from the background value
	    diffval = gsl_matrix_get(obs->grism, ix, iy) - gsl_matrix_get(corr_img, ix, iy);
	  else
	    // use the naked background value
	    diffval = gsl_matrix_get(obs->grism, ix, iy);

	  // set the grism image pixel as dependent variable
	  //fbck_data->y_values[index] = gsl_matrix_get(obs->grism, ix, iy);
	  // set the difference as dependent variable
	  fbck_data->y_values[index] = diffval;

	  // set the x- and y-positions
	  fbck_data->x_pos[index] = ix;
	  fbck_data->y_pos[index] = iy;

	  // check whether the pixel was masked out
	  // or is part of an object
	  if (isnan(gsl_matrix_get(obs->grism, ix, iy))
		    || gsl_matrix_get(msk_img, ix, iy) < -9.0e+05)
	    // mask the pixel in the array
	    fbck_data->e_values[index] = 0.0;
	  else
	    {
	      // give the pixel full weight
	      //x,y,value_back,value_grism_imag
	      fbck_data->e_values[index] = 1.0;
	      sprintf (Buffer, "%i %i %e %e\n",ix, iy, gsl_matrix_get(obs->grism, ix, iy), gsl_matrix_get(mbck_img, ix, iy));
	      fputs (Buffer, fout);
	    }

	  // check whether the pixel is in the
	  // allowed are
	  if (ix < px_min || ix > px_max
	      || iy < py_min || iy > py_max)
	    // mask the pixel in the array
	    fbck_data->e_values[index] = 0.0;
	  // enhance the counter
	  index++;
	}
    }

  fclose (fout);
}

/*
 * Function: free_fitbck_data
 * The function releases the memory allocated in
 * a fit-data structure
 *
 * Parameters:
 * @param fbck_data - the allocated structure
 *
 * Returns:
 * @return -
 */
void
free_fitbck_data(fitbck_data *fbck_data)
{
  free(fbck_data->x_values);
  free(fbck_data->y_values);
  free(fbck_data->e_values);
  free(fbck_data->x_pos);
  free(fbck_data->y_pos);

  free(fbck_data);

  fbck_data = NULL;
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/* integration/qag.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
 *
 * 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 2 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>

#include "initialise.c"
#include "set_initial.c"
#include "qpsrt.c"
#include "util.c"

static int
qag (const gsl_function *f,
     const double a, const double b,
     const double epsabs, const double epsrel,
     const size_t limit,
     gsl_integration_workspace * workspace,
     double * result, double * abserr,
     gsl_integration_rule * q) ;

int
gsl_integration_qag (const gsl_function *f,
		     double a, double b,
		     double epsabs, double epsrel, size_t limit,
		     int key,
		     gsl_integration_workspace * workspace,
		     double * result, double * abserr)
{
  int status ;
  gsl_integration_rule * integration_rule = gsl_integration_qk15 ;

  if (key < GSL_INTEG_GAUSS15)
    {
      key = GSL_INTEG_GAUSS15 ;
    }
  else if (key > GSL_INTEG_GAUSS61)
    {
      key = GSL_INTEG_GAUSS61 ;
    }

  switch (key)
    {
    case GSL_INTEG_GAUSS15:
      integration_rule = gsl_integration_qk15 ;
      break ;
    case GSL_INTEG_GAUSS21:
      integration_rule = gsl_integration_qk21 ;
      break ;
    case GSL_INTEG_GAUSS31:
      integration_rule = gsl_integration_qk31 ;
      break ;
    case GSL_INTEG_GAUSS41:
      integration_rule = gsl_integration_qk41 ;
      break ;
    case GSL_INTEG_GAUSS51:
      integration_rule = gsl_integration_qk51 ;
      break ;
    case GSL_INTEG_GAUSS61:
      integration_rule = gsl_integration_qk61 ;
      break ;
    default:
      GSL_ERROR("value of key does specify a known integration rule",
		GSL_EINVAL) ;
    }

  status = qag (f, a, b, epsabs, epsrel, limit,
                workspace,
                result, abserr,
                integration_rule) ;

  return status ;
}

static int
qag (const gsl_function * f,
     const double a, const double b,
     const double epsabs, const double epsrel,
     const size_t limit,
     gsl_integration_workspace * workspace,
     double *result, double *abserr,
     gsl_integration_rule * q)
{
  double area, errsum;
  double result0, abserr0, resabs0, resasc0;
  double tolerance;
  size_t iteration = 0;
  int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;

  double round_off;

  /* Initialize results */

  initialise (workspace, a, b);

  *result = 0;
  *abserr = 0;

  if (limit > workspace->limit)
    {
      GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
    }

  if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
    {
      GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
		 GSL_EBADTOL);
    }

  /* perform the first integration */

  q (f, a, b, &result0, &abserr0, &resabs0, &resasc0);

  set_initial_result (workspace, result0, abserr0);

  /* Test on accuracy */

  tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));

  /* need IEEE rounding here to match original quadpack behavior */

  round_off = GSL_COERCE_DBL (50 * GSL_DBL_EPSILON * resabs0);

  if (abserr0 <= round_off && abserr0 > tolerance)
    {
      *result = result0;
      *abserr = abserr0;

      GSL_ERROR ("cannot reach tolerance because of roundoff error "
		 "on first attempt", GSL_EROUND);
    }
  else if ((abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0.0)
    {
      *result = result0;
      *abserr = abserr0;

      return GSL_SUCCESS;
    }
  else if (limit == 1)
    {
      *result = result0;
      *abserr = abserr0;

      GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);
    }

  area = result0;
  errsum = abserr0;

  iteration = 1;

  do
    {
      double a1, b1, a2, b2;
      double a_i, b_i, r_i, e_i;
      double area1 = 0, area2 = 0, area12 = 0;
      double error1 = 0, error2 = 0, error12 = 0;
      double resasc1, resasc2;
      double resabs1, resabs2;

      /* Bisect the subinterval with the largest error estimate */

      retrieve (workspace, &a_i, &b_i, &r_i, &e_i);

      a1 = a_i;
      b1 = 0.5 * (a_i + b_i);
      a2 = b1;
      b2 = b_i;

      q (f, a1, b1, &area1, &error1, &resabs1, &resasc1);
      q (f, a2, b2, &area2, &error2, &resabs2, &resasc2);

      area12 = area1 + area2;
      error12 = error1 + error2;

      errsum += (error12 - e_i);
      area += area12 - r_i;

      if (resasc1 != error1 && resasc2 != error2)
	{
	  double delta = r_i - area12;

	  if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)
	    {
	      roundoff_type1++;
	    }
	  if (iteration >= 10 && error12 > e_i)
	    {
	      roundoff_type2++;
	    }
	}

      tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));

      if (errsum > tolerance)
	{
	  if (roundoff_type1 >= 6 || roundoff_type2 >= 20)
	    {
	      error_type = 2;	/* round off error */
	    }

	  /* set error flag in the case of bad integrand behaviour at
	     a point of the integration range */

	  if (subinterval_too_small (a1, a2, b2))
	    {
	      error_type = 3;
	    }
	}

      update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);

      retrieve (workspace, &a_i, &b_i, &r_i, &e_i);

      iteration++;

    }
  while (iteration < limit && !error_type && errsum > tolerance);

  *result = sum_results (workspace);
  *abserr = errsum;

  if (errsum <= tolerance)
    {
      return GSL_SUCCESS;
    }
  else if (error_type == 2)
    {
      GSL_ERROR ("roundoff error prevents tolerance from being achieved",
		 GSL_EROUND);
    }
  else if (error_type == 3)
    {
      GSL_ERROR ("bad integrand behavior found in the integration interval",
		 GSL_ESING);
    }
  else if (iteration == limit)
    {
      GSL_ERROR ("maximum number of subdivisions reached", GSL_EMAXITER);
    }
  else
    {
      GSL_ERROR ("could not integrate function", GSL_EFAILED);
    }
}

ext: c


text:
/**
 *
 * @file testing_zhegv.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Azzam Haidar
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_zmain.h"

#undef REAL
#define COMPLEX

static int check_orthogonality(int itype, int uplo, int N,
                               PLASMA_Complex64_t *Z, int LDZ,
                               PLASMA_Complex64_t *B, PLASMA_Complex64_t *CHOLB, int LDB,
                               double eps);
static int check_reduction(int itype, int uplo, int N, double *D,
                           PLASMA_Complex64_t *A, int LDA,
                           PLASMA_Complex64_t *B, int LDB,
                           PLASMA_Complex64_t *Z, int LDZ,
                           double eps );
static int check_solution(int N, double *E1, double *E2, double eps);

int testing_zhegv(int argc, char **argv)
{
    /* Check for number of arguments*/
    if (argc != 3) {
        USAGE("HEGV", "N LDA LDB",
              "   - N    : size of the matrices A and B\n"
              "   - LDA  : leading dimension of the matrix A\n"
              "   - LDB  : leading dimension of the matrix B\n");
        return -1;
    }

    double      eps = LAPACKE_dlamch_work('e');
    PLASMA_enum vec = PlasmaVec;
    int    N        = atoi(argv[0]);
    int    LDA      = atoi(argv[1]);
    int    LDB      = atoi(argv[2]);
    int    LDQ      = LDA;
    int    LDAxN    = LDA*N;
    int    LDBxN    = LDB*N;
    int    LDQxN    = LDQ*N;

    int info_ortho     = 0;
    int info_solution  = 0;
    int info_reduction = 0;
    int i, u;

    PLASMA_Complex64_t *A1    = (PLASMA_Complex64_t *)malloc(LDAxN*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *A2    = (PLASMA_Complex64_t *)malloc(LDAxN*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *B1    = (PLASMA_Complex64_t *)malloc(LDBxN*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *B2    = (PLASMA_Complex64_t *)malloc(LDBxN*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *Q     = (PLASMA_Complex64_t *)malloc(LDQxN*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *Ainit = (PLASMA_Complex64_t *)malloc(LDAxN*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *Binit = (PLASMA_Complex64_t *)malloc(LDBxN*sizeof(PLASMA_Complex64_t));
    double *W1 = (double *)malloc(N*sizeof(double));
    double *W2 = (double *)malloc(N*sizeof(double));
    PLASMA_Complex64_t *work = (PLASMA_Complex64_t *)malloc(3*N* sizeof(PLASMA_Complex64_t));
    PLASMA_desc *T;

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)||(!B1)||(!B2)||(!Q)||(!Ainit)||(!Binit)){
        printf("Out of Memory \n ");
        return -2;
    }

    PLASMA_Disable(PLASMA_AUTOTUNING);
    PLASMA_Set(PLASMA_TILE_SIZE, 120);
    PLASMA_Set(PLASMA_INNER_BLOCK_SIZE, 20);

    PLASMA_Enable(PLASMA_WARNINGS);
    PLASMA_Enable(PLASMA_ERRORS);
    PLASMA_Alloc_Workspace_zhegv(N, N, &T);

    /*----------------------------------------------------------
    *  TESTING ZHEGV
    */
    /* Initialize A1 and Ainit */
    PLASMA_zplghe(0., N, A1, LDA, 5198);
    PLASMA_zlacpy(PlasmaUpperLower, N, N, A1, LDA, Ainit, LDA);

    /* Initialize B1 and Binit */
    PLASMA_zplghe((double)N, N, B1, LDB, 4321 );
    PLASMA_zlacpy(PlasmaUpperLower, N, N, B1, LDB, Binit, LDB);

    printf("\n");
    printf("------ TESTS FOR PLASMA ZHEGV ROUTINE -------  \n");
    printf("        Size of the Matrix %d by %d\n", N, N);
    printf("\n");
    printf(" The matrix A is randomly generated for each test.\n");
    printf("============\n");
    printf(" The relative machine precision (eps) is to be %e \n",eps);
    printf(" Computational tests pass if scaled residuals are less than 60.\n");

    /*----------------------------------------------------------
     *  TESTING ZHEGV
     */

    for (i=0; i<3; i++) {
        for (u=0; u<2; u++) {
            memcpy(A2, Ainit, LDAxN*sizeof(PLASMA_Complex64_t));
            memcpy(B2, Binit, LDBxN*sizeof(PLASMA_Complex64_t));
            /* CALL ZHEGV with itype= 1, 2, 3 and uplo= L,U */
            PLASMA_zhegv(itype[i], vec, uplo[u], N, A2, LDA, B2, LDB, W2, T, Q, LDQ);

            /* CALL LAPACK to compute eigenvalues. */
            memcpy(A1, Ainit, LDAxN*sizeof(PLASMA_Complex64_t));
            memcpy(B1, Binit, LDBxN*sizeof(PLASMA_Complex64_t));
            LAPACKE_zhegv( LAPACK_COL_MAJOR,
                     itype[i], lapack_const(vec), lapack_const(uplo[u]),
                     N, A1, LDA, B1, LDB, W1);

            /* CHECK */
            info_solution = check_solution(N, W1, W2, eps);
            /* Check the orthogonality, reduction and the eigen solutions */
            if (vec == PlasmaVec){
                info_ortho = check_orthogonality(itype[i], uplo[u], N, Q, LDQ, Binit, B1, LDB, eps);
                info_reduction = check_reduction(itype[i], uplo[u], N, W2, Ainit, LDA, Binit, LDB, Q, LDQ, eps );
            }

            if ( (info_ortho == 0) & (info_reduction == 0) & (info_solution == 0)) {
                printf("***************************************************\n");
                printf(" ---- TESTING ZHEGV (%s, %s) ...................... PASSED !\n", itypestr[i], uplostr[u]);
                printf("***************************************************\n");
            }
            else {
                printf("************************************************\n");
                printf(" - TESTING ZHEGV (%s, %s) ... FAILED !\n", itypestr[i], uplostr[u]);
                printf("************************************************\n");
            }
        }
    }

    PLASMA_Dealloc_Handle_Tile(&T);
    free(A1);
    free(A2);
    free(B1);
    free(B2);
    free(Q);
    free(Ainit);
    free(Binit);
    free(W1);
    free(W2);
    free(work);

    return 0;
}

/*-------------------------------------------------------------------
 * Check the orthogonality of Z
 */
static int check_orthogonality(int itype, int uplo, int N,
                               PLASMA_Complex64_t *Z, int LDZ,
                               PLASMA_Complex64_t *B, PLASMA_Complex64_t *CHOLB, int LDB,
                               double eps)
{
    static PLASMA_Complex64_t zone  =  1.0;
    static PLASMA_Complex64_t mzone = -1.0;
    static PLASMA_Complex64_t zzero =  0.0;
    static double             done  =  1.0;
    static double             mdone = -1.0;
    PLASMA_enum trans;
    double  normQ, result;
    int     info_ortho;
    double *work = (double *)malloc(N*sizeof(double));
    PLASMA_Complex64_t *TEMP    = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *Id = (PLASMA_Complex64_t *) malloc(N*N*sizeof(PLASMA_Complex64_t));
    char *str;

    /* Build the idendity matrix */
    LAPACKE_zlaset_work(LAPACK_COL_MAJOR, 'A', N, N, 0., 1., Id, N);

    /* Perform orth test*/
    if ( itype == 3 ) {
        /*
         * if ITYPE = 3, Z**H*inv(B)*Z = I.
         * inv(B) = inv(L)**H * inv(L) or inv(U)*inv(U)**H
         */
        str = "Z**H*inv(B)*Z";

        if( uplo==PlasmaUpper ) {
            trans = PlasmaConjTrans;
        } else {
            trans = PlasmaNoTrans;
        }

        /* Compute inv(L)*Z or inv(U)**H*Z */
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaUpperLower),
                            N, N, Z, LDZ, TEMP, N);

        cblas_ztrsm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, (CBLAS_TRANSPOSE)trans,
                    CblasNonUnit, N, N,
                    CBLAS_SADDR(zone), CHOLB, LDB,
                                       TEMP,  N);

        /*
         * Compute Z**H*inv(B)*Z-Id = (Z**H*inv(L)**H) * (inv(L)*Z) - Id
         *
         * Note: Z**H*inv(L)**H is the ConjTranspose of the previous result, so we use ZHERK
         */
        cblas_zherk(CblasColMajor, (CBLAS_UPLO)uplo, CblasConjTrans,
                    N, N, done, TEMP, N, mdone, Id, N);
    } else {
        /*
         * if ITYPE = 1 or 2, Z**H*B*Z = I;
         */
        str = "Z**H*B*Z";
        cblas_zhemm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, N, N,
                    CBLAS_SADDR(zone),  B,    LDB,
                                        Z,    LDZ,
                    CBLAS_SADDR(zzero), TEMP, N);
        cblas_zgemm(CblasColMajor, CblasConjTrans, CblasNoTrans, N, N, N,
                    CBLAS_SADDR(zone),  Z,    LDZ,
                                        TEMP, N,
                    CBLAS_SADDR(mzone), Id  , N);
    }

    normQ = LAPACKE_zlanhe_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), 'U', N, Id, N, work);
    result = normQ / (N * eps);
    printf(" ======================================================\n");
    printf(" | Id-%s |_oo / (N*eps)          : %15.3E \n",str, result );
    printf(" ======================================================\n");

    if ( isnan(result) || isinf(result) || (result > 60.0) ) {
        printf("-- Orthogonality is suspicious ! \n");
        info_ortho=1;
    }
    else {
        printf("-- Orthogonality is CORRECT ! \n");
        info_ortho=0;
    }

    free(work); free(Id); free(TEMP);
    return info_ortho;
}

/*------------------------------------------------------------
 *  Check the reduction
 */
static int check_reduction(int itype, int uplo, int N, double *D,
                           PLASMA_Complex64_t *A, int LDA,
                           PLASMA_Complex64_t *B, int LDB,
                           PLASMA_Complex64_t *Z, int LDZ,
                           double eps )
{
    PLASMA_Complex64_t zone  =  1.0;
    PLASMA_Complex64_t zzero =  0.0;
    PLASMA_Complex64_t mzone = -1.0;
    double Anorm, Znorm, Rnorm, result;
    int info_reduction;
    int i;
    char *str;
    PLASMA_Complex64_t *TEMP     = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *Residual = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
    double *work = (double *)malloc(N*sizeof(double));

    if( itype == 1 ) {
        /*
         * | A Z - B Z D | / ( |A| |Z| n ulp )
         */
        str = " A Z - B Z D ";

        /* Compute TEMP = Z*D */
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaUpperLower),
                            N, N, Z, LDZ, TEMP, N);

        for (i=0; i<N; i++) {
            cblas_zdscal(N, D[i], TEMP + i*N, 1);
        }

        /* Compute Residual = B*Z*D = B*TEMP */
        cblas_zhemm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, N, N,
                    CBLAS_SADDR(zone),  B,        LDB,
                                        TEMP,     N,
                    CBLAS_SADDR(zzero), Residual, N);

        /* Compute A*Z - B*Z*D */
        cblas_zhemm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, N, N,
                    CBLAS_SADDR(zone),  A,        LDA,
                                        Z,        LDZ,
                    CBLAS_SADDR(mzone), Residual, N);
    }
    else if( itype == 2 ) {
        /*
         * | A B Z - Z D | / ( |A| |Z| n ulp )
         */
        str = " A B Z - Z D ";

        /* Compute Residual = Z*D */
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaUpperLower),
                            N, N, Z, LDZ, Residual, N);

        for (i=0; i<N; i++) {
            cblas_zdscal(N, D[i], Residual + i*N, 1);
        }

        /* Compute TEMP = B*Z */
        cblas_zhemm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, N, N,
                    CBLAS_SADDR(zone),  B,    LDB,
                                        Z,    LDZ,
                    CBLAS_SADDR(zzero), TEMP, N);

        /* Compute A*B*Z - Z*D = A*TEMP-Residual */
        cblas_zhemm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, N, N,
                    CBLAS_SADDR(zone),  A,        LDA,
                                        TEMP,     N,
                    CBLAS_SADDR(mzone), Residual, N);
    } else {
        /*
         * | B A Z - Z D | / ( |A| |Z| n ulp )
         */
        str = " B A Z - Z D ";

        /* Compute Residual = Z*D */
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaUpperLower),
                            N, N, Z, LDZ, Residual, N);

        for (i=0; i<N; i++) {
            cblas_zdscal(N, D[i], Residual + i*N, 1);
        }

        /* Compute TEMP = A*Z */
        cblas_zhemm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, N, N,
                    CBLAS_SADDR(zone),  A,    LDA,
                                        Z,    LDZ,
                    CBLAS_SADDR(zzero), TEMP, N);

        /* Compute B*A*Z - Z*D = B*TEMP-Residual */
        cblas_zhemm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, N, N,
                    CBLAS_SADDR(zone),  B,        LDB,
                                        TEMP,     N,
                    CBLAS_SADDR(mzone), Residual, N);
    }

    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaOneNorm), N, N, Residual, N,   work);
    Anorm = LAPACKE_zlanhe_work(LAPACK_COL_MAJOR, lapack_const(PlasmaOneNorm), lapack_const(uplo), N, A,        LDA, work);
    Znorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaOneNorm), N, N, Z,        LDZ, work);

    result = Rnorm / ( Anorm * Znorm * N * eps);
    printf(" ======================================================\n");
    printf(" | %s |_oo/(|A|_oo.|Z|_oo.N.eps) : %15.3E \n", str, result );
    printf(" ======================================================\n");

    if ( isnan(result) || isinf(result) || (result > 60.0) ) {
        printf("-- Reduction is suspicious ! \n");
        info_reduction = 1;
    }
    else {
        printf("-- Reduction is CORRECT ! \n");
        info_reduction = 0;
    }

    free(TEMP); free(Residual); free(work);
    return info_reduction;
}
/*--------------------------------------------------------------
 * Check the solution
 */
static int check_solution(int N, double *E1, double *E2, double eps)
{
    int info_solution, i;
    double resid;
    double maxtmp;
    double maxel = fabs( fabs(E1[0]) - fabs(E2[0]) );
    double maxeig = max( fabs(E1[0]), fabs(E2[0]) );

    for (i=1; i<N; i++) {
        resid  = fabs(fabs(E1[i])-fabs(E2[i]));
        maxtmp = max(fabs(E1[i]), fabs(E2[i]));

        /* Update */
        maxeig = max(maxtmp, maxeig);
        maxel  = max(resid,  maxel );
    }

    maxel = maxel / (maxeig * N * eps);
    printf(" ======================================================\n");
    printf(" | D - eigcomputed | / (|D| * N * eps) : %15.3E \n",  maxel );
    printf(" ======================================================\n");

    printf("============\n");
    printf("Checking the eigenvalues of A\n");
    if ( isnan(maxel) || isinf(maxel) || (maxel > 100) ) {
        printf("-- The eigenvalues are suspicious ! \n");
        info_solution = 1;
    }
    else{
        printf("-- The eigenvalues are CORRECT ! \n");
        info_solution = 0;
    }
    return info_solution;
}