text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#include <petsc.h>
#include <petscmat.h>
#include <petscvec.h>
#include <petscksp.h>
#include <petscpc.h>
#include <StGermain/StGermain.h>
#include <StgDomain/StgDomain.h>
#include <StgFEM/StgFEM.h>
#include <PICellerator/PICellerator.h>
#include <Underworld/Underworld.h>
#include "Solvers/KSPSolvers/KSPSolvers.h"

#include "common-driver-utils.h"
#include "BSSCR.h" /* includes StokesBlockKSPInterface.h */

#define KSPBSSCR         "bsscr"

void bsscr_summary(KSP_BSSCR * bsscrp_self, KSP ksp_S, KSP ksp_inner,
		   Mat K,Mat K2,Mat D,Mat G,Mat C,Vec u,Vec p,Vec f,Vec h,Vec t,
		   double penaltyNumber,PetscTruth KisJustK,
			 double mgSetupTime,
			 double RHSSolveTime,
			 double scrSolveTime,
			 double a11SingleSolveTime){
    PetscTruth flg, found;
    PetscInt   uSize, pSize, lmax, lmin, iterations;
    PetscReal  rNorm, fNorm, uNorm, uNormInf, pNorm, pNormInf, p_sum, min, max;
    Vec q, qq, t2, t3;
    double solutionAnalysisTime;

      PetscPrintf( PETSC_COMM_WORLD,  "\n\nSCR Solver Summary:\n\n");
      if(bsscrp_self->mg)
	    PetscPrintf( PETSC_COMM_WORLD,     "  Multigrid setup:        = %.4g secs \n", mgSetupTime);

			iterations = bsscrp_self->solver->stats.velocity_presolve_its;
      PetscPrintf( PETSC_COMM_WORLD,     "  RHS V Solve:            = %.4g secs / %d its\n", RHSSolveTime, iterations);

      KSPGetIterationNumber( ksp_S, &iterations);
      bsscrp_self->solver->stats.pressure_its = iterations;

      PetscPrintf( PETSC_COMM_WORLD,     "  Pressure Solve:         = %.4g secs / %d its\n", scrSolveTime, iterations);
      KSPGetIterationNumber( ksp_inner, &iterations);

      bsscrp_self->solver->stats.velocity_backsolve_its = iterations;
      PetscPrintf( PETSC_COMM_WORLD,     "  Final V Solve:          = %.4g secs / %d its\n\n", a11SingleSolveTime, iterations);

      /***************************************************************************************************************/

      flg = PETSC_FALSE; /* Off by default */
	    PetscOptionsGetTruth( PETSC_NULL, "-scr_ksp_solution_summary", &flg, &found );

      if(flg) {
	    PetscScalar KuNorm;
	    solutionAnalysisTime = MPI_Wtime();
	    VecGetSize( u, &uSize );
	    VecGetSize( p, &pSize );
	    VecDuplicate( u, &t2 );
	    VecDuplicate( u, &t3 );
	    MatMult( K, u, t3); VecNorm( t3, NORM_2, &KuNorm );
	    double angle, kdot;
	    if(penaltyNumber > 1e-10){/* should change this to ifK2built maybe */
          MatMult( K2, u, t2); VecNorm( t2, NORM_2, &rNorm );
          VecDot(t2,t3,&kdot);
          angle = (kdot/(rNorm*KuNorm));
          PetscPrintf( PETSC_COMM_WORLD,  "  <K u, K2 u>/(|K u| |K2 u|)    = %.6e\n", angle);
	    }
	    VecNorm( t, NORM_2, &rNorm ); /* t = f- G p  should be the formal residual vector, calculated on line 267 in auglag-driver-DGTGD.c */
	    VecDot(t3,t,&kdot);
	    angle = (kdot/(rNorm*KuNorm));
            PetscPrintf( PETSC_COMM_WORLD,  "  <K u, (f-G p)>/(|K u| |f- G p|)    = %.6e\n\n", angle);

	    MatMult( K, u, t2); VecNorm(t2, NORM_2, &KuNorm);
	    VecAYPX( t2, -1.0, t ); /* t2 <- -t2 + t  : t = f- G p  should be the formal residual vector, calculated on line 267 in auglag-driver-DGTGD.c*/
	    VecNorm( t2, NORM_2, &rNorm );
	    VecNorm( f,  NORM_2, &fNorm );
	    if(KisJustK){
          PetscPrintf( PETSC_COMM_WORLD,"Velocity back-solve with original K matrix\n");
          PetscPrintf( PETSC_COMM_WORLD,"Solved    K u = G p -f\n");
          PetscPrintf( PETSC_COMM_WORLD,"Residual with original K matrix\n");
          PetscPrintf( PETSC_COMM_WORLD,  "  |f - K u - G p|                       = %.12e\n", rNorm);
          PetscPrintf( PETSC_COMM_WORLD,  "  |f - K u - G p|/|f|                   = %.12e\n", rNorm/fNorm);
          if(penaltyNumber > 1e-10){/* means the restore_K flag was used */
            //if(K2 && f2){
            MatAXPY(K,penaltyNumber,K2,DIFFERENT_NONZERO_PATTERN);/* Computes K = penaltyNumber*K2 + K */
            //VecAXPY(f,penaltyNumber,f2); /* f = penaltyNumber*f2 + f */
            KisJustK=PETSC_FALSE;
            MatMult( K, u, t2);
            MatMult( G, p, t);
            VecAYPX( t, -1.0, f ); /* t <- -t + f */
            VecAYPX( t2, -1.0, t ); /* t2 <- -t2 + t */
            VecNorm( t2, NORM_2, &rNorm );
            PetscPrintf( PETSC_COMM_WORLD,"Residual with K+K2 matrix and f rhs vector\n");
            PetscPrintf( PETSC_COMM_WORLD,  "  |(f) - (K + K2) u - G p|         = %.12e\n", rNorm);
            //}
          }
	    }
	    else{
          PetscPrintf( PETSC_COMM_WORLD,"Velocity back-solve with K+K2 matrix\n");
          PetscPrintf( PETSC_COMM_WORLD,"Solved    (K + K2) u = G p - (f)\n");
          PetscPrintf( PETSC_COMM_WORLD,"Residual with K+K2 matrix and f rhs vector\n");
          PetscPrintf( PETSC_COMM_WORLD,  "  |(f) - (K + K2) u - G p|         = %.12e\n", rNorm);
          PetscReal KK2Norm,KK2Normf;
          MatNorm(K,NORM_1,&KK2Norm);
          MatNorm(K,NORM_FROBENIUS,&KK2Normf);
          penaltyNumber = -penaltyNumber;
          MatAXPY(K,penaltyNumber,K2,DIFFERENT_NONZERO_PATTERN);/* Computes K = penaltyNumber*K2 + K */
          //VecAXPY(f,penaltyNumber,f2); /* f = penaltyNumber*f2 + f */
          KisJustK=PETSC_FALSE;
          MatMult( K, u, t2);    /* t2 = K*u  */
          MatMult( G, p, t);     /* t  = G*p  */
          VecAYPX( t, -1.0, f ); /* t <- f - t ; t = f - G*p  */
          VecAYPX( t2, -1.0, t ); /* t2 <- t - t2; t2 = f - G*p - K*u  */
          VecNorm( t2, NORM_2, &rNorm );
          PetscPrintf( PETSC_COMM_WORLD,"Residual with original K matrix\n");
          PetscPrintf( PETSC_COMM_WORLD,  "  |f - K u - G p|                       = %.12e\n", rNorm);
          PetscPrintf( PETSC_COMM_WORLD,  "  |f - K u - G p|/|f|                   = %.12e\n", rNorm/fNorm);
          PetscReal KNorm, K2Norm;
          MatNorm(K,NORM_1,&KNorm);	  MatNorm(K2,NORM_1,&K2Norm);
          PetscPrintf( PETSC_COMM_WORLD,"K and K2 norm_1    %.12e %.12e   ratio %.12e\n",KNorm,K2Norm,K2Norm/KNorm);
          MatNorm(K,NORM_INFINITY,&KNorm);  MatNorm(K2,NORM_INFINITY,&K2Norm);
          PetscPrintf( PETSC_COMM_WORLD,"K and K2 norm_inf  %.12e %.12e   ratio %.12e\n",KNorm,K2Norm,K2Norm/KNorm);
          MatNorm(K,NORM_FROBENIUS,&KNorm); MatNorm(K2,NORM_FROBENIUS,&K2Norm);
          PetscPrintf( PETSC_COMM_WORLD,"K and K2 norm_frob %.12e %.12e   ratio %.12e\n",KNorm,K2Norm,K2Norm/KNorm);
          PetscPrintf( PETSC_COMM_WORLD,"K+r*K2 norm_1    %.12e\n",KK2Norm);
          PetscPrintf( PETSC_COMM_WORLD,"K+r*K2 norm_frob %.12e\n",KK2Normf);
          penaltyNumber = -penaltyNumber;
          MatAXPY(K,penaltyNumber,K2,DIFFERENT_NONZERO_PATTERN);/* Computes K = penaltyNumber*K2 + K */
        }
        PetscPrintf( PETSC_COMM_WORLD,"\n");
        PetscPrintf( PETSC_COMM_WORLD,  "  |K u|    = %.12e\n", KuNorm);
        if(penaltyNumber > 1e-10){
          MatMult( K2, u, t2); VecNorm( t2, NORM_2, &rNorm );
          PetscPrintf( PETSC_COMM_WORLD,  "  |K2 u|   = %.12e\n", rNorm);
          PetscPrintf( PETSC_COMM_WORLD,"\n");
	    }

	    VecDuplicate( p, &q );
	    MatMult( D, u, q );   /* q = G'*u = D*u */
	    VecNorm( u, NORM_2, &uNorm );
	    VecNorm( q, NORM_2, &rNorm );

	    PetscPrintf( PETSC_COMM_WORLD,  "  |G^T u|_2               = %.6e\n", rNorm );
	    PetscPrintf( PETSC_COMM_WORLD,  "  |G^T u|_2/|u|_2         = %.6e\n", sqrt( (double) uSize / (double) pSize ) * rNorm / uNorm);

	    VecDuplicate( p, &qq );
	    MatMultTranspose( G, u, qq );
	    VecNorm( qq, NORM_2, &rNorm );
	    PetscPrintf( PETSC_COMM_WORLD,  "  |G^T u|/|u|             = %.8e\n", rNorm/uNorm ); /* to compare directly with Uzawa */

	    VecNorm( q, NORM_INFINITY, &rNorm );
	    PetscPrintf( PETSC_COMM_WORLD,  "  |G^T u|_infty/|u|_2     = %.6e\n", sqrt( (double) uSize ) * rNorm / uNorm);
	    /* create G'*u+C*p-h to check on this constraint */
	    /* already have q = D*u */
	    VecZeroEntries(qq);
	    if(C){
          MatMult( C, p, qq );
	    }
	    VecAYPX( q, 1.0, qq ); /* q = q+qq; G'*u + C*p*/
	    VecAXPY( q, -1.0, h ); /* q = q-h;  G'*u + C*p - h  */
	    VecNorm( q, NORM_2, &rNorm );
	    PetscPrintf( PETSC_COMM_WORLD,  "  |G^T u + C p - h|        = %.8e  :constraint\n", rNorm );

	    VecNorm( u, NORM_INFINITY, &uNormInf );
	    VecNorm( u, NORM_2,        &uNorm );
	    VecGetSize( u, &uSize );

	    VecNorm( p, NORM_INFINITY, &pNormInf );
	    VecNorm( p, NORM_2,        &pNorm );

	    PetscPrintf( PETSC_COMM_WORLD,  "  |u|_{\\infty}  = %.6e , u_rms = %.6e\n",
	                 uNormInf, uNorm / sqrt( (double) uSize ) );

	    PetscPrintf( PETSC_COMM_WORLD,  "  |p|_{\\infty}  = %.6e , p_rms = %.6e\n",
	                 pNormInf, pNorm / sqrt( (double) pSize ) );

	    VecMax( u, &lmax, &max );
	    VecMin( u, &lmin, &min );
	    PetscPrintf( PETSC_COMM_WORLD,  "  min/max(u)    = %.6e [%d] / %.6e [%d]\n",min,lmin,max,lmax);
        bsscrp_self->solver->stats.vmin = min;
        bsscrp_self->solver->stats.vmax = max;

	    VecMax( p, &lmax, &max );
	    VecMin( p, &lmin, &min );
	    PetscPrintf( PETSC_COMM_WORLD,  "  min/max(p)    = %.6e [%d] / %.6e [%d]\n",min,lmin,max,lmax);
        bsscrp_self->solver->stats.pmin = min;
        bsscrp_self->solver->stats.pmax = max;

	    VecSum( p, &p_sum );
	    PetscPrintf( PETSC_COMM_WORLD,  "  \\sum_i p_i    = %.6e \n", p_sum );
        bsscrp_self->solver->stats.p_sum=p_sum;

	    solutionAnalysisTime = MPI_Wtime() - solutionAnalysisTime;

	    PetscPrintf( PETSC_COMM_WORLD,  "\n  Time for this analysis  = %.4g secs\n\n",solutionAnalysisTime);

	    Stg_VecDestroy(&t2 );
	    Stg_VecDestroy(&t3 );
	    Stg_VecDestroy(&q );
	    Stg_VecDestroy(&qq );
      }

}


text:
#include <gsl/gsl_spline.h>
#include "eos_pres.h"
#include "global.h"

// "forward" EOS. input: rho, T. output: P
double
eos_pres (double rho, void *params)
{
  struct param *eos_pres_params = (struct param *) params;
  return (gsl_spline_eval (spline, rho, acc));
  // return (pow(rho, 1.0/params.Gamma)); // polytrope
}

ext: c


text:
/**
 *
 * @file dsungesv.c
 *
 *  PLASMA computational routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Emmanuel Agullo
 * @date 2010-11-15
 * @generated ds Tue Jan  7 11:45:09 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <lapacke.h>
#include "common.h"

#define PLASMA_dlag2s(_descA, _descSB)                \
  plasma_parallel_call_4(plasma_pdlag2s,              \
                         PLASMA_desc,      (_descA),  \
                         PLASMA_desc,      (_descSB), \
                         PLASMA_sequence*, sequence,  \
                         PLASMA_request*,  request)

#define PLASMA_slag2d(_descSA, _descB)                \
  plasma_parallel_call_4(plasma_pslag2d,              \
                         PLASMA_desc,      (_descSA), \
                         PLASMA_desc,      (_descB),  \
                         PLASMA_sequence*, sequence,  \
                         PLASMA_request*,  request)

#define PLASMA_dlange(_norm, _descA, _result, _work)   \
  _result = 0;                                         \
  plasma_parallel_call_6(plasma_pdlange,               \
                         PLASMA_enum,      (_norm),    \
                         PLASMA_desc,      (_descA),   \
                         double*,          (_work),    \
                         double*,          &(_result), \
                         PLASMA_sequence*, sequence,   \
                         PLASMA_request*,  request);

#define PLASMA_dlacpy(_descA, _descB)                        \
  plasma_parallel_call_5(plasma_pdlacpy,                     \
                         PLASMA_enum,      PlasmaUpperLower, \
                         PLASMA_desc,      (_descA),         \
                         PLASMA_desc,      (_descB),         \
                         PLASMA_sequence*, sequence,         \
                         PLASMA_request*,  request)

#define PLASMA_dgeadd(_alpha, _descA, _descB)                 \
  plasma_parallel_call_5(plasma_pdgeadd,                      \
                         double, (_alpha),        \
                         PLASMA_desc,        (_descA),        \
                         PLASMA_desc,        (_descB),        \
                         PLASMA_sequence*,   sequence,        \
                         PLASMA_request*,    request)

/***************************************************************************//**
 *
 * @ingroup double
 *
 *  PLASMA_dsungesv - Solves overdetermined or underdetermined linear systems involving an M-by-N
 *  matrix A using the QR or the LQ factorization of A.  It is assumed that A has full rank.
 *  The following options are provided:
 *
 *  # trans = PlasmaNoTrans and M >= N: find the least squares solution of an overdetermined
 *    system, i.e., solve the least squares problem: minimize || B - A*X ||.
 *
 *  # trans = PlasmaNoTrans and M < N:  find the minimum norm solution of an underdetermined
 *    system A * X = B.
 *
 *  Several right hand side vectors B and solution vectors X can be handled in a single call;
 *  they are stored as the columns of the M-by-NRHS right hand side matrix B and the N-by-NRHS
 *  solution matrix X.
 *
 *  PLASMA_dsungesv first attempts to factorize the matrix in COMPLEX and use this
 *  factorization within an iterative refinement procedure to produce a
 *  solution with COMPLEX*16 normwise backward error quality (see below).
 *  If the approach fails the method switches to a COMPLEX*16
 *  factorization and solve.
 *
 *  The iterative refinement is not going to be a winning strategy if
 *  the ratio COMPLEX performance over COMPLEX*16 performance is too
 *  small. A reasonable strategy should take the number of right-hand
 *  sides and the size of the matrix into account. This might be done
 *  with a call to ILAENV in the future. Up to now, we always try
 *  iterative refinement.
 *
 *  The iterative refinement process is stopped if ITER > ITERMAX or
 *  for all the RHS we have: RNRM < N*XNRM*ANRM*EPS*BWDMAX
 *  where:
 *
 *  - ITER is the number of the current iteration in the iterative refinement process
 *  - RNRM is the infinity-norm of the residual
 *  - XNRM is the infinity-norm of the solution
 *  - ANRM is the infinity-operator-norm of the matrix A
 *  - EPS is the machine epsilon returned by DLAMCH('Epsilon').
 *
 *  Actually, in its current state (PLASMA 2.1.0), the test is slightly relaxed.
 *
 *  The values ITERMAX and BWDMAX are fixed to 30 and 1.0D+00 respectively.
 *
 *  We follow Bjorck's algorithm proposed in "Iterative Refinement of Linear
 *  Least Squares solutions I", BIT, 7:257-278, 1967.4
 *
 *******************************************************************************
 *
 * @param[in] trans
 *          Intended usage:
 *          = PlasmaNoTrans:   the linear system involves A;
 *          = PlasmaTrans: the linear system involves A**H.
 *          Currently only PlasmaNoTrans is supported.
 *
 * @param[in] N
 *          The number of columns of the matrix A. N >= 0.
 *
 * @param[in] NRHS
 *          The number of right hand sides, i.e., the number of columns of the matrices B and X.
 *          NRHS >= 0.
 *
 * @param[in] A
  *          The M-by-N matrix A. This matrix is not modified.
 *
 * @param[in] LDA
 *          The leading dimension of the array A. LDA >= max(1,M).
 *
 * @param[in] B
 *          The M-by-NRHS matrix B of right hand side vectors, stored columnwise. Not modified.
 *
 * @param[in] LDB
 *          The leading dimension of the array B. LDB >= MAX(1,M,N).
 *
 * @param[out] X
 *          If return value = 0, the solution vectors, stored columnwise.
 *          if M >= N, rows 1 to N of B contain the least squares solution vectors; the residual
 *          sum of squares for the solution in each column is given by the sum of squares of the
 *          modulus of elements N+1 to M in that column;
 *          if M < N, rows 1 to N of B contain the minimum norm solution vectors;
 *
 * @param[in] LDX
 *          The leading dimension of the array B. LDB >= MAX(1,M,N).
 *
 * @param[out] ITER
 *          The number of the current iteration in the iterative refinement process
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 *******************************************************************************
 *
 * @sa PLASMA_dsungesv_Tile
 * @sa PLASMA_dsungesv_Tile_Async
 * @sa PLASMA_dsungesv
 * @sa PLASMA_dgels
 *
 ******************************************************************************/
int PLASMA_dsungesv(PLASMA_enum trans, int N, int NRHS,
                    double *A, int LDA,
                    double *B, int LDB,
                    double *X, int LDX, int *ITER)
{
    int NB;
    int status;
    PLASMA_desc  descA;
    PLASMA_desc  descB;
    PLASMA_desc *descT;
    PLASMA_desc  descX;
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_dsungesv", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }

    *ITER = 0;

    /* Check input arguments */
    if (trans != PlasmaNoTrans   &&
        trans != PlasmaTrans &&
        trans != PlasmaTrans )
    {
        plasma_error("PLASMA_dsungesv", "illegal value of trans");
        return -1;
    }
    if (trans != PlasmaNoTrans) {
        plasma_error("PLASMA_dsungesv", "only PlasmaNoTrans supported");
        return PLASMA_ERR_NOT_SUPPORTED;
    }
    if (N < 0) {
        plasma_error("PLASMA_dsungesv", "illegal value of N");
        return -2;
    }
    if (NRHS < 0) {
        plasma_error("PLASMA_dsungesv", "illegal value of NRHS");
        return -3;
    }
    if (LDA < max(1, N)) {
        plasma_error("PLASMA_dsungesv", "illegal value of LDA");
        return -5;
    }
    if (LDB < max(1, N)) {
        plasma_error("PLASMA_dsungesv", "illegal value of LDB");
        return -7;
    }
    if (LDX < max(1, N)) {
        plasma_error("PLASMA_dsungesv", "illegal value of LDX");
        return -9;
    }

    /* Quick return */
    if ( N == 0 )
        return PLASMA_SUCCESS;

    /* Tune NB & IB depending on M, N & NRHS; Set NBNB */
    status = plasma_tune(PLASMA_FUNC_DSGELS, N, N, NRHS);
    if (status != PLASMA_SUCCESS) {
        plasma_error("PLASMA_dsungesv", "plasma_tune() failed");
        return status;
    }

    NB = PLASMA_NB;

    plasma_sequence_create(plasma, &sequence);

    /* DOUBLE PRECISION INITIALIZATION */
    if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
        plasma_dooplap2tile( descA, A, NB, NB, LDA, N,    0, 0, N, N,    sequence, &request,
                             plasma_desc_mat_free(&(descA)) );
        plasma_dooplap2tile( descB, B, NB, NB, LDB, NRHS, 0, 0, N, NRHS, sequence, &request,
                             plasma_desc_mat_free(&(descA)); plasma_desc_mat_free(&(descB)) );
        plasma_ddesc_alloc(  descX, NB, NB, N, NRHS, 0, 0, N, NRHS, plasma_desc_mat_free(&(descA)); plasma_desc_mat_free(&(descB)); plasma_desc_mat_free(&(descX)) );
    } else {
        plasma_diplap2tile( descA, A, NB, NB, LDA, N,    0, 0, N, N,
                            sequence, &request);
        plasma_diplap2tile( descB, B, NB, NB, LDB, NRHS, 0, 0, N, NRHS,
                            sequence, &request);

        descX = plasma_desc_init(
            PlasmaRealDouble, NB, NB, (NB*NB),
            LDX, NRHS, 0, 0, N, NRHS);
        descX.mat = X;
    }

    /* Allocate workspace */
    PLASMA_Alloc_Workspace_dgels_Tile(N, N, &descT);

    /* Call the native interface */
    status = PLASMA_dsungesv_Tile_Async(PlasmaNoTrans, &descA, descT, &descB, &descX, ITER,
                                        sequence, &request);

    if (status == PLASMA_SUCCESS) {
        if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
            plasma_dooptile2lap( descX, X, NB, NB, LDX, NRHS,  sequence, &request);
            plasma_dynamic_sync();
            plasma_desc_mat_free(&descA);
            plasma_desc_mat_free(&descB);
            plasma_desc_mat_free(&descX);
        } else {
            plasma_diptile2lap( descA, A, NB, NB, LDA, N,     sequence, &request);
            plasma_diptile2lap( descB, B, NB, NB, LDB, NRHS,  sequence, &request);
            plasma_diptile2lap( descX, X, NB, NB, LDX, NRHS,  sequence, &request);
            plasma_dynamic_sync();
        }
    }

    PLASMA_Dealloc_Handle_Tile(&descT);
    plasma_sequence_destroy(plasma, sequence);
    return status;
}

/***************************************************************************//**
 *
 * @ingroup double_Tile
 *
 *  PLASMA_dsungesv_Tile - Solves symmetric linear system of equations using the tile QR
 *  or the tile LQ factorization and mixed-precision iterative refinement.
 *  Tile equivalent of PLASMA_dsungesv().
 *  Operates on matrices stored by tiles.
 *  All matrices are passed through descriptors.
 *  All dimensions are taken from the descriptors.
 *
 *******************************************************************************
 *
 * @param[in] trans
 *          Intended usage:
 *          = PlasmaNoTrans:   the linear system involves A;
 *          = PlasmaTrans: the linear system involves A**H.
 *          Currently only PlasmaNoTrans is supported.
 *
 * @param[in,out] A
 *          - If the iterative refinement converged, A is not modified;
 *          - otherwise, it fell back to double precision solution, and
 *          on exit the M-by-N matrix A contains:
 *          if M >= N, A is overwritten by details of its QR factorization as returned by
 *                     PLASMA_dgeqrf;
 *          if M < N, A is overwritten by details of its LQ factorization as returned by
 *                      PLASMA_dgelqf.
 *
 * @param[out] T
 *          On exit:
 *          - if the iterative refinement converged, T is not modified;
 *          - otherwise, it fell back to double precision solution,
 *          and then T is an auxiliary factorization data.
 *
 * @param[in,out] B
 *          On entry, the M-by-NRHS matrix B of right hand side vectors, stored columnwise;
 * @param[in] B
 *          The N-by-NRHS matrix of right hand side matrix B.
 *
 * @param[out] X
 *          If return value = 0, X is the solution vectors, stored columnwise:
 *          if M >= N, rows 1 to N of X contain the least squares solution vectors; the residual
 *          sum of squares for the solution in each column is given by the sum of squares of the
 *          modulus of elements N+1 to M in that column;
 *          if M < N, rows 1 to N of X contain the minimum norm solution vectors;
 *
 * @param[out] ITER
 *          The number of the current iteration in the iterative refinement process
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *
 *******************************************************************************
 *
 * @sa PLASMA_dsungesv
 * @sa PLASMA_dsungesv_Tile_Async
 * @sa PLASMA_dsungesv_Tile
 * @sa PLASMA_dgels_Tile
 *
 ******************************************************************************/
int PLASMA_dsungesv_Tile(PLASMA_enum trans, PLASMA_desc *A, PLASMA_desc *T,
                         PLASMA_desc *B, PLASMA_desc *X, int *ITER)
{
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;
    int status;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_dsungesv_Tile", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    plasma_sequence_create(plasma, &sequence);
    status = PLASMA_dsungesv_Tile_Async(trans, A, T, B, X, ITER, sequence, &request);
    if (status != PLASMA_SUCCESS)
        return status;
    plasma_dynamic_sync();
    status = sequence->status;
    plasma_sequence_destroy(plasma, sequence);
    return status;
}

/***************************************************************************//**
 *
 * @ingroup double_Tile_Async
 *
 *  PLASMA_dsungesv_Tile_Async - Solves symmetric linear system of equations using
 *  the tile QR or the tile LQ factorization and mixed-precision iterative refinement.
 *  Non-blocking equivalent of PLASMA_dsungesv_Tile().
 *  May return before the computation is finished.
 *  Allows for pipelining of operations at runtime.
 *
 *******************************************************************************
 *
 * @param[in] sequence
 *          Identifies the sequence of function calls that this call belongs to
 *          (for completion checks and exception handling purposes).
 *
 * @param[out] request
 *          Identifies this function call (for exception handling purposes).
 *
 *******************************************************************************
 *
 * @sa PLASMA_dsungesv
 * @sa PLASMA_dsungesv_Tile
 * @sa PLASMA_dsungesv_Tile_Async
 * @sa PLASMA_dgels_Tile_Async
 *
 ******************************************************************************/
int PLASMA_dsungesv_Tile_Async(PLASMA_enum trans, PLASMA_desc *A, PLASMA_desc *T,
                               PLASMA_desc *B, PLASMA_desc *X, int *ITER,
                               PLASMA_sequence *sequence, PLASMA_request *request)
{
    int N, NB, IB;
    PLASMA_desc descA;
    PLASMA_desc descT;
    PLASMA_desc descB;
    PLASMA_desc descX;
    PLASMA_desc descR, descSA, descST, descSX;
    plasma_context_t *plasma;
    double *work;

    const int    itermax = 30;
    const double bwdmax  = 1.0;
    const double negone = -1.0;
    const double one = 1.0;
    int iiter;
    double Anorm, cte, eps, Rnorm, Xnorm;
    *ITER=0;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_dsungesv_Tile", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    if (sequence == NULL) {
        plasma_fatal_error("PLASMA_dsungesv_Tile", "NULL sequence");
        return PLASMA_ERR_UNALLOCATED;
    }
    if (request == NULL) {
        plasma_fatal_error("PLASMA_dsungesv_Tile", "NULL request");
        return PLASMA_ERR_UNALLOCATED;
    }
    /* Check sequence status */
    if (sequence->status == PLASMA_SUCCESS)
        request->status = PLASMA_SUCCESS;
    else
        return plasma_request_fail(sequence, request, PLASMA_ERR_SEQUENCE_FLUSHED);
    /* Check descriptors for correctness */
    if (plasma_desc_check(A) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_dsungesv_Tile", "invalid first descriptor");
        return PLASMA_ERR_ILLEGAL_VALUE;
    } else {
        descA = *A;
    }
    if (plasma_desc_check(T) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_dsungesv_Tile", "invalid second descriptor");
        return PLASMA_ERR_ILLEGAL_VALUE;
    } else {
        descT = *T;
    }
    if (plasma_desc_check(B) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_dsungesv_Tile", "invalid third descriptor");
        return PLASMA_ERR_ILLEGAL_VALUE;
    } else {
        descB = *B;
    }
    if (plasma_desc_check(X) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_dsungesv_Tile", "invalid fourth descriptor");
        return PLASMA_ERR_ILLEGAL_VALUE;
    } else {
        descX = *X;
    }
    /* Check input arguments */
    if ( (descA.nb != descA.mb) || (descB.nb != descB.mb) || (descX.nb != descX.mb) ||
         (descA.mb != descB.mb) || (descB.mb != descX.mb) ) {
        plasma_error("PLASMA_dsungesv_Tile", "only square tiles of same size are supported");
        return PLASMA_ERR_ILLEGAL_VALUE;
    }
    if (trans != PlasmaNoTrans) {
        plasma_error("PLASMA_dsungesv_Tile", "only PlasmaNoTrans supported");
        return PLASMA_ERR_NOT_SUPPORTED;
    }

    /* Set N, NRHS, NB */
    N  = descA.m;
    NB = descA.nb;
    IB = descT.mb;

    work = (double *)plasma_shared_alloc(plasma, PLASMA_SIZE, PlasmaRealDouble);
    if (work == NULL) {
        plasma_error("PLASMA_dsungesv", "plasma_shared_alloc() failed");
        plasma_shared_free(plasma, work);
        return PLASMA_ERR_OUT_OF_RESOURCES;
    }

    plasma_ddesc_alloc( descR,  NB, NB, descB.m, descB.n, 0, 0, descB.m, descB.n, plasma_shared_free( plasma, work ); plasma_desc_mat_free(&descR) );
    plasma_sdesc_alloc( descSA, NB, NB, descA.m, descA.n, 0, 0, descA.m, descA.n, plasma_shared_free( plasma, work ); plasma_desc_mat_free(&descR); plasma_desc_mat_free(&descSA) );
    plasma_sdesc_alloc( descST, IB, NB, descT.m, descT.n, 0, 0, descT.m, descT.n, plasma_shared_free( plasma, work ); plasma_desc_mat_free(&descR); plasma_desc_mat_free(&descSA); plasma_desc_mat_free(&descST) );
    plasma_sdesc_alloc( descSX, NB, NB, descX.m, descX.n, 0, 0, descX.m, descX.n, plasma_shared_free( plasma, work ); plasma_desc_mat_free(&descR); plasma_desc_mat_free(&descSA); plasma_desc_mat_free(&descST); plasma_desc_mat_free(&descSX) );

    /* Compute some constants */
    PLASMA_dlange(PlasmaInfNorm, descA, Anorm, work);
    eps = LAPACKE_dlamch_work('e');

    /* Convert B from double precision to single precision and store
       the result in SX. */
    PLASMA_dlag2s(descB, descSX);
    if (sequence->status != PLASMA_SUCCESS)
        return plasma_request_fail(sequence, request, PLASMA_ERR_SEQUENCE_FLUSHED);

    /* Convert A from double precision to single precision and store
       the result in SA. */
    PLASMA_dlag2s(descA, descSA);
    if (sequence->status != PLASMA_SUCCESS)
        return plasma_request_fail(sequence, request, PLASMA_ERR_SEQUENCE_FLUSHED);

    /* Compute the QR factorization of SA */
    plasma_parallel_call_4(plasma_psgeqrf,
        PLASMA_desc, descSA,
        PLASMA_desc, descST,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    /* Compute the solve in simple */
    plasma_parallel_call_7(plasma_psormqr,
        PLASMA_enum, PlasmaLeft,
        PLASMA_enum, PlasmaTrans,
        PLASMA_desc, descSA,
        PLASMA_desc, descSX,
        PLASMA_desc, descST,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    plasma_parallel_call_9(plasma_pstrsm,
        PLASMA_enum, PlasmaLeft,
        PLASMA_enum, PlasmaUpper,
        PLASMA_enum, PlasmaNoTrans,
        PLASMA_enum, PlasmaNonUnit,
        float, 1.0,
        PLASMA_desc, descSA,
        PLASMA_desc, descSX,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    /* Convert SX back to double precision */
    PLASMA_slag2d(descSX, descX);

    /* Compute R = B - AX. */
    PLASMA_dlacpy(descB, descR);

    plasma_parallel_call_9(plasma_pdgemm,
        PLASMA_enum, PlasmaNoTrans,
        PLASMA_enum, PlasmaNoTrans,
        double, negone,
        PLASMA_desc, descA,
        PLASMA_desc, descX,
        double, one,
        PLASMA_desc, descR,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    /* Check whether the NRHS normwise backward error satisfies the
       stopping criterion. If yes return. Note that ITER=0 (already set). */
    PLASMA_dlange(PlasmaInfNorm, descX, Xnorm, work);
    PLASMA_dlange(PlasmaInfNorm, descR, Rnorm, work);

    /* Wait for the end of Anorm, Xnorm and Bnorm computations */
    plasma_dynamic_sync();

    cte = Anorm*eps*((double) N)*bwdmax;
    if (Rnorm < Xnorm * cte){
        /* The NRHS normwise backward errors satisfy the
           stopping criterion. We are good to exit. */
        plasma_desc_mat_free(&descSA);
        plasma_desc_mat_free(&descST);
        plasma_desc_mat_free(&descSX);
        plasma_desc_mat_free(&descR);
        plasma_shared_free(plasma, work);
        return PLASMA_SUCCESS;
    }

    /* Iterative refinement */
    for (iiter = 0; iiter < itermax; iiter++){

        /* Convert R from double precision to single precision
           and store the result in SX. */
        PLASMA_dlag2s(descR, descSX);

        plasma_parallel_call_7(plasma_psormqr,
            PLASMA_enum, PlasmaLeft,
            PLASMA_enum, PlasmaTrans,
            PLASMA_desc, descSA,
            PLASMA_desc, descSX,
            PLASMA_desc, descST,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        plasma_parallel_call_9(plasma_pstrsm,
            PLASMA_enum, PlasmaLeft,
            PLASMA_enum, PlasmaUpper,
            PLASMA_enum, PlasmaNoTrans,
            PLASMA_enum, PlasmaNonUnit,
            float, (float)1.0,
            PLASMA_desc, descSA,
            PLASMA_desc, descSX,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        /* Convert SX back to double precision and update the current
           iterate. */
        PLASMA_slag2d(descSX, descR);
        PLASMA_dgeadd(one, descR, descX);

        /* Compute R = B - AX. */
        PLASMA_dlacpy(descB,descR);
        plasma_parallel_call_9(plasma_pdgemm,
            PLASMA_enum, PlasmaNoTrans,
            PLASMA_enum, PlasmaNoTrans,
            double, negone,
            PLASMA_desc, descA,
            PLASMA_desc, descX,
            double, one,
            PLASMA_desc, descR,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        /* Check whether the NRHS normwise backward errors satisfy the
           stopping criterion. If yes, set ITER=IITER>0 and return. */
        PLASMA_dlange(PlasmaInfNorm, descX, Xnorm, work);
        PLASMA_dlange(PlasmaInfNorm, descR, Rnorm, work);

        /* Wait for the end of Xnorm and Bnorm computations */
        plasma_dynamic_sync();

        if (Rnorm < Xnorm * cte){
            /* The NRHS normwise backward errors satisfy the
               stopping criterion. We are good to exit. */
            *ITER = iiter;

            plasma_desc_mat_free(&descSA);
            plasma_desc_mat_free(&descST);
            plasma_desc_mat_free(&descSX);
            plasma_desc_mat_free(&descR);
            plasma_shared_free(plasma, work);
            return PLASMA_SUCCESS;
        }
    }

    /* We have performed ITER=itermax iterations and never satisified
       the stopping criterion, set up the ITER flag accordingly and
       follow up on double precision routine. */
    *ITER = -itermax - 1;

    plasma_desc_mat_free(&descSA);
    plasma_desc_mat_free(&descST);
    plasma_desc_mat_free(&descSX);
    plasma_desc_mat_free(&descR);
    plasma_shared_free(plasma, work);

    /* Single-precision iterative refinement failed to converge to a
       satisfactory solution, so we restart to double precision. */
    PLASMA_dlacpy(descB, descX);

    plasma_parallel_call_4(plasma_pdgeqrf,
        PLASMA_desc, descA,
        PLASMA_desc, descT,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    plasma_parallel_call_7(plasma_pdormqr,
        PLASMA_enum, PlasmaLeft,
        PLASMA_enum, PlasmaTrans,
        PLASMA_desc, descA,
        PLASMA_desc, descX,
        PLASMA_desc, descT,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    plasma_parallel_call_9(plasma_pdtrsm,
        PLASMA_enum, PlasmaLeft,
        PLASMA_enum, PlasmaUpper,
        PLASMA_enum, PlasmaNoTrans,
        PLASMA_enum, PlasmaNonUnit,
        double, (double)1.0,
        PLASMA_desc, descA,
        PLASMA_desc, descX,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    return PLASMA_SUCCESS;
}


text:
/* bst/gsl_bst_avl.h
 *
 * Copyright (C) 2018 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_BST_AVL_H__
#define __GSL_BST_AVL_H__

#include <gsl/gsl_math.h>

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

__BEGIN_DECLS

#ifndef GSL_BST_AVL_MAX_HEIGHT
#define GSL_BST_AVL_MAX_HEIGHT 32
#endif

/* AVL node */
struct gsl_bst_avl_node
{
  struct gsl_bst_avl_node *avl_link[2]; /* subtrees */
  void *avl_data;                       /* pointer to data */
  signed char avl_balance;              /* balance factor */
};

/* tree data structure */
typedef struct
{
  struct gsl_bst_avl_node *avl_root;          /* tree's root */
  gsl_bst_cmp_function *avl_compare;          /* comparison function */
  void *avl_param;                            /* extra argument to |avl_compare| */
  const gsl_bst_allocator *avl_alloc;         /* memory allocator */
  size_t avl_count;                           /* number of items in tree */
  unsigned long avl_generation;               /* generation number */
} gsl_bst_avl_table;

/* AVL traverser structure */
typedef struct
{
  const gsl_bst_avl_table *avl_table;                         /* tree being traversed */
  struct gsl_bst_avl_node *avl_node;                          /* current node in tree */
  struct gsl_bst_avl_node *avl_stack[GSL_BST_AVL_MAX_HEIGHT]; /* all the nodes above |avl_node| */
  size_t avl_height;                                          /* number of nodes in |avl_parent| */
  unsigned long avl_generation;                               /* generation number */
} gsl_bst_avl_traverser;

__END_DECLS

#endif /* __GSL_BST_AVL_H__ */

ext: h

max_stars_repo_name: vinej/sml

max_issues_repo_name: vinej/sml

max_forks_repo_name: vinej/sml


text:
#ifndef COORDINATE_H_Y7FJINQO
#define COORDINATE_H_Y7FJINQO

#include <gsl/gsl>
#include <sens_loc/math/eigen_types.h>
#include <type_traits>

namespace sens_loc::math {

/// This enum allows to discrimate in which coordinate system or frame of
/// reference a vector lives. It is solely for code readability and correctness
/// evaluation.
enum class frame {
    world,   ///< 3D, Global, World Coordinate System. \f$(U, V, W)\f$.
    camera,  ///< 3D, Camera coordinate system, but the camera center is the
             ///< frame reference. \f$(X, Y, Z)\f$
    sphere,  ///< 3D, projection from image into space. \f$(X_s, Y_s, Z_s)\f$.
             ///< Equivalent to the direction of a lightray for a pixel.
    image,   ///< 2D, coordinate in an image as floating point. \f$(x, y)\f$
    pixel,   ///< 2D, usually a discrete coordinate system that refers to pixels
             ///< in an image. \f$(u, v)\f$
};

namespace detail {
/// Helper type-traits to examine if a reference frame is 2D. If its is not
/// 2D, the coordinate system is 3D.
/// \sa math::frame
template <frame Frame>
struct is_2d
    : std::bool_constant<Frame == frame::image || Frame == frame::pixel> {};

/// Small helper class that holds the data for 2 dimensional coordinates.
/// \sa math::coordinate
template <typename Real>
struct coordinate2d {
    static_assert(std::is_arithmetic_v<Real>);

    coordinate2d()
        : d(Real(0.), Real(0.)) {}
    coordinate2d(Real x1, Real x2)
        : d(x1, x2) {}

    vec<Real, 2> d;
};

/// Small helper class that holds the data for 3 dimensional coordinates.
/// \sa math::coordinate
template <typename Real>
struct coordinate3d {
    static_assert(std::is_arithmetic_v<Real>);

    coordinate3d()
        : d(Real(0.), Real(0.), Real(0.)) {}
    coordinate3d(Real x1, Real x2, Real x3)
        : d(x1, x2, x3) {}

    vec<Real, 3> d;
};
}  // namespace detail

/// This class is a strongly typed coordinate class to ensure the implemented
/// math does not mix up reference frames.
///
/// The overall goal of this slightly complicated construct is to increase
/// readability of the code **USING** this class and to improve correctness.
///
/// This class allows only access with the coordinate conventions for each
/// frame.
/// It prevents mixing different coordinates while assignment or other
/// operations and therefore forces conversion functions. Such conversions
/// are well defined through camera intrinsics for example.
///
/// There are helper typedefs, for easier writing:
/// \c pixel_coord<Real>, \c image_coord<Real>, \c sphere_coord<Real>,
/// \c camera_coord<Real> and \c world_coord<Real>.
///
/// \code
///   const pixel_coord<int> pixel(42, 32);
///   pixel.u(); // Accessing the u-coordinate for that pixel.
///   pixel.x(); // DOES NOT COMPILE - WRONG CONVENTION!
///
///   // Changing coordinate systems requires an explicit transformation, for
///   // example with the camera calibration as in this case.
///   const camera_models::pinhole<float> i;
///   const image_coord<float> transformed = i.transform_to_image(pixel);
///   transformed.x() // Accessing x-coordinate for that pixel.
///                   // Same point but different coordinate system.
///   transformed.u() // DOES NOT COMPILE - WRONG CONVENTION!
///
///   // It is not possible to mix coordinates up.
///   pixel_coord<float> sub_pixel(42.42, 42.42);
///   image_coord<float> transformed = sub_pixel; // DOES NOT COMPILE
///                                               // TYPE MISSMATCH!!
/// \endcode
///
/// \tparam Real underlying type of the coordinates
/// \tparam Frame type of reference frame, results in either 2 or 3 dimensional
/// vectors. The vector data is inherited.
template <typename Real, frame Frame>
class coordinate
    : private std::conditional_t<detail::is_2d<Frame>::value,
                                 detail::coordinate2d<Real>,
                                 detail::coordinate3d<Real>> {
    static_assert(std::is_arithmetic_v<Real>,
                  "Arithmetic type for Real required!");

    using base_class = std::conditional_t<detail::is_2d<Frame>::value,
                                          detail::coordinate2d<Real>,
                                          detail::coordinate3d<Real>>;

  public:
    using real = Real;

    // Reuse constructors of baseclass
    using base_class::base_class;

    /// Calculate the euclidean norm for the coordinate.
    /// \note This function is only enabled if the \p Real type is a floating
    /// point number.
    template <typename T = Real,
              typename   = std::enable_if_t<!std::is_integral_v<T>>>
    [[nodiscard]] Real norm() const noexcept {
        return base_class::d.norm();
    }

    /// Return the normalized vector of this coordinate.
    /// \note This function is only enabled if the \p Real type is a floating
    /// point number and the vector is in 3 dimensional space.
    template <typename T = Real,
              typename U = frame,
              typename   = std::enable_if_t<!std::is_integral_v<T> &&
                                          !detail::is_2d<Frame>::value>>
    [[nodiscard]] coordinate<Real, Frame> normalized() const noexcept {
        const auto res = base_class::d.normalized();
        return {res[0], res[1], res[2]};
    }

    /// Calculate the dot product with \p other.
    [[nodiscard]] Real dot(const coordinate<Real, Frame>& other) const
        noexcept {
        return base_class::d.dot(other.d);
    }
    /// Calculate the cross product with \p other.
    template <
        typename T = frame,
        typename = std::enable_if_t<!(Frame == T::pixel || Frame == T::image)>>
    [[nodiscard]] coordinate<Real, Frame>
    cross(const coordinate<Real, Frame>& other) const noexcept {
        const auto res = base_class::d.cross(other.d);
        return {res[0], res[1], res[2]};
    }

    // clang-format off
    template <typename T = frame,
              typename   = std::enable_if_t<Frame == T::pixel>>
    [[nodiscard]] Real u() const noexcept { return base_class::d[0]; }
    template <typename T = frame,
              typename   = typename std::enable_if<Frame == T::pixel>>
    [[nodiscard]] Real v() const noexcept { return base_class::d[1]; }

    template <typename T = frame,
              typename   = typename std::enable_if_t<Frame == T::image>>
    [[nodiscard]] Real x() const noexcept { return base_class::d[0]; }
    template <typename T = frame,
              typename   = typename std::enable_if<Frame == T::image>>
    [[nodiscard]] Real y() const noexcept { return base_class::d[1]; }

    template <typename T = frame,
              typename   = typename std::enable_if_t<Frame == T::sphere>>
    [[nodiscard]] Real Xs() const noexcept { return base_class::d[0]; }
    template <typename T = frame,
              typename   = typename std::enable_if<Frame == T::sphere>>
    [[nodiscard]] Real Ys() const noexcept { return base_class::d[1]; }
    template <typename T = frame,
              typename   = typename std::enable_if<Frame == T::sphere>>
    [[nodiscard]] Real Zs() const noexcept { return base_class::d[2]; }

    template <typename T = frame,
              typename   = typename std::enable_if_t<Frame == T::camera>>
    [[nodiscard]] Real X() const noexcept { return base_class::d[0]; }
    template <typename T = frame,
              typename   = typename std::enable_if<Frame == T::camera>>
    [[nodiscard]] Real Y() const noexcept { return base_class::d[1]; }
    template <typename T = frame,
              typename   = typename std::enable_if<Frame == T::camera>>
    [[nodiscard]] Real Z() const noexcept { return base_class::d[2]; }

    template <typename T = frame,
              typename   = typename std::enable_if_t<Frame == T::world>>
    [[nodiscard]] Real U() const noexcept { return base_class::d[0]; }
    template <typename T = frame,
              typename   = typename std::enable_if<Frame == T::world>>
    [[nodiscard]] Real V() const noexcept { return base_class::d[1]; }
    template <typename T = frame,
              typename   = typename std::enable_if<Frame == T::world>>
    [[nodiscard]] Real W() const noexcept { return base_class::d[2]; }
    // clang-format on

    coordinate<Real, Frame> operator-(const coordinate<Real, Frame>& p) const
        noexcept {
        coordinate<Real, Frame> res;
        auto                    difference = base_class::d - p.d;

        res.d = difference;
        return res;
    }
};

template <typename Real>
coordinate<Real, frame::camera>
operator*(Real factor, const coordinate<Real, frame::camera>& pt) noexcept {
    static_assert(std::is_floating_point_v<Real>,
                  "Real must be a floating point type for scaling operations!");

    coordinate<Real, frame::camera> result{factor * pt.X(), factor * pt.Y(),
                                           factor * pt.Z()};
    Ensures(std::abs(factor * pt.norm() - result.norm()) < Real(0.0001));
    return result;
}

/// Scaling spherical coordinates results in camera coordinates.
template <typename Real>
coordinate<Real, frame::camera>
operator*(Real factor, const coordinate<Real, frame::sphere>& pt) noexcept {
    static_assert(std::is_floating_point_v<Real>,
                  "Real must be a floating point type for scaling operations!");
    coordinate<Real, frame::camera> result{factor * pt.Xs(), factor * pt.Ys(),
                                           factor * pt.Zs()};
    Ensures(std::abs(factor * pt.norm() - result.norm()) < Real(0.0001));
    return result;
}

template <typename Real = int>
using pixel_coord = coordinate<Real, frame::pixel>;
template <typename Real = float>
using image_coord = coordinate<Real, frame::image>;
template <typename Real = float>
using sphere_coord = coordinate<Real, frame::sphere>;
template <typename Real = float>
using camera_coord = coordinate<Real, frame::camera>;
template <typename Real = float>
using world_coord = coordinate<Real, frame::world>;

}  // namespace sens_loc::math

#endif /* end of include guard: COORDINATE_H_Y7FJINQO */

max_stars_repo_name: JonasToth/depth-conversions

max_issues_repo_name: JonasToth/depth-conversions

max_forks_repo_name: JonasToth/depth-conversions


text:
/* bspline/test.c
 *
 * Copyright (C) 2006, 2007, 2009 Brian Gough
 * Copyright (C) 2008, 2011 Rhys Ulerich
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_nan.h>

void
test_bspline(gsl_bspline_workspace * bw, gsl_bspline_deriv_workspace * dbw)
{
  gsl_vector *B;
  gsl_matrix *dB;
  size_t i, j;
  size_t n = 100;
  size_t ncoeffs = gsl_bspline_ncoeffs(bw);
  size_t order = gsl_bspline_order(bw);
  size_t nbreak = gsl_bspline_nbreak(bw);
  double a = gsl_bspline_breakpoint(0, bw);
  double b = gsl_bspline_breakpoint(nbreak - 1, bw);

  B  = gsl_vector_alloc(ncoeffs);
  dB = gsl_matrix_alloc(ncoeffs, 1);

  /* Ensure B-splines form a partition of unity */
  for (i = 0; i < n; i++)
    {
      double xi = a + (b - a) * (i / (n - 1.0));
      double sum = 0;
      gsl_bspline_eval(xi, B, bw);

      for (j = 0; j < ncoeffs; j++)
        {
          double Bj = gsl_vector_get(B, j);
          int s = (Bj < 0 || Bj > 1);
          gsl_test(s,
                   "basis-spline coefficient %u is in range [0,1] for x=%g",
                   j, xi);
          sum += Bj;
        }

      gsl_test_rel(sum, 1.0, order * GSL_DBL_EPSILON,
                   "basis-spline order %u is normalized for x=%g", order,
                   xi);
    }

  /* Ensure B-splines 0th derivatives agree with regular evaluation */
  for (i = 0; i < n; i++)
    {
      double xi = a + (b - a) * (i / (n - 1.0));
      gsl_bspline_eval(xi, B, bw);
      gsl_bspline_deriv_eval(xi, 0, dB, bw, dbw);

      for (j = 0; j < ncoeffs; j++)
        {
          gsl_test_abs(gsl_matrix_get(dB, j, 0), gsl_vector_get(B, j),
                       GSL_DBL_EPSILON,
                       "b-spline order %d basis #%d evaluation and 0th derivative consistent for x=%g",
                       order, j, xi);
        }

    }

  gsl_vector_free(B);
  gsl_matrix_free(dB);
}

int
main(int argc, char **argv)
{
  size_t order, breakpoints, i;

  gsl_ieee_env_setup();

  argc = 0;                     /* prevent warnings about unused parameters */
  argv = 0;

  for (order = 1; order < 10; order++)
    {
      for (breakpoints = 2; breakpoints < 100; breakpoints++)
        {
          double a = -1.23 * order, b = 45.6 * order;
          gsl_bspline_workspace *bw = gsl_bspline_alloc(order, breakpoints);
          gsl_bspline_deriv_workspace *dbw = gsl_bspline_deriv_alloc(order);
          gsl_bspline_knots_uniform(a, b, bw);
          test_bspline(bw, dbw);
          gsl_bspline_deriv_free(dbw);
          gsl_bspline_free(bw);
        }
    }

  for (order = 1; order < 10; order++)
    {
      for (breakpoints = 2; breakpoints < 100; breakpoints++)
        {
          double a = -1.23 * order, b = 45.6 * order;
          gsl_bspline_workspace *bw = gsl_bspline_alloc(order, breakpoints);
          gsl_bspline_deriv_workspace *dbw = gsl_bspline_deriv_alloc(order);
          gsl_vector *k = gsl_vector_alloc(breakpoints);
          for (i = 0; i < breakpoints; i++)
            {
              double f, x;
              f = sqrt(i / (breakpoints - 1.0));
              x = (1 - f) * a + f * b;
              gsl_vector_set(k, i, x);
            };
          gsl_bspline_knots(k, bw);
          test_bspline(bw, dbw);
          gsl_vector_free(k);
          gsl_bspline_deriv_free(dbw);
          gsl_bspline_free(bw);
        }
    }

  /* Spot check known 0th, 1st, 2nd derivative
     evaluations for a particular k = 2 case.  */
  {
    size_t i, j; /* looping */

    const double xloc[4]     =  { 0.0,  1.0,  6.0,  7.0};
    const double deriv[4][3] =
    {
      { -1.0/2.0,  1.0/2.0, 0.0     },
      { -1.0/2.0,  1.0/2.0, 0.0     },
      {      0.0, -1.0/5.0, 1.0/5.0 },
      {      0.0, -1.0/5.0, 1.0/5.0 }
    };

    gsl_bspline_workspace *bw = gsl_bspline_alloc(2, 3);
    gsl_bspline_deriv_workspace *dbw = gsl_bspline_deriv_alloc(2);
    gsl_matrix *dB = gsl_matrix_alloc(gsl_bspline_ncoeffs(bw),
                                      gsl_bspline_order(bw) + 1);

    gsl_vector *breakpts = gsl_vector_alloc(3);
    gsl_vector_set(breakpts, 0, 0.0);
    gsl_vector_set(breakpts, 1, 2.0);
    gsl_vector_set(breakpts, 2, 7.0);
    gsl_bspline_knots(breakpts, bw);

    for (i = 0; i < 4; ++i)  /* at each location */
      {
        /* Initialize dB with poison to ensure we overwrite it */
        gsl_matrix_set_all(dB, GSL_NAN);

        gsl_bspline_deriv_eval(xloc[i], gsl_bspline_order(bw), dB, bw, dbw);
        for (j = 0; j < gsl_bspline_ncoeffs(bw) ; ++j)
          {
            /* check basis function 1st deriv */
            gsl_test_abs(gsl_matrix_get(dB, j, 1), deriv[i][j], GSL_DBL_EPSILON,
                         "b-spline k=%d basis #%d derivative %d at x = %f",
                         gsl_bspline_order(bw), j, 1, xloc[i]);
          }
        for (j = 0; j < gsl_bspline_ncoeffs(bw); ++j)
          {
            /* check k order basis function has k-th deriv equal to 0 */
            gsl_test_abs(gsl_matrix_get(dB, j, gsl_bspline_order(bw)), 0.0,
                         GSL_DBL_EPSILON,
                         "b-spline k=%d basis #%d derivative %d at x = %f",
                         gsl_bspline_order(bw), j, gsl_bspline_order(bw),
                         xloc[i]);
          }
      }

    gsl_matrix_free(dB);
    gsl_bspline_deriv_free(dbw);
    gsl_bspline_free(bw);
    gsl_vector_free(breakpts);
  }

  /* Spot check known 0th, 1st, 2nd derivative
     evaluations for a particular k = 3 case.  */
  {
    size_t i, j; /* looping */
    const double xloc[5]     =  { 0.0, 5.0, 9.0, 12.0, 15.0 };
    const double eval[5][6]  =
    {
      { 4./25.,  69./100.,   3./ 20. ,  0.    , 0.   , 0.    },
      { 0.     ,  4./21. , 143./210. ,  9./70., 0.   , 0.    },
      { 0.     ,  0.     ,   3./ 10. ,  7./10., 0.   , 0.    },
      { 0.     ,  0.     ,   0.      ,  3./4. , 1./4., 0.    },
      { 0.     ,  0.     ,   0.      ,  1./3. , 5./9., 1./9. }
    };
    const double deriv[5][6] =
    {
      { -4./25.,  3./50.,   1./ 10.,  0.    , 0.    , 0.      },
      {  0.    , -2./21.,   1./105.,  3./35., 0.    , 0.      },
      {  0.    ,  0.    ,  -1./5.  ,  1./ 5., 0.    , 0.      },
      {  0.    ,  0.    ,   0.     , -1./ 6., 1./6. , 0.      },
      {  0.    ,  0.    ,   0.     , -1./ 9., 1./27., 2./27. }
    };
    const double deriv2[5][6] =
    {
      { 2./25., -17./150.,   1.0/30.0 ,  0.0     ,  0.     , 0.     },
      { 0.    ,   1./ 42., -11.0/210.0,  1.0/35.0,  0.     , 0.     },
      { 0.    ,   0.     ,   1.0/15.0 ,-11.0/90.0,  1./18. , 0.     },
      { 0.    ,   0.     ,   0.0      ,  1.0/54.0, -7./162., 2./81. },
      { 0.    ,   0.     ,   0.0      ,  1.0/54.0, -7./162., 2./81. }
    };

    gsl_bspline_workspace *bw = gsl_bspline_alloc(3, 5);
    gsl_bspline_deriv_workspace *dbw = gsl_bspline_deriv_alloc(3);

    gsl_matrix *dB = gsl_matrix_alloc(gsl_bspline_ncoeffs(bw),
                                      gsl_bspline_order(bw) + 1);

    gsl_vector *breakpts = gsl_vector_alloc(5);
    gsl_vector_set(breakpts, 0, -3.0);
    gsl_vector_set(breakpts, 1,  2.0);
    gsl_vector_set(breakpts, 2,  9.0);
    gsl_vector_set(breakpts, 3, 12.0);
    gsl_vector_set(breakpts, 4, 21.0);
    gsl_bspline_knots(breakpts, bw);

    for (i = 0; i < 5; ++i)  /* at each location */
      {
        /* Initialize dB with poison to ensure we overwrite it */
        gsl_matrix_set_all(dB, GSL_NAN);
        gsl_bspline_deriv_eval(xloc[i], gsl_bspline_order(bw), dB, bw, dbw);

        /* check basis function evaluation */
        for (j = 0; j < gsl_bspline_ncoeffs(bw); ++j)
          {
            gsl_test_abs(gsl_matrix_get(dB, j, 0), eval[i][j], GSL_DBL_EPSILON,
                         "b-spline k=%d basis #%d derivative %d at x = %f",
                         gsl_bspline_order(bw), j, 0, xloc[i]);
          }
        /* check 1st derivative evaluation */
        for (j = 0; j < gsl_bspline_ncoeffs(bw); ++j)
          {
            gsl_test_abs(gsl_matrix_get(dB, j, 1), deriv[i][j], GSL_DBL_EPSILON,
                         "b-spline k=%d basis #%d derivative %d at x = %f",
                         gsl_bspline_order(bw), j, 1, xloc[i]);
          }
        /* check 2nd derivative evaluation */
        for (j = 0; j < gsl_bspline_ncoeffs(bw); ++j)
          {
            gsl_test_abs(gsl_matrix_get(dB, j, 2), deriv2[i][j], GSL_DBL_EPSILON,
                         "b-spline k=%d basis #%d derivative %d at x = %f",
                         gsl_bspline_order(bw), j, 2, xloc[i]);
          }
      }

    gsl_matrix_free(dB);
    gsl_bspline_deriv_free(dbw);
    gsl_bspline_free(bw);
    gsl_vector_free(breakpts);
  }

  /* Check Greville abscissae functionality on a non-uniform k=1 */
  {
    size_t i; /* looping */

    /* Test parameters */
    const size_t k = 1;
    const double bpoint_data[]    = { 0.0, 0.2, 0.5, 0.75, 1.0 };
    const size_t nbreak           = sizeof(bpoint_data)/sizeof(bpoint_data[0]);

    /* Expected results */
    const double abscissae_data[] = { 0.1, 0.35, 0.625, 0.875 };
    const size_t nabscissae       = sizeof(abscissae_data)/sizeof(abscissae_data[0]);

    gsl_vector_const_view bpoints = gsl_vector_const_view_array(bpoint_data, nbreak);
    gsl_bspline_workspace *w = gsl_bspline_alloc(k, nbreak);
    gsl_bspline_knots((const gsl_vector *) &bpoints, w);

    gsl_test_int(nabscissae, gsl_bspline_ncoeffs(w),
        "b-spline k=%d number of abscissae", k);
    for (i = 0; i < nabscissae; ++i)
      {
        gsl_test_abs(gsl_bspline_greville_abscissa(i, w), abscissae_data[i], 2*k*GSL_DBL_EPSILON,
            "b-spline k=%d Greville abscissa #%d at x = %f", k, i, abscissae_data[i]);
      }

    gsl_bspline_free(w);
  }

  /* Check Greville abscissae functionality on a non-uniform k=2 */
  {
    size_t i; /* looping */

    /* Test parameters */
    const size_t k = 2;
    const double bpoint_data[]    = { 0.0, 0.2, 0.5, 0.75, 1.0 };
    const size_t nbreak           = sizeof(bpoint_data)/sizeof(bpoint_data[0]);

    /* Expected results */
    const double abscissae_data[] = { 0.0, 0.2, 0.5, 0.75, 1.0 };
    const size_t nabscissae       = sizeof(abscissae_data)/sizeof(abscissae_data[0]);

    gsl_vector_const_view bpoints = gsl_vector_const_view_array(bpoint_data, nbreak);
    gsl_bspline_workspace *w = gsl_bspline_alloc(k, nbreak);
    gsl_bspline_knots((const gsl_vector *) &bpoints, w);

    gsl_test_int(nabscissae, gsl_bspline_ncoeffs(w),
        "b-spline k=%d number of abscissae", k);
    for (i = 0; i < nabscissae; ++i)
      {
        gsl_test_abs(gsl_bspline_greville_abscissa(i, w), abscissae_data[i], 2*k*GSL_DBL_EPSILON,
            "b-spline k=%d Greville abscissa #%d at x = %f", k, i, abscissae_data[i]);
      }

    gsl_bspline_free(w);
  }

  /* Check Greville abscissae functionality on non-uniform k=3 */
  {
    size_t i; /* looping */

    /* Test parameters */
    const size_t k = 3;
    const double bpoint_data[]    = { 0.0, 0.2, 0.5, 0.75, 1.0 };
    const size_t nbreak           = sizeof(bpoint_data)/sizeof(bpoint_data[0]);

    /* Expected results */
    const double abscissae_data[] = {      0.0, 1.0/10.0, 7.0/20.0,
                                      5.0/ 8.0, 7.0/ 8.0,      1.0 };
    const size_t nabscissae       = sizeof(abscissae_data)/sizeof(abscissae_data[0]);

    gsl_vector_const_view bpoints = gsl_vector_const_view_array(bpoint_data, nbreak);
    gsl_bspline_workspace *w = gsl_bspline_alloc(k, nbreak);
    gsl_bspline_knots((const gsl_vector *) &bpoints, w);

    gsl_test_int(nabscissae, gsl_bspline_ncoeffs(w),
        "b-spline k=%d number of abscissae", k);
    for (i = 0; i < nabscissae; ++i)
      {
        gsl_test_abs(gsl_bspline_greville_abscissa(i, w), abscissae_data[i], 2*k*GSL_DBL_EPSILON,
            "b-spline k=%d Greville abscissa #%d at x = %f", k, i, abscissae_data[i]);
      }

    gsl_bspline_free(w);
  }

  /* Check Greville abscissae functionality on non-uniform k=4 */
  {
    size_t i; /* looping */

    /* Test parameters */
    const size_t k = 4;
    const double bpoint_data[]    = { 0.0, 0.2, 0.5, 0.75, 1.0 };
    const size_t nbreak           = sizeof(bpoint_data)/sizeof(bpoint_data[0]);

    /* Expected results */
    const double abscissae_data[] = { 0.0,  1.0/15.0,  7.0/30.0,  29.0/60.0,
                                            3.0/ 4.0, 11.0/12.0,        1.0 };
    const size_t nabscissae       = sizeof(abscissae_data)/sizeof(abscissae_data[0]);

    gsl_vector_const_view bpoints = gsl_vector_const_view_array(bpoint_data, nbreak);
    gsl_bspline_workspace *w = gsl_bspline_alloc(k, nbreak);
    gsl_bspline_knots((const gsl_vector *) &bpoints, w);

    gsl_test_int(nabscissae, gsl_bspline_ncoeffs(w),
        "b-spline k=%d number of abscissae", k);
    for (i = 0; i < nabscissae; ++i)
      {
        gsl_test_abs(gsl_bspline_greville_abscissa(i, w), abscissae_data[i], 2*k*GSL_DBL_EPSILON,
            "b-spline k=%d Greville abscissa #%d at x = %f", k, i, abscissae_data[i]);
      }

    gsl_bspline_free(w);
  }

  /* Knots computed from prescribed Greville abscissae for k = 4 */
  {
    size_t i; /* looping */

    /* Test parameters */
    const size_t k = 4;
    const double abscissae_data[] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 };
    const size_t nabscissae       = sizeof(abscissae_data)/sizeof(abscissae_data[0]);

    /* Expected results */
    const double bpoint_data[]    = { 1.0, 4.0, 4.0, 4.0, 7.0 };
    const size_t nbreak           = sizeof(bpoint_data)/sizeof(bpoint_data[0]);

    /* Compute knots from Greville abscissae */
    double abserr;
    gsl_vector_const_view abscissae
        = gsl_vector_const_view_array(abscissae_data, nabscissae);
    gsl_bspline_workspace *w = gsl_bspline_alloc(k, nbreak);
    gsl_bspline_knots_greville(&abscissae.vector, w, &abserr);

    for (i = 0; i < nbreak; ++i)
      {
        gsl_test_abs(gsl_bspline_breakpoint(i,w), bpoint_data[i], GSL_DBL_EPSILON*50,
            "b-spline k=%d knots_greville breakpoint #%d", k, i);
      }
    gsl_test_abs(abserr, 0.0, GSL_DBL_EPSILON*15,
        "b-spline k=%d nbreak=%d knots_greville abserr", k, nbreak);

    gsl_bspline_free(w);
  }

  /* Knots computed from prescribed Greville abscissae for k = 8 */
  {
    size_t i; /* looping */

    /* Test parameters */
    const size_t k = 8;

    const double abscissae_data[] = { 1.0, 10.0/7, 13.0/7, 16.0/7, 22.0/7,
                                      4.0, 34.0/7, 40.0/7, 43.0/7, 46.0/7, 7.0 };
    const size_t nabscissae       = sizeof(abscissae_data)/sizeof(abscissae_data[0]);

    /* Expected results */
    const double bpoint_data[]    = { 1.0, 4.0, 4.0, 4.0, 7.0 };
    const size_t nbreak           = sizeof(bpoint_data)/sizeof(bpoint_data[0]);

    /* Compute knots from Greville abscissae */
    double abserr;
    gsl_vector_const_view abscissae
        = gsl_vector_const_view_array(abscissae_data, nabscissae);
    gsl_bspline_workspace *w = gsl_bspline_alloc(k, nbreak);
    gsl_bspline_knots_greville(&abscissae.vector, w, &abserr);

    for (i = 0; i < nbreak; ++i)
      {
        gsl_test_abs(gsl_bspline_breakpoint(i,w), bpoint_data[i], GSL_DBL_EPSILON*50,
            "b-spline k=%d knots_greville breakpoint #%d", k, i);
      }
    gsl_test_abs(abserr, 0.0, GSL_DBL_EPSILON*15,
        "b-spline k=%d nbreak=%d knots_greville abserr", k, nbreak);

    gsl_bspline_free(w);
  }

  /* Knots computed from prescribed Greville abscissae for k = 2 */
  /* Not an interesting calculation but checks the k = 2 edge case */
  {
    size_t i; /* looping */

    /* Test parameters */
    const size_t k = 2;
    const double abscissae_data[] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 };
    const size_t nabscissae       = sizeof(abscissae_data)/sizeof(abscissae_data[0]);

    /* Expected results */
    const double bpoint_data[]    = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 };
    const size_t nbreak           = sizeof(bpoint_data)/sizeof(bpoint_data[0]);

    /* Compute knots from Greville abscissae */
    double abserr;
    gsl_vector_const_view abscissae
        = gsl_vector_const_view_array(abscissae_data, nabscissae);
    gsl_bspline_workspace *w = gsl_bspline_alloc(k, nbreak);
    gsl_bspline_knots_greville(&abscissae.vector, w, &abserr);

    for (i = 0; i < nbreak; ++i)
      {
        gsl_test_abs(gsl_bspline_breakpoint(i,w), bpoint_data[i], GSL_DBL_EPSILON,
            "b-spline k=%d knots_greville breakpoint #%d", k, i);
      }
    gsl_test_abs(abserr, 0.0, GSL_DBL_EPSILON,
        "b-spline k=%d nbreak=%d knots_greville abserr", k, nbreak);

    gsl_bspline_free(w);
  }

  /* Knots computed from prescribed abscissae for edge case when nbreak = 2 */
  {
    size_t i; /* looping */

    /* Test parameters */
    const size_t k = 4;
    double abscissae_data[] = { 1.0, 3.0, 5.0, 7.0 };
    const size_t nabscissae = sizeof(abscissae_data)/sizeof(abscissae_data[0]);

    /* Expected results */
    const double bpoint_data[] = { 1.0, 7.0 };
    const size_t nbreak        = sizeof(bpoint_data)/sizeof(bpoint_data[0]);

    /* Compute knots from Greville abscissae where abscissae are recoverable */
    double abserr;
    gsl_vector_view abscissae
        = gsl_vector_view_array(abscissae_data, nabscissae);
    gsl_bspline_workspace *w = gsl_bspline_alloc(k, nbreak);
    gsl_bspline_knots_greville(&abscissae.vector, w, &abserr);

    /* Check recovery of breakpoints and abscissae */
    for (i = 0; i < nbreak; ++i)
      {
        gsl_test_abs(gsl_bspline_breakpoint(i,w), bpoint_data[i], GSL_DBL_EPSILON,
            "b-spline k=%d knots_greville breakpoint #%d", k, i);
      }
    gsl_test_abs(abserr, 0.0, GSL_DBL_EPSILON,
        "b-spline k=%d nbreak=%d knots_greville abserr", k, nbreak);

    /* Modify interior abscissae so they cannot be recovered with nbreak = 2 */
    /* Then recompute breakpoints and check that abserr is as expected */
    abscissae_data[1] -= 1;
    abscissae_data[2] += 1;
    gsl_bspline_knots_greville(&abscissae.vector, w, &abserr);
    for (i = 0; i < nbreak; ++i)
      {
        gsl_test_abs(gsl_bspline_breakpoint(i,w), bpoint_data[i], GSL_DBL_EPSILON,
            "b-spline k=%d knots_greville breakpoint #%d", k, i);
      }
    gsl_test_abs(abserr, /* deliberate error */ 2.0, GSL_DBL_EPSILON,
        "b-spline k=%d nbreak=%d knots_greville abserr large", k, nbreak);

    gsl_bspline_free(w);
  }

  exit(gsl_test_summary());
}

max_stars_repo_name: ruslankuzmin/julia

max_issues_repo_name: ruslankuzmin/julia

max_forks_repo_name: ruslankuzmin/julia


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header for functions windowing and computing FFT/IFFT of time/frequency series.
 *
 */

#ifndef __GENERATEWAVEFORM_H__
#define __GENERATEWAVEFORM_H__ 1

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#define _XOPEN_SOURCE 500

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "struct.h"
#include "waveform.h"
#include "EOBNRv2HMROMstruct.h"
#include "EOBNRv2HMROM.h"
#include "fft.h"

/********************************** Structures ******************************************/

typedef enum GenWavetag {
  hlm,
  h22TD,
  hphcFD,
  hphcTD
} GenWavetag;

/* Parameters for the generation of a LISA waveform (in the form of a list of modes) */
typedef struct tagGenWaveParams {
  double tRef;               /* reference time (s) - GPS time at the frequency representing coalescence */
  double phiRef;             /* reference phase (rad) - phase at the frequency representing coalescence (or at fRef if specified) */
  double fRef;               /* reference frequency at which phiRef is set (Hz, default 0 which is interpreted as Mf=0.14) */
  double m1;                 /* mass of companion 1 (solar masses, default 2e6) */
  double m2;                 /* mass of companion 2 (solar masses, default 1e6) */
  double distance;           /* distance of source (Mpc, default 1e3) */
  double inclination;        /* inclination of source (rad, default pi/3) */
  double minf;               /* Minimal frequency, ignore if 0 (Hz, default=0) - will use first frequency covered by the ROM if higher */
  double maxf;               /* Maximal frequency, ignore if 0 (Hz, default=0) - will use last frequency covered by the ROM if lower */
  double deltatobs;          /* Observation duration (years, default=2) */
  int tagextpn;              /* Tag to allow PN extension of the waveform at low frequencies */
  double Mfmatch;            /* When PN extension allowed, geometric matching frequency: will use ROM above this value. If <=0, use ROM down to the lowest covered frequency */
  int nbmode;                /* Number of modes to generate (starting with 22) - defaults to 1 (22 only) */
  int taggenwave;            /* Tag selecting the desired output format */
  double f1windowbeg;        /* If generating h22TD/hphcTD, start frequency for windowing at the beginning - set to 0 to ignore and use max(fstartobs, fLowROM, minf), where fLowROM is either the lowest frequency covered by the ROM or simply minf if PN extension is used (Hz, default=0) */
  double f2windowbeg;        /* If generating h22TD/hphcTD, stop frequency for windowing at the beginning - set to 0 to ignore and use 1.1*f1windowbeg (Hz, default=0) */
  double f1windowend;        /* If generating h22TD/hphcTD, start frequency for windowing at the end - set to 0 to ignore and use 0.995*f2windowend (Hz, default=0) */
  double f2windowend;        /* If generating h22TD/hphcTD, stop frequency for windowing at the end - set to 0 to ignore and use min(maxf, fHighROM), where fHighROM is the highest frequency covered by the ROM (Hz, default=0) */
  int tagh22fromfile;        /* Tag choosing wether to load h22 FD downsampled Amp/Phase from file (default 0) */
  int nsamplesinfile;        /* Number of lines of inputs file */
  int binaryin;              /* Tag for loading the data in gsl binary form instead of text (default false) */
  char indir[256];           /* Input directory */
  char infile[256];          /* Input file name */
  int binaryout;             /* Tag for outputting the data in gsl binary form instead of text (default 0) */
  char outdir[256];          /* Path for the output directory */
  char outfile[256];         /* Path for the output file */
} GenWaveParams;

#if 0
{ /* so that editors will match succeeding brace */
#elif defined(__cplusplus)
}
#endif

#endif /* _GENERATEWAVEFORM_H */

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
#include <Python.h>
#include <arrayobject.h>
#include <time.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_sf_erf.h>

static PyObject* g_engineError;
static gsl_rng* g_rng;

/** low level helpers **/

/** numpy helpers **/

/**
 *  Check that PyArrayObject is a double (Float) type and a vector
 * 	return 1 if an error and raise exception.
 */
static int not_doublevector(PyArrayObject *vec)  {
	if (vec->descr->type_num != NPY_DOUBLE || vec->nd != 1)  {
		PyErr_SetString(PyExc_ValueError,
			"In not_doublevector: array must be of type Float and 1 dimensional (n).");
		return 1;
	}
	return 0;
}

static double lognormal(double m, double stddev) {
/*
m: mean of x
stddev: std of x
*/
	return gsl_ran_lognormal(g_rng,
			log(m) - 0.5 * log(1.0 + (stddev*stddev)/(m*m)),
			sqrt(log((stddev*stddev)/(m*m) + 1)) );
}

/** python callables **/

static PyObject* engine_bivariate_gaussian(PyObject *self, PyObject *args) {
	float sx, sy;
	if (!PyArg_ParseTuple(args, "ff", &sx, &sy))
        return NULL;
	double x, y;
	gsl_ran_bivariate_gaussian(g_rng, sx, sy, 0.0, &x, &y);
	return Py_BuildValue("(ff)", (float)x, (float)y);
}

static PyObject* engine_beta(PyObject *self, PyObject *args) {
	float a = 0;
	float b = 0;
	if (!PyArg_ParseTuple(args, "ff", &a, &b))
        return NULL;
	double sample = (float)gsl_ran_beta(g_rng, a, b);
	return Py_BuildValue("f", (float)sample);
}

static PyObject* engine_lognormal(PyObject *self, PyObject *args) {
	float m = 0;
	float stddev = 0;
	if (!PyArg_ParseTuple(args, "ff", &m, &stddev))
        return NULL;
	return Py_BuildValue("f", (float)lognormal(m, stddev));
}

static PyObject* engine_lognormal_cdf(PyObject *self, PyObject *args) {
	float m = 0;
	float stddev = 0;
	float x = 0;
	if (!PyArg_ParseTuple(args, "fff", &x, &m, &stddev))
        return NULL;
	double m2 = log(m) - 0.5 * log(1.0 + (stddev*stddev)/(m*m));
	double s2 = sqrt(log((stddev*stddev)/(m*m) + 1));
	return Py_BuildValue("f", (float)gsl_cdf_lognormal_P(x, m2, s2));
}

static PyObject* engine_lognormal_cdf_inv(PyObject *self, PyObject *args) {
	float m = 0;
	float stddev = 0;
	float p = 0;
	if (!PyArg_ParseTuple(args, "fff", &p, &m, &stddev))
		return NULL;
	double m2 = log(m) - 0.5 * log(1.0 + (stddev*stddev)/(m*m));
	double s2 = sqrt(log((stddev*stddev)/(m*m) + 1));
	return Py_BuildValue("f", (float)gsl_cdf_lognormal_Pinv(p, m2, s2));
}

static PyObject* engine_percentileofscore(PyObject *self, PyObject *args) {
	// python: engine.percentileofscore([scores], score)
	double score = 0;
	PyObject* listObj;
	if (!PyArg_ParseTuple(args, "O!d", &PyList_Type, &listObj, &score))
        return NULL;

	// loop through python list of doubles, counting scores less than or equal to score
	int arrsz = PyList_Size(listObj);
	if (arrsz == 0)
		return NULL;
	int i;
	int count=0;
	for (i=0; i<arrsz; i++) {
		double aScore = PyFloat_AsDouble(PyList_GetItem(listObj, i));
		if (aScore <= score)
			count++;
	}
	float percentile = count / (float)arrsz;
	return Py_BuildValue("f", percentile);
}

static PyObject* engine_percentileofscore_array(PyObject *self, PyObject *args) {
	// python: engine.percentileofscore(numpy.array[scores], score)
	double score = 0;
	PyArrayObject* arrayObj;
	if (!PyArg_ParseTuple(args, "O!d", &PyArray_Type, &arrayObj, &score))
        return NULL;
	if (not_doublevector(arrayObj)) {
		printf("invalid type!\n");
		return NULL;
	}

	// loop through python list of doubles, counting scores less than or equal to score
	double* cin = (double *)arrayObj->data;
	int arrsz = arrayObj->dimensions[0];
	if (arrsz == 0)
		return NULL;
	int i;
	int count=0;
	for (i=0; i<arrsz; i++) {
		double aScore = cin[i];
		printf("score: %lf\n", aScore);
		if (aScore <= score)
			count++;
	}
	float percentile = count / (float)arrsz;
	return Py_BuildValue("f", percentile);
}

static PyObject* engine_erf(PyObject *self, PyObject *args) {
	float x = 0;
	if (!PyArg_ParseTuple(args, "f", &x))
        return NULL;
	double sample = gsl_sf_erf(x);
	return Py_BuildValue("f", (float)sample);
}

static PyObject* engine_poisson(PyObject *self, PyObject *args) {
	float mean = 0;
	if (!PyArg_ParseTuple(args, "f", &mean))
        return NULL;
	unsigned int sample = gsl_ran_poisson(g_rng, mean);
	return Py_BuildValue("i", (int)sample);
}

static PyObject* engine_seed(PyObject *self, PyObject *args) {
	long seed = 0;
	if (!PyArg_ParseTuple(args, "l", &seed))
        return NULL;
	gsl_rng_set(g_rng, seed);
	Py_INCREF(Py_None);
	return Py_None;
}

/** debris module **/

static PyObject* debris_generate_item(PyObject* self, PyObject* args) {
	int typeid;
	float V, xord, yord, flighttime_mean, flighttime_stddev, cdav, mass_mean, mass_stddev, fa_mean, fa_stddev;
	if (!PyArg_ParseTuple(args, "fffffifffff",
			&V,
			&xord, &yord,
			&flighttime_mean,
			&flighttime_stddev,
			&typeid,
			&cdav,
			&mass_mean, &mass_stddev,
			&fa_mean, &fa_stddev))
		return NULL;

	double mass = lognormal(mass_mean, mass_stddev);
	double fa = lognormal(fa_mean, fa_stddev);
	double flight_time = lognormal(flighttime_mean, flighttime_stddev);

	double t = 9.81 * flight_time / V;
	double k = (1.2 * V * V) / (2 * 9.81 * (mass/fa) );
	double flight_distance=0;

	switch (typeid) {
	case 0:
		flight_distance = (V*V/9.81) * (1.0/k) * (0.2060*pow(k*t, 2) + 0.011*k*t);
		break;
	case 1:
		flight_distance = (V*V/9.81) * (1.0/k) * (0.072*pow(k*t, 2) + 0.3456*k*t);
		break;
	case 2:
		flight_distance = (V*V/9.81) * (1.0/k) * (0.0723*pow(k*t, 2) + 0.2376*k*t);
		break;
	}

	// Sample Impact Location
	double sigma_x = flight_distance / 3.0;
	double sigma_y = flight_distance / 12.0;
	double x, y;
	gsl_ran_bivariate_gaussian(g_rng, sigma_x, sigma_y, 1.0, &x, &y);

	// translate to footprint coords
	double item_X = xord - flight_distance + x;
	double item_Y = yord + y;

	// calculate momentum
	double pa = 1.2;
	double b = sqrt( (pa * cdav * fa) / (mass) );
	double beta_mean = 1 - exp(-b * sqrt(flight_distance));

    double dispersion = 0;
    if (beta_mean == 1.0) {
    	dispersion = 1/beta_mean + 3.0;
    	beta_mean -= 0.001;
    }
    else if (beta_mean == 0) {
    	dispersion = 4.0;
    }
    else {
        dispersion = fmax(1.0/beta_mean, 1.0/(1.0 - beta_mean)) + 3.0;
    }

    double beta_a = beta_mean * dispersion;
    double beta_b = dispersion * (1 - beta_mean);

    double sampled = gsl_ran_beta(g_rng, beta_a, beta_b);
    float item_momentum = (float)(sampled * V * mass);

	/* output: (item_momentum, item_X, item_Y)
	 */
	return Py_BuildValue("(fff)", item_momentum, item_X, item_Y);
}

/** Module setup **/

static PyMethodDef EngineMethods[] = {
	{"beta", engine_beta, METH_VARARGS, "Draw from Beta using given a,b"},
	{"lognormal", engine_lognormal, METH_VARARGS, "Draw from Lognormal using given mean, cov"},
	{"lognormal_cdf", engine_lognormal_cdf, METH_VARARGS, "lognormal cdf"},
	{"lognormal_cdf_inv", engine_lognormal_cdf_inv, METH_VARARGS, "lognormal inv cdf"},
	{"erf", engine_erf, METH_VARARGS, "erf"},
	{"bivariate_gaussian", engine_bivariate_gaussian, METH_VARARGS, "Bivariate normal distribution"},
	{"poisson", engine_poisson, METH_VARARGS, "Draw from Poisson using given mean"},
	{"percentileofscore", engine_percentileofscore, METH_VARARGS, "Return a quantile value of sorted_data"},
	{"percentileofscore_array", engine_percentileofscore_array, METH_VARARGS, "Return a quantile value of sorted_data"},
	{"seed", engine_seed, METH_VARARGS, "Set random number seed"},
	{"debris_generate_item", debris_generate_item, METH_VARARGS, "Generate debris item for source"},
	{NULL, NULL, 0, NULL}
};

PyMODINIT_FUNC initengine(void) {
	PyObject* m = Py_InitModule("engine", EngineMethods);
	if (m == NULL)
		return;
	import_array();

	g_engineError = PyErr_NewException("engine.error", NULL, NULL);
	Py_INCREF(g_engineError);
	PyModule_AddObject(m, "error", g_engineError);

	// setup GSL random number generators
	const gsl_rng_type* T;
	gsl_rng_env_setup();
	T = gsl_rng_default;
	g_rng = gsl_rng_alloc(T);
	gsl_rng_set(g_rng, (long)time(NULL));
}

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#ifndef REMESHER_GSL_HELPERS_HEADER
#define REMESHER_GSL_HELPERS_HEADER

#include <iostream>
#include <iomanip>
#include <string>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

namespace GSL
{
  void
  pretty_print_matrix (gsl_matrix* m, std::string const& name)
  {
    std::cout << name << " = [";
    for (std::size_t i = 0; i < m->size2; ++i)
    {
      if (i != 0) std::cout << ";" << std::endl << "     ";
      for (std::size_t j = 0; j < m->size1; ++j)
      {
        if (j != 0) std::cout << ", ";
        std::cout << std::setw(6) << std::fixed << std::setprecision(2)
            << gsl_matrix_get(m, i, j);
      }
    }
    std::cout << "];" << std::endl;
  }

  void
  pretty_print_vector (gsl_vector* v, std::string const& name)
  {
    std::cout << name << " = [";
    for (std::size_t i = 0; i < v->size; ++i)
    {
      if (i != 0) std::cout << "; ";
      std::cout << std::setw(6) << std::fixed << std::setprecision(2)
            << gsl_vector_get(v, i);
    }
    std::cout << "];" << std::endl;
  }

}

#endif /* REMESHER_GSL_HELPERS_HEADER */


text:
/* gmres.c
 *
 * Copyright (C) 2014 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_spmatrix.h>
#include <gsl/gsl_spblas.h>
#include <gsl/gsl_splinalg.h>

/*
 * The code in this module is based on the Householder GMRES
 * algorithm described in
 *
 * [1] H. F. Walker, Implementation of the GMRES method using
 *     Householder transformations, SIAM J. Sci. Stat. Comput.
 *     9(1), 1988.
 *
 * [2] Y. Saad, Iterative methods for sparse linear systems,
 *     2nd edition, SIAM, 2003.
 */

typedef struct
{
  size_t n;        /* size of linear system */
  size_t m;        /* dimension of Krylov subspace K_m */
  gsl_vector *r;   /* residual vector r = b - A*x */
  gsl_matrix *H;   /* Hessenberg matrix n-by-(m+1) */
  gsl_vector *tau; /* householder scalars */
  gsl_vector *y;   /* least squares rhs and solution vector */

  double *c;       /* Givens rotations */
  double *s;

  double normr;    /* residual norm ||r|| */
} gmres_state_t;

static void gmres_free(void *vstate);
static int gmres_iterate(const gsl_spmatrix *A, const gsl_vector *b,
                         const double tol, gsl_vector *x, void *vstate);

/*
gmres_alloc()
  Allocate a GMRES workspace for solving an n-by-n system A x = b

Inputs: n        - size of system
        krylov_m - size of Krylov subspace (ie: number of inner iterations)
                   if this parameter is 0, the value GSL_MIN(n,10) is
                   used

Return: pointer to workspace
*/

static void *
gmres_alloc(const size_t n, const size_t m)
{
  gmres_state_t *state;

  if (n == 0)
    {
      GSL_ERROR_NULL("matrix dimension n must be a positive integer",
                     GSL_EINVAL);
    }

  state = calloc(1, sizeof(gmres_state_t));
  if (!state)
    {
      GSL_ERROR_NULL("failed to allocate gmres state", GSL_ENOMEM);
    }

  state->n = n;

  /* compute size of Krylov subspace */
  if (m == 0)
    state->m = GSL_MIN(n, 10);
  else
    state->m = GSL_MIN(n, m);

  state->r = gsl_vector_alloc(n);
  if (!state->r)
    {
      gmres_free(state);
      GSL_ERROR_NULL("failed to allocate r vector", GSL_ENOMEM);
    }

  state->H = gsl_matrix_alloc(n, state->m + 1);
  if (!state->H)
    {
      gmres_free(state);
      GSL_ERROR_NULL("failed to allocate H matrix", GSL_ENOMEM);
    }

  state->tau = gsl_vector_alloc(state->m + 1);
  if (!state->tau)
    {
      gmres_free(state);
      GSL_ERROR_NULL("failed to allocate tau vector", GSL_ENOMEM);
    }

  state->y = gsl_vector_alloc(state->m + 1);
  if (!state->y)
    {
      gmres_free(state);
      GSL_ERROR_NULL("failed to allocate y vector", GSL_ENOMEM);
    }

  state->c = malloc(state->m * sizeof(double));
  state->s = malloc(state->m * sizeof(double));
  if (!state->c || !state->s)
    {
      gmres_free(state);
      GSL_ERROR_NULL("failed to allocate Givens vectors", GSL_ENOMEM);
    }

  state->normr = 0.0;

  return state;
} /* gmres_alloc() */

static void
gmres_free(void *vstate)
{
  gmres_state_t *state = (gmres_state_t *) vstate;

  if (state->r)
    gsl_vector_free(state->r);

  if (state->H)
    gsl_matrix_free(state->H);

  if (state->tau)
    gsl_vector_free(state->tau);

  if (state->y)
    gsl_vector_free(state->y);

  if (state->c)
    free(state->c);

  if (state->s)
    free(state->s);

  free(state);
} /* gmres_free() */

/*
gmres_iterate()
  Solve A*x = b using GMRES algorithm

Inputs: A    - sparse square matrix
        b    - right hand side vector
        tol  - stopping tolerance (see below)
        x    - (input/output) on input, initial estimate x_0;
               on output, solution vector
        work - workspace

Return:
GSL_SUCCESS if converged to solution (solution stored in x). In
this case the following will be true:

||b - A*x|| <= tol * ||b||

GSL_CONTINUE if not yet converged; in this case x contains the
most recent solution vector and calling this function more times
with the input x could result in convergence (ie: restarted GMRES)

Notes:
1) Based on algorithm 2.2 of (Walker, 1998 [1]) and algorithm 6.10 of
(Saad, 2003 [2])

2) On output, work->normr contains ||b - A*x||
*/

static int
gmres_iterate(const gsl_spmatrix *A, const gsl_vector *b,
              const double tol, gsl_vector *x,
              void *vstate)
{
  const size_t N = A->size1;
  gmres_state_t *state = (gmres_state_t *) vstate;

  if (N != A->size2)
    {
      GSL_ERROR("matrix must be square", GSL_ENOTSQR);
    }
  else if (N != b->size)
    {
      GSL_ERROR("matrix does not match right hand side", GSL_EBADLEN);
    }
  else if (N != x->size)
    {
      GSL_ERROR("matrix does not match solution vector", GSL_EBADLEN);
    }
  else if (N != state->n)
    {
      GSL_ERROR("matrix does not match workspace", GSL_EBADLEN);
    }
  else
    {
      int status = GSL_SUCCESS;
      const size_t maxit = state->m;
      const double normb = gsl_blas_dnrm2(b); /* ||b|| */
      const double reltol = tol * normb;      /* tol*||b|| */
      double normr;                           /* ||r|| */
      size_t m, k;
      double tau;                             /* householder scalar */
      gsl_matrix *H = state->H;               /* Hessenberg matrix */
      gsl_vector *r = state->r;               /* residual vector */
      gsl_vector *w = state->y;               /* least squares RHS */
      gsl_matrix_view Rm;                     /* R_m = H(1:m,2:m+1) */
      gsl_vector_view ym;                     /* y(1:m) */
      gsl_vector_view h0 = gsl_matrix_column(H, 0);

      /*
       * The Hessenberg matrix will have the following structure:
       *
       * H = [ ||r_0|| | v_1 v_2 ... v_m     ]
       *     [   u_1   | u_2 u_3 ... u_{m+1} ]
       *
       * where v_j are the orthonormal vectors spanning the Krylov
       * subpsace of length j + 1 and u_{j+1} are the householder
       * vectors of length n - j - 1.
       * In fact, u_{j+1} has length n - j since u_{j+1}[0] = 1,
       * but this 1 is not stored.
       */
      gsl_matrix_set_zero(H);

      /* Step 1a: compute r = b - A*x_0 */
      gsl_vector_memcpy(r, b);
      gsl_spblas_dgemv(CblasNoTrans, -1.0, A, x, 1.0, r);

      /* Step 1b */
      gsl_vector_memcpy(&h0.vector, r);
      tau = gsl_linalg_householder_transform(&h0.vector);

      /* store tau_1 */
      gsl_vector_set(state->tau, 0, tau);

      /* initialize w (stored in state->y) */
      gsl_vector_set_zero(w);
      gsl_vector_set(w, 0, gsl_vector_get(&h0.vector, 0));

      for (m = 1; m <= maxit; ++m)
        {
          size_t j = m - 1; /* C indexing */
          double c, s;      /* Givens rotation */

          /* v_m */
          gsl_vector_view vm = gsl_matrix_column(H, m);

          /* v_m(m:end) */
          gsl_vector_view vv = gsl_vector_subvector(&vm.vector, j, N - j);

          /* householder vector u_m for projection P_m */
          gsl_vector_view um = gsl_matrix_subcolumn(H, j, j, N - j);

          /* Step 2a: form v_m = P_m e_m = e_m - tau_m w_m */
          gsl_vector_set_zero(&vm.vector);
          gsl_vector_memcpy(&vv.vector, &um.vector);
          tau = gsl_vector_get(state->tau, j); /* tau_m */
          gsl_vector_scale(&vv.vector, -tau);
          gsl_vector_set(&vv.vector, 0, 1.0 - tau);

          /* Step 2a: v_m <- P_1 P_2 ... P_{m-1} v_m */
          for (k = j; k > 0 && k--; )
            {
              gsl_vector_view uk =
                gsl_matrix_subcolumn(H, k, k, N - k);
              gsl_vector_view vk =
                gsl_vector_subvector(&vm.vector, k, N - k);
              tau = gsl_vector_get(state->tau, k);
              gsl_linalg_householder_hv(tau, &uk.vector, &vk.vector);
            }

          /* Step 2a: v_m <- A*v_m */
          gsl_spblas_dgemv(CblasNoTrans, 1.0, A, &vm.vector, 0.0, r);
          gsl_vector_memcpy(&vm.vector, r);

          /* Step 2a: v_m <- P_m ... P_1 v_m */
          for (k = 0; k <= j; ++k)
            {
              gsl_vector_view uk = gsl_matrix_subcolumn(H, k, k, N - k);
              gsl_vector_view vk = gsl_vector_subvector(&vm.vector, k, N - k);
              tau = gsl_vector_get(state->tau, k);
              gsl_linalg_householder_hv(tau, &uk.vector, &vk.vector);
            }

          /* Steps 2c,2d: find P_{m+1} and set v_m <- P_{m+1} v_m */
          if (m < N)
            {
              /* householder vector u_{m+1} for projection P_{m+1} */
              gsl_vector_view ump1 = gsl_matrix_subcolumn(H, m, m, N - m);

              tau = gsl_linalg_householder_transform(&ump1.vector);
              gsl_vector_set(state->tau, j + 1, tau);
            }

          /* Step 2e: v_m <- J_{m-1} ... J_1 v_m */
          for (k = 0; k < j; ++k)
            {
              gsl_linalg_givens_gv(&vm.vector, k, k + 1,
                                   state->c[k], state->s[k]);
            }

          if (m < N)
            {
              /* Step 2g: find givens rotation J_m for v_m(m:m+1) */
              gsl_linalg_givens(gsl_vector_get(&vm.vector, j),
                                gsl_vector_get(&vm.vector, j + 1),
                                &c, &s);

              /* store givens rotation for later use */
              state->c[j] = c;
              state->s[j] = s;

              /* Step 2h: v_m <- J_m v_m */
              gsl_linalg_givens_gv(&vm.vector, j, j + 1, c, s);

              /* Step 2h: w <- J_m w */
              gsl_linalg_givens_gv(w, j, j + 1, c, s);
            }

          /*
           * Step 2i: R_m = [ R_{m-1}, v_m ] - already taken care
           * of due to our memory storage scheme
           */

          /* Step 2j: check residual w_{m+1} for convergence */
          normr = fabs(gsl_vector_get(w, j + 1));
          if (normr <= reltol)
            {
              /*
               * method has converged, break out of loop to compute
               * update to solution vector x
               */
              break;
            }
        }

      /*
       * At this point, we have either converged to a solution or
       * completed all maxit iterations. In either case, compute
       * an update to the solution vector x and test again for
       * convergence.
       */

      /* rewind m if we exceeded maxit iterations */
      if (m > maxit)
        m--;

      /* Step 3a: solve triangular system R_m y_m = w, in place */
      Rm = gsl_matrix_submatrix(H, 0, 1, m, m);
      ym = gsl_vector_subvector(w, 0, m);
      gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit,
                     &Rm.matrix, &ym.vector);

      /*
       * Step 3b: update solution vector x; the loop below
       * uses a different but equivalent formulation from
       * Saad, algorithm 6.10, step 14; store Krylov projection
       * V_m y_m in 'r'
       */
      gsl_vector_set_zero(r);
      for (k = m; k > 0 && k--; )
        {
          double ymk = gsl_vector_get(&ym.vector, k);
          gsl_vector_view uk = gsl_matrix_subcolumn(H, k, k, N - k);
          gsl_vector_view rk = gsl_vector_subvector(r, k, N - k);

          /* r <- n_k e_k + r */
          gsl_vector_set(r, k, gsl_vector_get(r, k) + ymk);

          /* r <- P_k r */
          tau = gsl_vector_get(state->tau, k);
          gsl_linalg_householder_hv(tau, &uk.vector, &rk.vector);
        }

      /* x <- x + V_m y_m */
      gsl_vector_add(x, r);

      /* compute new residual r = b - A*x */
      gsl_vector_memcpy(r, b);
      gsl_spblas_dgemv(CblasNoTrans, -1.0, A, x, 1.0, r);
      normr = gsl_blas_dnrm2(r);

      if (normr <= reltol)
        status = GSL_SUCCESS;  /* converged */
      else
        status = GSL_CONTINUE; /* not yet converged */

      /* store residual norm */
      state->normr = normr;

      return status;
    }
} /* gmres_iterate() */

static double
gmres_normr(const void *vstate)
{
  const gmres_state_t *state = (const gmres_state_t *) vstate;
  return state->normr;
} /* gmres_normr() */

static const gsl_splinalg_itersolve_type gmres_type =
{
  "gmres",
  &gmres_alloc,
  &gmres_iterate,
  &gmres_normr,
  &gmres_free
};

const gsl_splinalg_itersolve_type * gsl_splinalg_itersolve_gmres =
  &gmres_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <stdio.h>
#include <gsl/gsl_rstat.h>

int
main(void)
{
  double data[5] = {17.2, 18.1, 16.5, 18.3, 12.6};
  double mean, variance, largest, smallest, sd,
         rms, sd_mean, median, skew, kurtosis;
  gsl_rstat_workspace *rstat_p = gsl_rstat_alloc();
  size_t i, n;

  /* add data to rstat accumulator */
  for (i = 0; i < 5; ++i)
    gsl_rstat_add(data[i], rstat_p);

  mean     = gsl_rstat_mean(rstat_p);
  variance = gsl_rstat_variance(rstat_p);
  largest  = gsl_rstat_max(rstat_p);
  smallest = gsl_rstat_min(rstat_p);
  median   = gsl_rstat_median(rstat_p);
  sd       = gsl_rstat_sd(rstat_p);
  sd_mean  = gsl_rstat_sd_mean(rstat_p);
  skew     = gsl_rstat_skew(rstat_p);
  rms      = gsl_rstat_rms(rstat_p);
  kurtosis = gsl_rstat_kurtosis(rstat_p);
  n        = gsl_rstat_n(rstat_p);

  printf ("The dataset is %g, %g, %g, %g, %g\n",
         data[0], data[1], data[2], data[3], data[4]);

  printf ("The sample mean is %g\n", mean);
  printf ("The estimated variance is %g\n", variance);
  printf ("The largest value is %g\n", largest);
  printf ("The smallest value is %g\n", smallest);
  printf( "The median is %g\n", median);
  printf( "The standard deviation is %g\n", sd);
  printf( "The root mean square is %g\n", rms);
  printf( "The standard devation of the mean is %g\n", sd_mean);
  printf( "The skew is %g\n", skew);
  printf( "The kurtosis %g\n", kurtosis);
  printf( "There are %zu items in the accumulator\n", n);

  gsl_rstat_reset(rstat_p);
  n = gsl_rstat_n(rstat_p);
  printf( "There are %zu items in the accumulator\n", n);

  gsl_rstat_free(rstat_p);

  return 0;
}


text:
#ifndef LOGGER_INCLUDE
#define LOGGER_INCLUDE

#include <gsl/string_span>

#include "log/log.h"
BOOST_LOG_GLOBAL_LOGGER(exec_helper_core_logger, execHelper::log::LoggerType);

static const gsl::czstring<> LOG_CHANNEL = "core";
#define LOG(x)                                                                 \
    BOOST_LOG_STREAM_CHANNEL_SEV(exec_helper_core_logger::get(), LOG_CHANNEL,  \
                                 execHelper::log::x)                           \
        << boost::log::add_value(fileLog, __FILE__)                            \
        << boost::log::add_value(lineLog, __LINE__)

#endif /* LOGGER_INCLUDE */

max_stars_repo_name: exec-helper/source

max_issues_repo_name: exec-helper/source

max_forks_repo_name: exec-helper/source