text:
/* specfunc/coupling.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_gamma.h"
#include "gsl_sf_coupling.h"

#include "error.h"

inline
static
int locMax3(const int a, const int b, const int c)
{
  int d = GSL_MAX(a, b);
  return GSL_MAX(d, c);
}

inline
static
int locMin3(const int a, const int b, const int c)
{
  int d = GSL_MIN(a, b);
  return GSL_MIN(d, c);
}

inline
static
int locMin5(const int a, const int b, const int c, const int d, const int e)
{
  int f = GSL_MIN(a, b);
  int g = GSL_MIN(c, d);
  int h = GSL_MIN(f, g);
  return GSL_MIN(e, h);
}

/* See: [Thompson, Atlas for Computing Mathematical Functions] */

static
int
delta(int ta, int tb, int tc, gsl_sf_result * d)
{
  gsl_sf_result f1, f2, f3, f4;
  int status = 0;
  status += gsl_sf_fact_e((ta + tb - tc)/2, &f1);
  status += gsl_sf_fact_e((ta + tc - tb)/2, &f2);
  status += gsl_sf_fact_e((tb + tc - ta)/2, &f3);
  status += gsl_sf_fact_e((ta + tb + tc)/2 + 1, &f4);
  if(status != 0) {
    OVERFLOW_ERROR(d);
  }
  d->val = f1.val * f2.val * f3.val / f4.val;
  d->err = 4.0 * GSL_DBL_EPSILON * fabs(d->val);
  return GSL_SUCCESS;
}

static
int
triangle_selection_fails(int two_ja, int two_jb, int two_jc)
{
  return ((two_jb < abs(two_ja - two_jc)) || (two_jb > two_ja + two_jc));
}

static
int
m_selection_fails(int two_ja, int two_jb, int two_jc,
                  int two_ma, int two_mb, int two_mc)
{
  return (   abs(two_ma) > two_ja
          || abs(two_mb) > two_jb
	  || abs(two_mc) > two_jc
	  || GSL_IS_ODD(two_ja + two_ma)
	  || GSL_IS_ODD(two_jb + two_mb)
	  || GSL_IS_ODD(two_jc + two_mc)
          || (two_ma + two_mb + two_mc) != 0
	  );
}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_coupling_3j_e(int two_ja, int two_jb, int two_jc,
                        int two_ma, int two_mb, int two_mc,
			gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(two_ja < 0 || two_jb < 0 || two_jc < 0) {
    DOMAIN_ERROR(result);
  }
  else if(   triangle_selection_fails(two_ja, two_jb, two_jc)
          || m_selection_fails(two_ja, two_jb, two_jc, two_ma, two_mb, two_mc)
     ) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result n1_a, n1_b, n3_a, n3_b;
    gsl_sf_result d1_a, d1_b, d2_a, d2_b, d3_a, d3_b;
    gsl_sf_result n1, n2, n3;
    gsl_sf_result d1, d2, d3;
    double norm;
    double sign = (GSL_IS_ODD((two_ja - two_jb - two_mc)/2) ? -1.0 : 1.0);
    int tk, tkmin, tkmax;
    double sum_pos = 0.0;
    double sum_neg = 0.0;
    double phase;
    int status = 0;
    status += gsl_sf_fact_e((two_jc + two_ja - two_jb)/2, &n1_a);
    status += gsl_sf_fact_e((two_jc - two_ja + two_jb)/2, &n1_b);
    status += gsl_sf_fact_e((two_ja + two_jb - two_jc)/2, &n2);
    status += gsl_sf_fact_e((two_jc - two_mc)/2, &n3_a);
    status += gsl_sf_fact_e((two_jc + two_mc)/2, &n3_b);
    status += gsl_sf_fact_e((two_ja + two_jb + two_jc)/2 + 1, &d1);
    status += gsl_sf_fact_e((two_ja - two_ma)/2, &d2_a);
    status += gsl_sf_fact_e((two_ja + two_ma)/2, &d2_b);
    status += gsl_sf_fact_e((two_jb - two_mb)/2, &d3_a);
    status += gsl_sf_fact_e((two_jb + two_mb)/2, &d3_b);

    if(status != 0) {
      OVERFLOW_ERROR(result);
    }

    n1.val = n1_a.val * n1_b.val;
    n3.val = n3_a.val * n3_b.val;
    d2.val = d2_a.val * d2_b.val;
    d3.val = d3_a.val * d3_b.val;

    norm = sign * sqrt(n1.val*n2.val*n3.val)/sqrt(d1.val*d2.val*d3.val);

    tkmin = GSL_MAX(0, two_jb - two_ja - two_mc);
    tkmax = GSL_MIN(two_jc - two_ja + two_jb, two_jc - two_mc);

    phase = GSL_IS_ODD((tkmin + two_jb + two_mb)/2) ? -1.0 : 1.0;

    for(tk=tkmin; tk<=tkmax; tk += 2) {
      double term;

      status = 0;
      status += gsl_sf_fact_e((two_jb + two_jc + two_ma - tk)/2, &n1);
      status += gsl_sf_fact_e((two_ja - two_ma + tk)/2, &n2);
      status += gsl_sf_fact_e(tk/2, &d1_a);
      status += gsl_sf_fact_e((two_jc - two_ja + two_jb - tk)/2, &d1_b);
      status += gsl_sf_fact_e((two_jc - two_mc - tk)/2, &d2);
      status += gsl_sf_fact_e((two_ja - two_jb + two_mc + tk)/2, &d3);

      if(status != 0) {
        OVERFLOW_ERROR(result);
      }

      d1.val = d1_a.val * d1_b.val;

      term = phase * n1.val * n2.val / (d1.val * d2.val * d3.val);
      phase = -phase;

      if(norm*term >= 0.0) {
        sum_pos += norm*term;
      }
      else {
        sum_neg -= norm*term;
      }
    }

    result->val  = sum_pos - sum_neg;
    result->err  = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
    result->err += 2.0 * GSL_DBL_EPSILON * (tkmax - tkmin) * fabs(result->val);

    return GSL_SUCCESS;
  }
}

int
gsl_sf_coupling_6j_e(int two_ja, int two_jb, int two_jc,
                        int two_jd, int two_je, int two_jf,
			gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(   two_ja < 0 || two_jb < 0 || two_jc < 0
     || two_jd < 0 || two_je < 0 || two_je < 0
     ) {
    DOMAIN_ERROR(result);
  }
  else if(   triangle_selection_fails(two_ja, two_jb, two_je)
          || triangle_selection_fails(two_ja, two_jc, two_jf)
          || triangle_selection_fails(two_jb, two_jd, two_jf)
          || triangle_selection_fails(two_jc, two_jd, two_je)
     ) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result n1;
    gsl_sf_result d1, d2, d3, d4, d5, d6;
    double norm;
    int tk, tkmin, tkmax;
    double phase;
    double sum_pos = 0.0;
    double sum_neg = 0.0;
    double sumsq_err = 0.0;
    int status = 0;
    status += delta(two_ja, two_jb, two_je, &d1);
    status += delta(two_ja, two_jc, two_jf, &d2);
    status += delta(two_jb, two_jd, two_jf, &d3);
    status += delta(two_jc, two_jd, two_je, &d4);
    if(status != GSL_SUCCESS) {
      OVERFLOW_ERROR(result);
    }
    norm = sqrt(d1.val) * sqrt(d2.val) * sqrt(d3.val) * sqrt(d4.val);

    tkmin = locMax3(0,
                   two_ja + two_jd - two_je - two_jf,
                   two_jb + two_jc - two_je - two_jf);

    tkmax = locMin5(two_ja + two_jb + two_jc + two_jd + 2,
                    two_ja + two_jb - two_je,
		    two_jc + two_jd - two_je,
		    two_ja + two_jc - two_jf,
		    two_jb + two_jd - two_jf);

    phase = GSL_IS_ODD((two_ja + two_jb + two_jc + two_jd + tkmin)/2)
            ? -1.0
	    :  1.0;

    for(tk=tkmin; tk<=tkmax; tk += 2) {
      double term;
      double term_err;
      gsl_sf_result den_1, den_2;
      gsl_sf_result d1_a, d1_b;
      status = 0;

      status += gsl_sf_fact_e((two_ja + two_jb + two_jc + two_jd - tk)/2 + 1, &n1);
      status += gsl_sf_fact_e(tk/2, &d1_a);
      status += gsl_sf_fact_e((two_je + two_jf - two_ja - two_jd + tk)/2, &d1_b);
      status += gsl_sf_fact_e((two_je + two_jf - two_jb - two_jc + tk)/2, &d2);
      status += gsl_sf_fact_e((two_ja + two_jb - two_je - tk)/2, &d3);
      status += gsl_sf_fact_e((two_jc + two_jd - two_je - tk)/2, &d4);
      status += gsl_sf_fact_e((two_ja + two_jc - two_jf - tk)/2, &d5);
      status += gsl_sf_fact_e((two_jb + two_jd - two_jf - tk)/2, &d6);

      if(status != GSL_SUCCESS) {
        OVERFLOW_ERROR(result);
      }

      d1.val = d1_a.val * d1_b.val;
      d1.err = d1_a.err * fabs(d1_b.val) + fabs(d1_a.val) * d1_b.err;

      den_1.val  = d1.val*d2.val*d3.val;
      den_1.err  = d1.err * fabs(d2.val*d3.val);
      den_1.err += d2.err * fabs(d1.val*d3.val);
      den_1.err += d3.err * fabs(d1.val*d2.val);

      den_2.val  = d4.val*d5.val*d6.val;
      den_2.err  = d4.err * fabs(d5.val*d6.val);
      den_2.err += d5.err * fabs(d4.val*d6.val);
      den_2.err += d6.err * fabs(d4.val*d5.val);

      term  = phase * n1.val / den_1.val / den_2.val;
      phase = -phase;
      term_err  = n1.err / fabs(den_1.val) / fabs(den_2.val);
      term_err += fabs(term / den_1.val) * den_1.err;
      term_err += fabs(term / den_2.val) * den_2.err;

      if(term >= 0.0) {
        sum_pos += norm*term;
      }
      else {
        sum_neg -= norm*term;
      }

      sumsq_err += norm*norm * term_err*term_err;
    }

    result->val  = sum_pos - sum_neg;
    result->err  = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
    result->err += sqrt(sumsq_err / (0.5*(tkmax-tkmin)+1.0));
    result->err += 2.0 * GSL_DBL_EPSILON * (tkmax - tkmin + 2.0) * fabs(result->val);

    return GSL_SUCCESS;
  }
}

int
gsl_sf_coupling_9j_e(int two_ja, int two_jb, int two_jc,
                        int two_jd, int two_je, int two_jf,
			int two_jg, int two_jh, int two_ji,
			gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(   two_ja < 0 || two_jb < 0 || two_jc < 0
     || two_jd < 0 || two_je < 0 || two_jf < 0
     || two_jg < 0 || two_jh < 0 || two_ji < 0
     ) {
    DOMAIN_ERROR(result);
  }
  else if(   triangle_selection_fails(two_ja, two_jb, two_jc)
          || triangle_selection_fails(two_jd, two_je, two_jf)
          || triangle_selection_fails(two_jg, two_jh, two_ji)
          || triangle_selection_fails(two_ja, two_jd, two_jg)
          || triangle_selection_fails(two_jb, two_je, two_jh)
          || triangle_selection_fails(two_jc, two_jf, two_ji)
     ) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    int tk;
    int tkmin = locMax3(abs(two_ja-two_ji), abs(two_jh-two_jd), abs(two_jb-two_jf));
    int tkmax = locMin3(two_ja + two_ji, two_jh + two_jd, two_jb + two_jf);
    double sum_pos = 0.0;
    double sum_neg = 0.0;
    double sumsq_err = 0.0;
    double phase;
    for(tk=tkmin; tk<=tkmax; tk += 2) {
      gsl_sf_result s1, s2, s3;
      double term;
      double term_err;
      int status = 0;
      status += gsl_sf_coupling_6j_e(two_ja, two_ji, two_jd,  two_jh, tk, two_jg,  &s1);
      status += gsl_sf_coupling_6j_e(two_jb, two_jf, two_jh,  two_jd, tk, two_je,  &s2);
      status += gsl_sf_coupling_6j_e(two_ja, two_ji, two_jb,  two_jf, tk, two_jc,  &s3);
      if(status != GSL_SUCCESS) {
        OVERFLOW_ERROR(result);
      }
      term = s1.val * s2.val * s3.val;
      term_err  = s1.err * fabs(s2.val*s3.val);
      term_err += s2.err * fabs(s1.val*s3.val);
      term_err += s3.err * fabs(s1.val*s2.val);

      if(term >= 0.0) {
        sum_pos += (tk + 1) * term;
      }
      else {
        sum_neg -= (tk + 1) * term;
      }

      sumsq_err += ((tk+1) * term_err) * ((tk+1) * term_err);
    }

    phase = GSL_IS_ODD(tkmin) ? -1.0 : 1.0;

    result->val  = phase * (sum_pos - sum_neg);
    result->err  = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
    result->err += sqrt(sumsq_err / (0.5*(tkmax-tkmin)+1.0));
    result->err += 2.0 * GSL_DBL_EPSILON * (tkmax-tkmin + 2.0) * fabs(result->val);

    return GSL_SUCCESS;
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_coupling_3j(int two_ja, int two_jb, int two_jc,
                          int two_ma, int two_mb, int two_mc)
{
  EVAL_RESULT(gsl_sf_coupling_3j_e(two_ja, two_jb, two_jc,
                                   two_ma, two_mb, two_mc,
                                   &result));
}

double gsl_sf_coupling_6j(int two_ja, int two_jb, int two_jc,
                          int two_jd, int two_je, int two_jf)
{
  EVAL_RESULT(gsl_sf_coupling_6j_e(two_ja, two_jb, two_jc,
                                   two_jd, two_je, two_jf,
                                   &result));
}

double gsl_sf_coupling_9j(int two_ja, int two_jb, int two_jc,
                          int two_jd, int two_je, int two_jf,
                          int two_jg, int two_jh, int two_ji)
{
  EVAL_RESULT(gsl_sf_coupling_9j_e(two_ja, two_jb, two_jc,
                                   two_jd, two_je, two_jf,
                                   two_jg, two_jh, two_ji,
                                   &result));
}


text:
/* Copyright Koby Hayashi 2018 */

#ifndef DIMTREE_DDTTENSOR_H_
#define DIMTREE_DDTTENSOR_H_

#include <cblas.h>
#include <math.h>
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <time.h>

/**
This is a header file that contains defiition for:
        1) KTENSOR
        2) TENSOR
        3) Hadamard Product - Hadamard()
        4) Multi_KR() - function do to multiple KR products
        5) printM() - function to print a matrix in column major odering
        6) compareM() - function for comparing two matrices using a certain
tolerance 7) TransposeM() 8) MHada()

Move function to its own c file with MTTKRP
make inputs to KR just a ktensor struct

NOTES and QUESTIONS:
-should we make functions that act like constructors that take points to
        various structs and fill them with given argument values;
        EXAMPLE: ktensor* ktensor_const(double ** factors, long int * dims, long
int nmodes, long int rank){//FILL VALUES}

*/

/**
KTENSOR definition
Should Include:
        1) Array of Factor Matrices     (double **)
        2) Array of Dimensions      (long int *)
        3) Number of Modes      (long int)
        4) Rank             (long int)
        5) Weights or Lambda values
*/

typedef struct {
  double **factors;
  long int *dims;
  long int nmodes;
  long int rank;
  long int dims_product;
  double *lambdas;
} ktensor;

/**
TENSOR definition
Should Include:
        1) Data         (double *)
        2) Dims         (long int *)
        3) Number of Modes  (long int)

        - the Data is stored column wise in a the first mode matricization of
the tensor.
*/

typedef struct {
  double *data;
  long int *dims;
  long int nmodes;
  long int dims_product;
} tensor;

/**
        For processing command line inputs
*/
typedef struct {
  long int rank;
  long int nmodes;
  double noise;
  long int num_threads;
  long int max_iters;
  double tolerance;
  long int *dims;
} tensor_inputs;

/**
        Left or Right direction type
        Used as an input for the MTTKRP 2Step function
*/
typedef enum { left, right } direction;

/**
        Function Prototypes
*/

/**
        Utility fuctions

        printM_ColMajor()       - prints a marix in column major order
        printM_RowMajor()       - prints a Matrix in row-major order
        print_Ktensor_RowMajor()    - prints a Ktensor object
        TransposeM()            - transposes the matrix A into A_T
        CompareM()          - compares two matrices of the same size
   element by element compute_Iteration_Space() - divides up an iteration space
   amongst threads Gen_Tensor()         - generates a rank R tensor
        clean_Up_Gen_Tensor()       - cleans up a tensor and ktensor
        print_dgemm_inputs()        - prints inputs to a cblas_dgemm call
        print_dgemv_inputs()        - prints inputs to a cblas_dgemv call
        Full_nMode_Matricization_RowMajor() - creates a matricization of a
   tensor from a tensor model
*/
void printM_ColMajor(double *M, long int num_cols, long int num_rows);
void printM_RowMajor(double *M, long int num_cols, long int num_rows);
void print_Ktensor_RowMajor(ktensor *Y);
void TransposeM(double *A, double *A_T, long int rowsA, long int colsA);
long int CompareM(double *A, double *B, long int nRows, long int nCols,
                  double eps);
void compute_Iteration_Space(long int thread_id, long int num_threads,
                             long int iter_space, long int *start,
                             long int *end);
void Gen_Tensor(ktensor *Y, tensor *T, long int R, long int N, long int *D,
                double noise);
void clean_Up_Gen_Tensor(ktensor *Y, tensor *T);
void print_tensor(tensor *T, long int show_data);
void print_dgemm_inputs(CBLAS_ORDER dgemm_layout, CBLAS_TRANSPOSE transA,
                        CBLAS_TRANSPOSE transB, long int m, long int n,
                        long int k, double alpha, long int strideA,
                        long int strideB, double beta, long int strideC);
void print_dgemv_inputs(CBLAS_ORDER dgemv_layout, CBLAS_TRANSPOSE transA,
                        long int m, long int n, double alpha, long int T_offset,
                        long int tensor_stride, long int A_offset,
                        long int A_stride, double beta,
                        long int output_col_stride, long int output_stride);
void Full_nMode_Matricization_RowMajor(tensor *T, ktensor *Y, long int n);

/**
        Sequential KRP functions

        Multi_KR_RowMajor()     - Does the rev-KRP of all factor
   matrices of Y, skipping the nth mode, has reuse
        KR_RowMajor()           - Does the rev-KRP of 2 factor matrices
        update_Partial_Hadamards()  - Updates partial hadamards for
   Multi_KR_RowMajor()

*/
void Multi_KR_RowMajor(ktensor *Y, double *C, long int n);
void KR_RowMajor(ktensor *Y, double *C, long int n);
void update_Partial_Hadamards(ktensor *Y, long int *indexers,
                              double *partial_Hadamards,
                              double **reordered_Factors);

/**
        Parallel KRP functions

        parallel_Multi_KR_RowMajor()        - computes the pieces of a KRP
   of all given factor matrices in Y
        parallel_KR_RowMajor()          - handles the Y->nmodes == 2
   case fo above parallel_update_Partial_Hadamards()    - updates the
   Partial_Hadamards for the above
        wrapper_Parallel_Multi_revKRP()     - wrapper function for doing a
   parallel revKRP, uses the 3 above functions
*/
void parallel_Multi_KR_RowMajor(ktensor *Y, double *C, long int start,
                                long int end);
void parallel_KR_RowMajor(ktensor *Y, double *C, long int start, long int end);
void parallel_update_Partial_Hadamards(ktensor *Y, long int *indexers,
                                       double *partial_Hadamards);
void wrapper_Parallel_Multi_revKRP(ktensor *Y, long int num_threads,
                                   double *KRP);

/**
        CP_ALS functions

        MTTKRP_rowMajor()           - performs a slicing(no
   permuting needed) MTTKRP in row major Upper_Hadamard_RowMajor()
   - Hadamard product of upper part of symmetric matrices in row major order
        normalize_Factor_Matrix_RowMajor()  - normalizes the nth factor
   matrix for a row major ktensor normalize_Ktensor_RowMajor()      -
   normalizes all the factor matrices for a row major ktensor
        do_SYRKs_RowMajor()         - does the SYRKS of all factor
   matrices in a given Y and stores them approximation_Error()
   - comptues the error of a given CP model against the original tensor
        MHada_RowMajor()            - does the hadamard product of
   upper triangular matrices
*/
void MTTKRP_RowMajor(tensor *T, double *K, double *C, long int rank,
                     long int n);
void Upper_Hadamard_RowMajor(long int nRows, long int nCols, double *A,
                             double *B, double *C);
void normalize_Factor_Matrix_RowMajor(ktensor *Y, long int n);
void normalize_Ktensor_RowMajor(ktensor *Y);
void do_SYRKs_RowMajor(ktensor *Y, double **SYRKs);
double approximation_Error(double *X, double *KR, ktensor *Y, long int n);
double CP_ALS_efficient_error_computation(ktensor *Y, long int n,
                                          double *MTTKRP, double *V, double *S,
                                          double tensor_norm);
double CP_ALS_naive_error_computation(tensor *T, ktensor *Y, double *X,
                                      double *KRP, long int n);
void MHada_RowMajor(ktensor *Y, double **SYRKs, double *V, long int n);

void reorder_Factors(ktensor *Y, double **reordered_Factors,
                     long int *reordered_Dims, long int n);
void reorder_Ktensor(ktensor *Y, ktensor *nY, long int n);
void compute_KRP_Indices(long int j, ktensor *Y, long int *indeces);

void update_Partial_Hadamards_2(ktensor *Y, long int *indexers,
                                double *partial_Hadamards);

/**
        Input and object validity functions

        dims_check()        - checks an array of dimensions with the lenght
   to verify the validity process_inputs()  - processes general inputs for
   generating a rank R tensor
*/
void dims_check(long int *dims, long int length);
void process_inputs(long int argc, char *argv[], tensor_inputs *inputs);
void destroy_inputs(tensor_inputs *inputs);  // Parallel MTTKRP function

// Functions for manipulating and manageing ktensor and tensor objects
void LR_Ktensor_Reordering_newY(ktensor *Y, ktensor *nY, long int n,
                                direction D);
void LR_Ktensor_Reordering_existingY(ktensor *Y, ktensor *nY, long int n,
                                     direction D);
void destruct_Ktensor(ktensor *Y, long int clear_factors);
void destruct_Tensor(tensor *T);
void LR_tensor_Reduction(tensor *T, tensor *nT, long int n, direction D);
void tensor_data_swap(tensor *t1, tensor *t2);
void ktensor_copy_constructor(ktensor *Y, ktensor *nY);
direction opposite_direction(direction D);
void remove_mode_Ktensor(ktensor *Y, long int n);

// mkl unavailable functions
void vdMul(long int, double *, double *, double *);

// define dgemm
extern "C" void dgemm_(char *transa, char *transb, int *m, int *n, int *k, double *alpha,
           double *a, int *lda, double *b, int *ldb, double *beta, double *c,
           int *ldc);

#endif  // DIMTREE_DDTTENSOR_H_

ext: h

max_stars_repo_name: lanl/DnMFkCPP

max_issues_repo_name: rvangara/DnMFk

max_forks_repo_name: rvangara/DnMFk


text:
#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"

int calc_units(fb_hier_t hier, fb_units_t *units)
{
  /* mass unit is total system mass */
  units->m = hier.hier[hier.hi[1]+0].m + hier.hier[hier.hi[1]+1].m +	\
    hier.hier[hier.hi[1]+2].m + hier.hier[hier.hi[1]+3].m +		\
    hier.hier[hier.hi[1]+4].m;
  /* length unit is one AU */
  units->l = FB_CONST_AU;
  /* everything else is derived */
  units->E = FB_CONST_G * fb_sqr(units->m) / units->l;
  units->v = sqrt(units->E/units->m);
  units->t = units->l / units->v;
  return(0);
}

int main(int argc, char *argv[])
{
  int j;
  double 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];

  /* set input parameters */
  input.ks = 0; /* turn K-S regularization off */
  input.tstop = 1.0e4; /* stopping time in units of units.t */
  input.Dflag = 0; /* don't output dynamical info to stdout */
  input.dt = 0.0; /* irrelevant when Dflag=0 */
  input.tcpustop = 120.0; /* stopping CPU time in seconds */
  input.absacc = 1.0e-9; /* integrator absolute accuracy */
  input.relacc = 1.0e-9; /* integrator relative accuracy */
  input.ncount = 500; /* number of integration steps between calls to fb_classify() */
  input.tidaltol = 1.0e-6; /* tidal perturbation required to force numerical */
                           /* integration of a binary node */
  input.fexp = 1.0; /* radius expansion factor of merger products */
  fb_debug = 0;

  /* initialize a few things */
  t = 0.0;
  hier.nstarinit = 5;
  hier.nstar = 5;
  fb_malloc_hier(&hier);
  fb_init_hier(&hier);

  /* set stellar properties */
  for (j=0; j<hier.nstar; j++) {
    hier.hier[hier.hi[1]+j].ncoll = 1;
    hier.hier[hier.hi[1]+j].id[0] = j;
    snprintf(hier.hier[hier.hi[1]+j].idstring, FB_MAX_STRING_LENGTH, "%d", j);
    hier.hier[hier.hi[1]+j].n = 1;
    hier.hier[hier.hi[1]+j].obj[0] = NULL;
    hier.hier[hier.hi[1]+j].obj[1] = NULL;
    hier.hier[hier.hi[1]+j].Eint = 0.0;
    hier.hier[hier.hi[1]+j].Lint[0] = 0.0;
    hier.hier[hier.hi[1]+j].Lint[1] = 0.0;
    hier.hier[hier.hi[1]+j].Lint[2] = 0.0;
    hier.hier[hier.hi[1]+j].R = 1.0e-3 * FB_CONST_RSUN;
    hier.hier[hier.hi[1]+j].m = FB_CONST_MSUN;
    hier.hier[hier.hi[1]+j].x[2] = 0.0; /* everything in x-y plane */
    hier.hier[hier.hi[1]+j].v[2] = 0.0;
  }

  /* set some positions and velocities */
  hier.hier[hier.hi[1]+0].x[0] = 10.0 * FB_CONST_AU;
  hier.hier[hier.hi[1]+0].x[1] = 0.0 * FB_CONST_AU;
  hier.hier[hier.hi[1]+0].v[0] = 10.0e5;
  hier.hier[hier.hi[1]+0].v[1] = 0.0;

  hier.hier[hier.hi[1]+1].x[0] = 10.0 * FB_CONST_AU;
  hier.hier[hier.hi[1]+1].x[1] = 10.0 * FB_CONST_AU;
  hier.hier[hier.hi[1]+0].v[0] = -10.0e5;
  hier.hier[hier.hi[1]+0].v[1] = 0.0;

  hier.hier[hier.hi[1]+2].x[0] = -30.0 * FB_CONST_AU;
  hier.hier[hier.hi[1]+2].x[1] = -5.0 * FB_CONST_AU;
  hier.hier[hier.hi[1]+0].v[0] = 0.0;
  hier.hier[hier.hi[1]+0].v[1] = 5.0e5;

  hier.hier[hier.hi[1]+3].x[0] = 0.0 * FB_CONST_AU;
  hier.hier[hier.hi[1]+3].x[1] = 27.0 * FB_CONST_AU;
  hier.hier[hier.hi[1]+0].v[0] = -5.0e5;
  hier.hier[hier.hi[1]+0].v[1] = -4.0e5;

  hier.hier[hier.hi[1]+4].x[0] = 10.0 * FB_CONST_AU;
  hier.hier[hier.hi[1]+4].x[1] = -32.0 * FB_CONST_AU;
  hier.hier[hier.hi[1]+0].v[0] = 5.0e5;
  hier.hier[hier.hi[1]+0].v[1] = -1.0e5;

  /* set units with our own routine, then normalize */
  calc_units(hier, &units);
  fb_normalize(&hier, units);

  /* call Fewbody to evolve system */
  retval = fewbody(input, &hier, &t);

  /* all the rest is parsing the output */
  if (retval.retval == 1) {
    fprintf(stderr, "encounter complete.\n");
  } else {
    fprintf(stderr, "encounter NOT complete.\n");
  }

  fprintf(stderr, "final configuration:  %s  (%s)\n",
	  fb_sprint_hier(hier, string1),
	  fb_sprint_hier_hr(hier, string2));

  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, "orbital parameters of outermost binaries:\n");
  for (j=0; j<hier.nobj; j++) {
    if (hier.obj[j]->n >= 2) {
      fprintf(stderr, "j=%d  a=%.6g AU  e=%.6g\n", j,
	      hier.obj[j]->a * units.l / FB_CONST_AU, hier.obj[j]->e);
    }
  }

  fb_free_hier(hier);
  return(0);
}

ext: c

max_stars_repo_name: gnodvi/cosmos

max_issues_repo_name: gnodvi/cosmos

max_forks_repo_name: gnodvi/cosmos


text:
/*
	Simbicon 1.5 Controller Editor Framework,
	Copyright 2009 Stelian Coros, Philippe Beaudoin and Michiel van de Panne.
	All rights reserved. Web: www.cs.ubc.ca/~van/simbicon_cef

	This file is part of the Simbicon 1.5 Controller Editor Framework.

	Simbicon 1.5 Controller Editor Framework 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.

	Simbicon 1.5 Controller Editor Framework 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 Simbicon 1.5 Controller Editor Framework.
	If not, see <http://www.gnu.org/licenses/>.
*/

#pragma once

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include <MathLib/MathLib.h>
#include <MathLib/ThreeTuple.h>

#include <gsl/matrix/gsl_matrix.h>
#include <gsl/vector/gsl_vector.h>

#include <MathLib/Matrix.h>

#define VECTOR_AT(v, i) (*((v->data + (i * v->tda ))))

/*====================================================================================================================================================================*
 | This class will be used to represent matrices of arbitrary sizes (m rows by n columns) that have elements of type double. The underlying data strucutre used by    |
 | this class is gsl's (Gnu Scientific Library) matrix class. This class also makes use of the ATLAS implementation of BLAS for some operations such as matrix-matrix |
 | multiplication. This class is meant to improve performance, not necessarily ease of use.                                                                           |
 *====================================================================================================================================================================*/

class Vector : public Matrix{
public:
	/**
		constructor - creates an n row vector that is not initialized to any particular values
	*/
	Vector(int n);

	/**
		default constructor
	*/
	Vector();

	/**
		copy constructor - performs a deep copy of the matrix passed in as a parameter.
	*/
	Vector(const Vector& other);

	/**
		destructor.
	*/
    virtual ~Vector();

	/**
		copy operator - performs a deep copy of the Vector passed in as a parameter.
	*/
	Vector& operator=(const Vector &other);

	/**
		copy operator - performs a deep copy of the Vector passed in as a parameter.
	*/
	Vector& operator=(const Matrix &other);

	/**
		this method performs a shallow copy of the Vector that is passed in as a parameter.
	*/
	void shallowCopy(const Matrix& other);

	/**
		this method performs a deep copy of the vector that is passed in as a paramerer.
	*/
	void deepCopy(const Matrix& other);

	/**
		This method sets the current vector to be equal to one of the products: A * b or A'*b.
		The value of transA indicates if A is transposed or not
	*/
	void setToProductOf(const Matrix& A, const Matrix& B, bool transA = false, bool transB = false);

	/**
		This method sets the current vector to be equal to one of the rows of A - shallow column only!
	*/
	void setToRow(const Matrix& A, int row, int start = 0, int howMany = -1);

	/**
		This method sets the current vector to be equal to one of the cols of A - shallow column only!
	*/
	void setToCol(const Matrix& A, int col, int start = 0, int howMany = -1);

	/**
		This method prints the contents of the matrix - testing purpose only.
	*/
	void printVector() const;

	/**
		This method returns a copy of the value of the matrix at (i,j)
	*/
	double get(int i) const;

	/**
		This method sets the value of the matrix at (i,j) to newVal.
	*/
	void set(int i, double newVal);

	/**
		Computes the 2-norm squared for the current vector.
	*/
	inline double normSquared(){
		int r = getRowCount();
		double result = 0;
		for (int i=0;i<r;i++){
			double n = MATRIX_AT(matrix, i, 0) ;
			result += n*n;
		}
		return result;
	}

};

void testVectorClass();


text:
#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <time.h>
#include <string.h>

/*
 * A benchmark to test performance of calculating determinants of
 * 6x6 matrices with the GNU Scientific Library (GSL) package.
 */

int main(int argc, char * argv[])
{
  if (argc != 3 || strcmp(argv[1], "-h") == 0)
  {
    printf("Usage: ./gslDeterm N_MATRIX(<=1000) ITERATIONS(~1000)\n");
    return 1;
  }

  srand(time(NULL));
  int MAT_DIM = 6;
  int N_MATRIX = atoi(argv[1]);
  int ITERATIONS = atoi(argv[2]);
  int i,j,k, count;
  double start, end, t_time = 0.0;

  printf("N_MATRIX: %d, ITERATIONS: %d\n", N_MATRIX, ITERATIONS);

  gsl_permutation *p;
  gsl_matrix *m_list[N_MATRIX];

  // Do this ITERATIONS times, timing only the loop with the determinant
  // finding
  for (count=0; count<ITERATIONS; count++)
  {
    // initialise each element of m_list as a gsl_matrix
    for (i=0; i<N_MATRIX; i++)
      m_list[i] = gsl_matrix_alloc(MAT_DIM,MAT_DIM);

    // initialise each element of each gsl_matrix as some random double
    // less than 100
    for (i=0; i<N_MATRIX; i++)
      for (j=0; j<MAT_DIM; j++)
        for (k=0; k<MAT_DIM; k++)
          gsl_matrix_set (m_list[i], j, k, (double) (rand()%100));

    volatile double det;
    int signum;
    p = gsl_permutation_alloc(m_list[0]->size1);

    //start time: find the determinant of each matrix in m_list
    start = (double)clock() /(double) CLOCKS_PER_SEC;

    for (i=0; i<N_MATRIX; i++)
    {
      gsl_linalg_LU_decomp(m_list[0], p, &signum);
      det = gsl_linalg_LU_det(m_list[0], signum);
    }

    end = (double)clock() / (double) CLOCKS_PER_SEC;
    t_time += end - start;
  }

  printf("GSL: %d matrices, %d iterations, cpu time: %fs\n",
                    N_MATRIX, ITERATIONS, t_time);

  gsl_permutation_free(p);
  for (i=0; i<N_MATRIX; i++)
    gsl_matrix_free(m_list[i]);

  return 0;
}

max_stars_repo_name: tcrundall/chronostar

max_issues_repo_name: tcrundall/chronostar

max_forks_repo_name: tcrundall/chronostar


text:
#ifndef boundary_h
#define boundary_h

#include <petsc.h>

typedef PetscErrorCode BCFunc(PetscInt, PetscReal, const PetscReal *, PetscInt,
                              PetscScalar *, void *);
BCFunc BCMMS, BCZero, BCClamp;

// -----------------------------------------------------------------------------
// Boundary Functions
// -----------------------------------------------------------------------------
// Note: If additional boundary conditions are added, an update is needed in
//         elasticity.h for the boundaryOptions variable.

// BCMMS - boundary function
// Values on all points of the mesh is set based on given solution below
// for u[0], u[1], u[2]
PetscErrorCode BCMMS(PetscInt dim, PetscReal load_increment,
                     const PetscReal coords[], PetscInt num_comp_u,
                     PetscScalar *u, void *ctx);

// BCClamp - fix boundary values with affine transformation at fraction of load
//   increment
PetscErrorCode BCClamp(PetscInt dim, PetscReal load_increment,
                       const PetscReal coords[], PetscInt num_comp_u,
                       PetscScalar *u, void *ctx);

#endif //boundary_h

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
/* interpolation/interp.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_interp.h>

#define DISCARD_STATUS(s) if ((s) != GSL_SUCCESS) { GSL_ERROR_VAL("interpolation error", (s),  GSL_NAN); }

gsl_interp *
gsl_interp_alloc (const gsl_interp_type * T, size_t size)
{
  gsl_interp * interp;

  if (size < T->min_size)
    {
      GSL_ERROR_NULL ("insufficient number of points for interpolation type",
                      GSL_EINVAL);
    }

  interp = (gsl_interp *) malloc (sizeof(gsl_interp));

  if (interp == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for interp struct",
                      GSL_ENOMEM);
    }

  interp->type = T;
  interp->size = size;

  if (interp->type->alloc == NULL)
    {
      interp->state = NULL;
      return interp;
    }

  interp->state = interp->type->alloc(size);

  if (interp->state == NULL)
    {
      free (interp);
      GSL_ERROR_NULL ("failed to allocate space for interp state", GSL_ENOMEM);
    };

  return interp;
}

int
gsl_interp_init (gsl_interp * interp, const double x_array[], const double y_array[], size_t size)
{
  if (size != interp->size)
    {
      GSL_ERROR ("data must match size of interpolation object", GSL_EINVAL);
    }

  interp->xmin = x_array[0];
  interp->xmax = x_array[size - 1];

  {
    int status = interp->type->init(interp->state, x_array, y_array, size);
    return status;
  }
}

const char *
gsl_interp_name(const gsl_interp * interp)
{
  return interp->type->name;
}

unsigned int
gsl_interp_min_size(const gsl_interp * interp)
{
  return interp->type->min_size;
}

void
gsl_interp_free (gsl_interp * interp)
{
  if (interp->type->free)
    interp->type->free (interp->state);
  free (interp);
}

int
gsl_interp_eval_e (const gsl_interp * interp,
                   const double xa[], const double ya[], double x,
                   gsl_interp_accel * a, double *y)
{
  if (x < interp->xmin)
    {
      *y = ya[0];
      return GSL_EDOM;
    }
  else if (x > interp->xmax)
    {
      *y = ya[interp->size - 1];
      return GSL_EDOM;
    }

  return interp->type->eval (interp->state, xa, ya, interp->size, x, a, y);
}

double
gsl_interp_eval (const gsl_interp * interp,
                 const double xa[], const double ya[], double x,
                 gsl_interp_accel * a)
{
  double y;
  int status = interp->type->eval (interp->state, xa, ya, interp->size, x, a, &y);

  DISCARD_STATUS(status);

  return y;
}

int
gsl_interp_eval_deriv_e (const gsl_interp * interp,
                         const double xa[], const double ya[], double x,
                         gsl_interp_accel * a,
                         double *dydx)
{
  if (x < interp->xmin)
    {
      *dydx = 0.0;
      return GSL_EDOM;
    }
  else if (x > interp->xmax)
    {
      *dydx = 0.0;
      return GSL_EDOM;
    }

  return interp->type->eval_deriv (interp->state, xa, ya, interp->size, x, a, dydx);
}

double
gsl_interp_eval_deriv (const gsl_interp * interp,
                       const double xa[], const double ya[], double x,
                       gsl_interp_accel * a)
{
  double dydx;
  int status = interp->type->eval_deriv (interp->state, xa, ya, interp->size, x, a, &dydx);

  DISCARD_STATUS(status);

  return dydx;
}

int
gsl_interp_eval_deriv2_e (const gsl_interp * interp,
                          const double xa[], const double ya[], double x,
                          gsl_interp_accel * a,
                          double * d2)
{
  if (x < interp->xmin)
    {
      *d2 = 0.0;
      return GSL_EDOM;
    }
  else if (x > interp->xmax)
    {
      *d2 = 0.0;
      return GSL_EDOM;
    }

  return interp->type->eval_deriv2 (interp->state, xa, ya, interp->size, x, a, d2);
}

double
gsl_interp_eval_deriv2 (const gsl_interp * interp,
                        const double xa[], const double ya[], double x,
                        gsl_interp_accel * a)
{
  double d2;
  int status = interp->type->eval_deriv2 (interp->state, xa, ya, interp->size, x, a, &d2);

  DISCARD_STATUS(status);

  return d2;
}

int
gsl_interp_eval_integ_e (const gsl_interp * interp,
                         const double xa[], const double ya[],
                         double a, double b,
                         gsl_interp_accel * acc,
                         double * result)
{
  if (a > b || a < interp->xmin || b > interp->xmax)
    {
      *result = 0.0;
      return GSL_EDOM;
    }
  else if(a == b)
    {
      *result = 0.0;
      return GSL_SUCCESS;
    }

  return interp->type->eval_integ (interp->state, xa, ya, interp->size, acc, a, b, result);
}

double
gsl_interp_eval_integ (const gsl_interp * interp,
                       const double xa[], const double ya[],
                       double a, double b,
                       gsl_interp_accel * acc)
{
  double result;
  int status = interp->type->eval_integ (interp->state, xa, ya, interp->size, acc, a, b, &result);

  DISCARD_STATUS(status);

  return result;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#ifndef _PYGSL_SOLVER_INTERN_H_
#define _PYGSL_SOLVER_INTERN_H_ 1

#define _PyGSL_SOLVER_API_MODULE 1
#include <pygsl/solver.h>
#undef _PyGSL_SOLVER_API_MODULE

#include <pygsl/general_helpers.h>
#include <pygsl/block_helpers.h>
#include <pygsl/function_helpers.h>
#include <setjmp.h>
#include <gsl/gsl_math.h>

/*
 * Collect all element accessor methods which just need a pointer to the struct
 * and return a value.
 */
struct _ElementAccessor{
     int_m_t method;
     const char * name;
};

static PyTypeObject PyGSL_solver_pytype;
#define PyGSL_solver_check(ob) ((ob)->ob_type == &PyGSL_solver_pytype)
#define _PyGSL_solver_n_init PyGSL_solver_n_init

static PyGSL_solver*
_PyGSL_solver_init(const struct _SolverStatic *mstatic);

static int
PyGSL_solver_set_called(PyGSL_solver *self);

#define PyGSL_SOLVER_SET_CALLED(ob) \
        (((ob)->set_called == 1) ? GSL_SUCCESS: PyGSL_solver_set_called((ob)))

static int
PyGSL_solver_func_set(PyGSL_solver *self, PyObject *args, PyObject *f,
		       PyObject *df, PyObject *fdf);

#if 0
static PyObject *
_PyGSL_solver_np_init(PyObject *self, PyObject *args, const solver_alloc_struct * alloc);

static PyObject *
_PyGSL_solver_i_vvdd(PyObject * self, PyObject * args, int_f_vvdd_t func);

#endif
#endif /* _PYGSL_SOLVER_INTERN_H_ */

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/*
 * Gsl_matrix: test gsl_matrix
 *
 * Copyright (c) 2012 Jérémie Decock
 *
 * Required: GSL library (libgsl0-dev)
 * Usage: gcc gsl_matrix.c -lgsl -lgslcblas -lm
 *    or: gcc gsl_matrix.c $(pkg-config --libs gsl)
 */

#include <gsl/gsl_matrix.h>

#define N 10

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

    gsl_matrix * m1 = gsl_matrix_calloc(N, N);
    gsl_matrix * m2 = gsl_matrix_calloc(N, N);

    gsl_matrix_set_all(m1, 2.0);
    gsl_matrix_set_all(m2, 6.0);

    gsl_matrix_mul_elements(m1, m2);

    gsl_matrix_fprintf(stdout, m1, "%.2f");

    gsl_matrix_free(m1);
    gsl_matrix_free(m2);

    return 0;

}

max_stars_repo_name: jeremiedecock/snippets

max_issues_repo_name: jeremiedecock/snippets

max_forks_repo_name: jeremiedecock/snippets


text:
#pragma once

#include <type_traits>

#include <cuda/define_specifiers.hpp>
#include <cuda/runtime_api.hpp>

#include <thrust/iterator/counting_iterator.h>
#include <thrust/iterator/transform_iterator.h>
#include <thrust/iterator/zip_iterator.h>

#include <gsl-lite/gsl-lite.hpp>

#include <cub/cub.cuh>

#include <thrustshift/constant.h>
#include <thrustshift/math.h>
#include <thrustshift/not-a-vector.h>

namespace thrustshift {

namespace device_function {

namespace implicit_unroll {

/*! \brief Write the selected values to the result if select_op is true.
 *
 *  All threads check if their value should be selected. Afterwards, it is
 *  determined within the warp how many values were selected. Then, the first
 *  in the warp determines the range in the result, which is reserved for the
 *  warp via an atomic increment on `result_pos`. Subsequently, each thread with
 *  a selected value and `valid_write` writes the result into the reserved range.
 *  Only full warps are allowed to enter this function.
 *
 *  \param x the value for each thread
 *  \param selected_values result range.
 *  \param selected_values_pos to an integer holding the ID of the first
 *      written result element by this function (usually initialized with
 *      zero). This counter is read and written atomically by the first lane of
 *      the warp. Usually this counter is shared among threads because it is
 *      incremented and read atomically.
 *  \param lane_id ID of the thread within the warp.
 *  \param valid_write some threads might be disabled with this flag (useful for
 *      handling range borders).
 *  \param select_op lambda to select the values.
 */
template <typename T, typename It, typename Bool, class F>
CUDA_FD void select_if_warp_aggregated(T x,
                                       It selected_values,
                                       int* selected_values_pos,
                                       int lane_id,
                                       Bool valid_write,
                                       F select_op) {

	const bool p = select_op(x) && valid_write;

	const unsigned umask = __ballot_sync(0xffffffff, p);
	const int num_selected_values = __popc(umask);

	constexpr int source_lane = 0;
	int pos;
	if (lane_id == source_lane) {
		pos = atomicAdd(selected_values_pos, num_selected_values);
	}
	const unsigned umask_before_me =
	    (~unsigned(0)) >> (sizeof(unsigned) * 8 - lane_id);
	const int num_selected_values_before_me = __popc(umask & umask_before_me);
	pos = __shfl_sync(0xffffffff, pos, source_lane);
	pos += num_selected_values_before_me;
	if (p) {
		selected_values[pos] = x;
	}
}

/*! \brief select_if as a block level primitive.
 *
 *  \param values iterator to range of length N
 *  \param N length of input values
 *  \param selected_values iterator to range of length `*selected_values_pos + N`
 *  \param selected_values_pos pointer to an integer holding the ID of the first
 *      written selected_values element by this function (usually initialized
 *      with zero). This counter must be shared among all threads, which
 *      enter this function because it is incremented and read atomically.
 *  \param tid thread ID
 *  \param num_threads number of threads which enter this function.
 *  \param select_op lambda to select the values.
 */
template <typename It, typename ItResult, typename I0, typename I1, class F>
CUDA_FD void select_if(It values,
                       I0 N,
                       ItResult selected_values,
                       int* selected_values_pos,
                       int tid,
                       I1 num_threads,
                       F select_op) {

	const int num_tiles = thrustshift::ceil_divide(N, num_threads);
	for (int tile_id = 0; tile_id < num_tiles - 1; ++tile_id) {
		const int j = tile_id * num_threads + tid;
		auto x = values[j];
		select_if_warp_aggregated(x,
		                          selected_values,
		                          selected_values_pos,
		                          tid % warp_size,
		                          std::bool_constant<true>(),
		                          select_op);
	}
	auto rest = N % num_threads;
	if (rest > 0) {
		const int j = (num_tiles - 1) * num_threads + tid;
		const bool valid_rw = tid < rest;
		using T = typename std::remove_const<
		    typename std::iterator_traits<It>::value_type>::type;
		const auto x = [&]() -> T {
			if (valid_rw) {
				return values[j];
			}
			return T{};
		}();
		select_if_warp_aggregated(x,
		                          selected_values,
		                          selected_values_pos,
		                          tid % warp_size,
		                          valid_rw,
		                          select_op);
	}
}

} // namespace implicit_unroll

} // namespace device_function

namespace async {

template <class ValuesRange,
          class SelectedRange,
          class NumSelectedPtr,
          class SelectOp,
          class MemoryResource>
void select_if(cuda::stream_t& stream,
               ValuesRange&& values,
               SelectedRange&& selected,
               NumSelectedPtr num_selected_ptr,
               SelectOp select_op,
               MemoryResource& delayed_memory_resource) {

	const std::size_t N = values.size();

	gsl_Expects(selected.size() == N);

	size_t tmp_bytes_size = 0;
	void* tmp_ptr = nullptr;

	using T = typename std::remove_reference<ValuesRange>::type::value_type;

	auto exec = [&] {
		cuda::throw_if_error(cub::DeviceSelect::If(tmp_ptr,
		                                           tmp_bytes_size,
		                                           values.data(),
		                                           selected.data(),
		                                           num_selected_ptr,
		                                           N,
		                                           select_op,
		                                           stream.handle()));
	};
	exec();
	auto tmp =
	    make_not_a_vector<uint8_t>(tmp_bytes_size, delayed_memory_resource);
	tmp_ptr = tmp.to_span().data();
	exec();
}

template <class ValuesIt,
          class SelectedIt,
          class NumSelectedPtr,
          class SelectOp,
          class MemoryResource>
void select_if(cuda::stream_t& stream,
               ValuesIt values,
               SelectedIt selected,
               NumSelectedPtr num_selected_ptr,
               size_t N,
               SelectOp select_op,
               MemoryResource& delayed_memory_resource) {

	size_t tmp_bytes_size = 0;
	void* tmp_ptr = nullptr;

	auto exec = [&] {
		cuda::throw_if_error(cub::DeviceSelect::If(tmp_ptr,
		                                           tmp_bytes_size,
		                                           values,
		                                           selected,
		                                           num_selected_ptr,
		                                           N,
		                                           select_op,
		                                           stream.handle()));
	};
	exec();
	auto tmp =
	    make_not_a_vector<uint8_t>(tmp_bytes_size, delayed_memory_resource);
	tmp_ptr = tmp.to_span().data();
	exec();
}

/*! \brief Select values based on predicate with their index.
 *
 *  \param values Range of length N with value_type `T`
 *  \param selected Range of length N with a value type which is assignable
 *      by a thrust::tuple<T, int>.
 *  \param select_op select predicate with signature
 *      ```
 *      select_op = [] __device__ (const thrust::tuple<T, int>& tup) { ... };
 *      ```
 */
template <class ValuesRange,
          class SelectedRange,
          class NumSelectedPtr,
          class SelectOp,
          class MemoryResource>
void select_if_with_index(cuda::stream_t& stream,
                          ValuesRange&& values,
                          SelectedRange&& selected,
                          NumSelectedPtr num_selected_ptr,
                          SelectOp select_op,
                          MemoryResource& delayed_memory_resource) {

	const std::size_t N = values.size();

	gsl_Expects(selected.size() == N);

	size_t tmp_bytes_size = 0;
	void* tmp_ptr = nullptr;

	using T = typename std::remove_reference<ValuesRange>::type::value_type;

	auto cit = thrust::make_counting_iterator(0);
	auto it = thrust::make_zip_iterator(thrust::make_tuple(values.data(), cit));

	auto exec = [&] {
		cuda::throw_if_error(cub::DeviceSelect::If(tmp_ptr,
		                                           tmp_bytes_size,
		                                           it,
		                                           selected.data(),
		                                           num_selected_ptr,
		                                           N,
		                                           select_op,
		                                           stream.handle()));
	};
	exec();
	auto tmp =
	    make_not_a_vector<uint8_t>(tmp_bytes_size, delayed_memory_resource);
	tmp_ptr = tmp.to_span().data();
	exec();
}

template <class ValuesRange,
          class SelectedRange,
          class SelectedIndexRange,
          class NumSelectedPtr,
          class SelectOp,
          class MemoryResource>
void select_if_with_index(cuda::stream_t& stream,
                          ValuesRange&& values,
                          SelectedRange&& selected,
                          SelectedIndexRange&& selected_indices,
                          NumSelectedPtr num_selected_ptr,
                          SelectOp select_op,
                          MemoryResource& delayed_memory_resource) {

	const std::size_t N = values.size();

	gsl_Expects(selected.size() == N);

	size_t tmp_bytes_size = 0;
	void* tmp_ptr = nullptr;

	using T = typename std::remove_reference<ValuesRange>::type::value_type;

	auto cit = thrust::make_counting_iterator(0);
	auto it = thrust::make_zip_iterator(thrust::make_tuple(values.data(), cit));
	auto result_it = thrust::make_zip_iterator(
	    thrust::make_tuple(selected.data(), selected_indices.data()));

	auto exec = [&] {
		cuda::throw_if_error(cub::DeviceSelect::If(tmp_ptr,
		                                           tmp_bytes_size,
		                                           it,
		                                           result_it,
		                                           num_selected_ptr,
		                                           N,
		                                           select_op,
		                                           stream.handle()));
	};
	exec();
	auto tmp =
	    make_not_a_vector<uint8_t>(tmp_bytes_size, delayed_memory_resource);
	tmp_ptr = tmp.to_span().data();
	exec();
}

} // namespace async

} // namespace thrustshift

max_stars_repo_name: pauleonix/thrustshift

max_issues_repo_name: pauleonix/thrustshift

max_forks_repo_name: pauleonix/thrustshift


text:
/*
 * testgen.c
 * Patrick Alken
 *
 * Compile: gcc -g -O2 -Wall -o testgen testgen.c -lm -lgsl -llapack -lf77blas -lcblas -latlas -lg2c
 *
 * Usage: testgen [options]
 *
 * -i             : incremental matrices
 * -z             : compute Schur vectors and test them
 * -n size        : size of matrices
 * -l lower-bound : lower bound for matrix elements
 * -u upper-bound : upper bound for matrix elements
 * -c num         : number of matrices to solve
 */

#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>

typedef struct
{
  gsl_eigen_gen_workspace *gen_p;
  gsl_matrix *A;
  gsl_matrix *B;
  gsl_vector_complex *alpha;
  gsl_vector *beta;
  gsl_vector_complex *evals;

  gsl_matrix *Q;
  gsl_matrix *Z;
  int compute_schur;

  size_t n_evals;
} gen_workspace;

gen_workspace *gen_alloc(size_t n, int compute_schur);
void gen_free(gen_workspace *w);
int gen_proc(gen_workspace *w);

typedef struct
{
  gsl_matrix *A;
  gsl_matrix *B;
  gsl_matrix *Q;
  gsl_matrix *Z;
  int N;

  char jobvsl;
  char jobvsr;
  char sort;
  int selctg;
  int lda;
  int ldb;
  int sdim;
  double *alphar;
  double *alphai;
  gsl_vector *beta;
  int ldvsr;
  int lwork;
  int info;
  double *work;

  gsl_vector_complex *evals;
  gsl_vector_complex *alpha;
  size_t n_evals;
} lapack_workspace;

lapack_workspace *lapack_alloc(const size_t n);
void lapack_free(lapack_workspace *w);
int lapack_proc(lapack_workspace *w);

void dgges_(char *jobvsl, char *jobvsr, char *sort, int *selctg, int *n,
            double *a, int *lda, double *b, int *ldb, int *sdim,
            double *alphar, double *alphai, double *beta, double *vsl,
            int *ldvsl, double *vsr, int *ldvsr, double *work, int *lwork,
            int *bwork, int *info);

/*
 * Global variables
 */
unsigned long count = 0;

/*
 * Prototypes
 */

void make_random_matrix(gsl_matrix *m, gsl_rng *r, int lower, int upper);
void make_random_integer_matrix(gsl_matrix *m, gsl_rng *r, int lower,
                                int upper);
void make_start_matrix(gsl_matrix *m, int lower);
int inc_matrix (gsl_matrix *m, int lower, int upper);
void output_matrix(gsl_matrix *m);
void print_matrix(gsl_matrix *m, const char *str);
int test_evals(gsl_vector_complex *obs, gsl_vector_complex *expected,
               gsl_matrix *A, gsl_matrix *B,
               const char *obsname, const char *expname);
int test_alpha(gsl_vector_complex *obs, gsl_vector_complex *expected,
               gsl_matrix *A, gsl_matrix *B, const char *obsname,
               const char *expname);
int test_beta(gsl_vector *obs, gsl_vector *expected,
              gsl_matrix *A, gsl_matrix *B, const char *obsname,
              const char *expname);
void test_schur(gsl_matrix *A, gsl_matrix *S, gsl_matrix *Q, gsl_matrix *Z);
void print_vector(gsl_vector_complex *eval, const char *str);
int cmp(double a, double b);
int compare(const void *a, const void *b);
void sort_complex_vector(gsl_vector_complex *v);

gen_workspace *
gen_alloc(size_t n, int compute_schur)
{
  gen_workspace *w;

  w = (gen_workspace *) calloc(1, sizeof(gen_workspace));

  w->gen_p = gsl_eigen_gen_alloc(n);

  w->A = gsl_matrix_alloc(n, n);
  w->B = gsl_matrix_alloc(n, n);
  w->alpha = gsl_vector_complex_alloc(n);
  w->beta = gsl_vector_alloc(n);
  w->evals = gsl_vector_complex_alloc(n);
  w->compute_schur = compute_schur;

  if (compute_schur)
    {
      w->Q = gsl_matrix_alloc(n, n);
      w->Z = gsl_matrix_alloc(n, n);
      gsl_eigen_gen_params(1, 1, 0, w->gen_p);
    }

  return (w);
} /* gen_alloc() */

void
gen_free(gen_workspace *w)
{
  if (w->gen_p)
    gsl_eigen_gen_free(w->gen_p);

  if (w->A)
    gsl_matrix_free(w->A);

  if (w->B)
    gsl_matrix_free(w->B);

  if (w->alpha)
    gsl_vector_complex_free(w->alpha);

  if (w->beta)
    gsl_vector_free(w->beta);

  if (w->evals)
    gsl_vector_complex_free(w->evals);

  if (w->Q)
    gsl_matrix_free(w->Q);

  if (w->Z)
    gsl_matrix_free(w->Z);

  free(w);
}

int
gen_proc(gen_workspace *w)
{
  int s;

  s = gsl_eigen_gen_QZ(w->A, w->B, w->alpha, w->beta, w->Q, w->Z, w->gen_p);

  w->n_evals = w->gen_p->n_evals;

  return s;
} /* gen_proc() */

lapack_workspace *
lapack_alloc(const size_t n)
{
  lapack_workspace *w;
  double work[1];

  w = (lapack_workspace *) calloc(1, sizeof(lapack_workspace));

  w->A = gsl_matrix_alloc(n, n);
  w->B = gsl_matrix_alloc(n, n);
  w->Q = gsl_matrix_alloc(n, n);
  w->Z = gsl_matrix_alloc(n, n);
  w->alphar = malloc(n * sizeof(double));
  w->alphai = malloc(n * sizeof(double));
  w->beta = gsl_vector_alloc(n);
  w->alpha = gsl_vector_complex_alloc(n);
  w->evals = gsl_vector_complex_alloc(n);

  w->N = (int) n;
  w->n_evals = 0;

  w->jobvsl = 'N';
  w->jobvsr = 'N';
  w->sort = 'N';
  w->info = 0;

  w->lwork = -1;
  dgges_(&w->jobvsl,
         &w->jobvsr,
         &w->sort,
         (int *) 0,
         &w->N,
         w->A->data,
         (int *) &w->A->tda,
         w->B->data,
         (int *) &w->B->tda,
         &w->sdim,
         w->alphar,
         w->alphai,
         w->beta->data,
         w->Q->data,
         (int *) &w->Q->tda,
         w->Z->data,
         (int *) &w->Z->tda,
         work,
         &w->lwork,
         (int *) 0,
         &w->info);

  w->lwork = (int) work[0];
  w->work = malloc(w->lwork * sizeof(double));

  return (w);
} /* lapack_alloc() */

void
lapack_free(lapack_workspace *w)
{
  if (w->A)
    gsl_matrix_free(w->A);

  if (w->B)
    gsl_matrix_free(w->B);

  if (w->Q)
    gsl_matrix_free(w->Q);

  if (w->Z)
    gsl_matrix_free(w->Z);

  if (w->work)
    free(w->work);

  if (w->alphar)
    free(w->alphar);

  if (w->alphai)
    free(w->alphai);

  if (w->beta)
    gsl_vector_free(w->beta);

  if (w->alpha)
    gsl_vector_complex_free(w->alpha);

  if (w->evals)
    gsl_vector_complex_free(w->evals);

  free(w);
} /* lapack_free() */

int
lapack_proc(lapack_workspace *w)
{
  dgges_(&w->jobvsl,
         &w->jobvsr,
         &w->sort,
         (int *) 0,
         &w->N,
         w->A->data,
         (int *) &w->A->tda,
         w->B->data,
         (int *) &w->B->tda,
         &w->sdim,
         w->alphar,
         w->alphai,
         w->beta->data,
         w->Q->data,
         (int *) &w->Q->tda,
         w->Z->data,
         (int *) &w->Z->tda,
         w->work,
         &w->lwork,
         (int *) 0,
         &w->info);

  return (w->info);
} /* lapack_proc() */

/**********************************************
 * General routines
 **********************************************/

void
make_random_matrix(gsl_matrix *m, gsl_rng *r, int lower, int upper)
{
  size_t i, j;
  size_t N = m->size1;

  for (i = 0; i < N; ++i)
  {
    for (j = 0; j < N; ++j)
    {
      gsl_matrix_set(m,
                     i,
                     j,
                     gsl_rng_uniform(r) * (upper - lower) + lower);
    }
  }
} /* make_random_matrix() */

void
make_random_integer_matrix(gsl_matrix *m, gsl_rng *r, int lower, int upper)
{
  size_t i, j;
  size_t N = m->size1;

  for (i = 0; i < N; ++i)
  {
    for (j = 0; j < N; ++j)
    {
      double a = gsl_rng_uniform(r) * (upper - lower) + lower;
      gsl_matrix_set(m, i, j, floor(a));
    }
  }
} /* make_random_integer_matrix() */

void
make_start_matrix(gsl_matrix *m, int lower)

{
  size_t i, j;
  size_t N = m->size1;

  for (i = 0; i < N; ++i)
    for (j = 0; j < N; ++j)
      gsl_matrix_set(m, i, j, (double)lower);
} /* make_start_matrix() */

int
inc_matrix (gsl_matrix *m, int lower, int upper)
{
  size_t i = 0;
  size_t N = m->size1 * m->size2;
  int carry = 1;

  for (i = 0; carry > 0 && i < N; i++)
    {
      double v = m->data[i] + carry;
      carry = (v > upper) ? 1 : 0;
      m->data[i] = (v > upper) ? lower : v;
    }

  return carry;
} /* inc_matrix() */

void
output_matrix(gsl_matrix *m)
{
  size_t i, j;
  size_t N = m->size1;
  size_t M = m->size2;

  for (i = 0; i < N; ++i)
  {
    for (j = 0; j < M; ++j)
      {
        printf("%10.18e%s",
        /*printf("%10.18e%s",*/
               gsl_matrix_get(m, i, j),
               (j < M - 1) ? "," : ";\n");
      }
  }
}

void
print_matrix(gsl_matrix *m, const char *str)

{
  size_t i, j;
  size_t N = m->size1;
  size_t M = m->size2;
  gsl_matrix_view v;
  size_t rows, cols;
  size_t r, c;
  char buf[100];

  /*print_octave(m, str);
  return;*/

  /*rows = GSL_MIN(15, N);*/
  rows = N;
  cols = GSL_MIN(15, N);
  /*cols = N;*/

  for (i = 0; i < N; i += rows)
  {
    for (j = 0; j < M; j += cols)
    {
      r = GSL_MIN(rows, N - i);
      c = GSL_MIN(cols, N - j);

      v = gsl_matrix_submatrix(m, i, j, r, c);

      sprintf(buf, "%s(%u:%u,%u:%u)",
              str,
              i + 1,
              i + r,
              j + 1,
              j + c);

      printf("%s = [\n", buf);

      output_matrix(&v.matrix);

      printf("]\n");
    }
  }
} /* print_matrix() */

int
test_evals(gsl_vector_complex *obs, gsl_vector_complex *expected,
           gsl_matrix *A, gsl_matrix *B, const char *obsname,
           const char *expname)
{
  size_t N = expected->size;
  size_t i, k;
  double max, max_abserr, max_relerr;

  max = 0.0;
  max_abserr = 0.0;
  max_relerr = 0.0;
  k = 0;

  for (i = 0; i < N; ++i)
    {
      gsl_complex z = gsl_vector_complex_get(expected, i);
      max = GSL_MAX_DBL(max, gsl_complex_abs(z));
    }

  for (i = 0; i < N; ++i)
    {
      gsl_complex z_obs = gsl_vector_complex_get(obs, i);
      gsl_complex z_exp = gsl_vector_complex_get(expected, i);

      double x_obs = GSL_REAL(z_obs);
      double y_obs = GSL_IMAG(z_obs);
      double x_exp = GSL_REAL(z_exp);
      double y_exp = GSL_IMAG(z_exp);

      double abserr_x = fabs(x_obs - x_exp);
      double abserr_y = fabs(y_obs - y_exp);
      double noise = max * GSL_DBL_EPSILON * N * N;

      max_abserr = GSL_MAX_DBL(max_abserr, abserr_x + abserr_y);

      if (abserr_x < noise && abserr_y < noise)
        continue;

      if (abserr_x > 1.0e-6 || abserr_y > 1.0e-6)
        ++k;
    }

    if (k)
      {
        printf("==== CASE %lu ===========================\n\n", count);

        print_matrix(A, "A");
        print_matrix(B, "B");

        printf("=== eval - %s ===\n", expname);
        print_vector(expected, expname);

        printf("=== eval - %s ===\n", obsname);
        print_vector(obs, obsname);

        printf("max abserr = %g  max relerr = %g\n", max_abserr, max_relerr);

        printf("=========================================\n\n");
      }

    return k;
} /* test_evals() */

int
test_alpha(gsl_vector_complex *obs, gsl_vector_complex *expected,
           gsl_matrix *A, gsl_matrix *B, const char *obsname,
           const char *expname)
{
  size_t N = expected->size;
  size_t i, k;
  double max, max_abserr, max_relerr;

  max = 0.0;
  max_abserr = 0.0;
  max_relerr = 0.0;
  k = 0;

  for (i = 0; i < N; ++i)
    {
      gsl_complex z = gsl_vector_complex_get(expected, i);
      max = GSL_MAX_DBL(max, gsl_complex_abs(z));
    }

  for (i = 0; i < N; ++i)
    {
      gsl_complex z_obs = gsl_vector_complex_get(obs, i);
      gsl_complex z_exp = gsl_vector_complex_get(expected, i);

      double x_obs = GSL_REAL(z_obs);
      double y_obs = fabs(GSL_IMAG(z_obs));
      double x_exp = GSL_REAL(z_exp);
      double y_exp = fabs(GSL_IMAG(z_exp));

      double abserr_x = fabs(x_obs - x_exp);
      double abserr_y = fabs(y_obs - y_exp);
      double noise = max * GSL_DBL_EPSILON * N * N;

      max_abserr = GSL_MAX_DBL(max_abserr, abserr_x + abserr_y);

      if (abserr_x < noise && abserr_y < noise)
        continue;

      if (abserr_x > 1.0e-6 || abserr_y > 1.0e-6)
        ++k;
    }

    if (k)
      {
        printf("==== CASE %lu ===========================\n\n", count);

        print_matrix(A, "A");
        print_matrix(B, "B");

        printf("=== alpha - %s ===\n", expname);
        print_vector(expected, expname);

        printf("=== alpha - %s ===\n", obsname);
        print_vector(obs, obsname);

        printf("max abserr = %g  max relerr = %g\n", max_abserr, max_relerr);

        printf("=========================================\n\n");
      }

    return k;
} /* test_alpha() */

int
test_beta(gsl_vector *obs, gsl_vector *expected,
          gsl_matrix *A, gsl_matrix *B, const char *obsname,
          const char *expname)
{
  size_t N = expected->size;
  size_t i, k;
  double max, max_abserr, max_relerr;

  max = 0.0;
  max_abserr = 0.0;
  max_relerr = 0.0;
  k = 0;

  for (i = 0; i < N; ++i)
    {
      double z = gsl_vector_get(expected, i);
      max = GSL_MAX_DBL(max, fabs(z));
    }

  for (i = 0; i < N; ++i)
    {
      double v_obs = gsl_vector_get(obs, i);
      double v_exp = gsl_vector_get(expected, i);
      double abserr = fabs(v_obs - v_exp);
      double noise = max * GSL_DBL_EPSILON * N * N;

      max_abserr = GSL_MAX_DBL(max_abserr, abserr);

      if (abserr < noise)
        continue;

      if (abserr > 1.0e-6)
        ++k;
    }

    if (k)
      {
        printf("==== CASE %lu ===========================\n\n", count);

        print_matrix(A, "A");
        print_matrix(B, "B");

        printf("=== beta - %s ===\n", expname);
        printf("%s = [\n", expname);
        gsl_vector_fprintf(stdout, expected, "%.12e");
        printf("]\n");

        printf("=== beta - %s ===\n", obsname);
        printf("%s = [\n", obsname);
        gsl_vector_fprintf(stdout, obs, "%.12e");
        printf("]\n");

        printf("max abserr = %g  max relerr = %g\n", max_abserr, max_relerr);

        printf("=========================================\n\n");
      }

    return k;
} /* test_beta() */

/* test if A = Q S Z^t */
void
test_schur(gsl_matrix *A, gsl_matrix *S, gsl_matrix *Q, gsl_matrix *Z)
{
  const size_t N = A->size1;
  gsl_matrix *T1, *T2;
  size_t i, j, k;
  double lhs, rhs;
  double abserr;

  T1 = gsl_matrix_alloc(N, N);
  T2 = gsl_matrix_alloc(N, N);

  /* compute T1 = S Z^t */
  gsl_blas_dgemm(CblasNoTrans,
                 CblasTrans,
                 1.0,
                 S,
                 Z,
                 0.0,
                 T1);

  /* compute T2 = Q T1 = Q S Z^t */
  gsl_blas_dgemm(CblasNoTrans,
                 CblasNoTrans,
                 1.0,
                 Q,
                 T1,
                 0.0,
                 T2);

  k = 0;
  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          lhs = gsl_matrix_get(A, i, j);
          rhs = gsl_matrix_get(T2, i, j);

          abserr = fabs(lhs - rhs);

          if (abserr > 1.0e-6)
            ++k;
        }
    }

  if (k)
    {
      printf("==== CASE %lu ===========================\n\n", count);

      print_matrix(A, "A");

      printf("=== Schur Form matrix ===\n");
      print_matrix(S, "S");

      printf("=== Left Schur matrix ===\n");
      print_matrix(Q, "Q");

      printf("=== Right Schur matrix ===\n");
      print_matrix(Z, "Z");

      printf("=== Q S Z^t ===\n");
      print_matrix(T1, "Q S Z^t");

      printf("=== A - Q S Z^t ===\n");
      gsl_matrix_sub(T2, A);
      print_matrix(T1, "A - Q S Z^t");

      printf("=========================================\n\n");
    }

  gsl_matrix_free(T1);
  gsl_matrix_free(T2);
} /* test_schur() */

void
print_vector(gsl_vector_complex *eval, const char *str)
{
  size_t N = eval->size;
  size_t i;
  gsl_complex z;

  printf("%s = [\n", str);

  for (i = 0; i < N; ++i)
    {
      z = gsl_vector_complex_get(eval, i);
      printf("%.18e %.18e;\n", GSL_REAL(z), GSL_IMAG(z));
    }

  printf("]\n");
} /* print_vector() */

int
cmp(double a, double b)
{
  return ((a > b) ? 1 : ((a < b) ? -1 : 0));
} /* cmp() */

int
compare(const void *a, const void *b)
{
  const double *x = a;
  const double *y = b;
  int r1 = cmp(y[0], x[0]);
  int r2 = cmp(y[1], x[1]);

  if (!gsl_finite(x[0]))
    return 1;
  if (!gsl_finite(y[0]))
    return -1;

  if (fabs(x[0] - y[0]) < 1.0e-8)
    {
      /* real parts are very close to each other */
      return r2;
    }
  else
    {
      return r1 ? r1 : r2;
    }
} /* compare() */

void
sort_complex_vector(gsl_vector_complex *v)
{
  qsort(v->data, v->size, 2 * sizeof(double), &compare);
} /* sort_complex_vector() */

int
main(int argc, char *argv[])
{
  gen_workspace *gen_workspace_p;
  lapack_workspace *lapack_workspace_p;
  size_t N;
  int c;
  int lower;
  int upper;
  int incremental;
  size_t nmat;
  gsl_matrix *A, *B;
  gsl_rng *r;
  int s;
  int compute_schur;
  size_t i;

  gsl_ieee_env_setup();
  gsl_rng_env_setup();

  N = 30;
  lower = -10;
  upper = 10;
  incremental = 0;
  nmat = 0;
  compute_schur = 0;

  while ((c = getopt(argc, argv, "ic:n:l:u:z")) != (-1))
    {
      switch (c)
        {
          case 'i':
            incremental = 1;
            break;

          case 'n':
            N = strtol(optarg, NULL, 0);
            break;

          case 'l':
            lower = strtol(optarg, NULL, 0);
            break;

          case 'u':
            upper = strtol(optarg, NULL, 0);
            break;

          case 'c':
            nmat = strtoul(optarg, NULL, 0);
            break;

          case 'z':
            compute_schur = 1;
            break;

          case '?':
          default:
            printf("usage: %s [-i] [-z] [-n size] [-l lower-bound] [-u upper-bound] [-c num]\n", argv[0]);
            exit(1);
            break;
        } /* switch (c) */
    }

  A = gsl_matrix_alloc(N, N);
  B = gsl_matrix_alloc(N, N);
  gen_workspace_p = gen_alloc(N, compute_schur);
  lapack_workspace_p = lapack_alloc(N);

  r = gsl_rng_alloc(gsl_rng_default);

  if (incremental)
    {
      make_start_matrix(A, lower);

      /* we need B to be non-singular */
      make_random_integer_matrix(B, r, lower, upper);
    }

  fprintf(stderr, "testing N = %d", N);
  if (incremental)
    fprintf(stderr, " incrementally");
  else
    fprintf(stderr, " randomly");

  fprintf(stderr, " on element range [%d, %d]", lower, upper);

  if (compute_schur)
    fprintf(stderr, ", with Schur vectors");

  fprintf(stderr, "\n");

  while (1)
    {
      if (nmat && (count >= nmat))
        break;

      ++count;

      if (!incremental)
        {
          make_random_matrix(A, r, lower, upper);
          make_random_matrix(B, r, lower, upper);
        }
      else
        {
          s = inc_matrix(A, lower, upper);
          if (s)
            break; /* all done */

          make_random_integer_matrix(B, r, lower, upper);
        }

      /*if (count != 89120)
        continue;*/

      /* make copies of matrices */
      gsl_matrix_memcpy(gen_workspace_p->A, A);
      gsl_matrix_memcpy(gen_workspace_p->B, B);
      gsl_matrix_transpose_memcpy(lapack_workspace_p->A, A);
      gsl_matrix_transpose_memcpy(lapack_workspace_p->B, B);

      /* compute eigenvalues with LAPACK */
      s = lapack_proc(lapack_workspace_p);

      if (s != GSL_SUCCESS)
        {
          printf("LAPACK failed, case %lu\n", count);
          exit(1);
        }

#if 0
      print_matrix(A, "A");
      print_matrix(B, "B");
      gsl_matrix_transpose(lapack_workspace_p->A);
      gsl_matrix_transpose(lapack_workspace_p->B);
      print_matrix(lapack_workspace_p->A, "S_lapack");
      print_matrix(lapack_workspace_p->B, "T_lapack");
#endif

      /* compute eigenvalues with GSL */
      s = gen_proc(gen_workspace_p);

      if (s != GSL_SUCCESS)
        {
          printf("=========== CASE %lu ============\n", count);
          printf("Failed to converge: found %u eigenvalues\n",
                 gen_workspace_p->n_evals);
          print_matrix(A, "A");
          print_matrix(B, "B");
          print_matrix(gen_workspace_p->A, "Af");
          print_matrix(gen_workspace_p->B, "Bf");
          print_matrix(lapack_workspace_p->A, "Ae");
          print_matrix(lapack_workspace_p->B, "Be");
          exit(1);
        }

#if 0
      print_matrix(gen_workspace_p->A, "S_gsl");
      print_matrix(gen_workspace_p->B, "T_gsl");
#endif

      /* compute alpha / beta vectors */
      for (i = 0; i < N; ++i)
        {
          double beta;
          gsl_complex alpha, z;

          beta = gsl_vector_get(gen_workspace_p->beta, i);
          if (beta == 0.0)
            GSL_SET_COMPLEX(&z, GSL_POSINF, GSL_POSINF);
          else
            {
              alpha = gsl_vector_complex_get(gen_workspace_p->alpha, i);
              z = gsl_complex_div_real(alpha, beta);
            }

          gsl_vector_complex_set(gen_workspace_p->evals, i, z);

          beta = gsl_vector_get(lapack_workspace_p->beta, i);
          GSL_SET_COMPLEX(&alpha,
                          lapack_workspace_p->alphar[i],
                          lapack_workspace_p->alphai[i]);

          if (beta == 0.0)
            GSL_SET_COMPLEX(&z, GSL_POSINF, GSL_POSINF);
          else
            z = gsl_complex_div_real(alpha, beta);

          gsl_vector_complex_set(lapack_workspace_p->evals, i, z);
          gsl_vector_complex_set(lapack_workspace_p->alpha, i, alpha);
        }

#if 0
      gsl_sort_vector(gen_workspace_p->beta);
      gsl_sort_vector(lapack_workspace_p->beta);
      sort_complex_vector(gen_workspace_p->alpha);
      sort_complex_vector(lapack_workspace_p->alpha);

      s = test_alpha(gen_workspace_p->alpha,
                     lapack_workspace_p->alpha,
                     A,
                     B,
                     "gen",
                     "lapack");
      s = test_beta(gen_workspace_p->beta,
                    lapack_workspace_p->beta,
                    A,
                    B,
                    "gen",
                    "lapack");
#endif
#if 1
      sort_complex_vector(gen_workspace_p->evals);
      sort_complex_vector(lapack_workspace_p->evals);

      s = test_evals(gen_workspace_p->evals,
                     lapack_workspace_p->evals,
                     A,
                     B,
                     "gen",
                     "lapack");
#endif

      if (compute_schur)
        {
          test_schur(A,
                     gen_workspace_p->A,
                     gen_workspace_p->Q,
                     gen_workspace_p->Z);
          test_schur(B,
                     gen_workspace_p->B,
                     gen_workspace_p->Q,
                     gen_workspace_p->Z);
        }
    }

  gsl_matrix_free(A);
  gsl_matrix_free(B);
  gen_free(gen_workspace_p);
  lapack_free(lapack_workspace_p);

  if (r)
    gsl_rng_free(r);

  return 0;
} /* main() */

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:
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* Program fwhm_by_modsq
* to compute the Full Width at Half Maximum of an equivalent long exposure
* with the power spectrum (square modulus of the FFT)
* (by fitting a Kolmogorof law)
*
* JLP
* Version 11-05-99
-------------------------------------------------------------------*/
/*
#define DEBUG
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <fcntl.h>
#include <math.h>
#include <fftw.h>
#include <jlp_ftoc.h>

static int fit_kolmo_center(float *in_image, int nx, int ny,
                        float *sigx, float *sigy);
static int aris_r0calc(float *modsq, int nx, int ny, float *r0);

main(argc,argv)
int argc;
char *argv[];
{
char in_name[61], comments[81];
float *in_image, *testi;
INT_PNTR pntr_image;
INT4 istatus, nx, ny;
float sigx, sigy;
register int i, j;

printf(" Program FWHM_BY_MODSQ to compute the Full Width at Half Maximum from the power spectrum)\n");
printf(" JLP Version 11-05-99 \n");

/* One parameters only is allowed to run the program: */
/* Carefull: 7 parameters always, using JLP "runs" */
if(argc == 7 && *argv[3]) argc = 4;
if(argc == 7 && *argv[2]) argc = 3;
if(argc == 7 && *argv[1]) argc = 2;
if(argc != 2 && argc != 1)
  {
  printf(" Syntax: fwhm_by_modsq filename\n");
  printf(" Fatal: Syntax error: argc=%d\n",argc);
  exit(-1);
  }

/* Interactive input of parameters: */
if (argc == 2 )
 {
  strcpy(in_name,argv[1]);
 }
else
 {
  printf(" Input file := ");scanf("%s",in_name);
 }

/**********************************************************/
  JLP_BEGIN();
  JLP_INQUIFMT();
  istatus = JLP_VM_READIMAG1(&pntr_image,&nx,&ny,in_name,comments);
  if(istatus != 0) exit(-1);

  in_image = (float *) pntr_image;
/* Computing the power spectrum (for debug purpose)
  testi = (float *)malloc(nx * ny * sizeof(float));
  for(i = 0; i < nx * ny; i++) testi[i] = 0.;

  fftw_float(in_image,testi,nx,ny,1);
  for(i = 0; i < nx * ny; i++) in_image[i] = in_image[i] * in_image[i]
                                             + testi[i] * testi[i];
*/

  fit_kolmo_center(in_image, nx, ny, &sigx, &sigy);

JLP_END();
}
/*****************************************************************
* To fit a Kolmogorof law to f1(xx,yy)
*
* f1: array of values f1(xx[i],yy[i])
* xx, yy: arrays of coordinates x, y
* npts: number of points
*
* Parameters:
* sigx, sigy, xc, yc, rho
* error[sigx,sigy,xc,yc,rho]
*
* ifail = 0 if correct
*         -3: all values are negative or null!
*****************************************************************/
static int fit_kolmo_center(float *in_image, INT4 nx, INT4 ny,
                               float *sigx, float *sigy)
{
float *ff2, *gg;
INT4 istatus, npts, ifail, nx1, ny1, istart, jstart, iw, jw;
float aa, bb, r0, fwhm;
register int i, j;

nx1 = 20; ny1 = 20;
istart = nx/2 - nx1/2;
jstart = ny/2 - ny1/2;
ff2 = (float *) malloc(nx1 * ny1 * sizeof(float));
gg = (float *) malloc(nx1 * ny1 * sizeof(float));
if(ff2 == NULL || gg == NULL)
  {
  printf("FWHM/Fatal error alocating memory space \n");
  exit(-1);
  }
#ifdef DEBUG
printf(" istart=%d jstart=%d nx1=%d ny1=%d\n",istart,jstart,nx1,ny1);
#endif
for(j = 0; j < ny1; j++)
  for(i = 0; i < nx1; i++)
   {
   iw = i + istart;
   jw = j + jstart;
   ff2[i + j * nx1] = (iw - nx/2) * (iw - nx/2) + (jw - ny/2) * (jw - ny/2);
   gg[i + j * nx1] = in_image[iw + jw * nx];
   }
npts = nx1 * ny1;

jlp_fit_kolmogorof(ff2, gg, &npts, &aa, &bb, &ifail);

/* JLP99: I notice that the ratio between long exposure determination
* and this program is about 15: */
fwhm = pow((double)(bb/6.88),-0.6)/15.;
/* fwhm = lambda /r0 */
r0 = 1./fwhm;
#ifdef DEBUG
printf(" Full Kolmogorov fit: aa=%f bb=%f r0=%f fwhm~%f (10 mm)\n",
         aa,bb,r0,fwhm);
#else
printf(" Full Kolmogorov fit: r0=%f fwhm~%f (10 mm)\n",r0,fwhm);
#endif

aris_r0calc(in_image, nx, ny, &r0);
/* To make it compatible with the previous fit: */
r0 *= 4.;
printf(" Eric Aristidi approximation: r0=%f fwhm=%f\n",r0,1./r0);

return(0);
}

/***************************************************************
* fit_kolmogorof
* To fit a Kolmogorof law to Log_intensity(|x|^2) as a function of
* in x^2 + y^2
*
*   g(f) = rho * exp [ -6.88 * (lambda * f / r_0)^5/3 ]
*
*   g(f) = rho * exp [ - b * f^5/3 ]    with b = 6.88 * (lambda / r_0)^5/3
*
*  Log_Intensity(x,y)
*    G(f) = a - b f^(5/3)    with f^2 = x^2 + y^2
*
*  where:
*  The problem is to find the coefficients (a, b)
*  which minimize the sums:
*   SUM on all the selected disk (x, y) of ( Log_intensity - G(x,y) ) **2
*
* The normal equation can be written as:
*
*    n a - b SUM f^5/3 = SUM Log_g
*    a SUM f^5/3 - b SUM f^10/3 = SUM f^5/3 Log_g
*
* Hence:
*    a = 1/n * (SUM Log_g + b SUM f^5/3)
*    b = (SUM f^5/3 Log_g  - 1/n SUM Log_g SUM f^5/3)
*          / ( 1/n * SUM f^5/3 * SUM f^5/3 - SUM f^10/3)
*
*   where n is the number of points (SUM 1)
*
*
* gg: array of values gg(xx[i],yy[i])
*
* ff2: arrays of x^2 + y^2
* npts: number of points
*
* ifail = 0 if correct
*         -3: all values are negative or null!
*****************************************************************/
int jlp_fit_kolmogorof(float *ff2, float *gg, INT4 *npts,
                       float *aa, float *bb, INT4 *ifail)
{
/* gg measured intensities */
/* log_g log of measured intensities */
double *log_g, *ff;
double sum_g, sum_f1, sum_f2, sum_gf, w1;
register int i, j, k;

*ifail = 0;

 if((log_g = (double *) malloc(*npts * sizeof(double))) == NULL ||
     (ff = (double *) malloc(*npts * sizeof(double))) == NULL )
  {
  printf("jlp_fit_gauss/Error allocating memory for array (npts=%d)\n",*npts);
  *ifail = -1;
  return(-1);
  }

/* Transfer to double precision arrays and conversion of intensity to Log : */
i = 0;
for(k = 0; k < *npts; k++)
  {
       if(gg[k] > 0)
         {
          log_g[i] = log((double)gg[k]);
          ff[i] = ff2[k];
          i++;
         }
  }
*npts = i;

if(*npts == 0)
  {
  printf("jlp_fit_gauss/All input values are null!\n");
  *ifail = -1;
  return(-1);
  }

/* Compute the sums:  5/3 = 1.667 (but 5/6=0.833 since x^2+y^2)*/
/* Compute the sums:  10/3 = 0.833 (but 5/3=1.667 since x^2+y^2)*/
sum_f1 = 0.; sum_f2 = 0.;
sum_g = 0.; sum_gf = 0.;
for(k = 0; k < *npts; k++)
  {
  w1 = pow(ff[k],0.833);
  sum_f1 += w1;
  sum_f2 += pow(ff[k],1.667);
  sum_g  += log_g[k];
  sum_gf += w1 * log_g[k];
  }

#ifdef DEBUG
printf("sum_f1=%f sum_f2=%f sum_g=%f sum_gf=%f npts=%d\n",
        sum_f1,sum_f2,sum_g,sum_gf,*npts);
#endif
/* Result:
*    a = 1/n * (SUM Log_g + b SUM f^5/3)
*    b = (SUM f^5/3 Log_g  - 1/n SUM Log_g SUM f^5/3)
*          / ( 1/n * SUM f^5/3 * SUM f^5/3 - SUM f^10/3)
*/
  *bb = (sum_gf - (1. / *npts) * sum_g * sum_f1)
           / ( (1. / *npts) * sum_f1 * sum_f1 - sum_f2);
  *aa = (1. / *npts) * (sum_g + *bb * sum_f1);

return(0);
}
/**********************************************************/
/* Eric Aristidi's fit (approximation only) */
static int aris_r0calc(float *modsq, int nx, int ny, float *r0)
{
int     f12_1, f12_2, dim, ixc, iyc;
float   pix,lambda,x,y1,y2,r01,r02;
dim = nx;
ixc = nx/2;
iyc = ny/2;

 x=modsq[ixc + iyc*dim];

 f12_1=4;
 y1=modsq[(ixc+f12_1) + iyc * dim] + modsq[ixc + (iyc+f12_1) * dim];
 y1/=2.;

 f12_2=3;
 y2=modsq[(ixc+f12_2) + iyc * dim] + modsq[ixc + (iyc+f12_2) * dim];
 y2/=2.;

//  Compute r0
//  ---------
 pix=0.0119/206265.;lambda=0.00000065;
 r01=pow(3.44/log(x/y1),0.6)*lambda*(float)f12_1/pix/dim;
 r02=pow(3.44/log(x/y2),0.6)*lambda*(float)f12_2/pix/dim;
     *r0=(r02+r01)/2.;
}

max_stars_repo_name: jlprieur/sourcc

max_issues_repo_name: jlprieur/sourcc

max_forks_repo_name: jlprieur/sourcc


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_eigen.h>

#include "chisq_dist.h"

int chisq_dist(int mode, int do_chisq, int nophotoerr, int ncalc, int ncol,
	       double *covmat, double *c, double *slope,
	       double *pivotmag, double *refmag, double *refmagerr, double *magerr,
	       double *color, double *lupcorr, double *dist, double sigint) {

  //int i,j,k;
  int i,j,k;
  int nmag = ncol+1;

  //gsl_matrix *mcovmat;
  gsl_matrix_view mvcovmat;
  gsl_matrix *cobs;
  gsl_matrix *cimat;
  gsl_matrix *cobsmat;
  gsl_matrix *cobstemp;
  gsl_matrix *mmetric;
  gsl_matrix *mrotmat;
  gsl_vector *vdcm;
  gsl_permutation *pp;
  gsl_vector *vdc;
  int s;
  double *covmat_temp = NULL;

  int covmat_stride;
  int test;
  double val;

  double chisq,norm;

  int use_refmagerr=0;

  norm = 1.0;

  covmat_stride = ncol*ncol; //*sizeof(double);

  mrotmat = gsl_matrix_alloc(ncol,nmag);
  gsl_matrix_set_zero(mrotmat);  // important!
  for (i=0;i<ncol;i++) {
    gsl_matrix_set(mrotmat, i, i, 1.0);
    gsl_matrix_set(mrotmat, i, i+1, -1.0);
  }
  cobs = gsl_matrix_alloc(nmag, nmag);
  cobsmat = gsl_matrix_alloc(ncol,ncol);
  cobstemp = gsl_matrix_alloc(ncol,nmag);
  pp=gsl_permutation_alloc(ncol);
  mmetric=gsl_matrix_alloc(ncol,ncol);
  vdc = gsl_vector_alloc(ncol);
  vdcm=gsl_vector_alloc(ncol);

  cimat = NULL;
  if (refmagerr != NULL) {
      use_refmagerr = 1;
      cimat = gsl_matrix_alloc(ncol,ncol);
  }

  if (mode == 0) {
    // mode = 0
    //   Many galaxies, one redshift
    //     - refmag/refmagerr is an array with ncalc (==ngal)
    //     - color is a matrix with ncol x ncalc (==ngal)
    //     - magerr is a matrix with nmag x ncalc (==ngal)
    //     - c is an array with ncol
    //     - slope is an array with ncol
    //     - pivotmag is a single value
    //     - lupcorr is a matrix with ncol x ncalc (==ngal)
    //     //- covmat is a matrix with ncol x ncol x ncalc (==ngal)
    //     - covmat is a matrix with ncol x ncol
    //
    // // covmat[ncol,ncol,ncalc]: (k*ncol+j)*ncol + i
    // covmat[ncol,ncol]: k*ncol + j
    // c[ncol]: j
    // slope[ncol]: j
    // pivotmag[0]
    // refmag[ncalc]: i
    // refmagerr[ncalc]: i
    // magerr[nmag,ncalc]: i*nmag + j
    // color[ncol,ncalc]: i*ncol + j
    // lupcorr[ncol,ncalc]: i*ncol + j

    // first make a buffer for the local copy of the covmat (which gets overwritten)
    if ((covmat_temp = (double *)calloc(ncol*ncol, sizeof(double))) == NULL) {
	return -1;
    }

    for (i=0;i<ncalc;i++) {

      memcpy(covmat_temp,covmat,sizeof(double)*ncol*ncol);

      gsl_matrix_set_identity(cobs);
      for (j=0;j<nmag;j++) {
	gsl_matrix_set(cobs, j, j, magerr[i*nmag+j]*magerr[i*nmag+j]);
      }

      gsl_matrix_set_zero(cobstemp);

      gsl_blas_dgemm(CblasNoTrans, CblasTrans,
		     1.0, mrotmat, cobs,
		     0.0, cobstemp);
      gsl_blas_dgemm(CblasNoTrans, CblasTrans,
		     1.0, mrotmat, cobstemp,
		     0.0, cobsmat);

      //mvcovmat = gsl_matrix_view_array(&covmat[covmat_stride*i], ncol, ncol);
      mvcovmat = gsl_matrix_view_array(covmat_temp, ncol, ncol);

      // and the ci matrix
      if (use_refmagerr) {
	  // make the C_i matrix for the refmag err
	  //  this is a matrix with ncol x ncol size
	  gsl_matrix_set_zero(cimat);
	  for (j=0;j<ncol;j++) {
	      for (k=j;k<ncol;k++) {
		  val = slope[j] * slope[k] * refmagerr[i] * refmagerr[i];
		  gsl_matrix_set(cimat, j, k, val);
		  if (k != j) {
		      gsl_matrix_set(cimat, k, j, val);
		  }
	      }
	  }
      }

      // check sigint
      test = 1;
      for (j=0;j<ncol;j++) {
	if (gsl_matrix_get(&mvcovmat.matrix, j, j) < sigint*sigint) {
	  test = 0;
	}
      }

      if (test == 0) {
	if (do_chisq) {
	  dist[i] = 1e11;
	} else {
	  dist[i] = -1e11;
	}
      } else {
	if (!nophotoerr) {
	  gsl_matrix_add(&mvcovmat.matrix, cobsmat);
	}

	if (use_refmagerr) {
	    gsl_matrix_add(&mvcovmat.matrix, cimat);
	}

	// check and fix the matrix if necessary
	check_and_fix_covmat(&mvcovmat.matrix);

	gsl_linalg_LU_decomp(&mvcovmat.matrix, pp, &s);
	gsl_linalg_LU_invert(&mvcovmat.matrix, pp, mmetric);

	if (!do_chisq) {
	  // need the determinant
	  norm = gsl_linalg_LU_det(&mvcovmat.matrix, s);
	}

	// now need the slope, etc.
	for (j=0;j<ncol;j++) {
	  gsl_vector_set(vdc,j,
			 (c[j]+slope[j]*(refmag[i]-pivotmag[0])) + lupcorr[i*ncol+j] -
			 color[i*ncol+j]);
	}

	gsl_blas_dgemv(CblasNoTrans, 1.0, mmetric, vdc, 0.0, vdcm);
	gsl_blas_ddot(vdcm, vdc, &chisq);

	if (do_chisq) {
	  dist[i] = chisq;
	} else {
	  dist[i]=-0.5*chisq-0.5*log(norm);
	}
      }

    }
    free(covmat_temp);

  } else if (mode==1) {
    // mode = 1
    //   One galaxy, many redshifts
    //     - refmag/refmagerr is a single value
    //     - color is an array with ncol values
    //     - magerr is an array with nmag=ncol+1 values
    //     - c is a matrix with ncol x ncalc (==nz)
    //     - slope is a matrix with ncol x ncalc (==nz)
    //     - pivotmag is an array with ncalc (==nz)
    //     - lupcorr is a matrix with ncol x ncalc (==nz)
    //     - covmat is a matrix with ncol x ncol x ncalc (==nz)
    //
    // covmat[ncol,ncol,ncalc] : (k*ncol+j)*ncol+i
    // c[ncol,ncalc]: i*ncol+j
    // slope[ncol,ncalc]: i*ncol+j
    // pivotmag[ncalc]: i
    // magerr[nmag]: j
    // color[ncol]: j
    // lupcorr[ncol,ncalc]: i*ncol + j
    // refmag[0]
    // refmagerr[0]

    // first make a buffer for the local copy of the covmat (which gets overwritten)
    if ((covmat_temp = (double *)calloc(ncol*ncol*ncalc, sizeof(double))) == NULL) {
	return -1;
    }
    memcpy(covmat_temp,covmat,sizeof(double)*ncol*ncol*ncalc);

    // We can generate the C_obs matrix (cobsmat) just once
    gsl_matrix_set_identity(cobs);
    for (j=0;j<nmag;j++) {
      gsl_matrix_set(cobs,j,j,magerr[j]*magerr[j]);
    }

    gsl_matrix_set_zero(cobstemp);

    gsl_blas_dgemm(CblasNoTrans, CblasTrans,
		   1.0, mrotmat, cobs,
		   0.0, cobstemp);
    gsl_blas_dgemm(CblasNoTrans, CblasTrans,
		   1.0, mrotmat, cobstemp,
		   0.0, cobsmat);

    for (i=0;i<ncalc;i++) {
      mvcovmat = gsl_matrix_view_array(&covmat_temp[covmat_stride*i], ncol, ncol);

      if (use_refmagerr) {
	  // make the C_i matrix for the refmag err
	  //  this is a matrix with ncol x ncol size...
	  gsl_matrix_set_zero(cimat);
	  for (j=0;j<ncol;j++) {
	      for (k=j;k<ncol;k++) {
		  val = slope[i*ncol+j] * slope[i*ncol+k] * refmagerr[0] * refmagerr[0];
		  gsl_matrix_set(cimat, j, k, val);
		  if (k != j) {
		      gsl_matrix_set(cimat, k, j, val);
		  }
	      }
	  }
      }

      // check sigint
      test = 1;
      for (j=0;j<ncol;j++) {
	if (gsl_matrix_get(&mvcovmat.matrix, j, j) < sigint*sigint) {
	  test = 0;
	}
      }

      if (test == 0) {
	if (do_chisq) {
	  dist[i] = 1e11;
	} else {
	  dist[i] = -1e11;
	}
      } else {
	if (!nophotoerr) {
	  gsl_matrix_add(&mvcovmat.matrix, cobsmat);
	}

	if (use_refmagerr) {
	    gsl_matrix_add(&mvcovmat.matrix, cimat);
	}

	// check and fix the matrix if necessary
	check_and_fix_covmat(&mvcovmat.matrix);

	gsl_linalg_LU_decomp(&mvcovmat.matrix, pp, &s);
	gsl_linalg_LU_invert(&mvcovmat.matrix, pp, mmetric);

	if (!do_chisq) {
	  norm = gsl_linalg_LU_det(&mvcovmat.matrix, s);
	}

	for (j=0;j<ncol;j++) {
	  gsl_vector_set(vdc,j,
			 (c[i*ncol+j] + slope[i*ncol+j]*(refmag[0] - pivotmag[i])) + lupcorr[i*ncol+j] -
			 color[j]);
	}

	gsl_blas_dgemv(CblasNoTrans, 1.0, mmetric, vdc, 0.0, vdcm);
	gsl_blas_ddot(vdcm,vdc, &chisq);

	if (do_chisq) {
	  dist[i] = chisq;
	} else {
	  dist[i]=-0.5*chisq-0.5*log(norm);
	}

      }
    }
    free(covmat_temp);
  } else {
      // mode = 2
      //  Many galaxies, many redshifts
      //     - refmag/refmagerr is an array with ncalc (==ngal)
      //     - color is a matrix with ncol x ncalc (==ngal)
      //     - magerr is a matrix with nmag x ncalc (==ngal)
      //     - c is a matrix with ncol x ncalc (==ngal)
      //     - slope is a matrix with ncol x ncalc (==ngal)
      //     - pivotmag is an array with ncalc (==ngal)
      //     - lupcorr is a matrix with ncol x ncalc (==ngal)
      //     - covmat is a matrix with ncol x ncol x ncalc (==ngal)
      //
      // covmat[ncol,ncol,ncalc] : (k*ncol+j)*ncol+i   -> mode 0
      // c[ncol,ncalc]: i*ncol+j                       -> mode 1
      // slope[ncol,ncalc]: i*ncol+j                   -> mode 1
      // pivotmag[ncalc]: i                              -> mode 1
      // refmag[ncalc]: i                                -> mode 0
      // refmagerr[ncalc]: i                             -> mode 0
      // magerr[nmag,ncalc]: i*nmag + j                -> mode 0
      // color[ncol,ncalc]: i*ncol + j                 -> mode 0
      // lupcorr[ncol,ncalc]: i*ncol + j               -> mode 0, 1

      if ((covmat_temp = (double *)calloc(ncol*ncol*ncalc, sizeof(double))) == NULL) {
	  return -1;
      }
      memcpy(covmat_temp,covmat,sizeof(double)*ncol*ncol*ncalc);

      for (i=0;i<ncalc;i++) {
	  // copy from mode 0
	  gsl_matrix_set_identity(cobs);
	  for (j=0;j<nmag;j++) {
	      // okay
	      gsl_matrix_set(cobs, j, j, magerr[i*nmag+j]*magerr[i*nmag+j]);
	  }

	  gsl_matrix_set_zero(cobstemp);

	  gsl_blas_dgemm(CblasNoTrans, CblasTrans,
			 1.0, mrotmat, cobs,
			 0.0, cobstemp);
	  gsl_blas_dgemm(CblasNoTrans, CblasTrans,
			 1.0, mrotmat, cobstemp,
			 0.0, cobsmat);

	  mvcovmat = gsl_matrix_view_array(&covmat_temp[covmat_stride*i], ncol, ncol);

	  // and the ci matrix
	  // copy from mode 0
	  if (use_refmagerr) {
	      // make the C_i matrix for the refmag err
	      //  this is a matrix with ncol x ncol size
	      gsl_matrix_set_zero(cimat);
	      for (j=0;j<ncol;j++) {
		  for (k=j;k<ncol;k++) {
		      val = slope[j] * slope[k] * refmagerr[i] * refmagerr[i];
		      gsl_matrix_set(cimat, j, k, val);
		      if (k != j) {
			  gsl_matrix_set(cimat, k, j, val);
		      }
		  }
	      }
	  }

	  // check sigint
	  // copy from mode 0
	  test = 1;
	  for (j=0;j<ncol;j++) {
	      if (gsl_matrix_get(&mvcovmat.matrix, j, j) < sigint*sigint) {
		  test = 0;
	      }
	  }

	  if (test == 0) {
	      if (do_chisq) {
		  dist[i] = 1e11;
	      } else {
		  dist[i] = -1e11;
	      }
	  } else {
	      // copy from mode 0
	      if (!nophotoerr) {
		  gsl_matrix_add(&mvcovmat.matrix, cobsmat);
	      }

	      // copy from mode 0
	      if (use_refmagerr) {
		  gsl_matrix_add(&mvcovmat.matrix, cimat);
	      }

	      // check and fix the matrix if necessary
	      // copy from mode 0
	      check_and_fix_covmat(&mvcovmat.matrix);

	      gsl_linalg_LU_decomp(&mvcovmat.matrix, pp, &s);
	      gsl_linalg_LU_invert(&mvcovmat.matrix, pp, mmetric);

	      if (!do_chisq) {
		  // need the determinant
		  norm = gsl_linalg_LU_det(&mvcovmat.matrix, s);
	      }

	      // vdc is a vector with ncol...
	      // copy from mode 1
	      for (j=0;j<ncol;j++) {
		  gsl_vector_set(vdc,j,
				 (c[i*ncol+j] + slope[i*ncol+j]*(refmag[i] - pivotmag[i])) + lupcorr[i*ncol+j] -
				 color[i*ncol+j]);
	      }

	      gsl_blas_dgemv(CblasNoTrans, 1.0, mmetric, vdc, 0.0, vdcm);
	      gsl_blas_ddot(vdcm,vdc, &chisq);

	      if (do_chisq) {
		  dist[i] = chisq;
	      } else {
		  dist[i]=-0.5*chisq-0.5*log(norm);
	      }
	  }

      }
      free(covmat_temp);
  }

  gsl_matrix_free(mrotmat);
  gsl_matrix_free(cobs);
  gsl_matrix_free(cobsmat);
  gsl_matrix_free(cobstemp);
  if (use_refmagerr) gsl_matrix_free(cimat);
  gsl_permutation_free(pp);
  gsl_matrix_free(mmetric);
  gsl_vector_free(vdcm);
  gsl_vector_free(vdc);

  return 0;
}

int check_and_fix_covmat(gsl_matrix *covmat) {
  int i, s, test;
  double eigenval_i;
  gsl_vector_view diag;

  int nelt;
  gsl_matrix *mat;
  gsl_vector *eigenval;
  gsl_vector *eigenval_temp;
  gsl_eigen_symm_workspace *wval;
  gsl_eigen_symmv_workspace *wvec;
  gsl_matrix *Q;
  gsl_matrix *Qinv;
  gsl_matrix *Lambda;
  gsl_matrix *temp;
  gsl_permutation *pp;

  nelt = covmat->size1;

  mat=gsl_matrix_alloc(nelt,nelt);
  wval = gsl_eigen_symm_alloc(nelt);
  wvec = gsl_eigen_symmv_alloc(nelt);
  eigenval = gsl_vector_alloc(nelt);

  // don't destroy the input matrix!
  gsl_matrix_memcpy(mat, covmat);

  // calculate eigenvalues...
  gsl_eigen_symm(mat, eigenval, wval);

  // test if the eigenvalues are negative...
  test = 0;
  for (i=0;i<nelt;i++) {
    eigenval_i = gsl_vector_get(eigenval,i);
    if (eigenval_i < MIN_EIGENVAL) {
      test = 1;
      gsl_vector_set(eigenval,i,MIN_EIGENVAL);
    }
  }
  if (test) {
    // initialize
    Q=gsl_matrix_alloc(nelt,nelt);
    Qinv=gsl_matrix_alloc(nelt,nelt);
    Lambda=gsl_matrix_alloc(nelt,nelt);
    temp=gsl_matrix_alloc(nelt,nelt);
    eigenval_temp=gsl_vector_alloc(nelt);
    pp=gsl_permutation_alloc(nelt);

    // reset the matrix
    gsl_matrix_memcpy(mat, covmat);

    // calculate eigenvalues and eigenvector matrix Q
    gsl_eigen_symmv(mat, eigenval_temp, Q, wvec);

    // invert eigenvector matrix Q-> Qinv (leaving Q in place)
    gsl_matrix_memcpy(temp,Q);
    gsl_linalg_LU_decomp(temp, pp, &s);
    gsl_linalg_LU_invert(temp, pp, Qinv);

    // create a diagonal matrix Lambda
    diag = gsl_matrix_diagonal(Lambda);
    gsl_matrix_set_zero(Lambda);
    gsl_vector_memcpy(&diag.vector, eigenval);

    // do the multiplication
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Q, Lambda, 0.0, temp);
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, temp, Qinv, 0.0, covmat);

    gsl_matrix_free(Q);
    gsl_matrix_free(Qinv);
    gsl_matrix_free(Lambda);
    gsl_matrix_free(temp);
    gsl_vector_free(eigenval_temp);
    gsl_permutation_free(pp);

  }

  gsl_matrix_free(mat);
  gsl_vector_free(eigenval);
  gsl_eigen_symm_free(wval);
  gsl_eigen_symmv_free(wvec);

  return 0;

}

max_stars_repo_name: jacobic/redmapper

max_issues_repo_name: jacobic/redmapper

max_forks_repo_name: jacobic/redmapper


text:
#ifndef utils_h
#define utils_h

#include <ceed.h>
#include <petsc.h>

// Translate PetscMemType to CeedMemType
static inline CeedMemType MemTypeP2C(PetscMemType mem_type) {
  return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST;
}

#endif // utils_h