text:
/////////////////////////////////////////////////////////////////////
// = NMatrix
//
// A linear algebra library for scientific computation in Ruby.
// NMatrix is part of SciRuby.
//
// NMatrix was originally inspired by and derived from NArray, by
// Masahiro Tanaka: http://narray.rubyforge.org
//
// == Copyright Information
//
// SciRuby is Copyright (c) 2010 - 2014, Ruby Science Foundation
// NMatrix is Copyright (c) 2012 - 2014, John Woods and the Ruby Science Foundation
//
// Please see LICENSE.txt for additional copyright notices.
//
// == Contributing
//
// By contributing source code to SciRuby, you agree to be bound by
// our Contributor Agreement:
//
// * https://github.com/SciRuby/sciruby/wiki/Contributor-Agreement
//
// == getrs.h
//
// getrs function in native C++.
//

/*
 *             Automatically Tuned Linear Algebra Software v3.8.4
 *                    (C) Copyright 1999 R. Clint Whaley
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *   1. Redistributions of source code must retain the above copyright
 *      notice, this list of conditions and the following disclaimer.
 *   2. Redistributions in binary form must reproduce the above copyright
 *      notice, this list of conditions, and the following disclaimer in the
 *      documentation and/or other materials provided with the distribution.
 *   3. The name of the ATLAS group or the names of its contributers may
 *      not be used to endorse or promote products derived from this
 *      software without specific written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE ATLAS GROUP OR ITS CONTRIBUTORS
 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */

#ifndef POTRS_H
#define POTRS_H

extern "C" {
#if defined HAVE_CBLAS_H
  #include <cblas.h>
#elif defined HAVE_ATLAS_CBLAS_H
  #include <atlas/cblas.h>
#endif
}

namespace nm { namespace math {

/*
 * Solves a system of linear equations A*X = B with a symmetric positive definite matrix A using the Cholesky factorization computed by POTRF.
 *
 * From ATLAS 3.8.0.
 */
template <typename DType, bool is_complex>
int potrs(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const int N, const int NRHS, const DType* A,
           const int lda, DType* B, const int ldb)
{
  // enum CBLAS_DIAG Lunit, Uunit; // These aren't used. Not sure why they're declared in ATLAS' src.

  CBLAS_TRANSPOSE MyTrans = is_complex ? CblasConjTrans : CblasTrans;

  if (!N || !NRHS) return 0;

  const DType ONE = 1;

  if (Order == CblasColMajor) {
    if (Uplo == CblasUpper) {
      nm::math::trsm<DType>(Order, CblasLeft, CblasUpper, MyTrans, CblasNonUnit, N, NRHS, ONE, A, lda, B, ldb);
      nm::math::trsm<DType>(Order, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, N, NRHS, ONE, A, lda, B, ldb);
    } else {
      nm::math::trsm<DType>(Order, CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, N, NRHS, ONE, A, lda, B, ldb);
      nm::math::trsm<DType>(Order, CblasLeft, CblasLower, MyTrans, CblasNonUnit, N, NRHS, ONE, A, lda, B, ldb);
    }
  } else {
    // There's some kind of scaling operation that normally happens here in ATLAS. Not sure what it does, so we'll only
    // worry if something breaks. It probably has to do with their non-templated code and doesn't apply to us.

    if (Uplo == CblasUpper) {
      nm::math::trsm<DType>(Order, CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, NRHS, N, ONE, A, lda, B, ldb);
      nm::math::trsm<DType>(Order, CblasRight, CblasUpper, MyTrans, CblasNonUnit, NRHS, N, ONE, A, lda, B, ldb);
    } else {
      nm::math::trsm<DType>(Order, CblasRight, CblasLower, MyTrans, CblasNonUnit, NRHS, N, ONE, A, lda, B, ldb);
      nm::math::trsm<DType>(Order, CblasRight, CblasLower, CblasNoTrans, CblasNonUnit, NRHS, N, ONE, A, lda, B, ldb);
    }
  }
  return 0;
}

/*
* Function signature conversion for calling LAPACK's potrs functions as directly as possible.
*
* For documentation: http://www.netlib.org/lapack/double/dpotrs.f
*
* This function should normally go in math.cpp, but we need it to be available to nmatrix.cpp.
*/
template <typename DType, bool is_complex>
inline int clapack_potrs(const enum CBLAS_ORDER order, const enum CBLAS_UPLO uplo, const int n, const int nrhs,
                         const void* a, const int lda, void* b, const int ldb) {
  return potrs<DType,is_complex>(order, uplo, n, nrhs, reinterpret_cast<const DType*>(a), lda, reinterpret_cast<DType*>(b), ldb);
}

} } // end nm::math

#endif // POTRS_H

max_stars_repo_name: blackwinter-attic/nmatrix

max_issues_repo_name: blackwinter-attic/nmatrix

max_forks_repo_name: blackwinter-attic/nmatrix


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 <petscversion.h>
#if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR >=3) )
  #if (PETSC_VERSION_MINOR >=6)
     #include <petsc/private/kspimpl.h>
  #else
     #include <petsc-private/kspimpl.h>
  #endif
#else
  #include <private/kspimpl.h>
#endif

#include "common-driver-utils.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 "petsccompat.h"
#include "BSSCR.h"
#include "stokes_block_scaling.h"
#include "stokes_mvblock_scaling.h"
#include "mg.h"
#include "ksp_pressure_nullspace.h"

typedef struct {
    int num_nsp_vecs;
    Vec *nsp_vecs;
} nsp_data_t;

typedef struct {
    KSP kspA11;
    SLE_Solver *sleSolver;
} WendyContext;

//PetscErrorCode BSSCR_NSPRemoveAll(Vec v, void *_data);

PetscErrorCode BSSCR_KSPNormInfMonitor( KSP ksp, PetscInt iteration, PetscReal residualNorm, void *dummy);
PetscErrorCode BSSCR_KSPNormInfToNorm2Monitor( KSP ksp, PetscInt iteration, PetscReal residualNorm, void *dummy);

PetscErrorCode BSSCR_DRIVER_flex( KSP ksp, Mat stokes_A, Vec stokes_x, Vec stokes_b, Mat approxS, KSP ksp_K,
                                  MatStokesBlockScaling BA, PetscTruth sym, KSP_BSSCR * bsscrp_self )
{
    char name[PETSC_MAX_PATH_LEN];
    char ubefore[100];
    char uafter[100];
    char pbefore[100];
    char pafter[100];
    PetscTruth flg, flg2, truth, useAcceleratingSmoothingMG, useFancySmoothingMG;
    PetscTruth usePreviousGuess, useNormInfStoppingConditions, useNormInfMonitor, found, extractMats;
    Mat K,G,D,C;
    Vec u,p,f,h;
    Mat S;
    Vec h_hat,t,t2,q,v;

    KSP ksp_inner;
    KSP ksp_S;
    KSP ksp_cleaner;
    KSPType ksp_inner_type;

    PetscTruth has_cnst_nullspace;
    PC pc_S, pc_MG, pcInner;
    PetscInt monitor_index,max_it,min_it;
    Vec nsp_vec = PETSC_NULL;

    PetscReal scr_rtol;
    PetscReal inner_rtol;
    PetscReal vSolver_rtol;

    PetscScalar uNormInf, pNormInf;
    PetscScalar uNorm, pNorm, rNorm, fNorm;
    PetscInt  uSize, pSize;
    PetscInt    lmin,lmax;
    PetscInt  iterations, rhs_iterations;
    PetscReal   min,max;
    PetscReal p_sum;

    MGContext mgCtx;
    PC shellPC;

    double t0, t1;
    double mgSetupTime, rhsSolveTime, problemBuildTime, scrSolveTime, a11SingleSolveTime, solutionAnalysisTime;
    Index nx,ny,nz;
    PetscInt j,start,end;

    static int been_here = 0;  /* Ha Ha Ha !! */

    /* get sub matrix / vector objects */
    MatNestGetSubMat( stokes_A, 0,0, &K );
    MatNestGetSubMat( stokes_A, 0,1, &G );
    MatNestGetSubMat( stokes_A, 1,0, &D );
    MatNestGetSubMat( stokes_A, 1,1, &C );

    VecNestGetSubVec( stokes_x, 0, &u );
    VecNestGetSubVec( stokes_x, 1, &p );

    VecNestGetSubVec( stokes_b, 0, &f );
    VecNestGetSubVec( stokes_b, 1, &h );

    /* PetscPrintf( PETSC_COMM_WORLD,  "\t Adress of stokes_x is %p\n", stokes_x); */
    /* VecNorm( u, NORM_2, &uNorm ); */
    /* PetscPrintf( PETSC_COMM_WORLD,  "\t u Norm is %.6e in %s: address is %p\n",uNorm,__func__,u); */
    /* VecNorm( p, NORM_2, &pNorm ); */
    /* PetscPrintf( PETSC_COMM_WORLD,  "\t p Norm is %.6e in %s: addres is %p\n",pNorm,__func__,p); */
    /* Create Schur complement matrix */

    problemBuildTime = MPI_Wtime();

    //MatCreateSchurFromBlock( stokes_A, 0.0, "MatSchur_A11", &S );
    MatCreateSchurComplement(K,K,G,D,C, &S);
    /* configure inner solver */
    if (ksp_K!=PETSC_NULL) {
        MatSchurComplementSetKSP( S, ksp_K );
        MatSchurComplementGetKSP( S, &ksp_inner );
    }
    else {
        abort();
        MatSchurComplementGetKSP( S, &ksp_inner );
        KSPSetType( ksp_inner, "cg" );
    }
    KSPGetPC( ksp_inner, &pcInner );

    /* If we're using multigrid, replace the preconditioner here
       so we get the same options prefix. */

    if(bsscrp_self->mg) {
        mgSetupTime=setupMG( bsscrp_self, ksp_inner, pcInner, K, &mgCtx );
    }

    /* SETFROMOPTIONS MIGHT FUCK MG UP */
    KSPSetOptionsPrefix( ksp_inner, "A11_" );
    KSPSetFromOptions( ksp_inner );

    useNormInfStoppingConditions = PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL ,"-A11_use_norm_inf_stopping_condition", &useNormInfStoppingConditions, &found );
    if(useNormInfStoppingConditions)
        BSSCR_KSPSetNormInfConvergenceTest( ksp_inner );

    useNormInfMonitor = PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL, "-A11_ksp_norm_inf_monitor", &useNormInfMonitor, &found );
    if(useNormInfMonitor)
        KSPMonitorSet( ksp_inner, BSSCR_KSPNormInfMonitor, PETSC_NULL, PETSC_NULL );

    useNormInfMonitor = PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL, "-A11_ksp_norm_inf_to_norm_2_monitor", &useNormInfMonitor, &found );
    if(useNormInfMonitor)
        KSPMonitorSet( ksp_inner, BSSCR_KSPNormInfToNorm2Monitor, PETSC_NULL, PETSC_NULL );

    /* create right hand side */
    /* h_hat = G'*inv(K)*f - h */

    MatGetVecs(K,PETSC_NULL,&t);
    MatGetVecs( S, PETSC_NULL, &h_hat );

    KSPSetOptionsPrefix( ksp_inner, "A11_" );
    KSPSetFromOptions( ksp_inner );

    problemBuildTime = MPI_Wtime() - problemBuildTime;

    rhsSolveTime = MPI_Wtime();
    KSPSolve(ksp_inner,f,t);/* t=f/K */
    rhsSolveTime = MPI_Wtime() - rhsSolveTime;
    KSPGetIterationNumber( ksp_inner, &rhs_iterations);

    //bsscr_writeVec( t, "ts", "Writing t vector");
    MatMult(D,t,h_hat);/* G'*t */
    VecAXPY(h_hat, -1.0, h);/* h_hat = h_hat - h */
    Stg_VecDestroy(&t);
    //bsscr_writeVec( h_hat, "h_hat", "Writing h_hat Vector in Solver");
    //MatSchurApplyReductionToVecFromBlock( S, stokes_b, h_hat );

    /* create solver for S p = h_hat */

    KSPCreate( PETSC_COMM_WORLD, &ksp_S );
    KSPSetOptionsPrefix( ksp_S, "scr_");
    Stg_KSPSetOperators( ksp_S, S,S, SAME_NONZERO_PATTERN );
    KSPSetType( ksp_S, "cg" );

    /* Build preconditioner for S */
    KSPGetPC( ksp_S, &pc_S );
    BSSCR_BSSCR_StokesCreatePCSchur2( K,G,D,C,approxS,pc_S,sym, bsscrp_self );

    KSPSetFromOptions(ksp_S);

    /* Set specific monitor test */
    KSPGetTolerances( ksp_S, PETSC_NULL, PETSC_NULL, PETSC_NULL, &max_it );
    //BSSCR_KSPLogSetMonitor( ksp_S, max_it, &monitor_index );

    /* Pressure / Velocity Solve */
    scrSolveTime = MPI_Wtime();
    // PetscPrintf( PETSC_COMM_WORLD, "\t* Pressure / Velocity Solve \n");

    usePreviousGuess = PETSC_FALSE;
    if(been_here)
        PetscOptionsGetTruth( PETSC_NULL, "-scr_use_previous_guess", &usePreviousGuess, &found );

    if(usePreviousGuess) {   /* Note this should actually look at checkpoint information */
        KSPSetInitialGuessNonzero( ksp_S, PETSC_TRUE );
    }
    else {
        KSPSetInitialGuessNonzero( ksp_S, PETSC_FALSE );
    }

    //KSPSetRelativeRhsConvergenceTest( ksp_S );

    useNormInfStoppingConditions = PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL ,"-scr_use_norm_inf_stopping_condition", &useNormInfStoppingConditions, &found );
    if(useNormInfStoppingConditions)
        BSSCR_KSPSetNormInfConvergenceTest(ksp_S);

    useNormInfMonitor = PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL, "-scr_ksp_norm_inf_monitor", &useNormInfMonitor, &found );
    if(useNormInfMonitor)
        KSPMonitorSet( ksp_S, BSSCR_KSPNormInfToNorm2Monitor, PETSC_NULL, PETSC_NULL );

    // PetscPrintf( PETSC_COMM_WORLD, "\t* KSPSolve( ksp_S, h_hat, p )\n");
    /* if h_hat needs to be fixed up ..take out any nullspace vectors here */
    /* we want to check that there is no "noise" in the null-space in the h vector */
    /* this causes problems when we are trying to solve a Jacobian system when the Residual is almost converged */

    if(bsscrp_self->check_pressureNS){
        bsscrp_self->buildPNS(ksp);/* build and set nullspace vectors on bsscr - which is on ksp (function pointer is set in KSPSetUp_BSSCR */
    }

    PetscScalar norm, a, a1, a2, hnorm, pnorm, gnorm;
    MatNorm(G,NORM_INFINITY,&gnorm);
    VecNorm(h_hat, NORM_2, &hnorm);
    hnorm=hnorm/gnorm;

    if((hnorm < 1e-6) && (hnorm > 1e-20)){
        VecScale(h_hat,1.0/hnorm);
    }
    /* test to see if v or t are in nullspace of G and orthogonalize wrt h_hat if needed */
    KSPRemovePressureNullspace_BSSCR(ksp, h_hat);
    /***************************************/
    /* set convergence test to use min_it */
    found = PETSC_FALSE;
    min_it = 0;
    PetscOptionsGetInt( PETSC_NULL,"-scr_ksp_set_min_it_converge", &min_it, &found);
    if(found && min_it > 0){
        BSSCR_KSPSetConvergenceMinIts(ksp_S, min_it, bsscrp_self);
    }
    KSPSolve( ksp_S, h_hat, p );
    sprintf(pafter,"psafter_%d",been_here);
    // bsscr_writeVec( p, pafter, "Writing p Vector in Solver");
    /***************************************/
    if((hnorm < 1e-6) && (hnorm > 1e-20)){
        VecScale(h_hat,hnorm);
        VecScale(p,hnorm);
    }
    KSPRemovePressureNullspace_BSSCR(ksp, p);

    scrSolveTime =  MPI_Wtime() - scrSolveTime;
    // PetscPrintf( PETSC_COMM_WORLD, "  SCR Solve  Finished in time: %lf seconds\n\n", scrSolveTime);
    /* Resolve with this pressure to obtain solution for u */

    /* obtain solution for u */
    VecDuplicate( u, &t );
    MatMult( G, p, t);
    VecAYPX( t, -1.0, f ); /* t <- -t + f */
    MatSchurComplementGetKSP( S, &ksp_inner );

    a11SingleSolveTime = MPI_Wtime();         /* ----------------------------------  Final V Solve */

    if(usePreviousGuess)
        KSPSetInitialGuessNonzero( ksp_inner, PETSC_TRUE );

    KSPSetOptionsPrefix( ksp_inner, "backsolveA11_" );
    KSPSetFromOptions( ksp_inner );
    KSPSolve( ksp_inner, t, u );       /* Solve, then restore default tolerance and initial guess */

    a11SingleSolveTime = MPI_Wtime() - a11SingleSolveTime;            /* ------------------ Final V Solve */

    PetscPrintf( PETSC_COMM_WORLD,  "SCR Solver Summary:\n\n");

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

    PetscPrintf( PETSC_COMM_WORLD,     "  RHS / Pre Solve:        = %.4g secs / %d its\n", rhsSolveTime, rhs_iterations);
    KSPGetIterationNumber( ksp_S, &iterations);
    PetscPrintf( PETSC_COMM_WORLD,     "  Pressure Solve:         = %.4g secs / %d its\n", scrSolveTime, iterations);
    KSPGetIterationNumber( ksp_inner, &iterations);
    PetscPrintf( PETSC_COMM_WORLD,     "  Final V Solve:          = %.4g secs / %d its\n\n", a11SingleSolveTime, iterations);

    /* Analysis of solution:
       This can be somewhat time consuming as it requires allocation / de-allocation,
       computing vector norms etc. So we make it optional..
       This should be put into a proper KSP  monitor now?
    */
    flg = PETSC_TRUE;
    PetscOptionsGetTruth( PETSC_NULL, "-scr_ksp_solution_summary", &flg, &found );

    if(flg) {

        solutionAnalysisTime = MPI_Wtime();

        VecGetSize( u, &uSize );
        VecGetSize( p, &pSize );

        VecDuplicate( u, &t2 );
        MatMult( K, u, t2);
        VecAYPX( t2, -1.0, t ); /* t2 <- -t2 + t  ... should be the formal residual vector */

        VecNorm( t2, NORM_2, &rNorm );
        VecNorm( f,  NORM_2, &fNorm );

        PetscPrintf( PETSC_COMM_WORLD,  "  |f - K u - G p|/|f|     = %.6e\n", rNorm/fNorm );

        VecDuplicate( p, &q );
        MatMult( D, u, q );
        VecNorm( u, NORM_2, &uNorm );
        VecNorm( q, NORM_2, &rNorm );

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

        VecNorm( q, NORM_INFINITY, &rNorm );

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

        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);

        VecMax( p, &lmax, &max );
        VecMin( p, &lmin, &min );
        PetscPrintf( PETSC_COMM_WORLD,  "  min/max(p)    = %.6e [%d] / %.6e [%d]\n",min,lmin,max,lmax);

        VecSum( p, &p_sum );
        PetscPrintf( PETSC_COMM_WORLD,  "  \\sum_i p_i    = %.6e \n", 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(&q );

    }

    if(bsscrp_self->mg) {
        //MG_inner_solver_pcmg_shutdown( pcInner );
    }

    Stg_VecDestroy(&t );

//      KSPLogDestroyMonitor( ksp_S );

    Stg_KSPDestroy(&ksp_S );
    //Stg_KSPDestroy(&ksp_inner );
    Stg_VecDestroy(&h_hat );
    Stg_MatDestroy(&S );

    /* Destroy nullspace vector if it exists. */
    if(nsp_vec)
        Stg_VecDestroy(&nsp_vec);

    //been_here = 1;
    been_here++;
    PetscFunctionReturn(0);
}

/*
  A pointwise stopping condition (infinity norm)

  To activate this stopping condition, I add this piece of code somewhere before the call to KSPSolve().

  {
  PetscTruth pw,flg;
  pw = PETSC_FALSE;
  PetscOptionsGetTruth(0,"-pointwise_stopping_condition", &pw,&flg);
  if(pw==PETSC_TRUE) { KSPSetPointwiseConvergenceTest(ksp_T); }
  }

  You can activate a monitor for the max pointwise residual by using the cmd line option
  -XXX_ksp_pointwise_monitor
  where XXX is the solver prefix specified by the call KSPSetOptionsPrefix().

  DAM
*/

typedef struct {
    PetscTruth initialrtol;    /* default relative residual decrease is computing from initial residual, not rhs */
    PetscTruth mininitialrtol; /* default relative residual decrease is computing from min of initial residual and rhs */
    Vec        work;
    PetscReal  pointwise_max;
} KSPPWConvergedCtx;

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPPWConvergedCreate"
PetscErrorCode BSSCR_KSPPWConvergedCreate(void **ctx)
{
    PetscErrorCode            ierr;
    KSPPWConvergedCtx         *cctx;

    PetscFunctionBegin;
    ierr = Stg_PetscNew(KSPPWConvergedCtx,&cctx);CHKERRQ(ierr);
    *ctx = cctx;
    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPPWConvergedDestroy"
PetscErrorCode BSSCR_KSPPWConvergedDestroy(void *ctx)
{
    PetscErrorCode         ierr;
    KSPPWConvergedCtx *cctx = (KSPPWConvergedCtx*) ctx;

    PetscFunctionBegin;
    if (cctx->work) {ierr = Stg_VecDestroy(&cctx->work);CHKERRQ(ierr);}
    ierr = PetscFree(cctx);CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPNormInfConverged"
PetscErrorCode BSSCR_KSPNormInfConverged(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
{
    PetscErrorCode         ierr;
    KSPPWConvergedCtx     *cctx = (KSPPWConvergedCtx*)ctx;
    KSPNormType            normtype;
    PetscReal              min, max, R_max, R_min, R_Ninf;
    Vec                    R, work, w1,w2;

    PetscFunctionBegin;
    PetscValidHeaderSpecific(ksp,KSP_COOKIE,1);
    PetscValidPointer(reason,4);
    *reason = KSP_CONVERGED_ITERATING;

    ierr = VecDuplicate(ksp->vec_rhs,&work);CHKERRQ(ierr);
    ierr = VecDuplicate(ksp->vec_rhs,&w1);CHKERRQ(ierr);
    ierr = VecDuplicate(ksp->vec_rhs,&w2);CHKERRQ(ierr);

    KSPBuildResidual( ksp, w1,w2, &R );
    VecNorm( R, NORM_INFINITY, &R_Ninf );

    //PetscPrintf( PETSC_COMM_WORLD, "Norm inf convergence %s\n ", ksp->prefix);

    cctx->pointwise_max = R_Ninf;

    ierr = KSPGetNormType(ksp,&normtype);

    CHKERRQ(ierr);
    if (normtype == KSP_NORM_NO)
        Stg_SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Use BSSCR_KSPSkipConverged() with KSPNormType of KSP_NORM_NO");

    if (!cctx)
        Stg_SETERRQ(PETSC_ERR_ARG_NULL,"Convergence context must have been created with BSSCR_KSPDefaultConvergedCreate()");

    if (!n) {
        /* if user gives initial guess need to compute norm of b */
        if (!ksp->guess_zero && !cctx->initialrtol) {
            PetscReal      snorm;
            if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->pc_side == PC_RIGHT) {
                ierr = PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of RHS\n");
                CHKERRQ(ierr);
                ierr = VecNorm(ksp->vec_rhs,NORM_INFINITY,&snorm);CHKERRQ(ierr);        /*     <- b'*b */

                PetscPrintf( PETSC_COMM_WORLD,  "Non Zero Guess; RHS - %g\n", snorm);

            }
            else {
                Vec z;
                if (!cctx->work) {
                    ierr = VecDuplicate(ksp->vec_rhs,&cctx->work);CHKERRQ(ierr);
                }
                z = cctx->work;
                ierr = KSP_PCApply(ksp,ksp->vec_rhs,z);CHKERRQ(ierr);
                if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
                    ierr = PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of preconditioned RHS\n");CHKERRQ(ierr);
                    ierr = VecNorm(z,NORM_INFINITY,&snorm);CHKERRQ(ierr);                 /*    dp <- b'*B'*B*b */

                }
                else if (ksp->normtype == KSP_NORM_NATURAL) {
                    PetscScalar norm;
                    Vec bz;
                    ierr = PetscInfo(ksp,"user has provided nonzero initial guess, computing natural norm of RHS\n");CHKERRQ(ierr);
                    //        ierr  = VecDot(ksp->vec_rhs,z,&norm);
                    //        snorm = sqrt(PetscAbsScalar(norm));                            /*    dp <- b'*B*b */
                    VecDuplicate( z, &bz );
                    VecPointwiseMult( bz, ksp->vec_rhs, z );
                    ierr = VecNorm(bz,NORM_INFINITY,&snorm);CHKERRQ(ierr);
                    Stg_VecDestroy(&bz);
                }
            }
            /* handle special case of zero RHS and nonzero guess */
            if (!snorm) {
                ierr = PetscInfo(ksp,"Special case, user has provided nonzero initial guess and zero RHS\n");CHKERRQ(ierr);
                snorm = rnorm;
            }
            if (cctx->mininitialrtol) {
                ksp->rnorm0 = PetscMin(snorm,rnorm);
            } else {
                ksp->rnorm0 = snorm;
            }
        } else {
            ksp->rnorm0 = rnorm;
        }
        ksp->ttol   = PetscMax(ksp->rtol*ksp->rnorm0,ksp->abstol);
    }

    // if (n <= ksp->chknorm) PetscFunctionReturn(0);

    if ( R_Ninf != R_Ninf ) {
        ierr = PetscInfo(ksp,"Linear solver has created a not a number (NaN) as the pointwise residual norm, declaring divergence \n");CHKERRQ(ierr);
        *reason = KSP_DIVERGED_NAN;
    }
    else
        if (R_Ninf <= ksp->ttol) {
            if (R_Ninf < ksp->abstol) {
                ierr = PetscInfo3(ksp,"Linear solver has converged. Pointwise residual %G is less than absolute tolerance %G at iteration %D\n",R_Ninf,ksp->abstol,n);
                CHKERRQ(ierr);
                *reason = KSP_CONVERGED_ATOL;
            }
            else {
                if (cctx->initialrtol) {
                    ierr = PetscInfo4(ksp,"Linear solver has converged. Norm_infinity %G is less than relative tolerance %G times initial Norm_infinity %G at iteration %D\n",R_Ninf,ksp->rtol,ksp->rnorm0,n);
                    CHKERRQ(ierr);
                }
                else {
                    ierr = PetscInfo4(ksp,"Linear solver has converged. Norm_infinity %G is less than relative tolerance %G times initial norm_infinity right hand side %G at iteration %D\n",R_Ninf,ksp->rtol,ksp->rnorm0,n);CHKERRQ(ierr);
                }
                *reason = KSP_CONVERGED_RTOL;
            }
        }
        else
            if (R_Ninf >= ksp->divtol*ksp->rnorm0) {
                ierr = PetscInfo3(ksp,"Linear solver is diverging. Initial right hand size Norm_infinity value %G, current residual norm %G at iteration %D\n",ksp->rnorm0,R_Ninf,n);CHKERRQ(ierr);
                *reason = KSP_DIVERGED_DTOL;
            }

    /* trash all work vectors here */

    Stg_VecDestroy(&work);
    Stg_VecDestroy(&w1);
    Stg_VecDestroy(&w2);

    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPNormInfToNorm2Monitor"
PetscErrorCode BSSCR_KSPNormInfToNorm2Monitor(KSP ksp,PetscInt n,PetscReal rnorm, void *dummy)
{
    PetscErrorCode          ierr;
    PetscViewerASCIIMonitor viewer  = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm);
    PetscReal              R_normInf, R_norm2;
    PetscInt                     R_size;
    Vec                    R, work, w1, w2;

    PetscFunctionBegin;

    ierr = VecDuplicate(ksp->vec_rhs,&work);CHKERRQ(ierr);
    ierr = VecDuplicate(ksp->vec_rhs,&w1);CHKERRQ(ierr);
    ierr = VecDuplicate(ksp->vec_rhs,&w2);CHKERRQ(ierr);

    KSPBuildResidual( ksp, w1,w2, &R );

    VecNorm( R, NORM_INFINITY, &R_normInf );
    VecNorm( R, NORM_2,        &R_norm2 );
    VecGetSize( R, &R_size );

    ierr = PetscViewerASCIIMonitorCreate(((PetscObject)ksp)->comm,"stdout",0,&viewer); CHKERRQ(ierr);

    if(R_norm2 == 0.0) {
        ierr = PetscViewerASCIIMonitorPrintf(viewer,"%3D KSP Residual Spikiness - INFINITY (%s) \n",
                                             n, ((PetscObject)ksp)->prefix ? ((PetscObject)ksp)->prefix: "" );
        CHKERRQ(ierr);
    }
    else{
        ierr = PetscViewerASCIIMonitorPrintf(viewer,"%3D KSP Residual Spikiness %14.12e (%s) Max / Rms \n",
                                             n,  sqrt((double) R_size) * R_normInf / R_norm2,
                                             ((PetscObject)ksp)->prefix ? ((PetscObject)ksp)->prefix: "");
        CHKERRQ(ierr);
    }

    ierr = PetscViewerASCIIMonitorDestroy(viewer);
    CHKERRQ(ierr);

    Stg_VecDestroy(&work);
    Stg_VecDestroy(&w1);
    Stg_VecDestroy(&w2);

    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPNormInfMonitor"
PetscErrorCode BSSCR_KSPNormInfMonitor(KSP ksp,PetscInt n,PetscReal rnorm, void *dummy)
{
    PetscErrorCode          ierr;
    PetscViewerASCIIMonitor viewer;
    PetscReal               R_normInf;
    Vec                     R, work, w1,w2;

    PetscFunctionBegin;

    ierr = VecDuplicate(ksp->vec_rhs,&work);CHKERRQ(ierr);
    ierr = VecDuplicate(ksp->vec_rhs,&w1);CHKERRQ(ierr);
    ierr = VecDuplicate(ksp->vec_rhs,&w2);CHKERRQ(ierr);

    KSPBuildResidual( ksp, w1,w2, &R );
    VecNorm( R, NORM_INFINITY, &R_normInf );

    ierr = PetscViewerASCIIMonitorCreate(((PetscObject)ksp)->comm,"stdout",0,&viewer);
    CHKERRQ(ierr);

    ierr = PetscViewerASCIIMonitorPrintf(viewer,"%3D KSP Residual Ninf %14.12e (%s) \n",
                                         n, R_normInf,
                                         ((PetscObject)ksp)->prefix ? ((PetscObject)ksp)->prefix: "");
    CHKERRQ(ierr);

    ierr = PetscViewerASCIIMonitorDestroy(viewer);
    CHKERRQ(ierr);

    Stg_VecDestroy(&work);
    Stg_VecDestroy(&w1);
    Stg_VecDestroy(&w2);

    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPNorm2RawMonitor"
PetscErrorCode BSSCR_KSPNorm2RawMonitor(KSP ksp,PetscInt n,PetscReal rnorm, void *dummy)
{
    PetscErrorCode          ierr;
    PetscViewerASCIIMonitor viewer;
    PetscReal               R_norm2;
    Vec                     R, work, w1,w2;

    PetscFunctionBegin;

    ierr = VecDuplicate(ksp->vec_rhs,&work);CHKERRQ(ierr);
    ierr = VecDuplicate(ksp->vec_rhs,&w1);CHKERRQ(ierr);
    ierr = VecDuplicate(ksp->vec_rhs,&w2);CHKERRQ(ierr);

    KSPBuildResidual( ksp, w1,w2, &R );
    VecNorm( R, NORM_2, &R_norm2 );

    ierr = PetscViewerASCIIMonitorCreate(((PetscObject)ksp)->comm,"stdout",0,&viewer);
    CHKERRQ(ierr);

    ierr = PetscViewerASCIIMonitorPrintf(viewer,"%3D KSP Residual Norm2raw %14.12e (%s) \n",
                                         n, R_norm2,
                                         ((PetscObject)ksp)->prefix ? ((PetscObject)ksp)->prefix: "");
    CHKERRQ(ierr);

    ierr = PetscViewerASCIIMonitorDestroy(viewer);
    CHKERRQ(ierr);

    Stg_VecDestroy(&work);
    Stg_VecDestroy(&w1);
    Stg_VecDestroy(&w2);

    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPSetNormInfConvergenceTest"
PetscErrorCode BSSCR_KSPSetNormInfConvergenceTest(KSP ksp)
{
    KSPPWConvergedCtx *ctx;
    PetscTruth monitor, flg;
    char *opt;
    const char *prefix;

    BSSCR_KSPPWConvergedCreate( (void**)&ctx );
    KSPGetOptionsPrefix( ksp, &prefix );

#if(PETSC_VERSION_MAJOR == 3)
    KSPSetConvergenceTest( ksp, BSSCR_KSPNormInfConverged, ctx, BSSCR_KSPPWConvergedDestroy );     // 3.0.0 only
#else
    KSPSetConvergenceTest( ksp, BSSCR_KSPNormInfConverged, ctx );                            // 2.3.3
#endif

    PetscFunctionReturn(0);
}

ext: c


text:
#include <stdio.h>
#include <gsl/gsl_blas.h>
#include <time.h>
#include <stdlib.h>
#include <limits.h>
//gcc program.c -lgsl -lgslcblas -lm -o program

void saveTime(FILE *fp, char *type, int vSize, clock_t start, clock_t end) {
	//save time and size to file
	fprintf(fp, "%s, %d, %f\n", type, vSize, ((double) (end - start)) / CLOCKS_PER_SEC);
}

void ddot(FILE *fp, int vSize, gsl_vector *A, gsl_vector *B) {
	clock_t start, end;
	double result = 0;

	start = clock();
	gsl_blas_ddot(A, B, &result);
	end = clock();

	//save time and size to file
	saveTime(fp, "ddot", vSize, start, end);
}

void dgemv(FILE *fp, int vSize, gsl_vector *A) {
	clock_t start, end;
	gsl_matrix *M = gsl_matrix_alloc(vSize, vSize);
	gsl_vector *B = gsl_vector_alloc(vSize);

	for(int i = 0; i < vSize; i++) {
		for(int j = 0; j < vSize; j++) {
			gsl_matrix_set(M, i, j, (double) rand() * i / (double) (rand() + 1));
		}
		gsl_vector_set(B, i, rand() * i / (rand() + 1));
	}

	start = clock();
	gsl_blas_dgemv(CblasNoTrans, 1.0, M, A, 0.0, B);
	end = clock();
	saveTime(fp, "dgemv", vSize, start, end);

	gsl_matrix_free(M);
	gsl_vector_free(B);

}

void measure(FILE *fp, int size) {
	gsl_vector *x = gsl_vector_alloc(size);
	gsl_vector *y = gsl_vector_alloc(size);
	for(int t = 0; t < 10; t++) {
		for(int i = 0; i < size; i++) {
			gsl_vector_set(x, i, rand() * i / (rand() + 1));
			gsl_vector_set(y, i, rand() * i / (rand() + 1));
		}
		ddot(fp, size, x, y);
		dgemv(fp, size, x);
	}

	gsl_vector_free(x);
	gsl_vector_free(y);
}

int main (void) {
	//uruchomić i zmierzyć czasy działania obydwu funkcji dla różnych rozmiarów wektorów -> 10 pomiarów dla każdego rozmiaru wektora
	FILE *fp;
	srand(time(NULL));
	fp = fopen("wyniki.csv", "w+");
	fprintf(fp, "typ, rozmiar, czas\n"); //make headline

	for(int i = 100; i < 1000; i += 100) {
		int n = i;
			measure(fp, n);
			n += i;
	}

	fclose(fp);
	return 0;
}

max_stars_repo_name: katabana/numerical-methods

max_issues_repo_name: katabana/numerical-methods

max_forks_repo_name: katabana/numerical-methods


text:
#pragma once

#include "Core/Pointers.h"
#include "Graphics/GraphicsResource.h"
#include "Graphics/TextureInfo.h"

#include <glad/gl.h>
#include <gsl/span>

#include <vector>

class Texture;
struct Viewport;

namespace Fb
{
   enum class DepthStencilType
   {
      None,
      Depth24Stencil8,
      Depth32FStencil8
   };

   enum class Target
   {
      Framebuffer = GL_FRAMEBUFFER,
      ReadFramebuffer = GL_READ_FRAMEBUFFER,
      DrawFramebuffer = GL_DRAW_FRAMEBUFFER
   };

   enum class CubeFace
   {
      Front,
      Back,
      Top,
      Bottom,
      Left,
      Right
   };

   struct Specification
   {
      GLsizei width = 0;
      GLsizei height = 0;

      GLsizei samples = 0;
      bool cubeMap = false;

      DepthStencilType depthStencilType = DepthStencilType::Depth24Stencil8;

      gsl::span<const Tex::InternalFormat> colorAttachmentFormats;

      bool operator==(const Specification& other) const;
   };

   struct Attachments
   {
      SPtr<Texture> depthStencilAttachment;
      std::vector<SPtr<Texture>> colorAttachments;
   };

   Attachments generateAttachments(const Specification& specification);
}

namespace std
{
   template<>
   struct hash<Fb::Specification>
   {
      size_t operator()(const Fb::Specification& specification) const;
   };
}

class Framebuffer : public GraphicsResource
{
public:
   Framebuffer();
   Framebuffer(const Framebuffer& other) = delete;
   Framebuffer(Framebuffer&& other);
   ~Framebuffer();
   Framebuffer& operator=(const Framebuffer& other) = delete;
   Framebuffer& operator=(Framebuffer&& other);

private:
   void move(Framebuffer&& other);
   void release();

public:
   using SpecificationType = Fb::Specification;
   static SPtr<Framebuffer> create(const Fb::Specification& specification);
   static const char* labelSuffix();

   static void blit(Framebuffer& source, Framebuffer& destination, GLenum readBuffer, GLenum drawBuffer, GLbitfield mask, GLenum filter);
   static void bindDefault(Fb::Target target = Fb::Target::Framebuffer);

   void bind(Fb::Target target = Fb::Target::Framebuffer);
   bool isBound(Fb::Target target = Fb::Target::Framebuffer) const;

   const Fb::Attachments& getAttachments() const
   {
      return attachments;
   }

   const SPtr<Texture>& getDepthStencilAttachment() const;
   SPtr<Texture> getColorAttachment(int index) const;

   void setAttachments(Fb::Attachments newAttachments);

   bool getViewport(Viewport& viewport) const;

   bool isCubeMap() const;
   void setActiveFace(Fb::CubeFace face);

private:
   const SPtr<Texture>* getFirstValidAttachment() const;

   Fb::Attachments attachments;
};

ext: h

max_stars_repo_name: aaronmjacobs/Swap

max_issues_repo_name: aaronmjacobs/Swap

max_forks_repo_name: aaronmjacobs/Swap


text:
#pragma once
#include <gsl/gsl_math.h>

double brent(gsl_function *F, double x_lo, double x_hi, int *status1, int *status2);


text:
#ifndef __BaseStar_h__
#define __BaseStar_h__

#include <gsl/gsl_roots.h>
#include <gsl/gsl_sf_erf.h>

#include "constants.h"
#include "typedefs.h"
#include "profiling.h"
#include "utils.h"
#include "vector3d.h"

#include "Options.h"
#include "Log.h"
#include "Errors.h"

class BaseStar {

public:

    BaseStar();

    BaseStar(const unsigned long int p_RandomSeed,
             const double            p_MZAMS,
             const double            p_Metallicity,
             const KickParameters    p_KickParameters,
             const double            p_RotationalVelocity = -1.0);

    virtual ~BaseStar() {}

    // object identifiers - all classes have these
    OBJECT_ID           ObjectId() const                                                        { return m_ObjectId; }
    OBJECT_TYPE         ObjectType() const                                                      { return m_ObjectType; }
    STELLAR_TYPE        InitialStellarType() const                                              { return m_InitialStellarType; }
    STELLAR_TYPE        StellarType() const                                                     { return m_StellarType; }
    STELLAR_TYPE        StellarTypePrev() const                                                 { return m_StellarTypePrev; }

    // getters - alphabetically
            double              Age() const                                                     { return m_Age; }
            double              AngularMomentum() const                                         { return CalculateGyrationRadius() * m_Radius * RSOL_TO_AU * m_Radius * RSOL_TO_AU * m_Omega; }
            double              BindingEnergy_Fixed() const                                     { return m_BindingEnergies.fixed; }
            double              BindingEnergy_Nanjing() const                                   { return m_BindingEnergies.nanjing; }
            double              BindingEnergy_Loveridge() const                                 { return m_BindingEnergies.loveridge; }
            double              BindingEnergy_LoveridgeWinds() const                            { return m_BindingEnergies.loveridgeWinds; }
            double              BindingEnergy_Kruckow() const                                   { return m_BindingEnergies.kruckow; }
            bool                CHonMS() const                                                  { return m_CHE; }
            double              COCoreMass() const                                              { return m_COCoreMass; }
            double              CoreMass() const                                                { return m_CoreMass; }
            int                 DominantMassLossRate() const                                    { return static_cast<int>(m_DominantMassLossRate); }
            double              Dt() const                                                      { return m_Dt; }
            double              DtPrev() const                                                  { return m_DtPrev; }
            ERROR               Error() const                                                   { return m_Error; }
            bool                ExperiencedCCSN() const                                         { return (m_SupernovaDetails.events.past & SN_EVENT::CCSN) == SN_EVENT::CCSN; }
            bool                ExperiencedECSN() const                                         { return (m_SupernovaDetails.events.past & SN_EVENT::ECSN) == SN_EVENT::ECSN; }
            bool                ExperiencedPISN() const                                         { return (m_SupernovaDetails.events.past & SN_EVENT::PISN) == SN_EVENT::PISN; }
            bool                ExperiencedPPISN() const                                        { return (m_SupernovaDetails.events.past & SN_EVENT::PPISN) == SN_EVENT::PPISN; }
            SN_EVENT            ExperiencedSN_Type() const                                      { return utils::SNEventType(m_SupernovaDetails.events.past); }
            bool                ExperiencedUSSN() const                                         { return (m_SupernovaDetails.events.past & SN_EVENT::USSN) == SN_EVENT::USSN; }
            double              HeCoreMass() const                                              { return m_HeCoreMass; }
            bool                IsCCSN() const                                                  { return (m_SupernovaDetails.events.current & SN_EVENT::CCSN) == SN_EVENT::CCSN; }
    virtual bool                IsDegenerate() const                                            { return false; }   // default is not degenerate - White Dwarfs, NS and BH are degenerate
            bool                IsECSN() const                                                  { return (m_SupernovaDetails.events.current & SN_EVENT::ECSN) == SN_EVENT::ECSN; }
    virtual bool                IsMassRatioUnstable(const double p_AccretorMass, const bool p_AccretorIsDegenerate) const { return false; } // default is stable
            bool                IsOneOf(const STELLAR_TYPE_LIST p_List) const;
            bool                IsPISN() const                                                  { return (m_SupernovaDetails.events.current & SN_EVENT::PISN) == SN_EVENT::PISN; }
            bool                IsPPISN() const                                                 { return (m_SupernovaDetails.events.current & SN_EVENT::PPISN) == SN_EVENT::PPISN; }
            bool                IsUSSN() const                                                  { return (m_SupernovaDetails.events.current & SN_EVENT::USSN) == SN_EVENT::USSN; }
            double              Lambda_Dewi() const                                             { return m_Lambdas.dewi; }
            double              Lambda_Fixed() const                                            { return m_Lambdas.fixed; }
            double              Lambda_Kruckow() const                                          { return m_Lambdas.kruckow; }
            double              Lambda_KruckowBottom() const                                    { return m_Lambdas.kruckowBottom; }
            double              Lambda_KruckowMiddle() const                                    { return m_Lambdas.kruckowMiddle; }
            double              Lambda_KruckowTop() const                                       { return m_Lambdas.kruckowTop; }
            double              Lambda_Loveridge() const                                        { return m_Lambdas.loveridge; }
            double              Lambda_LoveridgeWinds() const                                   { return m_Lambdas.loveridgeWinds; }
            double              Lambda_Nanjing() const                                          { return m_Lambdas.nanjing; }
            bool                LBV_Phase_Flag() const                                          { return m_LBVphaseFlag; }
            double              Luminosity() const                                              { return m_Luminosity; }
            double              Mass() const                                                    { return m_Mass; }
            double              Mass0() const                                                   { return m_Mass0; }
            double              MassPrev() const                                                { return m_MassPrev; }
            STYPE_VECTOR        MassTransferDonorHistory() const                                { return m_MassTransferDonorHistory; }
            std::string         MassTransferDonorHistoryString() const;
            double              Mdot() const                                                    { return m_Mdot; }
            double              Metallicity() const                                             { return m_Metallicity; }
            double              MZAMS() const                                                   { return m_MZAMS; }
            double              Omega() const                                                   { return m_Omega; }
            double              OmegaCHE() const                                                { return m_OmegaCHE; }
            double              OmegaBreak() const                                              { return CalculateOmegaBreak(); }
            double              OmegaPrev() const                                               { return m_OmegaPrev; }
            double              OmegaZAMS() const                                               { return m_OmegaZAMS; }
            COMPAS_VARIABLE     PropertyValue(const T_ANY_PROPERTY p_Property) const;
            double              Pulsar_MagneticField() const                                    { return m_PulsarDetails.magneticField; }
            double              Pulsar_SpinPeriod() const                                       { return m_PulsarDetails.spinPeriod; }
            double              Pulsar_SpinFrequency() const                                    { return m_PulsarDetails.spinFrequency; }
            double              Pulsar_SpinDownRate() const                                     { return m_PulsarDetails.spinDownRate; }
            double              Radius() const                                                  { return m_Radius; }
            double              RadiusPrev() const                                              { return m_RadiusPrev; }
            unsigned long int   RandomSeed() const                                              { return m_RandomSeed; }
            double              RZAMS() const                                                   { return m_RZAMS; }
            double              SN_CoreMassAtCOFormation() const                                { return m_SupernovaDetails.coreMassAtCOFormation; }
            double              SN_COCoreMassAtCOFormation() const                              { return m_SupernovaDetails.COCoreMassAtCOFormation; }
            SupernovaDetailsT   SN_Details() const                                              { return m_SupernovaDetails; }
            double              SN_DrawnKickMagnitude() const                                   { return m_SupernovaDetails.drawnKickMagnitude; }
            double              SN_EccentricAnomaly() const                                     { return m_SupernovaDetails.eccentricAnomaly; }
            double              SN_FallbackFraction() const                                     { return m_SupernovaDetails.fallbackFraction; }
            double              SN_HeCoreMassAtCOFormation() const                              { return m_SupernovaDetails.HeCoreMassAtCOFormation; }
            bool                SN_IsHydrogenPoor() const                                       { return m_SupernovaDetails.isHydrogenPoor; }
            double              SN_KickMagnitude() const                                        { return m_SupernovaDetails.kickMagnitude; }
            double              SN_MeanAnomaly() const                                          { return m_SupernovaDetails.meanAnomaly; }
            double              SN_Phi() const                                                  { return m_SupernovaDetails.phi; }
            double              SN_TotalMassAtCOFormation() const                               { return m_SupernovaDetails.totalMassAtCOFormation; }
            double              SN_TrueAnomaly() const                                          { return m_SupernovaDetails.trueAnomaly; }
            double              SN_Theta() const                                                { return m_SupernovaDetails.theta; }
            SN_EVENT            SN_Type() const                                                 { return utils::SNEventType(m_SupernovaDetails.events.current); }
            double              SN_KickMagnitudeRandom() const                                  { return m_SupernovaDetails.kickMagnitudeRandom; }
            double              Speed() const                                                   { return m_ComponentVelocity.Magnitude(); }
            COMPAS_VARIABLE     StellarPropertyValue(const T_ANY_PROPERTY p_Property) const;
            double              Tau() const                                                     { return m_Tau; }
            double              Temperature() const                                             { return m_Temperature; }
            double              Time() const                                                    { return m_Time; }
            double              Timescale(TIMESCALE p_Timescale) const                          { return m_Timescales[static_cast<int>(p_Timescale)]; }
            double              XExponent() const                                               { return m_XExponent; }

    // setters
            void                SetInitialType(STELLAR_TYPE p_InitialType)                      { m_InitialStellarType = p_InitialType; }                                           // JR Could do some sanity checks here
            void                SetOmega(double p_vRot)                                         { if (p_vRot >= 0.0) m_Omega = p_vRot; };                                           // Do nothing if sanity check fails (JR: I don't really like this, but I think unavoidable - at least for now)

            void                SetSNCurrentEvent(SN_EVENT p_SNEvent)                           { m_SupernovaDetails.events.current |= p_SNEvent; }                                 // Set supernova primary event/state for current timestep
            void                SetSNPastEvent(const SN_EVENT p_SNEvent)                        { m_SupernovaDetails.events.past |= p_SNEvent; }                                    // Set supernova primary event/state for any past timestep

            void                UpdateComponentVelocity(const Vector3d p_newVelocity);

            void                UpdateMassTransferDonorHistory();

    // member functions - alphabetically
            void            ApplyMassTransferRejuvenationFactor()                                               { m_Age *= CalculateMassTransferRejuvenationFactor(); }             // Apply age rejuvenation factor

            void            CalculateBindingEnergies(const double p_CoreMass, const double p_EnvMass, const double p_Radius);

            double          CalculateDynamicalTimescale() const                                                 { return CalculateDynamicalTimescale_Static(m_Mass, m_Radius); }         // Use class member variables

            double          CalculateEddyTurnoverTimescale();

    virtual void            CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams) { }                                                                                // Default is NO-OP
    virtual void            CalculateGBParams()                                                                 { CalculateGBParams(m_Mass0, m_GBParams); }                         // Use class member variables

    virtual double          CalculateGyrationRadius() const                                                     { return 0.0; }                                                     // Default is 0.0

            void            CalculateLambdas()                                                                  { CalculateLambdas(m_Mass - m_CoreMass); }                          // Use class member variables
            void            CalculateLambdas(const double p_EnvMass);

    virtual DBL_DBL         CalculateMassAcceptanceRate(const double p_DonorMassRate,
                                                        const double p_AccretorMassRate = 0.0);

            double          CalculateMassLossValues(const bool p_UpdateMDot = false, const bool p_UpdateMDt = false);                                                               // JR: todo: better name?

    virtual double          CalculateMomentOfInertia(const double p_RemnantRadius = 0.0) const                  { return 0.0; }                                                     // Use inheritance hierarchy
    virtual double          CalculateMomentOfInertiaAU(const double p_RemnantRadius = 0.0) const                { return 0.0; }                                                     // Use inheritance hierarchy

            double          CalculateNuclearTimescale() const                                                   { return CalculateNuclearTimescale_Static(m_Mass, m_Luminosity); }  // Use class member variables

            double          CalculateOmegaCHE(const double p_MZAMS, const double p_Metallicity) const;

            double          CalculateRadialChange() const                                                       { return (utils::Compare(m_RadiusPrev,0)<=0)? 0 : std::abs(m_Radius - m_RadiusPrev) / m_RadiusPrev; }                    // Return fractional radial change (if previous radius is negative or zero, return 0 to avoid NaN

            double          CalculateRadialExpansionTimescale() const                                           { return CalculateRadialExpansionTimescale_Static(m_StellarType, m_StellarTypePrev, m_Radius, m_RadiusPrev, m_DtPrev); } // Use class member variables

            void            CalculateSNAnomalies(const double p_Eccentricity);

            double          CalculateSNKickMagnitude(const double p_RemnantMass, const double p_EjectaMass, const STELLAR_TYPE p_StellarType);

            double          CalculateThermalMassAcceptanceRate(const double p_Radius) const;
            double          CalculateThermalMassAcceptanceRate() const                                          { return CalculateThermalMassAcceptanceRate(m_Radius); }

    virtual double          CalculateThermalMassLossRate() const                                                { return m_Mass / CalculateThermalTimescale(); }                    // Use class member variables - and inheritance hierarchy

    virtual double          CalculateThermalTimescale(const double p_Radius) const;                                                                                // Use inheritance hierarchy
    virtual double          CalculateThermalTimescale() const                                                   { return CalculateThermalTimescale(m_Radius); }                     // Use inheritance hierarchy

            double          CalculateTimestep();

    virtual double          CalculateZeta(ZETA_PRESCRIPTION p_ZetaPrescription)                                 { return 0.0; }                                                     // Use inheritance hierarchy

            void            ClearCurrentSNEvent()                                                               { m_SupernovaDetails.events.current = SN_EVENT::NONE; }             // Clear supernova event/state for current timestep

    virtual ENVELOPE        DetermineEnvelopeType() const                                                       { return ENVELOPE::REMNANT; }                                       // Default is REMNANT - but should never be called

            void            IncrementOmega(const double p_OmegaDelta)                                           { m_Omega += p_OmegaDelta; }                                        // Apply delta to current m_Omega

            void            ResolveAccretion(const double p_AccretionMass)                                      { m_Mass = std::max(0.0, m_Mass + p_AccretionMass); }               // Handles donation and accretion - won't let mass go negative

    virtual STELLAR_TYPE    ResolveEnvelopeLoss(bool p_NoCheck = false)                                         { return m_StellarType; }

    virtual void            ResolveMassLoss();

    virtual STELLAR_TYPE    ResolveRemnantAfterEnvelopeLoss()                                                   { return m_StellarType; }

            void            SetStellarTypePrev(const STELLAR_TYPE p_StellarTypePrev)                            { m_StellarTypePrev = p_StellarTypePrev; }

            void            StashSupernovaDetails(const STELLAR_TYPE p_StellarType)                             { LOGGING->StashSSESupernovaDetails(this, p_StellarType); }

    virtual void            UpdateAgeAfterMassLoss() { }                                                                                // Default is NO-OP

            STELLAR_TYPE    UpdateAttributesAndAgeOneTimestep(const double p_DeltaMass,
                                                              const double p_DeltaMass0,
                                                              const double p_DeltaTime,
                                                              const bool   p_ForceRecalculate);

    virtual void            UpdateInitialMass() { }                                                                                // Default is NO-OP

    virtual void            UpdateMagneticFieldAndSpin(const bool   p_CommonEnvelope,
                                                       const bool   p_RecyclesNS,
                                                       const double p_Stepsize,
                                                       const double p_MassGainPerTimeStep,
                                                       const double p_Epsilon) { }                                                                                // Default is NO-OP

    // printing functions
            bool            PrintDetailedOutput(const int p_Id) const                                           { return OPTIONS->DetailedOutput() ? LOGGING->LogSSEDetailedOutput(this, p_Id) : true; } // Write record to SSE Detailed Output log file
            bool            PrintSupernovaDetails() const                                                       { return LOGGING->LogSSESupernovaDetails(this); }                   // Write record to SSE Supernovae log file
            bool            PrintStashedSupernovaDetails()                                                      { return LOGGING->LogStashedSSESupernovaDetails(this); }            // Write record to SSE Supernovae log file
            bool            PrintSwitchLog() const                                                              { return OPTIONS->SwitchLog() ? LOGGING->LogSSESwitchLog(this) : true; } // Write record to SSE Switchlog log file
            bool            PrintSystemParameters() const                                                       { return LOGGING->LogSSESystemParameters(this); }                   // Write record to SSE System Parameters file

protected:

    OBJECT_ID               m_ObjectId;                                 // Instantiated object's unique object id
    OBJECT_TYPE             m_ObjectType;                               // Instantiated object's object type
    STELLAR_TYPE            m_InitialStellarType;                       // Stellar type at birth, defined in Hurley et al. 2000
    STELLAR_TYPE            m_StellarType;                              // Stellar type defined in Hurley et al. 2000

    ERROR                   m_Error;                                    // Records most recent error encountered for this star

    // member variables - alphabetical in groups

    bool                    m_CHE;                                      // CHE flag - true if the star spent entire MS as a CH star; false if evolved CH->MS

    // Stellar variables
    unsigned long int       m_RandomSeed;                               // Seeds the random number generator for this star

    // Zero Age Main Sequence
    double                  m_LZAMS;                                    // ZAMS Luminosity
    double                  m_MZAMS;                                    // ZAMS Mass
    double                  m_OmegaZAMS;                                // ZAMS Angular Frequency
    double                  m_OmegaCHE;                                 // Minimum angular frequency at which CHE will occur (calculated at ZAMS)
    double                  m_RZAMS;                                    // ZAMS Radius
    double                  m_TZAMS;                                    // ZAMS Temperature

    // Effective Zero Age Main Sequence
    double                  m_LZAMS0;                                   // Effective ZAMS Luminosity
    double                  m_RZAMS0;                                   // Effective ZAMS Radius
    double                  m_TZAMS0;                                   // Effective ZAMS Temperature

    // Current timestep variables
    double                  m_Age;                                      // Current effective age (changes with mass loss/gain)(myrs)
    double                  m_COCoreMass;                               // Current CO core mass (Msol)
    double                  m_CoreMass;                                 // Current core mass (Msol)
    double                  m_CoreRadius;                               // Current core radius (Rsol)                   JR: todo: I don't think this is used anywhere...
    double                  m_Dt;                                       // Current timestep (myrs)
    double                  m_HeCoreMass;                               // Current He core mass (Msol)
    bool                    m_LBVphaseFlag;                             // Flag to know if the star satisfied the conditions, at any point in its evolution, to be considered a Luminous Blue Variable (LBV)
    double                  m_Luminosity;                               // Current luminosity (Lsol)
    double                  m_Mass;                                     // Current mass (Msol)
    double                  m_Mass0;                                    // Current effective initial mass (Msol)        JR: todo: fix this one day - it is not always initial mass
    double                  m_MinimumLuminosityOnPhase;                 // JR: Only required for CHeB stars, but only needs to be calculated once per star
    double                  m_Mdot;                                     // Current mass loss rate (Msol per ?)
    MASS_LOSS_TYPE          m_DominantMassLossRate;                     // Current dominant mass loss rate
    double                  m_Mu;                                       // Current small envelope parameter mu
    double                  m_Omega;                                    // Current angular frequency (yr-1)
    double                  m_Radius;                                   // Current radius (Rsol)
    double                  m_Tau;                                      // Relative time
    double                  m_Temperature;                              // Current temperature (Tsol)
    double                  m_Time;                                     // Current physical time the star has been evolved(myrs)

    // Previous timestep variables
    double                  m_DtPrev;                                   // Previous timestep
    double                  m_MassPrev;                                 // Previous mass (Msol)
    double                  m_OmegaPrev;                                // Previous angular frequency (yr-1)
    double                  m_RadiusPrev;                               // Previous radius (Rsol)
    STELLAR_TYPE            m_StellarTypePrev;                          // Stellar type at previous timestep

    // Metallicity variables
    double                  m_LogMetallicityRho;                        // logMetallicityXi + 1.0       - called rho in Hurley et al 2000
    double                  m_LogMetallicitySigma;                      // log10(Metallicity)           - called sigma in Hurley et al 2000
    double                  m_LogMetallicityXi;                         // log10(Metallicity / Zsol)    - called xi in Hurley et al 2000
    double                  m_Metallicity;                              // Metallicity

    // Metallicity dependent constants
    double                  m_Alpha1;                                   // alpha1 in Hurly et al. 2000, just after eq 49
    double                  m_Alpha3;                                   // alpha4 in Hurley et al. 2000, just after eq 56
    double                  m_Alpha4;                                   // alpha4 in Hurley et al. 2000, just after eq 57
    double                  m_XExponent;                                // exponent to which R depends on M - 'x' in Hurley et al. 2000, eq 47

    // constants only calculated once
    double                  m_BaryonicMassOfMaximumNeutronStarMass;      // baryonic mass of MaximumNeutronStarMass

    // JR:
    // I initially implemented the following vectors as unordered_maps.  The code worked
    // quite well, except for one small problem - access times (presumably due to hashing)
    // were prohibitive when accessed hundreds of thousands, and in some cases, millions,
    // of times as we evolve the star.  So I used vectors instead - the code is not as
    // elegant, but performance is better by an order of magnitude

    // Timescales, Giant Branch parameters, mass cutoffs
    DBL_VECTOR              m_GBParams;                                 // Giant Branch Parameters
    DBL_VECTOR              m_MassCutoffs;                              // Mass cutoffs
    DBL_VECTOR              m_Timescales;                               // Timescales

    // Luminosity, Radius, a(n) and b(n) coefficients
    DBL_VECTOR              m_AnCoefficients;                           // a(n) coefficients
    DBL_VECTOR              m_BnCoefficients;                           // b(n) coefficients
    DBL_VECTOR              m_LCoefficients;                            // Luminosity coefficients
    DBL_VECTOR              m_RCoefficients;                            // Radius coefficients

    // Luminosity, Radius and Gamma constants
    // JR:
    // These are calculated in CalculateAnCoefficients()
    // Calculating the a(n) coefficients requires one of the R constants, and the L, R and Gamma
    // constants are calculated using the a(n) coefficients, so it seemed logical to calculate
    // them all in the same function
    DBL_VECTOR              m_GammaConstants;                           // Gamma constants
    DBL_VECTOR              m_LConstants;                               // Luminosity constants
    DBL_VECTOR              m_RConstants;                               // Radius constants

    // Binding energies, Lambdas and Zetas
    BindingEnergiesT        m_BindingEnergies;                          // Binding enery values
    LambdasT                m_Lambdas;                                  // Lambda values

    // Stellar details squirrelled away...
    SupernovaDetailsT       m_SupernovaDetails;                         // Supernova attributes
    PulsarDetailsT          m_PulsarDetails;                            // Pulsar attributes

    // Star vector velocity
    Vector3d                m_ComponentVelocity;                        // Isolated star velocity vector (binary's center-of-mass velocity for bound binary)

    // Star mass transfer history
    STYPE_VECTOR            m_MassTransferDonorHistory;                 // List of MT donor stellar types - mostly relevent for binary stars

    // member functions - alphabetically
            void                AgeOneTimestepPreamble(const double p_DeltaTime);

            double              ApplyBlackHoleKicks(const double p_vK, const double p_FallbackFraction, const double p_BlackHoleMass);

            double              CalculateAlpha1() const;
            double              CalculateAlpha3() const;
            double              CalculateAlpha4() const;

            void                CalculateAnCoefficients(DBL_VECTOR &p_AnCoefficients,
                                                        DBL_VECTOR &p_LConstants,
                                                        DBL_VECTOR &p_RConstants,
                                                        DBL_VECTOR &p_GammaConstants);

    virtual void                CalculateAndSetPulsarParameters() { }                                                                                // NO-OP for most stellar types

            double              CalculateBindingEnergy(const double p_CoreMass, const double p_EnvMass, const double p_Radius, const double p_Lambda) const;

            void                CalculateBnCoefficients(DBL_VECTOR &p_BnCoefficients);

    virtual double              CalculateCOCoreMassAtPhaseEnd() const                                                   { return m_COCoreMass; }                                                    // Default is NO-OP
    virtual double              CalculateCOCoreMassOnPhase() const                                                      { return m_COCoreMass; }                                                    // Default is NO-OP

    virtual double              CalculateCoreMassAtPhaseEnd() const                                                     { return m_CoreMass; }                                                      // Default is NO-OP
    static  double              CalculateCoreMassGivenLuminosity_Static(const double p_Luminosity, const DBL_VECTOR &p_GBParams);
    virtual double              CalculateCoreMassOnPhase() const                                                        { return m_CoreMass; }                                                      // Default is NO-OP

    static  double              CalculateDynamicalTimescale_Static(const double p_Mass, const double p_Radius);

    virtual double              CalculateEddingtonCriticalRate() const                                                  { return 2.08E-3 / 1.7 * m_Radius * MYR_TO_YEAR; }       // Hurley+, 2002, Eq. (67); should never be called...

            double              CalculateGBRadiusXExponent() const;

    virtual double              CalculateHeCoreMassAtPhaseEnd() const                                                   { return m_HeCoreMass; }                                                    // Default is NO-OP
    virtual double              CalculateHeCoreMassOnPhase() const                                                      { return m_HeCoreMass; }                                                    // Default is NO-OP

    static  double              CalculateHeRateConstant_Static()                                                        { return HE_RATE_CONSTANT; }                                                // Only >= CHeB stars need AHe, but no drama if other stars calculate (retrieve it) - it's only a constant (we could just use the constant inline...)
    static  double              CalculateHHeRateConstant_Static()                                                       { return HHE_RATE_CONSTANT; }                                               // Only TPAGB stars need AHHe, but no drama if other stars calculate (retrieve it) - it's only a constant (we could just use the constant inline...)

    static  double              CalculateInitialEnvelopeMass_Static(const double p_Mass);

    virtual double              CalculateLambdaDewi() const                                                             { SHOW_WARN(ERROR::NO_LAMBDA_DEWI, "Default used: 1.0"); return 1.0; }      // Not supported: show error
            double              CalculateLambdaKruckow(const double p_Radius, const double p_Alpha) const;
            double              CalculateLambdaLoveridgeEnergyFormalism(const double p_EnvMass, const double p_IsMassLoss = false) const;
    virtual double              CalculateLambdaNanjing() const                                                          { SHOW_WARN(ERROR::NO_LAMBDA_NANJING, "Default used: 1.0"); return 1.0; }   // Not supported: show error

            void                CalculateLCoefficients(const double p_LogMetallicityXi, DBL_VECTOR &p_LCoefficients);

            double              CalculateLifetimeToBAGB(const double p_tHeI, const double p_tHe) const;
            double              CalculateLifetimeToBGB(const double p_Mass) const;

            double              CalculateLogBindingEnergyLoveridge(bool p_IsMassLoss) const;

            double              CalculateLuminosityAtBAGB(double p_Mass) const;
    virtual double              CalculateLuminosityAtPhaseEnd() const                                                   { return m_Luminosity; }                                                    // Default is NO-OP
            double              CalculateLuminosityAtZAMS(const double p_MZAMS);
            double              CalculateLuminosityGivenCoreMass(const double p_CoreMass) const;
    virtual double              CalculateLuminosityOnPhase() const                                                      { return m_Luminosity; }                                                    // Default is NO-OP

            void                CalculateMassCutoffs(const double p_Metallicity, const double p_LogMetallicityXi, DBL_VECTOR &p_MassCutoffs);

    static  double              CalculateMassLoss_Static(const double p_Mass, const double p_Mdot, const double p_Dt);

            double              CalculateMassLossRate();
    virtual double              CalculateMassLossRateHurley();
            double              CalculateMassLossRateKudritzkiReimers() const;
            double              CalculateMassLossRateLBV(const LBV_PRESCRIPTION p_LBV_prescription);
            double              CalculateMassLossRateLBVHurley(const double p_HD_limit_fac) const;
            double              CalculateMassLossRateLBVBelczynski() const;
            double              CalculateMassLossRateNieuwenhuijzenDeJager() const;
            double              CalculateMassLossRateOB(const double p_Teff);
            double              CalculateMassLossRateVassiliadisWood() const;
    virtual double              CalculateMassLossRateVink();
            double              CalculateMassLossRateWolfRayetZDependent(const double p_Mu) const;
            double              CalculateMassLossRateWolfRayet3() const;                                                                                // JR: Never called - do we need it?
            double              CalculateMassLossRateWolfRayet(const double p_Mu) const;

    virtual double              CalculateMassTransferRejuvenationFactor() const;

            double              CalculateMaximumCoreMass(double p_Mass) const;

    static  double              CalculateNuclearTimescale_Static(const double p_Mass, const double p_Luminosity);

            double              CalculateOmegaBreak() const;

    static  double              CalculateOStarRotationalVelocityAnalyticCDF_Static(const double p_Ve);
    static  double              CalculateOStarRotationalVelocityAnalyticCDFInverse_Static(double p_Ve, void *p_Params);
    static  double              CalculateOStarRotationalVelocity_Static(const double p_Xmin, const double p_Xmax);

            double              CalculatePerturbationB(const double p_Mass) const;
            double              CalculatePerturbationC(double p_Mass) const;
    virtual double              CalculatePerturbationMu() const                                                         { return m_Mu; }                                                            // Default is NO-OP
    virtual double              CalculatePerturbationMuAtPhaseEnd() const                                               { return CalculatePerturbationMuOnPhase(); }                                // Same as on phase
    virtual double              CalculatePerturbationMuOnPhase() const                                                  { return CalculatePerturbationMu(); }
            double              CalculatePerturbationQ(const double p_Radius, const double p_Rc) const;
            double              CalculatePerturbationR(const double p_Mu, const double p_Mass, const double p_Radius, const double p_Rc) const;
            double              CalculatePerturbationS(const double p_Mu, const double p_Mass) const;

    static  double              CalculateRadialExpansionTimescale_Static(const STELLAR_TYPE p_StellarType,
                                                                     const STELLAR_TYPE p_StellarTypePrev,
                                                                     const double       p_Radius,
                                                                     const double       p_RadiusPrev,
                                                                     const double       p_DtPrev);

            virtual double      CalculateRadialExtentConvectiveEnvelope() const                                         { return m_Radius; }                                                        // default for stars with no convective envelope

    virtual double              CalculateRadiusAtPhaseEnd() const                                                       { return m_Radius; }                                                        // Default is NO-OP
            double              CalculateRadiusAtZAMS(const double p_MZAMS) const;
    virtual double              CalculateRadiusOnPhase() const                                                          { return m_Radius; }                                                        // Default is NO-OP
    virtual std::tuple <double, STELLAR_TYPE> CalculateRadiusAndStellarTypeOnPhase() const                              { return std::make_tuple(CalculateRadiusOnPhase(), m_StellarType); }

            void                CalculateRCoefficients(const double p_LogMetallicityXi, DBL_VECTOR &p_RCoefficients);

            double              CalculateRotationalVelocity(double p_MZAMS) const;

    virtual double              CalculateTauOnPhase() const                                                             { return m_Tau; }                                                           // Default is NO-OP
    virtual double              CalculateTauAtPhaseEnd() const                                                          { return m_Tau; }                                                           // Default is NO-OP

    virtual double              CalculateTemperatureAtPhaseEnd() const                                                  { return CalculateTemperatureAtPhaseEnd(m_Luminosity, m_Radius); }
    virtual double              CalculateTemperatureAtPhaseEnd(const double p_Luminosity, const double p_Radius) const  { return CalculateTemperatureOnPhase(p_Luminosity, p_Radius); }             // Same as on phase
            double              CalculateTemperatureKelvinOnPhase(const double p_Luminosity, const double p_Radius) const;
    virtual double              CalculateTemperatureOnPhase() const                                                     { return CalculateTemperatureOnPhase(m_Luminosity, m_Radius); }
    virtual double              CalculateTemperatureOnPhase(const double p_Luminosity, const double p_Radius) const     { return CalculateTemperatureOnPhase_Static(p_Luminosity, p_Radius); }
    static  double              CalculateTemperatureOnPhase_Static(const double p_Luminosity, const double p_Radius);

    virtual void                CalculateTimescales()                                                                   { CalculateTimescales(m_Mass0, m_Timescales); }                             // Use class member variables
    virtual void                CalculateTimescales(const double p_Mass, DBL_VECTOR &p_Timescales) { }                                                                                // Default is NO-OP

            double              CalculateZadiabaticHurley2002(const double p_CoreMass) const;
            double              CalculateZadiabaticSPH(const double p_CoreMass) const;
            double              CalculateZadiabatic(ZETA_PRESCRIPTION p_ZetaPrescription);

            double              CalculateZAMSAngularFrequency(const double p_MZAMS, const double p_RZAMS) const;

    virtual double              ChooseTimestep(const double p_Time) const                                               { return m_Dt; }

            double              DrawKickMagnitudeBrayEldridge(const double p_EjectaMass,
                                                              const double p_RemnantMass,
                                                              const double p_Alpha,
                                                              const double p_Beta) const;

            double              DrawKickMagnitudeDistributionFlat(const double p_MaxVK, const double p_Rand) const;
            double              DrawKickMagnitudeDistributionMaxwell(const double p_Sigma, const double p_Rand) const;

            double              DrawRemnantKickMuller(const double p_COCoreMass) const;

            double              DrawRemnantKickMullerMandel(const double p_COCoreMass,
                                                            const double p_Rand,
                                                            const double p_RemnantMass) const;

            double              DrawSNKickMagnitude(const double p_Sigma,
                                                    const double p_COCoreMass,
                                                    const double p_Rand,
                                                    const double p_EjectaMass,
                                                    const double p_RemnantMass) const;

    virtual void                EvolveOneTimestepPreamble() { };                                                                                // Default is NO-OP

            STELLAR_TYPE        EvolveOnPhase();

    virtual STELLAR_TYPE        EvolveToNextPhase()                                                                     { return m_StellarType; }

    virtual bool                IsEndOfPhase() const                                                                    { return false; }
    virtual bool                IsSupernova() const                                                                     { return false; }

            double              LimitTimestep(const double p_Dt);

    /*
     * Perturb Luminosity and Radius
     *
     * See Hurley et al. 2000, section 6.3
     *
     * The default is no perturbation - this function does nothing and is called
     * only if the stellar class doesn't define its own perturbation function.
     * See the stellar class perturbation functions for perturbation details specific
     * to the stellar class.
     *
     * Perturbation is disabled by default when DEBUG is enabled - except when
     * DEBUG_PERTURB is defined (see below).  The stellar class perturbation
     * functions are defined away if DEBUG is defined - so this generic Star
     * function is called (and does nothing).
     *
     * If DEBUG_PERTURB is defined then perturbation is not disabled while debbuging.
     * To enable perturbation while DEBUG is enabled, define DEBUG_PERTURB.
     */
    virtual void                PerturbLuminosityAndRadius() { }                                                                                // NO-OP
    virtual void                PerturbLuminosityAndRadiusAtPhaseEnd()                                                  { PerturbLuminosityAndRadiusOnPhase(); }                                    // Same as on phase
    virtual void                PerturbLuminosityAndRadiusOnPhase()                                                     { PerturbLuminosityAndRadius(); }

            STELLAR_TYPE        ResolveEndOfPhase();
    virtual void                ResolveHeliumFlash() { }
    virtual STELLAR_TYPE        ResolveSkippedPhase()                                                                   { return EvolveToNextPhase(); }                                             // Default is evolve to next phase
    virtual STELLAR_TYPE        ResolveSupernova()                                                                      { return m_StellarType; }                                                   // Default is NO-OP

    virtual void                SetSNHydrogenContent()                                                                  { m_SupernovaDetails.isHydrogenPoor = false; }                              // Default is false

            bool                ShouldBeMasslessRemnant() const                                                         { return (m_Mass <= 0.0 || m_StellarType==STELLAR_TYPE::MASSLESS_REMNANT); }
    virtual bool                ShouldEvolveOnPhase() const                                                             { return true; }
    virtual bool                ShouldSkipPhase() const                                                                 { return false; }                                                           // Default is false

            void                UpdateAttributesAndAgeOneTimestepPreamble(const double p_DeltaMass, const double p_DeltaMass0, const double p_DeltaTime);

};

#endif // __BaseStar_h__

max_stars_repo_name: LiekeVanSon/COMPAS

max_issues_repo_name: LiekeVanSon/COMPAS

max_forks_repo_name: LiekeVanSon/COMPAS


text:
#include <config.h>
#include <math.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_sort_vector.h>

static int compute_rank(gsl_vector *v);

#define BASE_LONG_DOUBLE
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef  BASE_LONG_DOUBLE

#define BASE_DOUBLE
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef  BASE_DOUBLE

#define BASE_FLOAT
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef  BASE_FLOAT

#define BASE_ULONG
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef  BASE_ULONG

#define BASE_LONG
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef  BASE_LONG

#define BASE_UINT
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef  BASE_UINT

#define BASE_INT
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef  BASE_INT

#define BASE_USHORT
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef  BASE_USHORT

#define BASE_SHORT
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef  BASE_SHORT

#define BASE_UCHAR
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef  BASE_UCHAR

#define BASE_CHAR
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef  BASE_CHAR

/*
compute_rank()
  Compute rank of a sorted vector

Inputs: v - sorted data vector on input; rank vector on output

Notes: ranks are always computed in double precision
*/

static int
compute_rank(gsl_vector *v)
{
  const size_t n = v->size;
  size_t i = 0;

  while (i < n - 1)
    {
      double vi = gsl_vector_get(v, i);

      if (vi == gsl_vector_get(v, i + 1))
        {
          size_t j = i + 2;
          size_t k;
          double rank = 0.0;

          /* we have detected a tie, find number of equal elements */
          while (j < n && vi == gsl_vector_get(v, j))
            ++j;

          /* compute rank */
          for (k = i; k < j; ++k)
            rank += k + 1.0;

          /* divide by number of ties */
          rank /= (double) (j - i);

          for (k = i; k < j; ++k)
            gsl_vector_set(v, k, rank);

          i = j;
        }
      else
        {
          /* no tie - set rank to natural ordered position */
          gsl_vector_set(v, i, i + 1.0);
          ++i;
        }
    }

  if (i == n - 1)
    gsl_vector_set(v, n - 1, (double) n);

  return GSL_SUCCESS;
} /* compute_rank() */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

#include <math.h>
#include <stdlib.h>
#include <core_blas.h>
#include <cblas.h>
#include "dplasma.h"
#include "dplasmatypes.h"
#include "dplasma/lib/dplasmaaux.h"
#include "parsec/private_mempool.h"

#include "dplasma/lib/zhetrf.h"
#include "dplasma/lib/ztrmdm.h"

/*
 * dplasma_zhetrf_New()
 */
parsec_taskpool_t*
dplasma_zhetrf_New( parsec_tiled_matrix_dc_t *A, int *INFO)
{
    int ldwork, lwork, ib;
    parsec_taskpool_t *parsec_zhetrf = NULL;
    parsec_memory_pool_t *pool_0, *pool_1;

    ib = A->mb;

    /* ldwork and lwork are necessary for the macros zhetrf_pool_0_SIZE and zhetrf_pool_1_SIZE */
    ldwork = (A->nb+1)*ib;
    lwork = (A->mb+1)*A->nb + ib*ib;

    pool_0 = (parsec_memory_pool_t*)malloc(sizeof(parsec_memory_pool_t));
    parsec_private_memory_init( pool_0, zhetrf_pool_0_SIZE );

    pool_1 = (parsec_memory_pool_t*)malloc(sizeof(parsec_memory_pool_t));
    parsec_private_memory_init( pool_1, zhetrf_pool_1_SIZE );

    parsec_zhetrf = (parsec_taskpool_t *)parsec_zhetrf_new(PlasmaLower, A, (parsec_data_collection_t *)A, ib, pool_1, pool_0, INFO);

    dplasma_add2arena_tile(((parsec_zhetrf_taskpool_t*)parsec_zhetrf)->arenas[PARSEC_zhetrf_DEFAULT_ARENA],
                           A->mb*A->nb*sizeof(parsec_complex64_t),
                           PARSEC_ARENA_ALIGNMENT_SSE,
                           parsec_datatype_double_complex_t, A->mb);

    return parsec_zhetrf;
}

void
dplasma_zhetrf_Destruct( parsec_taskpool_t *tp )
{
    parsec_zhetrf_taskpool_t *obut = (parsec_zhetrf_taskpool_t *)tp;

    parsec_matrix_del2arena( obut->arenas[PARSEC_zhetrf_DEFAULT_ARENA] );

    parsec_taskpool_free(tp);
}

/*
 * dplasma_ztrmdm_New()
 */
parsec_taskpool_t*
dplasma_ztrmdm_New( parsec_tiled_matrix_dc_t *A)
{
    parsec_taskpool_t *parsec_ztrmdm = NULL;

    parsec_ztrmdm = (parsec_taskpool_t *)parsec_ztrmdm_new(A);

    dplasma_add2arena_tile(((parsec_ztrmdm_taskpool_t*)parsec_ztrmdm)->arenas[PARSEC_ztrmdm_DEFAULT_ARENA],
                           A->mb*A->nb*sizeof(parsec_complex64_t),
                           PARSEC_ARENA_ALIGNMENT_SSE,
                           parsec_datatype_double_complex_t, A->mb);

    return parsec_ztrmdm;
}

void
dplasma_ztrmdm_Destruct( parsec_taskpool_t *tp )
{
    parsec_ztrmdm_taskpool_t *obut = (parsec_ztrmdm_taskpool_t *)tp;

    parsec_matrix_del2arena( obut->arenas[PARSEC_ztrmdm_DEFAULT_ARENA] );

    parsec_taskpool_free(tp);
}

/*
 * Blocking Interface
 */

int dplasma_zhetrf(parsec_context_t *parsec, parsec_tiled_matrix_dc_t *A)
{
    parsec_taskpool_t *parsec_zhetrf/*, *parsec_ztrmdm*/;
    int info = 0, ginfo = 0;

    parsec_zhetrf = dplasma_zhetrf_New(A, &info);
    parsec_context_add_taskpool(parsec, (parsec_taskpool_t *)parsec_zhetrf);
    dplasma_wait_until_completion(parsec);
    dplasma_zhetrf_Destruct(parsec_zhetrf);

    /*
    parsec_ztrmdm = dplasma_ztrmdm_New(A);
    parsec_context_add_taskpool(parsec, (parsec_taskpool_t *)parsec_ztrmdm);
    dplasma_wait_until_completion(parsec);
    dplasma_ztrmdm_Destruct(parsec_ztrmdm);
    */

    /* This covers both cases when we have not compiled with MPI, or we don't need to do the reduce */
    ginfo = info;
#if defined(PARSEC_HAVE_MPI)
    /* If we don't need to reduce, don't do it, this way we don't require MPI to be initialized */
    if( A->super.nodes > 1 )
        MPI_Allreduce( &info, &ginfo, 1, MPI_INT, MPI_MAX, *(MPI_Comm*)dplasma_pcomm);
#endif
    return ginfo;
}

max_stars_repo_name: NLAFET/ABFT

max_issues_repo_name: NLAFET/ABFT

max_forks_repo_name: NLAFET/ABFT


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

#include <config.h>

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

#include <gsl/gsl_errno.h>
#include <gsl/gsl_test.h>

#define CHECK(x) errors[n].number = x ; errors[n].name = #x ; n++ ;
#define MAX_ERRS 64

int verbose = 0 ;

int
main (void)
{
  int i, j, n = 0 ;

  struct {
    int number;
    const char * name;
  } errors[MAX_ERRS] ;

  CHECK(GSL_SUCCESS);
  CHECK(GSL_FAILURE);
  CHECK(GSL_CONTINUE);
  CHECK(GSL_EDOM);
  CHECK(GSL_ERANGE);
  CHECK(GSL_EFAULT);
  CHECK(GSL_EINVAL);
  CHECK(GSL_EFAILED);
  CHECK(GSL_EFACTOR);
  CHECK(GSL_ESANITY);
  CHECK(GSL_ENOMEM);
  CHECK(GSL_EBADFUNC);
  CHECK(GSL_ERUNAWAY);
  CHECK(GSL_EMAXITER);
  CHECK(GSL_EZERODIV);
  CHECK(GSL_EBADTOL);
  CHECK(GSL_ETOL);
  CHECK(GSL_EUNDRFLW);
  CHECK(GSL_EOVRFLW);
  CHECK(GSL_ELOSS);
  CHECK(GSL_EROUND);
  CHECK(GSL_EBADLEN);
  CHECK(GSL_ENOTSQR);
  CHECK(GSL_ESING);
  CHECK(GSL_EDIVERGE);
  CHECK(GSL_EUNSUP);
  CHECK(GSL_EUNIMPL);
  CHECK(GSL_ECACHE);
  CHECK(GSL_ETABLE);
  CHECK(GSL_ENOPROG);
  CHECK(GSL_ENOPROGJ);
  CHECK(GSL_ETOLF);
  CHECK(GSL_ETOLX);
  CHECK(GSL_ETOLG);
  CHECK(GSL_EOF);

  for (i = 0 ; i < n ; i++)
    {
      if (verbose) printf ("%s = %d\n", errors[i].name, errors[i].number) ;
    }

  for (i = 0; i < n; i++)
    {
      int status = 0;
      for (j = 0; j < n; j++)
	{
	  if (j != i)
	      status |= (errors[i].number == errors[j].number);
	}

      gsl_test (status, "%s is distinct from other error values",
		errors[i].name);
    }

  for (i = 0; i < n; i++)
    {
      int status = 0;
      int e1 = errors[i].number ;
      for (j = 0; j < n; j++)
	{
	  if (j != i)
	    {
	      int e2 = errors[j].number;
	      status |= (gsl_strerror(e1) == gsl_strerror(e2)) ;
	    }
	}
      gsl_test (status, "%s has a distinct error message",
		errors[i].name);
    }

  exit (gsl_test_summary ());
}