text:
/**
 *
 * @file pdpltmg_house.c
 *
 *  PLASMA auxiliary 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 Mathieu Faverge
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:45:13 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *  Parallel tile Householder matrix generation
 *  sequential call to either static or dynamic scheduling
 *
 * See http://www.mathworks.fr/fr/help/matlab/ref/gallery.html#f84-999993
 *
 * Householder matrix
 *
 * Generates a random column vector of size M, and returns the housholder matrix
 * H = eye(n,n) - beta*v*v' that satisfies the relationship
 *
 * H*x = -sign(x(1))*norm(x)*e1
 *
 * where e1 is the first column of eye(n,n). Note that if x is complex, then
 * sign(x) exp(i*arg(x)) (which equals x./abs(x) when x is nonzero).
 *
 */
void plasma_dpltmg_house( PLASMA_desc A, unsigned long long int seed,
                          PLASMA_sequence *sequence, PLASMA_request *request )
{
    plasma_context_t *plasma;
    Quark_Task_Flags task_flags = Quark_Task_Flags_Initializer;
    PLASMA_desc descV;
    double *V;
    double tau;

    plasma = plasma_context_self();
    if (sequence->status != PLASMA_SUCCESS)
        return;
    QUARK_Task_Flag_Set(&task_flags, TASK_SEQUENCE, (intptr_t)sequence->quark_sequence);

    /*
     * Submit A initialization, such that other thread have work while main
     * thread is initializing V
     */
    plasma_dynamic_call_6(plasma_pdlaset,
                          PLASMA_enum,        PlasmaUpperLower,
                          double, 0.,
                          double, 1.,
                          PLASMA_desc,        A,
                          PLASMA_sequence*,   sequence,
                          PLASMA_request*,    request);

    /*
     * Allocate the descriptor of the Householder reflector.
     * Exploit the fact that for one column, the PlasmaCM and PlasmaCCRB storage
     * are identical
     */
    V = (double*) plasma_shared_alloc( plasma, A.m, PlasmaRealDouble );
    descV = plasma_desc_init(
        PlasmaRealDouble, A.mb, A.nb, A.mb*A.nb,
        A.m, 1, 0, 0, A.m, 1 );
    descV.mat = V;

    /* Initialize Householder vector */
    {
        CORE_dplrnt( A.m, 1, V, A.m, A.m, 0, 0, seed );
        LAPACKE_dlarfg_work( A.m, V, V+1, 1, &tau );
        V[0] = 1.;
    }

    /* Computes A = A - tau * V * V' (ger == gemm) */
    plasma_parallel_call_9(plasma_pdgemm,
                           PLASMA_enum,        PlasmaNoTrans,
                           PLASMA_enum,        PlasmaTrans,
                           double, (-tau),
                           PLASMA_desc,        descV,
                           PLASMA_desc,        descV,
                           double, 1.,
                           PLASMA_desc,        A,
                           PLASMA_sequence*,   sequence,
                           PLASMA_request*,    request);

    plasma_dynamic_sync();
    plasma_shared_free( plasma, V );
}


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Sun Nov  7 20:44:21 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -hc2real 10 */

/*
 * This function contains 34 FP additions, 14 FP multiplications,
 * (or, 26 additions, 6 multiplications, 8 fused multiply/add),
 * 20 stack variables, and 20 memory accesses
 */
static const fftw_real K500000000 = FFTW_KONST(+0.500000000000000000000000000000000000000000000);
static const fftw_real K1_902113032 = FFTW_KONST(+1.902113032590307144232878666758764286811397268);
static const fftw_real K1_175570504 = FFTW_KONST(+1.175570504584946258337411909278145537195304875);
static const fftw_real K2_000000000 = FFTW_KONST(+2.000000000000000000000000000000000000000000000);
static const fftw_real K1_118033988 = FFTW_KONST(+1.118033988749894848204586834365638117720309180);

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.41 1999/05/26 15:44:14 fftw Exp $
 * $Id: fft.ml,v 1.43 1999/05/17 19:44:18 fftw Exp $
 * $Id: to_c.ml,v 1.25 1999/10/26 21:41:32 stevenj Exp $
 */

void fftw_hc2real_10(const fftw_real *real_input, const fftw_real *imag_input, fftw_real *output, int real_istride, int imag_istride, int ostride)
{
     fftw_real tmp3;
     fftw_real tmp11;
     fftw_real tmp23;
     fftw_real tmp31;
     fftw_real tmp20;
     fftw_real tmp30;
     fftw_real tmp10;
     fftw_real tmp28;
     fftw_real tmp14;
     fftw_real tmp16;
     fftw_real tmp18;
     fftw_real tmp19;
     ASSERT_ALIGNED_DOUBLE;
     {
	  fftw_real tmp1;
	  fftw_real tmp2;
	  fftw_real tmp21;
	  fftw_real tmp22;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp1 = real_input[0];
	  tmp2 = real_input[5 * real_istride];
	  tmp3 = tmp1 - tmp2;
	  tmp11 = tmp1 + tmp2;
	  tmp21 = imag_input[4 * imag_istride];
	  tmp22 = imag_input[imag_istride];
	  tmp23 = tmp21 - tmp22;
	  tmp31 = tmp21 + tmp22;
     }
     tmp18 = imag_input[2 * imag_istride];
     tmp19 = imag_input[3 * imag_istride];
     tmp20 = tmp18 - tmp19;
     tmp30 = tmp18 + tmp19;
     {
	  fftw_real tmp6;
	  fftw_real tmp12;
	  fftw_real tmp9;
	  fftw_real tmp13;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp4;
	       fftw_real tmp5;
	       fftw_real tmp7;
	       fftw_real tmp8;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp4 = real_input[2 * real_istride];
	       tmp5 = real_input[3 * real_istride];
	       tmp6 = tmp4 - tmp5;
	       tmp12 = tmp4 + tmp5;
	       tmp7 = real_input[4 * real_istride];
	       tmp8 = real_input[real_istride];
	       tmp9 = tmp7 - tmp8;
	       tmp13 = tmp7 + tmp8;
	  }
	  tmp10 = tmp6 + tmp9;
	  tmp28 = K1_118033988 * (tmp6 - tmp9);
	  tmp14 = tmp12 + tmp13;
	  tmp16 = K1_118033988 * (tmp12 - tmp13);
     }
     output[5 * ostride] = tmp3 + (K2_000000000 * tmp10);
     {
	  fftw_real tmp32;
	  fftw_real tmp34;
	  fftw_real tmp29;
	  fftw_real tmp33;
	  fftw_real tmp27;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp32 = (K1_175570504 * tmp30) - (K1_902113032 * tmp31);
	  tmp34 = (K1_902113032 * tmp30) + (K1_175570504 * tmp31);
	  tmp27 = tmp3 - (K500000000 * tmp10);
	  tmp29 = tmp27 - tmp28;
	  tmp33 = tmp28 + tmp27;
	  output[7 * ostride] = tmp29 - tmp32;
	  output[3 * ostride] = tmp29 + tmp32;
	  output[ostride] = tmp33 - tmp34;
	  output[9 * ostride] = tmp33 + tmp34;
     }
     output[0] = tmp11 + (K2_000000000 * tmp14);
     {
	  fftw_real tmp24;
	  fftw_real tmp26;
	  fftw_real tmp17;
	  fftw_real tmp25;
	  fftw_real tmp15;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp24 = (K1_175570504 * tmp20) - (K1_902113032 * tmp23);
	  tmp26 = (K1_902113032 * tmp20) + (K1_175570504 * tmp23);
	  tmp15 = tmp11 - (K500000000 * tmp14);
	  tmp17 = tmp15 - tmp16;
	  tmp25 = tmp16 + tmp15;
	  output[2 * ostride] = tmp17 - tmp24;
	  output[8 * ostride] = tmp17 + tmp24;
	  output[6 * ostride] = tmp25 - tmp26;
	  output[4 * ostride] = tmp25 + tmp26;
     }
}

fftw_codelet_desc fftw_hc2real_10_desc =
{
     "fftw_hc2real_10",
     (void (*)()) fftw_hc2real_10,
     10,
     FFTW_BACKWARD,
     FFTW_HC2REAL,
     235,
     0,
     (const int *) 0,
};

ext: c

max_stars_repo_name: albertsgrc/ftdock-opt

max_issues_repo_name: albertsgrc/ftdock-opt

max_forks_repo_name: albertsgrc/ftdock-opt


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

double square(double x, void* params){
  return x*x;
}

double square_int(double x){
  return (x*x*x)/3;
}

double root(double x, void * params){
  return sqrt(x);
}

double root_int(double x){
  return (2.0/3) * sqrt(x*x*x);
}

double osc(double x, void* params){
  double a = *(double *) params;
  return a * cos(a* x);
}

double osc_int(double x, double a){
  return sin(a* x);
}

double rev(double x, void* params){
  return 1.0/ (x*x);
}

double rev_int(double x){
  return (1.0) /x;
}

int
main (void)
{
  gsl_integration_workspace * w
    = gsl_integration_workspace_alloc (1000);

  double result;
  double error;
  gsl_function F;
  printf("method,funct,result,expected,est_error,intervals\n" );

  F.function = □
  // F.params = α
  size_t neval = 5;
  double desired_error;

  //////////////////////////////
  // QNG
  desired_error = 1.0;
  for(int i=0; i<6; i++){
    F.function = □

    gsl_integration_qng(&F, 0, 1, 0, desired_error, &result, &error, &neval);
    printf("qng,square,%.18f,%.18f,%.18f,%zd\n",result,square_int(1),error,neval );
    F.function = root;
    gsl_integration_qng(&F, 0.001, 1, 0, desired_error, &result, &error, &neval);
    printf("qng,root,%.18f,%.18f,%.18f,%zd\n",result,root_int(1),error,neval );
    desired_error = desired_error / 10;
  }

  /////////////////////

  ////////////////////////
  //QAG

  desired_error = 1.0;

  neval = 15;
  for(int i=0; i<6; i++){
    F.function = □
    gsl_integration_qag(&F, 0, 1, 0, desired_error, 1000, i, w, &result, &error);
    printf("qag,square,%.18f,%.18f,%.18f,%zd\n",result,square_int(1),error,w->size);
    F.function = root;
    gsl_integration_qag(&F, 0.001, 1, 0, desired_error, 1000, i, w, &result, &error);
    printf("qag,root,%.18f,%.18f,%.18f,%zd\n",result,root_int(1),error,w->size );
    desired_error = desired_error / 10;
    if(i==0) neval+=6;
    else neval +=10;
  }

  /////////////////////

  ////////////////////////
  //QAGS

  desired_error = 1.0;

  neval = 15;
  for(int i=0; i<6; i++){
    F.function = □
    gsl_integration_qags(&F, 0, 1, 0, desired_error, 1000, w, &result, &error);
    printf("qags,square,%.18f,%.18f,%.18f,%zd\n",result,square_int(1),error,w->size);
    F.function = root;
    gsl_integration_qags(&F, 0.001, 1, 0, desired_error, 1000, w, &result, &error);
    printf("qags,root,%.18f,%.18f,%.18f,%zd\n",result,root_int(1),error,w->size);
    desired_error = desired_error / 10;
    if(i==0) neval+=6;
    else neval +=10;
  }
  /////

  ////////////////////////
  //QAGP

  double pts[2];
  pts[0] = 0;
  pts[1] = 1;
  desired_error = 1.0;

  neval = 15;
  for(int i=0; i<6; i++){
    F.function = □
    gsl_integration_qagp(&F, pts, 2, 0, desired_error, 1000, w, &result, &error);
    printf("qags,square,%.18f,%.18f,%.18f,%zd\n",result,square_int(1),error,w->size);
    F.function = root;
    gsl_integration_qagp(&F,pts, 2, 0, desired_error, 1000, w, &result, &error);
    printf("qags,root,%.18f,%.18f,%.18f,%zd\n",result,root_int(1),error,w->size);
    desired_error = desired_error / 10;
    if(i==0) neval+=6;
    else neval +=10;
  }
  /////

  /////////////////
  // Oscillation
  //////////
  // QNG
  double stop;
  double alpha;
  alpha = 4.0;
  F.params = α
  F.function = &osc;
  desired_error = 1.0;
  double start=1.0;
  stop = (9.0 / (2.0 * M_PI));
  for(int i=0; i<6; i++){
    gsl_integration_qng(&F, start, stop, 0, desired_error, &result, &error, &neval);
    printf("qng,osc,%.18f,%.18f,%.18f,%zd\n",result,(osc_int(stop, alpha) - osc_int(start, alpha)),error,neval);
    desired_error = desired_error / 10;
  }

  // QAG
  alpha = 4.0;
  F.params = α
  F.function = &osc;
  desired_error = 1.0;
  stop = (9.0 / (2.0 * M_PI));
  neval=15;
  for(int i=0; i<6; i++){
    gsl_integration_qag(&F, start, stop, 0, desired_error, 1000, i, w, &result, &error);
    printf("qag,osc,%.18f,%.18f,%.18f,%zd\n",result,(osc_int(stop, alpha) - osc_int(start, alpha)),error,neval);
    desired_error = desired_error / 10;
    if(i==0) neval+=6;
    else neval +=10;
  }

  //QAWO
  gsl_integration_qawo_table * wf = gsl_integration_qawo_table_alloc(alpha, (stop - start), GSL_INTEG_SINE, 2);
  desired_error = 1.0;

  for(int i=0; i<6; i++){
    gsl_integration_qawo(&F, start, 0, desired_error, 1000, w, wf, &result, &error);
    printf("qawo,osc,%.18f,%.18f,%.18f,%zd\n",result,(osc_int(stop, alpha) - osc_int(start, alpha)),error, w->size);
    desired_error = desired_error / 10;
    if(i==0) neval+=6;
    else neval +=10;
  }
  gsl_integration_qawo_table_free(wf);

  ///QAGIU
  F.function = &rev;
  desired_error = 1.0;
  start = 2.5;
  for(int i=0; i<6; i++){
    gsl_integration_qagiu(&F, start, 0, desired_error, 1000, w, &result, &error);
    printf("qagiu,osc,%.18f,%.18f,%.18f,%zd\n",result,(rev_int(start)),error, w->size);
    desired_error = desired_error / 10;
    if(i==0) neval+=6;
    else neval +=10;
  }

  gsl_integration_workspace_free (w);
  return 0;
}


text:
#ifndef KEXPR_H
#define KEXPR_H
#pragma warning( disable : 991)

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <plplot/plplot.h>
#include <stdint.h>
#include <stdio.h>
#include "stack.h"
#include "utf8.h"
#include "gdate.h"
#include <setjmp.h>

#if defined(_MSC_VER) || defined(_WIN32)
#include <Windows.h>
#endif

#include "api.h"

#ifdef _DEBUG
void ke_validate_parameter_qte(sml_t* sml, token_t *p, int nb_param, char * function_name);
void ke_validate_parameter_vtype(sml_t* sml, token_t * p, int vtype, char * param_name, char * function_name);
void ke_validate_parameter_ttype(sml_t* sml, token_t * p, int ttype, char * function_name);
void ke_validate_parameter_not_null(sml_t* sml, token_t * p, void * ptr, char * param_name, char * function_name);
void ke_validate_parameter_int_gt_zero(sml_t* sml, token_t * p, char * param_name, char * function_name);
#endif // _DEBUG

void ke_set_jump(sml_t *sml, int info);
sml_t * ke_create_sml();
void ke_free_sml(sml_t *sml);
void ke_init_memory_count(sml_t *sml);
char *ke_mystrndup(sml_t *sml, char *src, size_t n);
utf8 *ke_mystr(sml_t *sml, utf8 *src, size_t n);
token_t * ke_get_tok(sml_t *sml);
token_t * ke_get_tokidx(sml_t *sml, int idx);
fncp ke_function(sml_t *sml, char * name);
typedef int(*DLLPROC)();
void ke_import(sml_t *sml, char * s);
void ke_load_dll(sml_t *sml, char * p);
void ke_free_val(sml_t *sml);
void ke_free_tokens(sml_t *sml);
int ke_fill_list(sml_t *sml, kexpr_t *ke);
void ke_set_int(sml_t* sml, token_t *tokp, int64_t y);
void ke_set_real(sml_t* sml, token_t *tokp, double x);
void ke_set_vector(sml_t * sml, token_t *tokp, gsl_vector * vecp);
void ke_set_vector_int(sml_t * sml, token_t *tokp, gsl_vector_int * vecp);
void ke_set_date(sml_t * sml, token_t *tokp, GDate_t * datep);
void ke_set_matrix(sml_t *sml, token_t *tokp, gsl_matrix * matp);
void ke_set_complex(sml_t *sml, token_t *tokp, gsl_complex complex);
static void ke_set_str_internal(sml_t* sml, token_t *tokp, char * tmp);
void ke_set_str(sml_t* sml, token_t *tokp, char *x);
void ke_set_str_direct(sml_t *sml, token_t * e, char *x);
void ke_set_record(sml_t *sml, token_t *source, token_t *dest);
void ke_set_image(sml_t *sml, token_t *source, token_t *dest);
void ke_set_file(sml_t *sml, token_t *source, token_t *dest);
void ke_set_buffer(sml_t *sml, token_t *source, token_t *dest);
void ke_print_one_stack(token_t * tokp);
void ke_print_one(sml_t *sml, token_t * tokp);
void ke_print_stack(sml_t* sml, token_t *tokp, int top);
void ke_print_range(kexpr_t *kexpr, int starti, int endi);
void ke_print_error_one(sml_t * sml, kexpr_t *kexpr, char * name, token_t * e, int itok);
void ke_print(sml_t* sml, kexpr_t *kexpr);
void ke_fill_hash(sml_t *sml);
token_t* ke_get_val_index(sml_t *sml, int i);
void ke_set_val_index(sml_t *sml, int i, token_t *tokp);
void ke_set_val(sml_t* sml, token_t* e, token_t *q);
void ke_push_stack(sml_t * sml, token_t * tokp, int *top);
void ke_free_dll(sml_t * sml);
void ke_destroy(sml_t *sml, kexpr_t *kexpr);
void ke_free_hash(sml_t *sml);
void ke_free(sml_t* sml, kexpr_t *kexpr);
int peekfor(sml_t* sml);
int popfor(sml_t* sml);
void pushfor(sml_t* sml, int val);

#endif

ext: h

max_stars_repo_name: vinej/sml

max_issues_repo_name: vinej/sml

max_forks_repo_name: vinej/sml


text:
/* -*- linux-c -*- */
/* binsingle.c

   Copyright (C) 2002-2004 John M. Fregeau

   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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <getopt.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_rng.h>
#include "fewbody.h"
#include "binsingle_simple.h"

/* calculate the units used */
int calc_units(fb_obj_t *single, fb_obj_t *binary, fb_units_t *units)
{
	units->v = sqrt(FB_CONST_G*(single->m + binary->m)/(single->m * binary->m) * \
			(binary->obj[0]->m * binary->obj[1]->m / binary->a));
	units->l = binary->a;
	units->t = units->l / units->v;
	units->m = units->l * fb_sqr(units->v) / FB_CONST_G;
	units->E = units->m * fb_sqr(units->v);

	return(0);
}

/* the main attraction */
int main(int argc, char *argv[])
{
	int i, n;
	unsigned long int seed;
	double m0, m10, m11, r0, r10, r11, a1, e1;
	double rtid, vinf, b, m1, M, mu, Ei, t;
	fb_hier_t hier;
	fb_input_t input;
	fb_ret_t retval;
	fb_units_t units;
	char string1[FB_MAX_STRING_LENGTH], string2[FB_MAX_STRING_LENGTH];
	gsl_rng *rng;
	const gsl_rng_type *rng_type=gsl_rng_mt19937;

	/* set parameters to default values */
	m0 = FB_M0;
	m10 = FB_M10;
	m11 = FB_M11;
	r0 = FB_R0;
	r10 = FB_R10;
	r11 = FB_R11;
	a1 = FB_A1;
	e1 = FB_E1;
	vinf = FB_VINF;
	b = FB_B;
	input.ks = FB_KS;
	input.tstop = FB_TSTOP;
	input.Dflag = 0;
	input.dt = FB_DT;
	input.tcpustop = FB_TCPUSTOP;
	input.absacc = FB_ABSACC;
	input.relacc = FB_RELACC;
	input.ncount = FB_NCOUNT;
	input.tidaltol = FB_TIDALTOL;
	input.fexp = FB_FEXP;
	seed = FB_SEED;
	fb_debug = FB_DEBUG;

	/* initialize a few things for integrator */
	t = 0.0;
	fbui_new_hier(&hier, 3);

	/* initialize GSL rng */
	gsl_rng_env_setup();
	rng = gsl_rng_alloc(rng_type);
	gsl_rng_set(rng, seed);

	/* create binary */
	fbui_make_pair(fbui_hierarchy_binary(&hier, 0),
		       fbui_hierarchy_single(&hier, 1), fbui_hierarchy_single(&hier, 2));
	fbui_obj_t_set(fbui_hierarchy_binary(&hier, 0), t);

	/* give the objects some properties */
	fbui_initialize_single(fbui_hierarchy_single(&hier, 0), 0, "0");
	fbui_initialize_single(fbui_hierarchy_single(&hier, 1), 1, "1");
	fbui_initialize_single(fbui_hierarchy_single(&hier, 2), 2, "2");

	fbui_obj_radius_set(fbui_hierarchy_single(&hier, 0), r0);
	fbui_obj_radius_set(fbui_hierarchy_single(&hier, 1), r10);
	fbui_obj_radius_set(fbui_hierarchy_single(&hier, 2), r11);

	fbui_obj_mass_set(fbui_hierarchy_single(&hier, 0), m0);
	fbui_obj_mass_set(fbui_hierarchy_single(&hier, 1), m10);
	fbui_obj_mass_set(fbui_hierarchy_single(&hier, 2), m11);

	fbui_obj_mass_set(fbui_hierarchy_binary(&hier, 0), m10+m11);

	fbui_obj_a_set(fbui_hierarchy_binary(&hier, 0), a1);
	fbui_obj_e_set(fbui_hierarchy_binary(&hier, 0), e1);

	/* get the units and normalize */
	calc_units(fbui_hierarchy_single(&hier, 0), fbui_hierarchy_binary(&hier, 0), &units);
	fb_normalize(&hier, units);

	/* move hierarchies analytically in from infinity along hyperbolic orbit */
	m0 = fbui_obj_mass_get(fbui_hierarchy_single(&hier, 0));
	m1 = fbui_obj_mass_get(fbui_hierarchy_binary(&hier, 0));
	M = m0 + m1;
	mu = m0 * m1 / M;

	Ei = 0.5 * mu * fb_sqr(vinf);

	a1 = fbui_obj_a_get(fbui_hierarchy_binary(&hier, 0));
	e1 = fbui_obj_e_get(fbui_hierarchy_binary(&hier, 0));
	m10 = fbui_obj_mass_get(fbui_obj_left_child_get(fbui_hierarchy_binary(&hier, 0)));
	m11 = fbui_obj_mass_get(fbui_obj_right_child_get(fbui_hierarchy_binary(&hier, 0)));

	rtid = pow(2.0*(m0+m1)/(m1*input.tidaltol), 1.0/3.0) * a1 * (1.0+e1);

	fb_init_scattering(fbui_hierarchy_single(&hier, 0), fbui_hierarchy_binary(&hier, 0), vinf, b, rtid);

	/* trickle down the binary properties, then back up */
	fb_randorient(fbui_hierarchy_binary(&hier, 0), rng);
	fb_downsync(fbui_hierarchy_binary(&hier, 0), t);

	/* call fewbody! */
	retval = fewbody(input, &hier, &t);

	/* print information to screen */
	fprintf(stderr, "OUTCOME:\n");
	if (retval.retval == 1) {
		fprintf(stderr, "  encounter complete:  t=%.6g (%.6g yr)  %s  (%s)\n\n",
			t, t * units.t/FB_CONST_YR,
			fb_sprint_hier(hier, string1),
			fb_sprint_hier_hr(hier, string2));
	} else {
		fprintf(stderr, "  encounter NOT complete:  t=%.6g (%.6g yr)  %s  (%s)\n\n",
			t, t * units.t/FB_CONST_YR,
			fb_sprint_hier(hier, string1),
			fb_sprint_hier_hr(hier, string2));
	}

	fb_dprintf("there were %ld integration steps\n", retval.count);
	fb_dprintf("fb_classify() was called %ld times\n", retval.iclassify);

	fprintf(stderr, "FINAL:\n");
	fprintf(stderr, "  t_final=%.6g (%.6g yr)  t_cpu=%.6g s\n", \
		t, t*units.t/FB_CONST_YR, retval.tcpu);

	fprintf(stderr, "  DeltaL/L0=%.6g  DeltaL=%.6g\n", retval.DeltaLfrac, retval.DeltaL);
	fprintf(stderr, "  DeltaE/E0=%.6g  DeltaE=%.6g\n", retval.DeltaEfrac, retval.DeltaE);
	fprintf(stderr, "  Rmin=%.6g (%.6g RSUN)  Rmin_i=%d  Rmin_j=%d\n", \
		retval.Rmin, retval.Rmin*units.l/FB_CONST_RSUN, retval.Rmin_i, retval.Rmin_j);
	fprintf(stderr, "  Nosc=%d (%s)\n", retval.Nosc, (retval.Nosc>=1?"resonance":"non-resonance"));

	fprintf(stderr, "Resulting objects:\n");
	for (i=0; i<hier.nobj; i++) {
		n = fbui_obj_n_get(fbui_tree(&hier, i));
		if (n == 1) {
			fprintf(stderr, "  single: m=%g idstring=%s\n",
				fbui_obj_mass_get(fbui_tree(&hier, i)) * units.m / FB_CONST_MSUN,
				fbui_obj_idstring_get(fbui_tree(&hier, i)));
		} else if (n == 2) {
			fprintf(stderr, "  binary: m_i=%g,%g a=%g e=%g idstring=%s\n",
				fbui_obj_mass_get(fbui_obj_left_child_get(fbui_tree(&hier, i))) * units.m / FB_CONST_MSUN,
				fbui_obj_mass_get(fbui_obj_right_child_get(fbui_tree(&hier, i))) * units.m / FB_CONST_MSUN,
				fbui_obj_a_get(fbui_tree(&hier, i)) * units.l / FB_CONST_AU,
				fbui_obj_e_get(fbui_tree(&hier, i)),
				fbui_obj_idstring_get(fbui_tree(&hier, i)));
		}
	}

	/* free GSL stuff */
	gsl_rng_free(rng);

	/* free our own stuff */
	fbui_delete_hier(&hier);

	/* done! */
	return(0);
}

max_stars_repo_name: gnodvi/cosmos

max_issues_repo_name: gnodvi/cosmos

max_forks_repo_name: gnodvi/cosmos


text:
/* matrix/gsl_matrix_uchar.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, 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 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.
 */

#ifndef __GSL_MATRIX_UCHAR_H__
#define __GSL_MATRIX_UCHAR_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_uchar.h>
#include <gsl/gsl_blas_types.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  unsigned char * data;
  gsl_block_uchar * block;
  int owner;
} gsl_matrix_uchar;

typedef struct
{
  gsl_matrix_uchar matrix;
} _gsl_matrix_uchar_view;

typedef _gsl_matrix_uchar_view gsl_matrix_uchar_view;

typedef struct
{
  gsl_matrix_uchar matrix;
} _gsl_matrix_uchar_const_view;

typedef const _gsl_matrix_uchar_const_view gsl_matrix_uchar_const_view;

/* Allocation */

GSL_FUN gsl_matrix_uchar *
gsl_matrix_uchar_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_uchar *
gsl_matrix_uchar_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_uchar *
gsl_matrix_uchar_alloc_from_block (gsl_block_uchar * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_FUN gsl_matrix_uchar *
gsl_matrix_uchar_alloc_from_matrix (gsl_matrix_uchar * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN gsl_vector_uchar *
gsl_vector_uchar_alloc_row_from_matrix (gsl_matrix_uchar * m,
                                        const size_t i);

GSL_FUN gsl_vector_uchar *
gsl_vector_uchar_alloc_col_from_matrix (gsl_matrix_uchar * m,
                                        const size_t j);

GSL_FUN void gsl_matrix_uchar_free (gsl_matrix_uchar * m);

/* Views */

GSL_FUN _gsl_matrix_uchar_view
gsl_matrix_uchar_submatrix (gsl_matrix_uchar * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_uchar_view
gsl_matrix_uchar_row (gsl_matrix_uchar * m, const size_t i);

GSL_FUN _gsl_vector_uchar_view
gsl_matrix_uchar_column (gsl_matrix_uchar * m, const size_t j);

GSL_FUN _gsl_vector_uchar_view
gsl_matrix_uchar_diagonal (gsl_matrix_uchar * m);

GSL_FUN _gsl_vector_uchar_view
gsl_matrix_uchar_subdiagonal (gsl_matrix_uchar * m, const size_t k);

GSL_FUN _gsl_vector_uchar_view
gsl_matrix_uchar_superdiagonal (gsl_matrix_uchar * m, const size_t k);

GSL_FUN _gsl_vector_uchar_view
gsl_matrix_uchar_subrow (gsl_matrix_uchar * m, const size_t i,
                         const size_t offset, const size_t n);

GSL_FUN _gsl_vector_uchar_view
gsl_matrix_uchar_subcolumn (gsl_matrix_uchar * m, const size_t j,
                            const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_uchar_view
gsl_matrix_uchar_view_array (unsigned char * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_uchar_view
gsl_matrix_uchar_view_array_with_tda (unsigned char * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_uchar_view
gsl_matrix_uchar_view_vector (gsl_vector_uchar * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_uchar_view
gsl_matrix_uchar_view_vector_with_tda (gsl_vector_uchar * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_submatrix (const gsl_matrix_uchar * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_uchar_const_view
gsl_matrix_uchar_const_row (const gsl_matrix_uchar * m,
                            const size_t i);

GSL_FUN _gsl_vector_uchar_const_view
gsl_matrix_uchar_const_column (const gsl_matrix_uchar * m,
                               const size_t j);

GSL_FUN _gsl_vector_uchar_const_view
gsl_matrix_uchar_const_diagonal (const gsl_matrix_uchar * m);

GSL_FUN _gsl_vector_uchar_const_view
gsl_matrix_uchar_const_subdiagonal (const gsl_matrix_uchar * m,
                                    const size_t k);

GSL_FUN _gsl_vector_uchar_const_view
gsl_matrix_uchar_const_superdiagonal (const gsl_matrix_uchar * m,
                                      const size_t k);

GSL_FUN _gsl_vector_uchar_const_view
gsl_matrix_uchar_const_subrow (const gsl_matrix_uchar * m, const size_t i,
                               const size_t offset, const size_t n);

GSL_FUN _gsl_vector_uchar_const_view
gsl_matrix_uchar_const_subcolumn (const gsl_matrix_uchar * m, const size_t j,
                                  const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_array (const unsigned char * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_array_with_tda (const unsigned char * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_vector (const gsl_vector_uchar * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_vector_with_tda (const gsl_vector_uchar * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_uchar_set_zero (gsl_matrix_uchar * m);
GSL_FUN void gsl_matrix_uchar_set_identity (gsl_matrix_uchar * m);
GSL_FUN void gsl_matrix_uchar_set_all (gsl_matrix_uchar * m, unsigned char x);

GSL_FUN int gsl_matrix_uchar_fread (FILE * stream, gsl_matrix_uchar * m) ;
GSL_FUN int gsl_matrix_uchar_fwrite (FILE * stream, const gsl_matrix_uchar * m) ;
GSL_FUN int gsl_matrix_uchar_fscanf (FILE * stream, gsl_matrix_uchar * m);
GSL_FUN int gsl_matrix_uchar_fprintf (FILE * stream, const gsl_matrix_uchar * m, const char * format);

GSL_FUN int gsl_matrix_uchar_memcpy(gsl_matrix_uchar * dest, const gsl_matrix_uchar * src);
GSL_FUN int gsl_matrix_uchar_swap(gsl_matrix_uchar * m1, gsl_matrix_uchar * m2);
GSL_FUN int gsl_matrix_uchar_tricpy(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_uchar * dest, const gsl_matrix_uchar * src);

GSL_FUN int gsl_matrix_uchar_swap_rows(gsl_matrix_uchar * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_uchar_swap_columns(gsl_matrix_uchar * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_uchar_swap_rowcol(gsl_matrix_uchar * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_uchar_transpose (gsl_matrix_uchar * m);
GSL_FUN int gsl_matrix_uchar_transpose_memcpy (gsl_matrix_uchar * dest, const gsl_matrix_uchar * src);
GSL_FUN int gsl_matrix_uchar_transpose_tricpy (CBLAS_UPLO_t Uplo_src, CBLAS_DIAG_t Diag, gsl_matrix_uchar * dest, const gsl_matrix_uchar * src);

GSL_FUN unsigned char gsl_matrix_uchar_max (const gsl_matrix_uchar * m);
GSL_FUN unsigned char gsl_matrix_uchar_min (const gsl_matrix_uchar * m);
GSL_FUN void gsl_matrix_uchar_minmax (const gsl_matrix_uchar * m, unsigned char * min_out, unsigned char * max_out);

GSL_FUN void gsl_matrix_uchar_max_index (const gsl_matrix_uchar * m, size_t * imax, size_t *jmax);
GSL_FUN void gsl_matrix_uchar_min_index (const gsl_matrix_uchar * m, size_t * imin, size_t *jmin);
GSL_FUN void gsl_matrix_uchar_minmax_index (const gsl_matrix_uchar * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_FUN int gsl_matrix_uchar_equal (const gsl_matrix_uchar * a, const gsl_matrix_uchar * b);

GSL_FUN int gsl_matrix_uchar_isnull (const gsl_matrix_uchar * m);
GSL_FUN int gsl_matrix_uchar_ispos (const gsl_matrix_uchar * m);
GSL_FUN int gsl_matrix_uchar_isneg (const gsl_matrix_uchar * m);
GSL_FUN int gsl_matrix_uchar_isnonneg (const gsl_matrix_uchar * m);

GSL_FUN unsigned char gsl_matrix_uchar_norm1 (const gsl_matrix_uchar * m);

GSL_FUN int gsl_matrix_uchar_add (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
GSL_FUN int gsl_matrix_uchar_sub (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
GSL_FUN int gsl_matrix_uchar_mul_elements (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
GSL_FUN int gsl_matrix_uchar_div_elements (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
GSL_FUN int gsl_matrix_uchar_scale (gsl_matrix_uchar * a, const double x);
GSL_FUN int gsl_matrix_uchar_scale_rows (gsl_matrix_uchar * a, const gsl_vector_uchar * x);
GSL_FUN int gsl_matrix_uchar_scale_columns (gsl_matrix_uchar * a, const gsl_vector_uchar * x);
GSL_FUN int gsl_matrix_uchar_add_constant (gsl_matrix_uchar * a, const double x);
GSL_FUN int gsl_matrix_uchar_add_diagonal (gsl_matrix_uchar * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_uchar_get_row(gsl_vector_uchar * v, const gsl_matrix_uchar * m, const size_t i);
GSL_FUN int gsl_matrix_uchar_get_col(gsl_vector_uchar * v, const gsl_matrix_uchar * m, const size_t j);
GSL_FUN int gsl_matrix_uchar_set_row(gsl_matrix_uchar * m, const size_t i, const gsl_vector_uchar * v);
GSL_FUN int gsl_matrix_uchar_set_col(gsl_matrix_uchar * m, const size_t j, const gsl_vector_uchar * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL unsigned char   gsl_matrix_uchar_get(const gsl_matrix_uchar * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void    gsl_matrix_uchar_set(gsl_matrix_uchar * m, const size_t i, const size_t j, const unsigned char x);
GSL_FUN INLINE_DECL unsigned char * gsl_matrix_uchar_ptr(gsl_matrix_uchar * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const unsigned char * gsl_matrix_uchar_const_ptr(const gsl_matrix_uchar * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE
INLINE_FUN
unsigned char
gsl_matrix_uchar_get(const gsl_matrix_uchar * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
        }
    }
#endif
  return m->data[i * m->tda + j] ;
}

INLINE_FUN
void
gsl_matrix_uchar_set(gsl_matrix_uchar * m, const size_t i, const size_t j, const unsigned char x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  m->data[i * m->tda + j] = x ;
}

INLINE_FUN
unsigned char *
gsl_matrix_uchar_ptr(gsl_matrix_uchar * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (unsigned char *) (m->data + (i * m->tda + j)) ;
}

INLINE_FUN
const unsigned char *
gsl_matrix_uchar_const_ptr(const gsl_matrix_uchar * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const unsigned char *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_UCHAR_H__ */

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#ifndef __GSL_SPMATRIX_H__
#define __GSL_SPMATRIX_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

enum
{
  GSL_SPMATRIX_COO = 0, /* coordinate/triplet representation */
  GSL_SPMATRIX_CSC = 1, /* compressed sparse column */
  GSL_SPMATRIX_CSR = 2, /* compressed sparse row */
  GSL_SPMATRIX_TRIPLET = GSL_SPMATRIX_COO,
  GSL_SPMATRIX_CCS = GSL_SPMATRIX_CSC,
  GSL_SPMATRIX_CRS = GSL_SPMATRIX_CSR
};

/* memory pool for binary tree node allocation */
struct gsl_spmatrix_pool_node
{
  struct gsl_spmatrix_pool_node * next;
  void * block_ptr;          /* pointer to memory block, of size n*tree_node_size */
  unsigned char * free_slot; /* pointer to next available slot */
};

typedef struct gsl_spmatrix_pool_node gsl_spmatrix_pool;

#define GSL_SPMATRIX_ISCOO(m)         ((m)->sptype == GSL_SPMATRIX_COO)
#define GSL_SPMATRIX_ISCSC(m)         ((m)->sptype == GSL_SPMATRIX_CSC)
#define GSL_SPMATRIX_ISCSR(m)         ((m)->sptype == GSL_SPMATRIX_CSR)

#define GSL_SPMATRIX_ISTRIPLET(m)     GSL_SPMATRIX_ISCOO(m)
#define GSL_SPMATRIX_ISCCS(m)         GSL_SPMATRIX_ISCSC(m)
#define GSL_SPMATRIX_ISCRS(m)         GSL_SPMATRIX_ISCSR(m)

#define GSL_SPMATRIX_FLG_GROW         (1 << 0) /* allow size of matrix to grow as elements are added */
#define GSL_SPMATRIX_FLG_FIXED        (1 << 1) /* sparsity pattern is fixed */

/* compare matrix entries (ia,ja) and (ib,jb) - sort by rows first, then by columns */
#define GSL_SPMATRIX_COMPARE_ROWCOL(m,ia,ja,ib,jb)   ((ia) < (ib) ? -1 : ((ia) > (ib) ? 1 : ((ja) < (jb) ? -1 : ((ja) > (jb)))))

/* common/utility functions */

GSL_FUN void gsl_spmatrix_cumsum(const size_t n, int * c);

#include <gsl/gsl_spmatrix_complex_long_double.h>
#include <gsl/gsl_spmatrix_complex_double.h>
#include <gsl/gsl_spmatrix_complex_float.h>

#include <gsl/gsl_spmatrix_long_double.h>
#include <gsl/gsl_spmatrix_double.h>
#include <gsl/gsl_spmatrix_float.h>

#include <gsl/gsl_spmatrix_ulong.h>
#include <gsl/gsl_spmatrix_long.h>

#include <gsl/gsl_spmatrix_uint.h>
#include <gsl/gsl_spmatrix_int.h>

#include <gsl/gsl_spmatrix_ushort.h>
#include <gsl/gsl_spmatrix_short.h>

#include <gsl/gsl_spmatrix_uchar.h>
#include <gsl/gsl_spmatrix_char.h>

#endif /* __GSL_SPMATRIX_H__ */


text:
/* $Id$ */
/*--------------------------------------------------------------------*/
/*;  Copyright (C) 2008                                               */
/*;  Associated Universities, Inc. Washington DC, USA.                */
/*;                                                                   */
/*;  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 Massachusetts Ave, Cambridge,     */
/*;  MA 02139, USA.                                                   */
/*;                                                                   */
/*; Correspondence about this software should be addressed as follows:*/
/*;         Internet email: bcotton@nrao.edu.                         */
/*;         Postal address: William Cotton                            */
/*;                         National Radio Astronomy Observatory      */
/*;                         520 Edgemont Road                         */
/*;                         Charlottesville, VA 22903-2475 USA        */
/*--------------------------------------------------------------------*/

#include "ObitUtil.h"
#include <gsl/gsl_multifit.h>

/*----------------Obit: Merx mollis mortibus nuper ------------------*/
/**
 * \file ObitUtil.c
 * ObitUtil utility function definitions.
 *
 */

/*----------------------Private functions---------------------------*/
/** Private: qsort ofloat comparison */
static int compare_gfloat  (const void* arg1,  const void* arg2);

/*----------------------Public functions---------------------------*/
/**
 * Get mean values of an array.
 * \param array array of values may be magic value blanked
 * \param n     dimension of array
 * \return      Mean value, possibly blanked
 */
ofloat meanValue (ofloat *array, olong incs, olong n)
{
  olong i, count;
  ofloat sum, fblank =ObitMagicF() ;

  if (n<=0) return fblank;
  sum = 0.0;
  count = 0;
    for (i=0; i<n; i++) {
      if (array[i*incs]!=fblank) {
	sum += array[i*incs];
	count++;
      }
    }
    if (count>0) sum /= count;
    else sum = fblank;

  return sum;
} /* end meanValue */

/**
 * Get median value of an array.
 * Does sort then returns value array[ngood/2]
 * \param array   array of values, on return will be in ascending order
 * \param n       dimension of array
 * \param inc     stride in array
 * \return        Median value, possibly blanked
 */
ofloat medianValue (ofloat *array, olong incs, olong n)
{
  olong i, ngood;
  ofloat out, fblank = ObitMagicF();
  gboolean blanked;

  if (n<=0) return fblank;

  /* Count good (non blank) points */
  ngood = 0;
  for (i=0; i<n; i++) if (array[i*incs]!=fblank) ngood++;
  if (ngood<1) return fblank;  /* Any good data? */
  blanked = ngood<n;   /* Anything blanked? */

  /* set blanked values to 2.0e20 for sort */
  for (i=0; i<n; i++) if (array[i*incs]==fblank) array[i*incs] = 2.0e20;

  /* Sort to ascending order */
  qsort ((void*)array, n, MAX(1,incs)*sizeof(ofloat), compare_gfloat);

  /* Set median */
  out = array[(ngood/2)*incs];

  /* reset blanked values */
  if (blanked) {
    for (i=0; i<n; i++) if (array[i*incs]>1.0e20) array[i*incs] = fblank;
  }

  return out;
} /* end medianValue */

/**
 * Return average of navg values around median (with magic value blanking)
 * Does sort then returns average of center of array
 * If there are fewer than navg valid data then the median is returned.
 * \param array array of values, on return will be in ascending order
 * \param incs  increment in data array, >1 => weights
 * \param navg  width of averaging about median
 * \param doWt  If True value after data is a weight
 * \param n     dimension of array
 * \return      Median/average, possibly blanked
 */
ofloat medianAvg (ofloat *array, olong incs, olong navg, gboolean doWt, olong n)
{
  olong i, cnt, wid, ngood, ind, hi, lo;
  ofloat temp, wt, out=0.0, fblank = ObitMagicF();
  ofloat center;
  gboolean blanked;

  if (n<=0) return fblank;

  /* Count good (non blank) points */
  ngood = 0;
  for (i=0; i<n; i++) if (array[i*incs]!=fblank) ngood++;
  if (ngood<1) return fblank;  /* Any good data? */
  blanked = ngood<n;   /* Anything blanked? */

  /* set blanked values to 2.0e20 for sort */
  for (i=0; i<n; i++) if (array[i*incs]==fblank) array[i*incs] = 2.0e20;

  /* Sort to ascending order */
  qsort ((void*)array, n, MAX(1,incs)*sizeof(ofloat), compare_gfloat);

  /* How many good values? */
  cnt = 0;
  for (i=0; i<n; i++) if (array[i*incs]<1.0e20) cnt++;
  wid = MIN (MAX (1, navg/2), ngood/2);  /* width of averaging */
  ngood = cnt;

  ind = (ngood-1)/2;  /* center good value */
  out = array[ind*incs];

  /* average central values */
  if (ngood>navg) {
    center = (ngood-1)/2.0;
    lo = MAX (1, (olong)(center - wid + 0.6));
    hi = MIN (ngood, (olong)(center + wid));
    /* Weighted? */
    if (doWt) {
      temp = 0.0; wt = 0.0;
      for (i=lo; i<=hi; i++) {
	if (array[i*incs]<1.0e20) {
	  wt += array[i*incs+1];
	  temp += array[i*incs]*array[i*incs+1];
	}
      }
      if (wt>0.0) out = temp/wt;
    } else { /* unweighted */
      temp = 0.0; cnt = 0;
      for (i=lo; i<=hi; i++) {
	if (array[i]<1.0e20) {
	  cnt++;
	  temp += array[i];
	}
	if (cnt>0) out = temp/cnt;
      }
    } /* end if weighting */
  }

  /* reset blanked values */
  if (blanked) {
    for (i=0; i<n; i++) if (array[i*incs]>1.0e20) array[i*incs] = fblank;
  }

  return out;
} /* end medianAvg */

/**
 * Return the running median of an array
 * \param n       Number of points
 * \param wind    Width of median window in cells
 * \param array   Array of values, fblank blanking supported
 * \param alpha   0 -> 1 = pure boxcar -> pure MWF (Alpha of the
 *                data samples in a window are discarded and
 *                the rest averaged).
 * \param rms     RMS of array, median average sigma from each wind of data
 * \param out     array of size of array to be with median values
 * \param work    work array of size of array
 */
void RunningMedian (olong n, olong wind, ofloat *array, ofloat alpha,
		    ofloat *RMS, ofloat *out, ofloat *work)
{
  ofloat *lwork=NULL;
  ofloat level, sigma, sigmaSum, sigmaCnt;
  ofloat fblank = ObitMagicF();
  olong i, j, k, op, ind, half, RMScnt=0;

  /* Create array */
  lwork = g_malloc0(wind*sizeof(ofloat));

  half = wind/2;
  ind  = 0;
  op   = 0;
  k    = 0;
  sigmaSum = 0.0;
  sigmaCnt = 1;

  /* First half wind filled with median of first wind points */
  for (j=ind; j<ind+wind; j++) lwork[k++] = array[j];
  ind++;
  level = MedianLevel (wind, lwork, alpha);
  sigma = MedianSigma (wind, lwork, level);
  if (sigma!=fblank) {
    sigmaSum = sigma;
    sigmaSum += sigma;
    work[RMScnt++] = sigma;
  }
  for (k=0; k<half; k++) out[op++] = level;

  /* Loop over middle of array */
  for (i=half; i<n-half; i++) {
    k = 0;
    for (j=ind; j<ind+wind; j++) lwork[k++] = array[j];
    ind++;
    level = MedianLevel (wind, lwork, alpha);
    sigma = MedianSigma (wind, lwork, level);
    if (sigma!=fblank) {
      sigmaSum += sigma;
      sigmaCnt++;
      work[RMScnt++] = sigma;
    }
    out[op++] = level;
  } /* end loop over array */

  /* Fill in bit at end */
  while (op<n) {
    out[op++] = level;
  }

  /* median Average sigmas */
  *RMS = MedianLevel (RMScnt, work, alpha);

  /* Cleanup */
  if (lwork) g_free(lwork);

} /* end RunningMedian */

/**
 * Determine alpha median value of a ofloat array
 * \param n       Number of points
 * \param value   Array of values
 * \param alpha   0 -> 1 = pure boxcar -> pure MWF (ALPHA of the
 *                data samples are discarded and the rest averaged).
 * \return alpha median value
 */
ofloat MedianLevel (olong n, ofloat *value, ofloat alpha)
{
  ofloat out=0.0;
  ofloat fblank = ObitMagicF();
  ofloat beta, sum;
  olong i, i1, i2, count;

  if (n<=0) return out;

  /* Sort to ascending order */
  qsort ((void*)value, n, sizeof(ofloat), compare_gfloat);

  out = value[n/2];

  beta = MAX (0.05, MIN (0.95, alpha)) / 2.0; /*  Average around median factor */

  /* Average around the center */
  i1 = MAX (0, (n/2)-(olong)(beta*n+0.5));
  i2 = MIN (n, (n/2)+(olong)(beta*n+0.5));

  if (i2>i1) {
    sum = 0.0;
    count = 0;
    for (i=i1; i<i2; i++) {
      if (value[i]!=fblank) {
	sum += value[i];
	count++;
      }
    }
    if (count>0) out = sum / count;
  }

  return out;
} /* end MedianLevel */

/**
 * Determine robust RMS value of a ofloat array about mean
 * Use center 90% of points, excluding at least one point from each end
 * \param n       Number of points, needs at least 4
 * \param value   Array of values assumed sorted
 * \param mean    Mean value of value
 * \return RMS value, fblank if cannot determine
 */
ofloat MedianSigma (olong n, ofloat *value, ofloat mean)
{
  ofloat fblank = ObitMagicF();
  ofloat out;
  ofloat sum;
  olong i, i1, i2, count;

  out = fblank;
  if (n<=4) return out;
  if (mean==fblank) return out;

  /* Get RMS around the center 90% */
  i1 = MAX (1,   (n/2)-(olong)(0.45*n+0.5));
  i2 = MIN (n-1, (n/2)+(olong)(0.45*n+0.5));

  if (i2>i1) {
    sum = 0.0;
    count = 0;
    for (i=i1; i<i2; i++) {
      if (value[i]!=fblank) {
	sum += (value[i]-mean)*(value[i]-mean);
	count++;
      }
    }
    if (count>1) out = sqrt(sum / (count-1));
  }

  return out;
} /* end MedianSigma */

/**
 * Fit polynomial y = f(poly, x) with magic value blanking
 * Use gsl package.
 * \param poly   [out] polynomial coef in order of increasing power of x
 * \param order  order of the polynomial
 * \param x      values at which y is sampled
 * \param y      values to be fitted
 * \param wt     weights for values
 * \param n      number of (x,y) value pairs
 */
void  FitPoly (ofloat *poly, olong order, ofloat *x, ofloat *y, ofloat *wt,
	       olong n)
{
  olong i, j, k, good, p=order+1;
  ofloat fgood=0.0;
  double xi, chisq;
  gsl_matrix *X, *cov;
  gsl_vector *yy, *w, *c;
  gsl_multifit_linear_workspace *work;
  ofloat fblank = ObitMagicF();

  /* Only use good data */
  good = 0;
  for (i=0; i<n; i++) if (y[i]!=fblank) {good++; fgood = y[i];}

  /* If only one good datum - use it */
  if (good==1) {
    poly[0] = fgood;
    poly[1] = 0.0;
  }

  if (good<(order+1)) {  /* Not enough good data */
    poly[0] = fblank;
    return;
  }

  /* If only one and order=0 use it */
  if ((good==1) && (order==0)) {
    poly[0] = y[0];
    return;
  }

  /* allocate arrays */
  X    = gsl_matrix_alloc(good, p);
  yy   = gsl_vector_alloc(good);
  w    = gsl_vector_alloc(good);
  c    = gsl_vector_alloc(p);
  cov  = gsl_matrix_alloc(p, p);
  work = gsl_multifit_linear_alloc (good, p);

  /* set data */
  k = 0;
  for (i=0; i<n; i++) {
    if (y[i]!=fblank) {
      gsl_vector_set(yy, k, (double)y[i]);
      gsl_vector_set(w, k, (double)wt[i]);
      xi = 1.0;
      for (j=0; j<p; j++) {
	gsl_matrix_set(X, k, j, xi);
	xi *= x[i];
      }
      k++;  /* good ones */
    }
  }

  /* Fit */
  gsl_multifit_wlinear (X, w, yy, c, cov, &chisq, work);

  /* get results */
  for (j=0; j<p; j++) poly[j] = (ofloat)gsl_vector_get(c, j);

  /* Deallocate arrays */
  gsl_matrix_free(X);
  gsl_vector_free(yy);
  gsl_vector_free(w);
  gsl_vector_free(c);
  gsl_matrix_free(cov);
  gsl_multifit_linear_free (work);
} /* end FitPoly */

/**
 * Evaluate polynomial y = f(poly, x)
 * \param order  order of the polynomial
 * \param poly   polynomial coef in order of increasing power of x
 * \param x      value at which polynomial is to be evaluated
 * \return evaluated polynomial
 */
ofloat  EvalPoly (olong order, ofloat *poly, ofloat x)
{
  olong i;
  ofloat sum = 0.0, aarg = 1.0;

  for (i=0; i<=order; i++) {
    sum  += poly[i]*aarg;
    aarg *= x;
  }
  return sum;
} /* end EvalPoly  */

/*----------------------Private functions---------------------------*/
/**
 * ofloat comparison of two arguments
 * \param arg1 first value to compare
 * \param arg2 second value to compare
 * \return negative if arg1 is less than arg2, zero if equal
 *  and positive if arg1 is greater than arg2.
 */
static int compare_gfloat  (const void* arg1,  const void* arg2)
{
  int out = 0;
  ofloat larg1, larg2;

  larg1 = *(ofloat*)arg1;
  larg2 = *(ofloat*)arg2;
  if (larg1<larg2) out = -1;
  else if (larg1>larg2) out = 1;
  return out;
} /* end compare_gfloat */

max_stars_repo_name: sarrvesh/Obit

max_issues_repo_name: sarrvesh/Obit

max_forks_repo_name: sarrvesh/Obit


text:
/*
 * Copyright (c) 2011-2018 The University of Tennessee and The University
 *                         of Tennessee Research Foundation.  All rights
 *                         reserved.
 * Copyright (c) 2013      Inria. All rights reserved.
 *
 * @precisions normal z -> c d s
 *
 */

#include "dplasma.h"
#include "dplasmatypes.h"
#include "dplasmaaux.h"
#include "parsec/data_dist/matrix/two_dim_rectangle_cyclic.h"
#include "parsec/data_dist/matrix/vector_two_dim_cyclic.h"

#include "zpltmg_chebvand.h"
#include "zpltmg_fiedler.h"
#include "zpltmg_hankel.h"
#include "zpltmg_toeppd.h"

#include <string.h>
#include <lapacke.h>

/**
 *******************************************************************************
 *
 *  Generic case
 *
 *******************************************************************************
 */
struct zpltmg_args_s {
    PLASMA_enum            mtxtype;
    unsigned long long int seed;
    parsec_complex64_t     *W;
};
typedef struct zpltmg_args_s zpltmg_args_t;

static int
dplasma_zpltmg_generic_operator( parsec_execution_stream_t *es,
                                 const parsec_tiled_matrix_dc_t *descA,
                                 void *_A,
                                 PLASMA_enum uplo, int m, int n,
                                 void *op_data )
{
    int tempmm, tempnn, ldam;
    zpltmg_args_t     *args = (zpltmg_args_t*)op_data;
    parsec_complex64_t *A    = (parsec_complex64_t*)_A;
    (void)es;
    (void)uplo;

    tempmm = (m == (descA->mt-1)) ? (descA->m - m * descA->mb) : descA->mb;
    tempnn = (n == (descA->nt-1)) ? (descA->n - n * descA->nb) : descA->nb;
    ldam   = BLKLDD( descA, m );

    if ( args->mtxtype == PlasmaMatrixCircul ) {
        return CORE_zpltmg_circul(
            tempmm, tempnn, A, ldam,
            descA->m, m*descA->mb, n*descA->nb, args->W );
    } else {
        return CORE_zpltmg(
            args->mtxtype, tempmm, tempnn, A, ldam,
            descA->m, descA->n, m*descA->mb, n*descA->nb, args->seed );
    }
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zpltmg_generic - Generic wrapper for cases that are based on the map
 * function. This is the default for many test matrices generation.
 *
 * See parsec_apply() for further information.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in] mtxtype
 *          Type of matrix to be generated.
 *
 * @param[in,out] A
 *          Descriptor of the distributed matrix A to generate. Any tiled matrix
 *          descriptor can be used.
 *
 * @param[out] W
 *          Workspace required by some generators.
 *
 * @param[in] seed
 *          The seed used in the random generation.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_cpltmg_genvect
 * @sa dplasma_dpltmg_genvect
 * @sa dplasma_spltmg_genvect
 *
 ******************************************************************************/
static inline int
dplasma_zpltmg_generic( parsec_context_t *parsec,
                        PLASMA_enum mtxtype,
                        parsec_tiled_matrix_dc_t *A,
                        parsec_complex64_t *W,
                        unsigned long long int seed)
{
    parsec_taskpool_t *parsec_zpltmg = NULL;
    zpltmg_args_t *params = (zpltmg_args_t*)malloc(sizeof(zpltmg_args_t));

    params->mtxtype = mtxtype;
    params->seed    = seed;
    params->W       = W;

    parsec_zpltmg = parsec_apply_New( PlasmaUpperLower, A, dplasma_zpltmg_generic_operator, params );
    if ( parsec_zpltmg != NULL )
    {
        parsec_context_add_taskpool(parsec, (parsec_taskpool_t*)parsec_zpltmg);
        dplasma_wait_until_completion(parsec);
        parsec_apply_Destruct( parsec_zpltmg );
        return 0;
    }
    return -101;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zpltmg_genvect - Generic wrapper for cases that are using two
 * datatypes: the default one, and one describing a vector.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in] mtxtype
 *          Type of matrix to be generated.
 *
 * @param[in,out] A
 *          Descriptor of the distributed matrix A to generate. Any tiled matrix
 *          descriptor can be used.
 *
 * @param[in] seed
 *          The seed used in the random generation.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_cpltmg_genvect
 * @sa dplasma_dpltmg_genvect
 * @sa dplasma_spltmg_genvect
 *
 ******************************************************************************/
static inline int
dplasma_zpltmg_genvect( parsec_context_t *parsec,
                        PLASMA_enum mtxtype,
                        parsec_tiled_matrix_dc_t *A,
                        unsigned long long int seed )
{
    size_t vectorsize = 0;
    parsec_taskpool_t* tp;

    switch( mtxtype ) {
    case PlasmaMatrixChebvand:
        tp = (parsec_taskpool_t*)parsec_zpltmg_chebvand_new( seed,
                                                            A );
        vectorsize = 2 * A->nb * sizeof(parsec_complex64_t);
        break;

    case PlasmaMatrixFiedler:
        tp = (parsec_taskpool_t*)parsec_zpltmg_fiedler_new( seed,
                                                            A );
        vectorsize = A->mb * sizeof(parsec_complex64_t);
        break;

    case PlasmaMatrixHankel:
        tp = (parsec_taskpool_t*)parsec_zpltmg_hankel_new( seed,
                                                           A );
        vectorsize = A->mb * sizeof(parsec_complex64_t);
        break;

    case PlasmaMatrixToeppd:
        tp = (parsec_taskpool_t*)parsec_zpltmg_toeppd_new( seed,
                                                           A );
        vectorsize = 2 * A->mb * sizeof(parsec_complex64_t);
        break;

    default:
        return -2;
    }

    if (tp != NULL) {
        parsec_zpltmg_hankel_taskpool_t *zpltmg_tp = (parsec_zpltmg_hankel_taskpool_t*)tp;

        /* Default type */
        dplasma_add2arena_tile( zpltmg_tp->arenas[PARSEC_zpltmg_hankel_DEFAULT_ARENA],
                                A->mb*A->nb*sizeof(parsec_complex64_t),
                                PARSEC_ARENA_ALIGNMENT_SSE,
                                parsec_datatype_double_complex_t, A->mb );

        /* Vector type */
        dplasma_add2arena_tile( zpltmg_tp->arenas[PARSEC_zpltmg_hankel_VECTOR_ARENA],
                                vectorsize,
                                PARSEC_ARENA_ALIGNMENT_SSE,
                                parsec_datatype_double_complex_t, A->mb );

        parsec_context_add_taskpool(parsec, tp);
        dplasma_wait_until_completion(parsec);

        parsec_matrix_del2arena( zpltmg_tp->arenas[PARSEC_zpltmg_hankel_DEFAULT_ARENA] );
        parsec_matrix_del2arena( zpltmg_tp->arenas[PARSEC_zpltmg_hankel_VECTOR_ARENA ] );
        parsec_taskpool_free(tp);
        return 0;
    }
    return -101;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zpltmg_circul - Generates a Circulant test matrix by tiles.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in,out] A
 *          Descriptor of the distributed matrix A to generate. Any tiled matrix
 *          descriptor can be used.
 *
 * @param[in] seed
 *          The seed used in the random generation.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_cpltmg_circul
 * @sa dplasma_dpltmg_circul
 * @sa dplasma_spltmg_circul
 *
 ******************************************************************************/
static inline int
dplasma_zpltmg_circul( parsec_context_t *parsec,
                       parsec_tiled_matrix_dc_t *A,
                       unsigned long long int seed )
{
    int info;
    parsec_complex64_t *V = (parsec_complex64_t*) malloc( A->m * sizeof(parsec_complex64_t) );

    CORE_zplrnt( A->m, 1, V, A->m, A->m, 0, 0, seed );

    info = dplasma_zpltmg_generic(parsec, PlasmaMatrixCircul, A, V, seed);

    free(V);
    return info;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_internal
 *
 * dplasma_zpltmg_condex - Generates a Condex test matrix by tiles.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in,out] A
 *          Descriptor of the distributed matrix A to generate. Any tiled matrix
 *          descriptor can be used.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_cpltmg_condex
 * @sa dplasma_dpltmg_condex
 * @sa dplasma_spltmg_condex
 *
 ******************************************************************************/
static inline int
dplasma_zpltmg_condex( parsec_context_t *parsec,
                       parsec_tiled_matrix_dc_t *A )
{
    /* gallery('condex', A->m, 4, 100.) */
    parsec_complex64_t theta = 100.;
    two_dim_block_cyclic_t *twodA = (two_dim_block_cyclic_t *)A;
    two_dim_block_cyclic_t Q;
    two_dim_block_cyclic_init( &Q, matrix_ComplexDouble, matrix_Tile,
                               1, A->super.myrank,
                               A->mb, A->nb, A->mb*A->mt, 3, 0, 0, A->m, 3, twodA->grid.strows, twodA->grid.stcols, 1 );
    Q.mat = parsec_data_allocate((size_t)Q.super.nb_local_tiles *
                                (size_t)Q.super.bsiz *
                                (size_t)parsec_datadist_getsizeoftype(Q.super.mtype));
    parsec_data_collection_set_key((parsec_data_collection_t*)&Q, "Q");

    if (A->super.myrank == 0) {
        parsec_complex64_t *Qmat;

        Qmat = (parsec_complex64_t*)(Q.mat);

        /* Initialize the Q matrix */
        CORE_zpltmg_condexq( A->m, A->n, Qmat, Q.super.lm );

        /*
         * Conversion to tile layout
         */
        {
            parsec_complex64_t *W = (parsec_complex64_t*) malloc (A->mb * sizeof(parsec_complex64_t) );
            int *leaders = NULL;
            int i, nleaders;

            /* Get all the cycles leaders and length
             * They are the same for each independent problem (each panel) */
            GKK_getLeaderNbr( Q.super.lmt, A->nb, &nleaders, &leaders );

            /* shift cycles. */
            for(i=0; i<nleaders; i++) {

                /* cycle #i belongs to this thread, so shift it */
                memcpy(W, Qmat + leaders[i*3] * A->mb, A->mb * sizeof(parsec_complex64_t) );
                CORE_zshiftw(leaders[i*3], leaders[i*3+1], A->mt, A->nb, A->mb, Qmat, W);
            }

            free(leaders); free(W);
        }
    }

    dplasma_zlaset( parsec, PlasmaUpperLower, 0., 1. + theta, A );
    dplasma_zgemm( parsec, PlasmaNoTrans, PlasmaConjTrans,
                   -theta, (parsec_tiled_matrix_dc_t*)&Q,
                           (parsec_tiled_matrix_dc_t*)&Q,
                   1.,     A );

    parsec_data_free(Q.mat);
    parsec_tiled_matrix_dc_destroy((parsec_tiled_matrix_dc_t*)&Q);
    return 0;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zpltmg_house - Generates a Householder test matrix by tiles.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in,out] A
 *          Descriptor of the distributed matrix A to generate. Any tiled matrix
 *          descriptor can be used.
 *
 * @param[in] seed
 *          The seed used in the random generation.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_cpltmg_house
 * @sa dplasma_dpltmg_house
 * @sa dplasma_spltmg_house
 *
 ******************************************************************************/
static inline int
dplasma_zpltmg_house( parsec_context_t *parsec,
                      parsec_tiled_matrix_dc_t *A,
                      unsigned long long int seed )
{
    /* gallery('house', random, 0 ) */
    vector_two_dim_cyclic_t V;
    parsec_complex64_t *Vmat, tau;

    vector_two_dim_cyclic_init( &V, matrix_ComplexDouble, PlasmaVectorDiag,
                                1, A->super.myrank,
                                A->mb, A->m, 0, A->m, 1 );
    V.mat = parsec_data_allocate((size_t)V.super.nb_local_tiles *
                                (size_t)V.super.bsiz *
                                (size_t)parsec_datadist_getsizeoftype(V.super.mtype));
    parsec_data_collection_set_key((parsec_data_collection_t*)&V, "V");
    Vmat = (parsec_complex64_t*)(V.mat);

    /* Initialize Householder vector */
    if (A->super.myrank == 0) {
        CORE_zplrnt( A->m, 1, Vmat, A->m, A->m, 0, 0, seed );
        LAPACKE_zlarfg_work( A->m, Vmat, Vmat+1, 1, &tau );
        Vmat[0] = 1.;
    }

#if defined(PARSEC_HAVE_MPI)
    /* If we don't need to broadcast, don't do it, this way we don't require MPI to be initialized */
    if( A->super.nodes > 1 )
        MPI_Bcast( &tau, 1, parsec_datatype_double_complex_t, 0, *(MPI_Comm*)dplasma_pcomm );
#endif

    /* Compute the Householder matrix I - tau v * v' */
    dplasma_zlaset( parsec, PlasmaUpperLower, 0., 1., A);
    dplasma_zgerc( parsec, -tau,
                   (parsec_tiled_matrix_dc_t*)&V,
                   (parsec_tiled_matrix_dc_t*)&V,
                   A );

    parsec_data_free(V.mat);
    parsec_tiled_matrix_dc_destroy((parsec_tiled_matrix_dc_t*)&V);

    return 0;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zpltmg - Generates a special test matrix by tiles.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in] mtxtype
 *          See PLASMA_zpltmg() for possible values and information on generated
 *          matrices.
 *
 * @param[in,out] A
 *          Descriptor of the distributed matrix A to generate. Any tiled matrix
 *          descriptor can be used.
 *          On exit, the matrix A generated.
 *
 * @param[in] seed
 *          The seed used in the random generation.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_cpltmg
 * @sa dplasma_dpltmg
 * @sa dplasma_spltmg
 *
 ******************************************************************************/
int
dplasma_zpltmg( parsec_context_t *parsec,
                PLASMA_enum mtxtype,
                parsec_tiled_matrix_dc_t *A,
                unsigned long long int seed)
{

    switch( mtxtype ) {
    case PlasmaMatrixCircul:
        return dplasma_zpltmg_circul(parsec, A, seed);
        break;

    case PlasmaMatrixChebvand:
        return dplasma_zpltmg_genvect(parsec, mtxtype,
                                      A, seed);
        break;

    case PlasmaMatrixCondex:
        return dplasma_zpltmg_condex(parsec, A);
        break;

    case PlasmaMatrixFiedler:
        return dplasma_zpltmg_genvect(parsec, mtxtype,
                                      A, seed);
        break;

    case PlasmaMatrixHankel:
        return dplasma_zpltmg_genvect(parsec, mtxtype,
                                      A, seed);
        break;

    case PlasmaMatrixHouse:
        return dplasma_zpltmg_house(parsec, A, seed);
        break;

    case PlasmaMatrixToeppd:
        return dplasma_zpltmg_genvect(parsec, mtxtype,
                                      A, seed);
        break;

    case PlasmaMatrixCauchy:
    case PlasmaMatrixCompan:
    case PlasmaMatrixDemmel:
    case PlasmaMatrixDorr:
    case PlasmaMatrixFoster:
    case PlasmaMatrixHadamard:
    case PlasmaMatrixHilb:
    case PlasmaMatrixInvhess:
    case PlasmaMatrixKms:
    case PlasmaMatrixLangou:
    case PlasmaMatrixLehmer:
    case PlasmaMatrixLotkin:
    case PlasmaMatrixMinij:
    case PlasmaMatrixMoler:
    case PlasmaMatrixOrthog:
    case PlasmaMatrixParter:
    case PlasmaMatrixRandom:
    case PlasmaMatrixRiemann:
    case PlasmaMatrixRis:
    case PlasmaMatrixWilkinson:
    case PlasmaMatrixWright:
        return dplasma_zpltmg_generic(parsec, mtxtype, A, NULL, seed);
        break;
    default:
        return -2;
    }

    return 0;
}

max_stars_repo_name: NLAFET/ABFT

max_issues_repo_name: NLAFET/ABFT

max_forks_repo_name: NLAFET/ABFT


text:
//#pragma once			/* guard */
#ifndef __UTILITIES_INCLUDED__
#define __UTILITIES_INCLUDED__
#define auxFloor(x) ((float)(long)(x))
#define auxCeil(x) ((float)(long)((x)+1))
/* fwd declarations */
/* includes */
#include <iostream>
#include <stdlib.h>
#include <sstream>
#include <fstream>
#include <math.h>
#include <vector>
#include <list>
//#include <unordered_map>
#include <map>
#include <set>
#include <algorithm>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
// declarations std lib
using namespace std;

//FUNCTION DECLARATIONS
void donowt(void);
unsigned int howmany_events(const gsl_rng *r, vector< double >& haz, double dt);
unsigned int howmany_events(const gsl_rng *r, double haz, double dt);
void gsl_ran_multinomial_vec (const gsl_rng * r,  unsigned int K,  unsigned int N, vector< double >& p, vector< unsigned int >& n);
void tridiagsolver( vector< double >& a, vector< double >& b, vector< double >& c, vector< double >& x, vector< double >& d );
void splinetable( vector< double >& X, vector< double >& Y, vector< double >& Z );
double splinefun(double x, vector< double >& X, vector< double >& Y, vector< double >& Z );
int bisect( double x, vector< double >& V, int lo, int hi);
unsigned int bisect_find (double x, vector< double >& V);
unsigned int spline_bisect_find( double x, vector< double >& V);
void find_events(const gsl_rng *r, list< unsigned int >& ES, unsigned int noeli, unsigned int noev, vector< double >& W );
int find_one_event(const gsl_rng *r, vector< double >& W );
int find_one_event(const gsl_rng *r, vector< int >& Wi );
int find_one_event(const gsl_rng *r, vector< unsigned int >& Wu );
void shuffle( const gsl_rng *r, vector< int >& w);
void MAsmooth( vector< double >& V, vector< double >& W, int L );
void MAsmooth( vector< double >& V,  int L );
int collapse_result_averager( vector< vector< double > >& ibmres, vector< vector< double > >& ibmres2 );
int collapse_result_averager( vector< vector< double > >& ibmres );

#endif

max_stars_repo_name: petedodd/ZAMSTAR

max_issues_repo_name: petedodd/ZAMSTAR

max_forks_repo_name: petedodd/ZAMSTAR


text:
#ifndef _H_LIMBER
#define _H_LIMBER

#include <gsl/gsl_spline.h>
#include <stdbool.h>
#include "interp2d.h"

// The integrated C_ell are in general allowed to be zero or negative if
// they describe cross-correlations. We use these statuses to describe
// errors where a log was requested also.
// LIMBER_STATUS_NEGATIVE is probably always an error
// but LIMBER_STATUS_ZERO is not necessarily
#define LIMBER_STATUS_OK 0
#define LIMBER_STATUS_ZERO 1
#define LIMBER_STATUS_NEGATIVE 2
#define LIMBER_STATUS_ERROR 3

// These are the options you can set for
// the Limber integrator.
typedef struct limber_config{
	bool xlog;  // The output spline will be in terms of log(ell) not ell
	bool ylog;  // The output spline will return log(C_ell) not C_ell
	int n_ell;  // Number of ell values you want in the spline
	double * ell;  // The chosen ell values you want
	double prefactor; //Scaling prefactor
    int status; // did everything go okay?
    double absolute_tolerance;
    double relative_tolerance;
} limber_config;

// Do a flat universe Limber approximation integral, of the form:
// C^{XY}(\ell)
// = A \int_0^{\chi_{\mathrm{max}}}
//       \frac{W^X(\chi) W^Y(\chi)}{\chi^2} P(\ell/\chi, \chi)

// The two splines and the matter power passed into the function
// be FUNCTIONS OF CHI NOT Z

// TODO - Fortran version

gsl_spline * limber_integral(limber_config * config,
	gsl_spline * WX, gsl_spline * WY, Interpolator2D * P);

#endif


text:
/*
 * Barcode
 * Copyright E.G.P. Bos and F.S. Kitaura
 *
 * Distributed under the terms of the MIT License.
 * The full license is in the file LICENSE, distributed with this software.
 */

#pragma once

// EGP: cleared out most includes, just left the two here:
#include <gsl/gsl_rng.h>  // For the gsl_rng type // EGP: added this one
#include "define_opt.h"  // For the ULONG and real_prec types
// #include "curses_funcs.h"
// EGP: toch maar weer terug gezet, gingen dingen mis
/* #include <math.h> */
/* #include <fstream> */
/* #include <iomanip> */
// #include <iostream>
/* #include <string>          */
/* #include <string.h> */
/* #include <cassert> */
/* #include <cfloat> */
/* #include <stdlib.h> */
/* #include <stdio.h> */
/* #include <fftw3.h> */
/* #include <omp.h> */
/* #include <sstream> */
/* #include <netinet/in.h> */

/* #include <gsl/gsl_errno.h> */
/* #include <gsl/gsl_spline.h> */
/* #include <gsl/gsl_randist.h> */
/* #include <gsl/gsl_integration.h> */
/* #include <gsl/gsl_sf.h> */

/* #include "fftw_array.h" */

/* #include "define_opt.h" */

#include <vector>
#include "HMC_models.h"
#include "struct_main.h"
#include <sstream>

using namespace std;

/* --- structures --- */

struct HAMIL_NUMERICAL {
  unsigned N1;  /* Number of cells per axis */
  unsigned N2;  /* Number of cells per axis */
  unsigned N3;  /* Number of cells per axis */
  ULONG N;
  ULONG Nhalf;

  int mk;

  real_prec L1;
  real_prec L2;
  real_prec L3;

  real_prec vol;

  real_prec d1;
  real_prec d2;
  real_prec d3;

  real_prec min1;
  real_prec min2;
  real_prec min3;

  // Redshift space distortions
  real_prec xobs;
  real_prec yobs;
  real_prec zobs;
  bool planepar;
  bool periodic;

  ULONG N_bin;

  ULONG iGibbs;
  ULONG rejections;
  ULONG total_steps_lim;

  // EGP: some extra parameters for statistics on performance:
  bool accepted;
  real_prec epsilon;
  ULONG Neps;
  real_prec dH, dK, dE, dprior, dlikeli, psi_prior, psi_likeli, psi_prior_i,
  psi_prior_f, psi_likeli_i, psi_likeli_f, H_kin_i, H_kin_f;

  int method;			/* which method */
  int filter;			/* which filter method */
  ULONG itmax;			/* maximal number of iterations */
  int i_reset;			/* number to reset */
  int INV_SUCCESS;		/* indicate inversion success */

  real_prec epsi;  /* precision */

  real_prec kth;

  // EGP: pseudo-timestep parameters
  real_prec N_eps_fac;
  real_prec eps_fac;
  real_prec eps_fac_target;
  real_prec eps_fac_power;
  ULONG s_eps_total;

  // mass type
  bool mass_fs;
  bool mass_rs;
  int mass_type;
  ULONG massnum_init;
  ULONG massnum_burn;

  // EGP: testing:
  real_prec mass_factor;

  /// TESTING ///
  real_prec grad_psi_prior_factor;
  real_prec grad_psi_likeli_factor;
  bool grad_psi_prior_conjugate;
  bool grad_psi_likeli_conjugate;
  bool grad_psi_prior_times_i;
  bool grad_psi_likeli_times_i;
  bool div_dH_by_N;
  int calc_h;
  real_prec deltaQ_factor;

  bool correct_delta;

  int particle_kernel;
  real_prec particle_kernel_h;

  // FFTW
  real_prec *in_r2c;
  real_prec *out_c2r;
  complex_prec *in_c2r;
  complex_prec *out_r2c;
  struct plan_pkg *R2Cplan;
  struct plan_pkg *C2Rplan;
};

struct HAMIL_DATA {
  struct HAMIL_NUMERICAL *numerical;  // numerical parameters

  // set arrays
  real_prec *mass_f;
  real_prec *mass_r;
  // real_prec *massi;
  real_prec *gradpsi;
  real_prec *x;
  // real_prec *growth;

  real_prec *posx;
  real_prec *posy;
  real_prec *posz;

  real_prec *signal_PS;
  real_prec *window;
  real_prec *func;
  real_prec *noise;
  real_prec *nobs;
  real_prec *deltaX;
  real_prec *corrf;

  // function pointers
  void (*partial_f_delta_x_log_like)(struct HAMIL_DATA*, real_prec*,
                                     real_prec*);
  real_prec (*log_like)(struct HAMIL_DATA*, real_prec*);
  void (*grad_f_delta_x_comp)(struct HAMIL_DATA *, real_prec *, real_prec *, unsigned int);

  real_prec (*log_prior)(struct HAMIL_DATA*, real_prec*);
  void (*grad_log_prior)(struct HAMIL_DATA*, real_prec*, real_prec*);

  // model flags
  int sfmodel;
  bool rsd_model;

  // set scalars
  real_prec Nmean_Gal;
  real_prec rho_c;
  real_prec code_norm;
  real_prec mu;

  real_prec sigma_fac;
  real_prec sigma_min;
  real_prec delta_min;

  real_prec biasP;
  real_prec biasE;

  real_prec z;
  real_prec ascale;

  real_prec D1;
  real_prec D2;

  real_prec OM;
  real_prec OL;
  real_prec h;

  // external inputs
  real_prec *M;

  // curses
  // struct CURSES_STRUCT *curses;

  // stuff for SPH_kernel_3D_cells
  const int N_cells;
  const std::vector<int> kernel_cells_i, kernel_cells_j, kernel_cells_k;

  // default constructor; initialize everything yourself!
  HAMIL_DATA() : N_cells(0), kernel_cells_i(0, 0), kernel_cells_j(0, 0),
  kernel_cells_k(0, 0) {
    numerical = new HAMIL_NUMERICAL;
  }

  // secondary constructor; used by call_hamil
  HAMIL_DATA(struct DATA *data, real_prec *A, real_prec *B_f, real_prec *B_r,
             real_prec *C, real_prec *D, real_prec *E, vector<int> cells_i,
             std::vector<int> cells_j, std::vector<int> cells_k,
             int N_cells_i) :
  N_cells(N_cells_i),
  kernel_cells_i(cells_i.begin(), cells_i.end()),
  kernel_cells_j(cells_j.begin(), cells_j.end()),
  kernel_cells_k(cells_k.begin(), cells_k.end()) {
    numerical = new HAMIL_NUMERICAL;

    numerical->N1 = data->numerical->N1;
    numerical->N2 = data->numerical->N2;
    numerical->N3 = data->numerical->N3;
    numerical->N = data->numerical->N;
    numerical->Nhalf = data->numerical->Nhalf;
    numerical->L1 = data->numerical->L1;
    numerical->L2 = data->numerical->L2;
    numerical->L3 = data->numerical->L3;
    numerical->vol = data->numerical->vol;
    numerical->d1 = data->numerical->d1;
    numerical->d2 = data->numerical->d2;
    numerical->d3 = data->numerical->d3;
    numerical->min1 = data->numerical->xllc;
    numerical->min2 = data->numerical->yllc;
    numerical->min3 = data->numerical->zllc;
    // Redshift space distortions
    numerical->xobs = data->numerical->xobs;
    numerical->yobs = data->numerical->yobs;
    numerical->zobs = data->numerical->zobs;
    numerical->planepar = data->numerical->planepar;
    numerical->periodic = data->numerical->periodic;

    numerical->N_bin = data->numerical->N_bin;
    numerical->mk = data->numerical->mk;

    numerical->kth = data->numerical->slength;

    numerical->iGibbs = data->numerical->iGibbs;

    numerical->rejections = 0; // EGP: added for rejection rate
    numerical->total_steps_lim = data->numerical->total_steps_lim;

    numerical->itmax = 2000;
    numerical->INV_SUCCESS = 0;/* indicate inversion success initialize with 0*/

    numerical->N_eps_fac = data->numerical->N_eps_fac;
    numerical->eps_fac = data->numerical->eps_fac;
    numerical->eps_fac_target = data->numerical->eps_fac_target;
    numerical->eps_fac_power = data->numerical->eps_fac_power;
    numerical->s_eps_total = data->numerical->s_eps_total;

    numerical->mass_type = data->numerical->mass_type;
    switch (numerical->mass_type) {
      case 0:  // 0: R: all ones (essentially: no mass term)
      numerical->mass_rs = true;
      numerical->mass_fs = false;
      break;
      case 1:  // 1: FS: inverse power spectrum
      numerical->mass_rs = false;
      numerical->mass_fs = true;
      break;
      case 2:  // 2: FS+FS: inverse power spectrum + likelihood force spectrum
      numerical->mass_rs = false;
      numerical->mass_fs = true;
      break;
      case 3:  // 3: FS+FS: inverse power spectrum + mean likelihood force (Wang+13)
      numerical->mass_rs = false;
      numerical->mass_fs = true;
      break;
      case 4:  // 4: FS: power spectrum
      numerical->mass_rs = false;
      numerical->mass_fs = true;
      break;
      case 5:  // 5: FS+R: inverse power spectrum + 1st order likelihood force expansion (Jasche+13)
      numerical->mass_rs = true;
      numerical->mass_fs = true;
      break;
      case 6:  // 1st order likelihood force expansion (Jasche+13) (R)
      numerical->mass_rs = true;
      numerical->mass_fs = false;
      break;
      case 60:  // type 0 until burn-in, type 6 afterwards
      numerical->mass_rs = true;
      numerical->mass_fs = false;
      break;
      default:
      stringstream message;
      message << "mass_type " << numerical->mass_type << " is not a valid value!";
      throw runtime_error(message.str());
    }

    numerical->massnum_init = data->numerical->massnum_init;
    numerical->massnum_burn = data->numerical->massnum_burn;

    gradpsi = A;
    mass_f = B_f;
    mass_r = B_r;
    posx = C;
    posy = D;
    posz = E;

    z = data->cosmology->z;
    ascale = data->cosmology->ascale;

    D1 = data->cosmology->D1;
    D2 = data->cosmology->D2;

    OM = data->cosmology->omega_m;
    OL = data->cosmology->omega_q;
    h = data->cosmology->h;

    biasP = data->observational->biasP;
    biasE = data->observational->biasE;

    x = data->observational->signal;
    nobs = data->observational->nobs;
    deltaX = data->observational->signalX;
    noise = data->observational->noise_sf;
    corrf = data->observational->corrf;

    Nmean_Gal = data->observational->Nmean_Gal;
    rho_c = data->observational->rho_c;
    code_norm = data->numerical->code_norm;
    mu = data->observational->mu;

    sigma_fac = data->observational->sigma_fac;
    sigma_min = data->observational->sigma_min;
    delta_min = data->observational->delta_min;

    signal_PS = data->observational->Power;

    window = data->observational->window;

    // function pointers
    log_like = data->observational->log_like;
    partial_f_delta_x_log_like = data->observational->partial_f_delta_x_log_like;
    grad_f_delta_x_comp = data->observational->grad_f_delta_x_comp;
    log_prior = data->observational->log_prior;
    grad_log_prior = data->observational->grad_log_prior;

    // model flags
    sfmodel = data->numerical->sfmodel;
    rsd_model = data->numerical->rsd_model;

    // Testing:
    numerical->mass_factor = data->numerical->mass_factor;
    numerical->grad_psi_prior_factor = data->numerical->grad_psi_prior_factor;
    numerical->grad_psi_likeli_factor = data->numerical->grad_psi_likeli_factor;
    numerical->grad_psi_prior_conjugate = data->numerical->grad_psi_prior_conjugate;
    numerical->grad_psi_likeli_conjugate = data->numerical->grad_psi_likeli_conjugate;
    numerical->grad_psi_prior_times_i = data->numerical->grad_psi_prior_times_i;
    numerical->grad_psi_likeli_times_i = data->numerical->grad_psi_likeli_times_i;
    numerical->div_dH_by_N = data->numerical->div_dH_by_N;
    numerical->calc_h = data->numerical->calc_h;
    numerical->deltaQ_factor = data->numerical->deltaQ_factor;

    numerical->particle_kernel = data->numerical->particle_kernel;
    numerical->particle_kernel_h = data->numerical->particle_kernel_h;

    // FFTW
    numerical->in_r2c = data->numerical->in_r2c;
    numerical->out_c2r = data->numerical->out_c2r;
    numerical->in_c2r = data->numerical->in_c2r;
    numerical->out_r2c = data->numerical->out_r2c;
    numerical->R2Cplan = data->numerical->R2Cplan;
    numerical->C2Rplan = data->numerical->C2Rplan;

    // curses
    // curses = data->curses;

    numerical->correct_delta = data->numerical->correct_delta;
  }

  ~HAMIL_DATA() {
    delete numerical;
  }
};

max_stars_repo_name: egpbos/barcode

max_issues_repo_name: egpbos/barcode

max_forks_repo_name: egpbos/barcode


text:
/* matrix/gsl_matrix_long_double.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_MATRIX_LONG_DOUBLE_H__
#define __GSL_MATRIX_LONG_DOUBLE_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_long_double.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  long double * data;
  gsl_block_long_double * block;
  int owner;
} gsl_matrix_long_double;

typedef struct
{
  gsl_matrix_long_double matrix;
} _gsl_matrix_long_double_view;

typedef _gsl_matrix_long_double_view gsl_matrix_long_double_view;

typedef struct
{
  gsl_matrix_long_double matrix;
} _gsl_matrix_long_double_const_view;

typedef const _gsl_matrix_long_double_const_view gsl_matrix_long_double_const_view;

/* Allocation */

gsl_matrix_long_double *
gsl_matrix_long_double_alloc (const size_t n1, const size_t n2);

gsl_matrix_long_double *
gsl_matrix_long_double_calloc (const size_t n1, const size_t n2);

gsl_matrix_long_double *
gsl_matrix_long_double_alloc_from_block (gsl_block_long_double * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

gsl_matrix_long_double *
gsl_matrix_long_double_alloc_from_matrix (gsl_matrix_long_double * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

gsl_vector_long_double *
gsl_vector_long_double_alloc_row_from_matrix (gsl_matrix_long_double * m,
                                        const size_t i);

gsl_vector_long_double *
gsl_vector_long_double_alloc_col_from_matrix (gsl_matrix_long_double * m,
                                        const size_t j);

void gsl_matrix_long_double_free (gsl_matrix_long_double * m);

/* Views */

_gsl_matrix_long_double_view
gsl_matrix_long_double_submatrix (gsl_matrix_long_double * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

_gsl_vector_long_double_view
gsl_matrix_long_double_row (gsl_matrix_long_double * m, const size_t i);

_gsl_vector_long_double_view
gsl_matrix_long_double_column (gsl_matrix_long_double * m, const size_t j);

_gsl_vector_long_double_view
gsl_matrix_long_double_diagonal (gsl_matrix_long_double * m);

_gsl_vector_long_double_view
gsl_matrix_long_double_subdiagonal (gsl_matrix_long_double * m, const size_t k);

_gsl_vector_long_double_view
gsl_matrix_long_double_superdiagonal (gsl_matrix_long_double * m, const size_t k);

_gsl_matrix_long_double_view
gsl_matrix_long_double_view_array (long double * base,
                             const size_t n1,
                             const size_t n2);

_gsl_matrix_long_double_view
gsl_matrix_long_double_view_array_with_tda (long double * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

_gsl_matrix_long_double_view
gsl_matrix_long_double_view_vector (gsl_vector_long_double * v,
                              const size_t n1,
                              const size_t n2);

_gsl_matrix_long_double_view
gsl_matrix_long_double_view_vector_with_tda (gsl_vector_long_double * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

_gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_submatrix (const gsl_matrix_long_double * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

_gsl_vector_long_double_const_view
gsl_matrix_long_double_const_row (const gsl_matrix_long_double * m,
                            const size_t i);

_gsl_vector_long_double_const_view
gsl_matrix_long_double_const_column (const gsl_matrix_long_double * m,
                               const size_t j);

_gsl_vector_long_double_const_view
gsl_matrix_long_double_const_diagonal (const gsl_matrix_long_double * m);

_gsl_vector_long_double_const_view
gsl_matrix_long_double_const_subdiagonal (const gsl_matrix_long_double * m,
                                    const size_t k);

_gsl_vector_long_double_const_view
gsl_matrix_long_double_const_superdiagonal (const gsl_matrix_long_double * m,
                                      const size_t k);

_gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_array (const long double * base,
                                   const size_t n1,
                                   const size_t n2);

_gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_array_with_tda (const long double * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

_gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_vector (const gsl_vector_long_double * v,
                                    const size_t n1,
                                    const size_t n2);

_gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_vector_with_tda (const gsl_vector_long_double * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

long double   gsl_matrix_long_double_get(const gsl_matrix_long_double * m, const size_t i, const size_t j);
void    gsl_matrix_long_double_set(gsl_matrix_long_double * m, const size_t i, const size_t j, const long double x);

long double * gsl_matrix_long_double_ptr(gsl_matrix_long_double * m, const size_t i, const size_t j);
const long double * gsl_matrix_long_double_const_ptr(const gsl_matrix_long_double * m, const size_t i, const size_t j);

void gsl_matrix_long_double_set_zero (gsl_matrix_long_double * m);
void gsl_matrix_long_double_set_identity (gsl_matrix_long_double * m);
void gsl_matrix_long_double_set_all (gsl_matrix_long_double * m, long double x);

int gsl_matrix_long_double_fread (FILE * stream, gsl_matrix_long_double * m) ;
int gsl_matrix_long_double_fwrite (FILE * stream, const gsl_matrix_long_double * m) ;
int gsl_matrix_long_double_fscanf (FILE * stream, gsl_matrix_long_double * m);
int gsl_matrix_long_double_fprintf (FILE * stream, const gsl_matrix_long_double * m, const char * format);

int gsl_matrix_long_double_memcpy(gsl_matrix_long_double * dest, const gsl_matrix_long_double * src);
int gsl_matrix_long_double_swap(gsl_matrix_long_double * m1, gsl_matrix_long_double * m2);

int gsl_matrix_long_double_swap_rows(gsl_matrix_long_double * m, const size_t i, const size_t j);
int gsl_matrix_long_double_swap_columns(gsl_matrix_long_double * m, const size_t i, const size_t j);
int gsl_matrix_long_double_swap_rowcol(gsl_matrix_long_double * m, const size_t i, const size_t j);
int gsl_matrix_long_double_transpose (gsl_matrix_long_double * m);
int gsl_matrix_long_double_transpose_memcpy (gsl_matrix_long_double * dest, const gsl_matrix_long_double * src);

long double gsl_matrix_long_double_max (const gsl_matrix_long_double * m);
long double gsl_matrix_long_double_min (const gsl_matrix_long_double * m);
void gsl_matrix_long_double_minmax (const gsl_matrix_long_double * m, long double * min_out, long double * max_out);

void gsl_matrix_long_double_max_index (const gsl_matrix_long_double * m, size_t * imax, size_t *jmax);
void gsl_matrix_long_double_min_index (const gsl_matrix_long_double * m, size_t * imin, size_t *jmin);
void gsl_matrix_long_double_minmax_index (const gsl_matrix_long_double * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

int gsl_matrix_long_double_isnull (const gsl_matrix_long_double * m);
int gsl_matrix_long_double_ispos (const gsl_matrix_long_double * m);
int gsl_matrix_long_double_isneg (const gsl_matrix_long_double * m);

int gsl_matrix_long_double_add (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
int gsl_matrix_long_double_sub (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
int gsl_matrix_long_double_mul_elements (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
int gsl_matrix_long_double_div_elements (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
int gsl_matrix_long_double_scale (gsl_matrix_long_double * a, const double x);
int gsl_matrix_long_double_add_constant (gsl_matrix_long_double * a, const double x);
int gsl_matrix_long_double_add_diagonal (gsl_matrix_long_double * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
int gsl_matrix_long_double_get_row(gsl_vector_long_double * v, const gsl_matrix_long_double * m, const size_t i);
int gsl_matrix_long_double_get_col(gsl_vector_long_double * v, const gsl_matrix_long_double * m, const size_t j);
int gsl_matrix_long_double_set_row(gsl_matrix_long_double * m, const size_t i, const gsl_vector_long_double * v);
int gsl_matrix_long_double_set_col(gsl_matrix_long_double * m, const size_t j, const gsl_vector_long_double * v);

/* inline functions if you are using GCC */

#ifdef HAVE_INLINE
extern inline
long double
gsl_matrix_long_double_get(const gsl_matrix_long_double * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
    }
#endif
  return m->data[i * m->tda + j] ;
}

extern inline
void
gsl_matrix_long_double_set(gsl_matrix_long_double * m, const size_t i, const size_t j, const long double x)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
    }
#endif
  m->data[i * m->tda + j] = x ;
}

extern inline
long double *
gsl_matrix_long_double_ptr(gsl_matrix_long_double * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (long double *) (m->data + (i * m->tda + j)) ;
}

extern inline
const long double *
gsl_matrix_long_double_const_ptr(const gsl_matrix_long_double * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (const long double *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_LONG_DOUBLE_H__ */

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified