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 <mpi.h>
#include <StGermain/StGermain.h>
#include <StgDomain/StgDomain.h>
#include <StgFEM/StgFEM.h>
#include <petsc.h>

#include "types.h"
#include "SemiLagrangianIntegrator.h"

#include <assert.h>

/** Textual name of this class */
const Type SemiLagrangianIntegrator_Type = "SemiLagrangianIntegrator";

SemiLagrangianIntegrator* _SemiLagrangianIntegrator_New(  SEMILAGRANGIANINTEGRATOR_DEFARGS  )
{
   SemiLagrangianIntegrator*		self;

   /* Allocate memory */
   assert( _sizeOfSelf >= sizeof(SemiLagrangianIntegrator) );
   /* The following terms are parameters that have been passed into this function but are being set before being passed onto the parent */
   /* This means that any values of these parameters that are passed into this function are not passed onto the parent function
      and so should be set to ZERO in any children of this class. */
   nameAllocationType = NON_GLOBAL;

   self = (SemiLagrangianIntegrator*) _Stg_Component_New(  STG_COMPONENT_PASSARGS  );

   /* General info */
   self->variableList = Stg_ObjectList_New();
   self->varStarList  = Stg_ObjectList_New();
   self->varOldList   = Stg_ObjectList_New();

   return self;
}

void* _SemiLagrangianIntegrator_Copy( void* slIntegrator, void* dest, Bool deep, Name nameExt, PtrMap* ptrMap ) {
   SemiLagrangianIntegrator*		self = (SemiLagrangianIntegrator*)slIntegrator;
   SemiLagrangianIntegrator*		newSemiLagrangianIntegrator;
   PtrMap*			map = ptrMap;
   Bool			ownMap = False;

   if( !map ) {
      map = PtrMap_New( 10 );
      ownMap = True;
   }

   newSemiLagrangianIntegrator = _Stg_Component_Copy( self, dest, deep, nameExt, map );

   if( deep ) {
      if( (newSemiLagrangianIntegrator->velocityField = PtrMap_Find( map, self->velocityField )) == NULL ) {
         newSemiLagrangianIntegrator->velocityField = Stg_Class_Copy( self->velocityField, NULL, deep, nameExt, map );
         PtrMap_Append( map, self->velocityField, newSemiLagrangianIntegrator->velocityField );
      }
   }
   else {
      newSemiLagrangianIntegrator->velocityField = Stg_Class_Copy( self->velocityField, NULL, deep, nameExt, map );
   }

   if( ownMap ) {
      Stg_Class_Delete( map );
   }

   return (void*)newSemiLagrangianIntegrator;
}

void _SemiLagrangianIntegrator_Delete( void* slIntegrator ) {
   SemiLagrangianIntegrator*		self = (SemiLagrangianIntegrator*)slIntegrator;
   Stg_Class_Delete( self->variableList );
   Stg_Class_Delete( self->varStarList );
   Stg_Class_Delete( self->varOldList );
}

void _SemiLagrangianIntegrator_Print( void* slIntegrator, Stream* stream ) {
   SemiLagrangianIntegrator*		self = (SemiLagrangianIntegrator*)slIntegrator;

   _Stg_Component_Print( self, stream );

   Journal_PrintPointer( stream, self->velocityField );
}

void* _SemiLagrangianIntegrator_DefaultNew( Name name ) {
   /* Variables set in this function */
   SizeT                                              _sizeOfSelf = sizeof(SemiLagrangianIntegrator);
   Type                                                      type = SemiLagrangianIntegrator_Type;
   Stg_Class_DeleteFunction*                              _delete = _SemiLagrangianIntegrator_Delete;
   Stg_Class_PrintFunction*                                _print = _SemiLagrangianIntegrator_Print;
   Stg_Class_CopyFunction*                                  _copy = _SemiLagrangianIntegrator_Copy;
   Stg_Component_DefaultConstructorFunction*  _defaultConstructor = _SemiLagrangianIntegrator_DefaultNew;
   Stg_Component_ConstructFunction*                    _construct = _SemiLagrangianIntegrator_AssignFromXML;
   Stg_Component_BuildFunction*                            _build = _SemiLagrangianIntegrator_Build;
   Stg_Component_InitialiseFunction*                  _initialise = _SemiLagrangianIntegrator_Initialise;
   Stg_Component_ExecuteFunction*                        _execute = _SemiLagrangianIntegrator_Execute;
   Stg_Component_DestroyFunction*                        _destroy = _SemiLagrangianIntegrator_Destroy;

   /* Variables that are set to ZERO are variables that will be set either by the current _New function or another parent _New function further up the hierachy */
   AllocationType  nameAllocationType = NON_GLOBAL /* default value NON_GLOBAL */;

   return (void*)_SemiLagrangianIntegrator_New(  SEMILAGRANGIANINTEGRATOR_PASSARGS  );
}

void _SemiLagrangianIntegrator_AssignFromXML( void* slIntegrator, Stg_ComponentFactory* cf, void* data ) {
   SemiLagrangianIntegrator*	self 		= (SemiLagrangianIntegrator*)slIntegrator;
   Dictionary*			dict;
   Dictionary_Entry_Value*		dev;
   unsigned			field_i;
   Name				fieldName;
   FeVariable*	   feVariable;

   Stg_Component_AssignFromXML( self, cf, data, False );

   self->context = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"Context", FiniteElementContext, False, data );
   if( !self->context  )
      self->context = Stg_ComponentFactory_ConstructByName( cf, (Name)"context", FiniteElementContext, True, data  );

   self->velocityField = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"VelocityField", FeVariable, False, NULL  );
   self->advectedField = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"AdvectedField", FeVariable, False, NULL  );

   dict = Dictionary_Entry_Value_AsDictionary( Dictionary_Get( cf->componentDict, (Dictionary_Entry_Key)self->name )  );
   dev  = Dictionary_Get( dict, (Dictionary_Entry_Key)"fields" );
   for( field_i = 0; field_i < Dictionary_Entry_Value_GetCount( dev ); field_i += 3  ) {
      fieldName = Dictionary_Entry_Value_AsString( Dictionary_Entry_Value_GetElement( dev, field_i ) );
      feVariable = Stg_ComponentFactory_ConstructByName( cf, (Name)fieldName, FeVariable, True, data  );
      Stg_ObjectList_Append( self->variableList, feVariable );

      /* the corresponding _* field */
      fieldName = Dictionary_Entry_Value_AsString( Dictionary_Entry_Value_GetElement( dev, field_i + 1 ) );
      feVariable = Stg_ComponentFactory_ConstructByName( cf, (Name)fieldName, FeVariable, True, data  );
      Stg_ObjectList_Append( self->varStarList, feVariable );

      /* the corresponding old field */
      fieldName = Dictionary_Entry_Value_AsString( Dictionary_Entry_Value_GetElement( dev, field_i + 2 ) );
      feVariable = Stg_ComponentFactory_ConstructByName( cf, (Name)fieldName, FeVariable, True, data  );
      Stg_ObjectList_Append( self->varOldList, feVariable );
   }

//   self->sle = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"SLE", Energy_SLE, False, NULL  );

   /* for problems with temporally evolving velocity */
   self->prevVelField = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"PreviousTimeStepVelocityField", FeVariable, False, data );
   if( self->prevVelField  ) {
      EP_AppendClassHook( Context_GetEntryPoint( self->context, AbstractContext_EP_UpdateClass ), SemiLagrangianIntegrator_UpdatePreviousVelocityField, self );
      EP_InsertClassHookAfter( Context_GetEntryPoint( self->context, AbstractContext_EP_UpdateClass ), "SemiLagrangianIntegrator_UpdatePreviousVelocityField", SemiLagrangianIntegrator_InitSolve, self );
   } else {
      EP_AppendClassHook( Context_GetEntryPoint( self->context, AbstractContext_EP_UpdateClass ), SemiLagrangianIntegrator_InitSolve, self );
   }

//   if( self->sle ) {
//      /** also set sle to run where required */
//      EP_InsertClassHookAfter( Context_GetEntryPoint( self->context, AbstractContext_EP_UpdateClass ), "SemiLagrangianIntegrator_InitSolve", SystemLinearEquations_GetRunEPFunction(), self->sle );
//      /** remember to disable the standard run at execute */
//      SystemLinearEquations_SetRunDuringExecutePhase( self->sle, False);
//   }

   self->isConstructed = True;
}

void _SemiLagrangianIntegrator_Build( void* slIntegrator, void* data ) {
   SemiLagrangianIntegrator*	self 		= (SemiLagrangianIntegrator*)slIntegrator;
   FeVariable*			feVariable;
   FeVariable*			feVarOld;
   FeVariable*			feVarStar;
   unsigned			   field_i;

   if(self->velocityField) Stg_Component_Build(self->velocityField, data, False);
   if(self->advectedField) Stg_Component_Build(self->advectedField, data, False);
   if(self->prevVelField)  Stg_Component_Build(self->prevVelField , data, False);

   for( field_i = 0; field_i < self->variableList->count; field_i++ ) {
      feVariable = (FeVariable*) self->variableList->data[field_i];
      feVarOld   = (FeVariable*) self->varOldList->data[field_i];
      feVarStar  = (FeVariable*) self->varStarList->data[field_i];
      if(feVariable) Stg_Component_Build(feVariable, data, False);
      if(feVarOld  ) Stg_Component_Build(feVarOld  , data, False);
      if(feVarStar ) Stg_Component_Build(feVarStar , data, False);
   }
}

void _SemiLagrangianIntegrator_Initialise( void* slIntegrator, void* data ) {
   SemiLagrangianIntegrator*	self 		= (SemiLagrangianIntegrator*)slIntegrator;
   FeVariable*			feVariable;
   FeVariable*			feVarOld;
   FeVariable*			feVarStar;
   unsigned			   field_i;

   if(self->velocityField) Stg_Component_Initialise(self->velocityField, data, False);
   if(self->advectedField) Stg_Component_Initialise(self->advectedField, data, False);
   if(self->prevVelField) {
      Stg_Component_Initialise(self->prevVelField , data, False);
      SemiLagrangianIntegrator_UpdatePreviousVelocityField( slIntegrator, NULL );
   }

   for( field_i = 0; field_i < self->variableList->count; field_i++ ) {
      feVariable = (FeVariable*) self->variableList->data[field_i];
      feVarOld   = (FeVariable*) self->varOldList->data[field_i];
      feVarStar  = (FeVariable*) self->varStarList->data[field_i];
      if(feVariable) Stg_Component_Initialise(feVariable, data, False);
      if(feVarOld  ) Stg_Component_Initialise(feVarOld  , data, False);
      if(feVarStar ) Stg_Component_Initialise(feVarStar , data, False);
   }
}

void _SemiLagrangianIntegrator_Execute( void* slIntegrator, void* data ) {
}

void _SemiLagrangianIntegrator_Destroy( void* slIntegrator, void* data ) {
   SemiLagrangianIntegrator*	self 	= (SemiLagrangianIntegrator*)slIntegrator;
   FeVariable*			feVariable;
   FeVariable*			feVarOld;
   FeVariable*			feVarStar;
   unsigned			   field_i;

   if(self->velocityField) Stg_Component_Destroy(self->velocityField, data, False);
   if(self->advectedField) Stg_Component_Destroy(self->advectedField, data, False);
   if(self->prevVelField)  Stg_Component_Destroy(self->prevVelField , data, False);

   for( field_i = 0; field_i < self->variableList->count; field_i++ ) {
      feVariable = (FeVariable*) self->variableList->data[field_i];
      feVarOld   = (FeVariable*) self->varOldList->data[field_i];
      feVarStar  = (FeVariable*) self->varStarList->data[field_i];
      if(feVariable) Stg_Component_Destroy(feVariable, data, False);
      if(feVarOld  ) Stg_Component_Destroy(feVarOld  , data, False);
      if(feVarStar ) Stg_Component_Destroy(feVarStar , data, False);
   }
}

void SemiLagrangianIntegrator_InitSolve( void* _self, void* _context ) {
   SemiLagrangianIntegrator*	self			= (SemiLagrangianIntegrator*) _self;
   unsigned			   field_i, node_i;
   FeVariable*			feVariable;
   FeVariable*			feVarOld;
   FeVariable*			feVarStar;
   double            phi[3];
   unsigned			   lMeshSize;
   FeMesh*				mesh;

   for( field_i = 0; field_i < self->variableList->count; field_i++ ) {
      feVariable = (FeVariable*) self->variableList->data[field_i];
      feVarStar  = (FeVariable*) self->varStarList->data[field_i];
      feVarOld   = (FeVariable*) self->varOldList->data[field_i];

      /* we're assuming that the solution vector has already been updated onto the FeVariable (in the SLE class) */
      mesh = feVariable->feMesh;
      lMeshSize = Mesh_GetLocalSize( mesh, MT_VERTEX );
      for( node_i = 0; node_i < lMeshSize; node_i++ ) {
         FeVariable_GetValueAtNode( feVariable, node_i, phi );
         FeVariable_SetValueAtNode( feVarOld, node_i, phi );
      }
      FeVariable_SyncShadowValues( feVarOld );

      /* generate the _* field */
      SemiLagrangianIntegrator_Solve( self, feVarOld, feVarStar );
   }
}

#define INV6 0.166666666667

void IntegrateRungeKutta( FeVariable* velocityField, double dt, double* origin, double* position ) {
   unsigned		ndims		     = Mesh_GetDimSize( velocityField->feMesh );
   unsigned		dim_i;
   double			min[3], max[3];
   double			k[4][3];
   double			coordPrime[3];
   unsigned*		periodic	     = ((CartesianGenerator*)velocityField->feMesh->generator)->periodic;

   Mesh_GetGlobalCoordRange( velocityField->feMesh, min, max );

   FieldVariable_InterpolateValueAt( velocityField, origin, k[0] );
   for( dim_i = 0; dim_i < ndims; dim_i++ ) {
      coordPrime[dim_i] = origin[dim_i] - 0.5 * dt * k[0][dim_i];
      PeriodicUpdate( coordPrime, min, max, dim_i, periodic[dim_i] );
   }
   FieldVariable_InterpolateValueAt( velocityField, coordPrime, k[1] );

   for( dim_i = 0; dim_i < ndims; dim_i++ ) {
      coordPrime[dim_i] = origin[dim_i] - 0.5 * dt * k[1][dim_i];
      PeriodicUpdate( coordPrime, min, max, dim_i, periodic[dim_i] );
   }
   FieldVariable_InterpolateValueAt( velocityField, coordPrime, k[2] );

   for( dim_i = 0; dim_i < ndims; dim_i++ ) {
      coordPrime[dim_i] = origin[dim_i] - dt * k[2][dim_i];
      PeriodicUpdate( coordPrime, min, max, dim_i, periodic[dim_i] );
   }
   FieldVariable_InterpolateValueAt( velocityField, coordPrime, k[3] );

   for( dim_i = 0; dim_i < ndims; dim_i++ ) {
      position[dim_i] = origin[dim_i] -
         INV6 * dt * ( k[0][dim_i] + 2.0 * k[1][dim_i] + 2.0 * k[2][dim_i] + k[3][dim_i] );
      PeriodicUpdate( position, min, max, dim_i, periodic[dim_i] );
   }
}

/* 2nd order acurate runge kutta algorithm for interpolating backwards in time through a temporally evolving velocity field
   Durran, D. "Numerical Methods for Wave Equations in Geophysical Fluid Dynamics" (1999), pages 310-313

   u(t^(n+0.5)) = 1.5u(t^n) - 0.5u(t^(n-1))
   x_* = x^(n+1) - 0.5*dt*u(x^(n+1),t^n)
   x_j^n = x^(n+1) - dt*u(x_*,t^(n+0.5))

   Force term:
      F = 0.5*[3*S(x_j^n,t^n) - S(x_j^(n-1),t^(n-1))]
 */
void IntegrateRungeKutta_StgVariableVelocity( FeVariable* currVelField, FeVariable* interVelField, double dt, double* origin, double* position ) {
   unsigned		nDims		     = Mesh_GetDimSize( currVelField->feMesh );
   unsigned		dim_i;
   double			min[3], max[3];
   double			midPoint[3];
   double			velCurr[3], velInter[3];
   CartesianGenerator*	gen		     = (CartesianGenerator*)currVelField->feMesh->generator;

   Mesh_GetGlobalCoordRange( currVelField->feMesh, min, max );

   FieldVariable_InterpolateValueAt( currVelField, origin, velCurr );
   for( dim_i = 0; dim_i < nDims; dim_i++ ) {
      midPoint[dim_i] = origin[dim_i] - 0.5 * dt * velCurr[dim_i];
      PeriodicUpdate( midPoint, min, max, dim_i, gen->periodic[dim_i] );
   }

   FieldVariable_InterpolateValueAt( interVelField, midPoint, velInter );

   /* 2nd order approximation of velocity at time current + dt/2 */
   for( dim_i = 0; dim_i < nDims; dim_i++ ) {
      position[dim_i] = origin[dim_i] - dt * velInter[dim_i];
      PeriodicUpdate( position, min, max, dim_i, gen->periodic[dim_i] );
   }
}

/* for case of temporally evolving velocity field, when we need to integrate backwards in time
 * throught this to find our take off point */
void SemiLagrangianIntegrator_UpdatePreviousVelocityField( void* _self, void* _context ) {
   SemiLagrangianIntegrator*	self		= (SemiLagrangianIntegrator*) _self;
   FeVariable*			currVelField	= self->velocityField;
   FeVariable*			prevVelField	= self->prevVelField;
   unsigned			node_i;
   double				vel[3];

   for( node_i = 0; node_i < Mesh_GetLocalSize( currVelField->feMesh, MT_VERTEX ); node_i++ ) {
      FeVariable_GetValueAtNode( currVelField, node_i, vel );
      FeVariable_SetValueAtNode( prevVelField, node_i, vel );
   }

   FeVariable_SyncShadowValues( currVelField );
   FeVariable_SyncShadowValues( prevVelField );
}

/* cubic Lagrangian interpoation in 1-D */
void InterpLagrange( double x, double* coords, double (*values)[3], unsigned numdofs, double* result ) {
   unsigned	node_i, dof_i;
   unsigned	otherIndices[3];
   unsigned	otherIndexCount, otherIndex_i;
   double		factor;

   for( dof_i = 0; dof_i < numdofs; dof_i++ )
      result[dof_i] = 0.0;

   for( node_i = 0; node_i < 4; node_i++ ) {
      otherIndexCount = 0;
      for( otherIndex_i = 0; otherIndex_i < 4; otherIndex_i++ )
         if( otherIndex_i != node_i )
            otherIndices[otherIndexCount++] = otherIndex_i;

      factor = 1.0;
      for( otherIndex_i = 0; otherIndex_i < 3; otherIndex_i++ )
         factor *= ( x - coords[otherIndices[otherIndex_i]] ) / ( coords[node_i] - coords[otherIndices[otherIndex_i]] );

      for( dof_i = 0; dof_i < numdofs; dof_i++ ) {
         result[dof_i] += ( values[node_i][dof_i] * factor );
      }
   }
}

Bool PeriodicUpdate( double* pos, double* min, double* max, unsigned dim, Bool isPeriodic ) {
  if( pos[dim] < min[dim] ) {
    pos[dim] = (isPeriodic) ? max[dim] - min[dim] + pos[dim] : min[dim];
    return True;
  }
  if( pos[dim] > max[dim] ) {
    pos[dim] = (isPeriodic) ? min[dim] - max[dim] + pos[dim] : max[dim];
    return True;
  }

  return False;
}

Bool BicubicInterpolatorNew( FeVariable* feVariable, FeVariable* stencilField, double* position, unsigned* sizes, double* result ) {
  /* Calculated the BicubicInterpolation of the feVariable at position
   *
   * Input Args:
   *   feVariable: the field to be interpolated at `position`.
   *   stencilField: the field of initial spline stencils for each node.
   *   position:   the position of interpolatation.
   *   sizes:      memory chunk of size = sizeof(double)*dim.
   *   results:    the interpolated value.
   *
   *
   * Returns Values:
   *  True: if interpolation successful
   *  False: position is not in 'domain' of local processor
   */

   FeMesh*	feMesh = feVariable->feMesh;
   unsigned	nInc ;
   int    ijk[3];
   int		x_i, y_i, z_i, *inc;
   double localMin[3], localMax[3], double_ijk[3];
   Index  gNode_I, lNode_I, elementIndex;
   double px[4], py[4], pz[4];
   unsigned	nodeIndex[4][4];
   unsigned	node_I3D[4][4][4];
   unsigned	nDims   = Mesh_GetDimSize( feMesh );
   unsigned	numdofs = feVariable->dofLayout->dofCounts[0];
   double   ptsX[4][3], ptsY[4][3], ptsZ[4][3];
   Bool     inDomain  = True;

   inDomain = Mesh_SearchElements( feMesh, position, &elementIndex ); // get the element id

   if (!inDomain) return False;

   FeMesh_GetElementNodes( feMesh, elementIndex, feVariable->inc ); // get the incidence graph (inc.) of nodes on the element
   nInc = IArray_GetSize( feVariable->inc );                        // from inc. get the number of nodes on the element
   inc = IArray_GetPtr( feVariable->inc );                          // get the node ids from inc.

   FeVariable_GetValueAtNode( stencilField, inc[0], &(double_ijk[0]) );
   ijk[0] = lround(double_ijk[0]);
   ijk[1] = lround(double_ijk[1]);
   ijk[2] = lround(double_ijk[2]);

   /* interpolate using Lagrange's formula */
   if( nDims == 2 ) {
      for( y_i = 0; y_i < 4; y_i++ )
         for( x_i = 0; x_i < 4; x_i++ ) {
            gNode_I = ijk[0] + x_i + ( ijk[1] + y_i ) * sizes[0];
            if( !Mesh_GlobalToDomain( feMesh, MT_VERTEX, gNode_I, &lNode_I ) ){
              printf("Error in %s, trying to build an interpolation to position (%g, %g) using node %d, a non domain node, in interpolator.\n", __func__, position[0], position[1], gNode_I);
              abort();
            }
            else
               nodeIndex[x_i][y_i] = lNode_I;
         }
   }
   else {
      for( z_i = 0; z_i < 4; z_i++ )
         for( y_i = 0; y_i < 4; y_i++ )
            for( x_i = 0; x_i < 4; x_i++ ) {
               gNode_I = ijk[0] + x_i + ( ijk[1] + y_i ) * sizes[0] + ( ijk[2] + z_i ) * sizes[0] * sizes[1];
               if( !Mesh_GlobalToDomain( feMesh, MT_VERTEX, gNode_I, &lNode_I ) ) {
                 printf("Error in %s, trying to build an interpolation to position (%g, %g, %g) using node %d, a non domain node, in interpolator.\n", __func__, position[0], position[1], position[2], gNode_I);
                 abort();
               }
               else
                  node_I3D[x_i][y_i][z_i] = lNode_I;
            }
   }

   if( nDims == 3 ) {
      for( x_i = 0; x_i < 4; x_i++ )
         px[x_i] = Mesh_GetVertex( feMesh, node_I3D[x_i][0][0] )[0];
      for( y_i = 0; y_i < 4; y_i++ )
         py[y_i] = Mesh_GetVertex( feMesh, node_I3D[0][y_i][0] )[1];
      for( z_i = 0; z_i < 4; z_i++ )
         pz[z_i] = Mesh_GetVertex( feMesh, node_I3D[0][0][z_i] )[2];

      for( z_i = 0; z_i < 4; z_i++ ) {
         for( y_i = 0; y_i < 4; y_i++ ) {
            for( x_i = 0; x_i < 4; x_i++ )
               FeVariable_GetValueAtNode( feVariable, node_I3D[x_i][y_i][z_i], ptsX[x_i] );

            InterpLagrange( position[0], px, ptsX, numdofs, ptsY[y_i] );
         }

         InterpLagrange( position[1], py, ptsY, numdofs, ptsZ[z_i] );
      }

      InterpLagrange( position[2], pz, ptsZ, numdofs, result );
   }
   else {
      for( x_i = 0; x_i < 4; x_i++ )
         px[x_i] = Mesh_GetVertex( feMesh, nodeIndex[x_i][0] )[0];
      for( y_i = 0; y_i < 4; y_i++ )
         py[y_i] = Mesh_GetVertex( feMesh, nodeIndex[0][y_i] )[1];

      for( y_i = 0; y_i < 4; y_i++ ) {
         for( x_i = 0; x_i < 4; x_i++ )
            FeVariable_GetValueAtNode( feVariable, nodeIndex[x_i][y_i], ptsX[x_i] );

         InterpLagrange( position[0], px, ptsX, numdofs, ptsY[y_i] );
      }

      InterpLagrange( position[1], py, ptsY, numdofs, result );
   }

   return True;
}

Bool SemiLagrangianIntegrator_PointsAreClose( double* p1, double* p2, int dim, double rtol, double atol ) {
  /* check if two points are within rtol (relative tolerance)
   * or atol (absolute tolerance)
   *
   * Intput Parameters:
   *  p1   : point 1
   *  p2   : point 2
   *  dim  : the dimensions of points
   *  rtol : the relative tolerance
   *  atol : the absolute tolerance
   *
   * Return Values:
   *  True:  if points are close
   *  False: if points are not close
   */

  double disp[3], p2_norm, length;

  StGermain_VectorSubtraction(disp, p1, p2, dim); // displacement vector
  length  = StGermain_VectorMagnitude(disp, dim); // length of vector
  p2_norm = StGermain_VectorMagnitude(p2, dim);   // original size of p2

  if (length < rtol*p2_norm + atol) return True;

  return False;
}

void SemiLagrangianIntegrator_BuildStaticStencils( FeVariable* stencilField ) {
  /* Function to build the node indices for cubic interpolation.
   * The idea is to find a record the starting node indices (ijk) for each interpolant stencil.
   *
   * **NOTE**: We never want to Sync the stencilField FeVariable shadow values, ie.
   *  FeVariable_SyncShadowValues( stencilField ), because the field contains processor
   *  specific ordering and Syncing with shadow space will produce erroroneous starting indicies
   *
   */

  FeMesh* feMesh = stencilField->feMesh;
  Grid **grid;
  int d_i,ijk[3];
  double double_ijk[3];
  unsigned *sizes,try,nDims,n_i, nNodes;
  Index gNode_I;

  nDims = Mesh_GetDimSize( feMesh );
  nNodes = Mesh_GetDomainSize( feMesh, MT_VERTEX );
  grid = (Grid**)ExtensionManager_Get( feMesh->info, feMesh, feMesh->vertGridId  );
  sizes = Grid_GetSizes(*grid);

  for( n_i=0; n_i<nNodes; n_i++) {
    /* For every domain node, build the stencil starting location, ijk */
    gNode_I = Mesh_DomainToGlobal( feMesh, MT_VERTEX, n_i );
    Grid_Lift( *grid, gNode_I, ijk );

    /* for every dim try go one node back */
    for(d_i=0; d_i<nDims; d_i++) {
      if( ijk[d_i] == 0 ) { continue; } // skip if at 0
      ijk[d_i]--;
      try = Grid_Project( *grid, ijk );
      // if not in domain space don't go back one
      if( !Mesh_GlobalToDomain( feMesh, MT_VERTEX, try, &try ) ) {
        ijk[d_i]++;
      }
    }

    /* for every dim try go 3 forward */
    for(d_i=0; d_i<nDims; d_i++) {
      ijk[d_i]+=3;

      // if we go off the edge, go back one
      if( ijk[d_i] >= sizes[d_i] ) {
        ijk[d_i] -= 4;
        continue;
      }
      try = Grid_Project( *grid, ijk );
      // if not in domain space go back one, else reset
      if( !Mesh_GlobalToDomain( feMesh, MT_VERTEX, try, &try ) ) {
        ijk[d_i]-=4;
      } else {
        ijk[d_i]-=3;
      }
    }

    double_ijk[0] = (double)ijk[0];
    double_ijk[1] = (double)ijk[1];
    double_ijk[2] = (double)ijk[2];

    // save ijk as the starting point
    FeVariable_SetValueAtNode( stencilField, n_i, &(double_ijk[0]) );
   }
}

void SemiLagrangianIntegrator_SolveNew( FeVariable* variableField, double dt, FeVariable* velocityField, FeVariable* varStarField, FeVariable* stencilField ) {
  /* Function evaluates varStarField - an interpolation of variableField taken at the departure points.
   * Departure points are positions taken from the nodes and advected backwards along the characteristic curves.
   * The interpolation method used is a cubic spline and the implementation is only compatible with orthogonal meshes.
   *
   * Input Args:
   *   feVariable:    the original field to be interpolated.
   *   dt:            the time step size to go backwards along the characteristic, it should NOT be > CFL condition.
   *   velocityField: the velocity field used to go backward.
   *   stencilField:  the field of initial spline stencils for each node.
   *   varStarField:  the resultant field of interpolated variableField taken at the departure points.
   *
   */

   FeMesh*   feMesh   = variableField->feMesh;
   unsigned  meshSize = Mesh_GetLocalSize( feMesh, MT_VERTEX );
   unsigned  nDims    = Mesh_GetDimSize( feMesh );
   Grid**    nodegrid = (Grid**) Mesh_GetExtension( feMesh, Grid*,  feMesh->vertGridId );
   unsigned* sizes    = Grid_GetSizes( *nodegrid );

   unsigned  node_I;
   double    var[3], x_i[3], delta[3], minLength, *x_0;
   Bool      result;

   Mesh_GetMinimumSeparation( feMesh, &minLength, delta );

   /* sync parallel field variables to get shadow values */
   FeVariable_SyncShadowValues( velocityField );
   FeVariable_SyncShadowValues( variableField );

   /* assume that the variable mesh is the same as the velocity mesh */
   for( node_I = 0; node_I < meshSize; node_I++ ) {
      /* find the position back in time (u*), x_i */
      x_0 = Mesh_GetVertex(feMesh, node_I);
      IntegrateRungeKutta( velocityField, dt, x_0, x_i );

      /* if the new x_i is "close" to original node, don't Bicubuic Interpolate, take original node value */
      if( SemiLagrangianIntegrator_PointsAreClose(x_i, x_0, nDims, 0, 1e-6*minLength) ) {
        FeVariable_GetValueAtNode( variableField, node_I, var );
      } else {
        result = BicubicInterpolatorNew(variableField, stencilField, x_i, sizes, var);

        /* if BicubicInerpolator returns false, x_i was not found in domain space.
         * Fallback to using the node value. */
        if(result == False) { FeVariable_GetValueAtNode( variableField, node_I, var ); }
      }

      FeVariable_SetValueAtNode( varStarField, node_I, var );
   }

   /* sync interpolated values */
   FeVariable_SyncShadowValues( varStarField );
}

Bool BicubicInterpolator( FeVariable* feVariable, double* position, double* delta, unsigned* nNodes, double* result ) {
  /* Calculated the BicubicInterpolation of the feVariable at position
   *
   * Input Args:
   *   feVariable: the field to be interpolated at `position`.
   *   position:   the position of interpolatation.
   *   delta:      memory chunk of size = sizeof(double)*dim.
   *   nNodes:     numbers of nodes in ijk representation.
   *   results:    the interpolated value.
   *
   *
   * Returns Values:
   *  True: if interpolation successful
   *  False: position is not in 'domain' of local processor
   */

   FeMesh*	feMesh			= feVariable->feMesh;
   unsigned	nInc ;
   int    ijk[3];
   int		x_i, y_i, z_i, *inc;
   double localMin[3], localMax[3];
   Index  gNode_I, lNode_I, elementIndex;
   double px[4], py[4], pz[4];
   unsigned	nodeIndex[4][4];
   unsigned	node_I3D[4][4][4];
   unsigned	nDims   = Mesh_GetDimSize( feMesh );
   unsigned	numdofs = feVariable->dofLayout->dofCounts[0];
   double   ptsX[4][3], ptsY[4][3], ptsZ[4][3];
   Bool     inDomain  = True;

   inDomain = Mesh_SearchElements( feMesh, position, &elementIndex ); // get the element id

   if (!inDomain) return False;

   FeMesh_GetElementNodes( feMesh, elementIndex, feVariable->inc ); // get the incidence graph (inc.) of nodes on the element
   nInc = IArray_GetSize( feVariable->inc );                        // from inc. get the number of nodes on the element
   inc = IArray_GetPtr( feVariable->inc );                          // get the node ids from inc.

   if( nInc % 3 == 0 ) /* quadratic elements */ {
      delta[0] = Mesh_GetVertex( feMesh, inc[1] )[0] - Mesh_GetVertex( feMesh, inc[0] )[0];
      delta[1] = Mesh_GetVertex( feMesh, inc[3] )[1] - Mesh_GetVertex( feMesh, inc[0] )[1];
      if( nDims == 3 )
         delta[2] = Mesh_GetVertex( feMesh, inc[9] )[2] - Mesh_GetVertex( feMesh, inc[0] )[2];
   }
   else {
      delta[0] = Mesh_GetVertex( feMesh, inc[1] )[0] - Mesh_GetVertex( feMesh, inc[0] )[0];
      delta[1] = Mesh_GetVertex( feMesh, inc[2] )[1] - Mesh_GetVertex( feMesh, inc[0] )[1];
      if( nDims == 3 )
         delta[2] = Mesh_GetVertex( feMesh, inc[4] )[2] - Mesh_GetVertex( feMesh, inc[0] )[2];
   }

   Mesh_GetLocalCoordRange( feMesh, localMin, localMax );
   gNode_I = Mesh_DomainToGlobal( feMesh, MT_VERTEX, inc[0] );

   {
    Grid **grid;
    int d_i;
    unsigned *sizes;
    grid = (Grid**)ExtensionManager_Get( feMesh->info, feMesh, feMesh->vertGridId  );
    Grid_Lift( *grid, gNode_I, ijk );
    sizes = Grid_GetSizes(*grid);

    if( nInc % 2 == 0 ) {
      unsigned try;
      /* for every dim try go one node back */
      for(d_i=0; d_i<nDims; d_i++) {
        if( ijk[d_i] == 0 ) { continue; } // skip if at 0
        ijk[d_i]--;
        try = Grid_Project( *grid, ijk );
        // if not in domain space don't go back one
        if( !Mesh_GlobalToDomain( feMesh, MT_VERTEX, try, &try ) ) {
          ijk[d_i]++;
        }
      }
      /* for every dim try go 3 forward */
      for(d_i=0; d_i<nDims; d_i++) {
        ijk[d_i]+=3;

        // if we go off the edge, go back one
        if( ijk[d_i] >= sizes[d_i] ) {
          ijk[d_i] -= 4;
          continue;
        }
        try = Grid_Project( *grid, ijk );
        // if not in domain space go back one, else reset
        if( !Mesh_GlobalToDomain( feMesh, MT_VERTEX, try, &try ) ) {
          ijk[d_i]-=4;
        } else {
          ijk[d_i]-=3;
        }
      }
    }
   }

   /* interpolate using Lagrange's formula */
   if( nDims == 2 ) {
      for( y_i = 0; y_i < 4; y_i++ )
         for( x_i = 0; x_i < 4; x_i++ ) {
            gNode_I = ijk[0] + x_i + ( ijk[1] + y_i ) * nNodes[0];
            if( !Mesh_GlobalToDomain( feMesh, MT_VERTEX, gNode_I, &lNode_I ) ){
              printf("Error in BicubicInterpolator(), trying to build an interpolation to position (%g, %g) using node %d, a non domain node, in interpolator.\n", position[0], position[1], gNode_I);
              abort();
            }
            else
               nodeIndex[x_i][y_i] = lNode_I;
         }
   }
   else {
      for( z_i = 0; z_i < 4; z_i++ )
         for( y_i = 0; y_i < 4; y_i++ )
            for( x_i = 0; x_i < 4; x_i++ ) {
               gNode_I = ijk[0] + x_i + ( ijk[1] + y_i ) * nNodes[0] + ( ijk[2] + z_i ) * nNodes[0] * nNodes[1];
               if( !Mesh_GlobalToDomain( feMesh, MT_VERTEX, gNode_I, &lNode_I ) ) {
                 printf("Error in BicubicInterpolator(), trying to build an interpolation to position (%g, %g, %g) using node %d, a non domain node, in interpolator.\n", position[0], position[1], position[2], gNode_I);
                 abort();
               }
               else
                  node_I3D[x_i][y_i][z_i] = lNode_I;
            }
   }

   if( nDims == 3 ) {
      for( x_i = 0; x_i < 4; x_i++ )
         px[x_i] = Mesh_GetVertex( feMesh, node_I3D[x_i][0][0] )[0];
      for( y_i = 0; y_i < 4; y_i++ )
         py[y_i] = Mesh_GetVertex( feMesh, node_I3D[0][y_i][0] )[1];
      for( z_i = 0; z_i < 4; z_i++ )
         pz[z_i] = Mesh_GetVertex( feMesh, node_I3D[0][0][z_i] )[2];

      for( z_i = 0; z_i < 4; z_i++ ) {
         for( y_i = 0; y_i < 4; y_i++ ) {
            for( x_i = 0; x_i < 4; x_i++ )
               FeVariable_GetValueAtNode( feVariable, node_I3D[x_i][y_i][z_i], ptsX[x_i] );

            InterpLagrange( position[0], px, ptsX, numdofs, ptsY[y_i] );
         }

         InterpLagrange( position[1], py, ptsY, numdofs, ptsZ[z_i] );
      }

      InterpLagrange( position[2], pz, ptsZ, numdofs, result );
   }
   else {
      for( x_i = 0; x_i < 4; x_i++ )
         px[x_i] = Mesh_GetVertex( feMesh, nodeIndex[x_i][0] )[0];
      for( y_i = 0; y_i < 4; y_i++ )
         py[y_i] = Mesh_GetVertex( feMesh, nodeIndex[0][y_i] )[1];

      for( y_i = 0; y_i < 4; y_i++ ) {
         for( x_i = 0; x_i < 4; x_i++ )
            FeVariable_GetValueAtNode( feVariable, nodeIndex[x_i][y_i], ptsX[x_i] );

         InterpLagrange( position[0], px, ptsX, numdofs, ptsY[y_i] );
      }

      InterpLagrange( position[1], py, ptsY, numdofs, result );
   }

   return True;
}

void SemiLagrangianIntegrator_Solve( void* slIntegrator, FeVariable* variableField, FeVariable* varStarField ) {
   SemiLagrangianIntegrator*	self 		     = (SemiLagrangianIntegrator*)slIntegrator;
   FiniteElementContext*		context		     = self->context;
   unsigned			node_I;
   FeMesh*			feMesh		     = variableField->feMesh;
   unsigned			meshSize	     = Mesh_GetLocalSize( feMesh, MT_VERTEX );
   FeVariable*			velocityField	     = self->velocityField;
   double				dt		     = AbstractContext_Dt( context );
   unsigned			dim_I;
   unsigned			nDims		     = Mesh_GetDimSize( feMesh );
   double				position[3];
   double				var[3];
   Grid**				grid		     = (Grid**) Mesh_GetExtension( feMesh, Grid*,  feMesh->elGridId );
   unsigned*    sizes		     = Grid_GetSizes( *grid );
   unsigned			nNodes[3];
   double				delta[3];
   unsigned			nInc;
   int*         inc;

   FeMesh_GetElementNodes( variableField->feMesh, 0, variableField->inc );
   nInc = IArray_GetSize( variableField->inc );
   inc = IArray_GetPtr( variableField->inc );

   delta[0] = Mesh_GetVertex( feMesh, inc[1] )[0] - Mesh_GetVertex( feMesh, inc[0] )[0];
   if( nInc % 3 == 0 ) /* quadratic elements */ {
      delta[1] = Mesh_GetVertex( feMesh, inc[3] )[1] - Mesh_GetVertex( feMesh, inc[0] )[1];
      if( nDims == 3 )
         delta[2] = Mesh_GetVertex( feMesh, inc[9] )[2] - Mesh_GetVertex( feMesh, inc[0] )[2];
      for( dim_I = 0; dim_I < nDims; dim_I++ )
         nNodes[dim_I] = 2 * sizes[dim_I] + 1;
   }
   else {
      delta[1] = Mesh_GetVertex( feMesh, inc[2] )[1] - Mesh_GetVertex( feMesh, inc[0] )[1];
      if( nDims == 3 )
         delta[2] = Mesh_GetVertex( feMesh, inc[4] )[2] - Mesh_GetVertex( feMesh, inc[0] )[2];

      for( dim_I = 0; dim_I < nDims; dim_I++ )
         nNodes[dim_I] = sizes[dim_I] + 1;
   }

   FeVariable_SyncShadowValues( velocityField );
   FeVariable_SyncShadowValues( variableField );

   /* assume that the variable mesh is the same as the velocity mesh */
   for( node_I = 0; node_I < meshSize; node_I++ ) {
      /* find the position back in time (u*) */
      IntegrateRungeKutta( velocityField, dt, Mesh_GetVertex(feMesh, node_I), position );

      /* create a bicubic interpolation of variableField at u* */
      BicubicInterpolator( variableField, position, delta, nNodes, var );

      FeVariable_SetValueAtNode( varStarField, node_I, var );
   }
   FeVariable_SyncShadowValues( varStarField );
}


text:
#pragma once

#include <bgfx/bgfx.h>
#include <napi/napi.h>
#include <gsl/gsl>
#include <map>
#include <optional>

namespace Babylon
{
    class VertexBuffer final
    {
    public:
        VertexBuffer(gsl::span<uint8_t> bytes, bool dynamic);
        ~VertexBuffer();

        void Dispose();

        void Update(Napi::Env env, gsl::span<uint8_t> bytes);
        bool CreateHandle(const bgfx::VertexLayout& layout);
        void PromoteToFloats(bgfx::AttribType::Enum attribType, uint32_t numElements, uint32_t byteOffset, uint32_t byteStride);
        void Set(bgfx::Encoder* encoder, uint8_t stream, uint32_t startVertex, uint32_t numVertices, bgfx::VertexLayoutHandle layoutHandle);

        struct InstanceVertexBufferRecord
        {
            VertexBuffer* Buffer{};
            uint32_t Offset{};
            uint32_t Stride{};
            uint32_t ElementSize{};
        };
        static void BuildInstanceDataBuffer(bgfx::InstanceDataBuffer& instanceDataBuffer, const std::map<bgfx::Attrib::Enum, InstanceVertexBufferRecord>& vertexBufferInstance);
    private:
        std::optional<std::vector<uint8_t>> m_bytes{};
        bool m_dynamic{};

        union
        {
            bgfx::VertexBufferHandle m_handle{bgfx::kInvalidHandle};
            bgfx::DynamicVertexBufferHandle m_dynamicHandle;
        };

        bool m_disposed{};
    };
}

ext: h

max_stars_repo_name: HarveyLijh/ReactNativeBabylon

max_issues_repo_name: HarveyLijh/ReactNativeBabylon

max_forks_repo_name: HarveyLijh/ReactNativeBabylon


text:
/**
 *
 * @file testing_sposv.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Bilel Hadri, Hatem Ltaief
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:45:18 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

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

enum blas_order_type {
            blas_rowmajor = 101,
            blas_colmajor = 102 };

enum blas_cmach_type {
            blas_base      = 151,
            blas_t         = 152,
            blas_rnd       = 153,
            blas_ieee      = 154,
            blas_emin      = 155,
            blas_emax      = 156,
            blas_eps       = 157,
            blas_prec      = 158,
            blas_underflow = 159,
            blas_overflow  = 160,
            blas_sfmin     = 161};

enum blas_norm_type {
            blas_one_norm       = 171,
            blas_real_one_norm  = 172,
            blas_two_norm       = 173,
            blas_frobenius_norm = 174,
            blas_inf_norm       = 175,
            blas_real_inf_norm  = 176,
            blas_max_norm       = 177,
            blas_real_max_norm  = 178 };

static void
BLAS_error(char *rname, int err, int val, int x) {
  fprintf( stderr, "%s %d %d %d\n", rname, err, val, x );
  abort();
}

static
void
BLAS_sge_norm(enum blas_order_type order, enum blas_norm_type norm,
  int m, int n, const float *a, int lda, float *res) {
  int i, j; float anorm, v;
  char rname[] = "BLAS_sge_norm";

  if (order != blas_colmajor) BLAS_error( rname, -1, order, 0 );

  if (norm == blas_frobenius_norm) {
    anorm = 0.0f;
    for (j = n; j; --j) {
      for (i = m; i; --i) {
        v = a[0];
        anorm += v * v;
        a++;
      }
      a += lda - m;
    }
    anorm = sqrt( anorm );
  } else if (norm == blas_inf_norm) {
    anorm = 0.0f;
    for (i = 0; i < m; ++i) {
      v = 0.0f;
      for (j = 0; j < n; ++j) {
        v += fabsf( a[i + j * lda] );
      }
      if (v > anorm)
        anorm = v;
    }
  } else {
    BLAS_error( rname, -2, norm, 0 );
    return;
  }

  if (res) *res = anorm;
}

static
float
BLAS_spow_di(float x, int n) {
  float rv = 1.0;

  if (n < 0) {
    n = -n;
    x = 1.0 / x;
  }

  for (; n; n >>= 1, x *= x) {
    if (n & 1)
      rv *= x;
  }

  return rv;
}

static
float
BLAS_sfpinfo(enum blas_cmach_type cmach) {
  float eps = 1.0, r = 1.0, o = 1.0, b = 2.0;
  int t = 53, l = 1024, m = -1021;
  char rname[] = "BLAS_sfpinfo";

  if ((sizeof eps) == sizeof(float)) {
    t = 24;
    l = 128;
    m = -125;
  } else {
    t = 53;
    l = 1024;
    m = -1021;
  }

  /* for (i = 0; i < t; ++i) eps *= half; */
  eps = BLAS_spow_di( b, -t );
  /* for (i = 0; i >= m; --i) r *= half; */
  r = BLAS_spow_di( b, m-1 );

  o -= eps;
  /* for (i = 0; i < l; ++i) o *= b; */
  o = (o * BLAS_spow_di( b, l-1 )) * b;

  switch (cmach) {
    case blas_eps: return eps;
    case blas_sfmin: return r;
    default:
      BLAS_error( rname, -1, cmach, 0 );
      break;
  }
  return 0.0;
}

static int check_factorization(int, float*, float*, int, int , float);
static int check_solution(int, int, float*, int, float*, float*, int, float);
static int check_estimator(PLASMA_enum, int, float *, int,
                           float *, float, float,
                           float);

int testing_sposv(int argc, char **argv)
{

    /* Check for number of arguments*/
    if (argc != 4){
        USAGE("POSV", "N LDA NRHS LDB",
              "   - N    : the size of the matrix\n"
              "   - LDA  : leading dimension of the matrix A\n"
              "   - NRHS : number of RHS\n"
              "   - LDB  : leading dimension of the RHS B\n");
        return -1;
    }

    int N     = atoi(argv[0]);
    int LDA   = atoi(argv[1]);
    int NRHS  = atoi(argv[2]);
    int LDB   = atoi(argv[3]);
    float eps;
    int info_solution, info_factorization;
    int u, trans1, trans2;

    float *A1   = (float *)malloc(LDA*N*sizeof(float));
    float *A2   = (float *)malloc(LDA*N*sizeof(float));
    float *B1   = (float *)malloc(LDB*NRHS*sizeof(float));
    float *B2   = (float *)malloc(LDB*NRHS*sizeof(float));

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

    eps = BLAS_sfpinfo( blas_eps );

    for(u=0; u<2; u++) {

        trans1 = uplo[u] == PlasmaUpper ? PlasmaTrans : PlasmaNoTrans;
        trans2 = uplo[u] == PlasmaUpper ? PlasmaNoTrans : PlasmaTrans;

        /*-------------------------------------------------------------
         *  TESTING SPOSV
         */

        /* Initialize A1 and A2 for Symmetric Positif Matrix */
        PLASMA_splgsy( (float)N, N, A1, LDA, 51 );
        PLASMA_slacpy( PlasmaUpperLower, N, N, A1, LDA, A2, LDA );

        /* Initialize B1 and B2 */
        PLASMA_splrnt( N, NRHS, B1, LDB, 371 );
        PLASMA_slacpy( PlasmaUpperLower, N, NRHS, B1, LDB, B2, LDB );

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

        /* PLASMA SPOSV */
        PLASMA_sposv(uplo[u], N, NRHS, A2, LDA, B2, LDB);

        /* Check the factorization and the solution */
        info_factorization = check_factorization( N, A1, A2, LDA, uplo[u], eps);
        info_solution = check_solution(N, NRHS, A1, LDA, B1, B2, LDB, eps);

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

        /*-------------------------------------------------------------
         *  TESTING SPOTRF + SPOTRS
         */

        /* Initialize A1 and A2 for Symmetric Positif Matrix */
        PLASMA_splgsy( (float)N, N, A1, LDA, 51 );
        PLASMA_slacpy( PlasmaUpperLower, N, N, A1, LDA, A2, LDA );

        /* Initialize B1 and B2 */
        PLASMA_splrnt( N, NRHS, B1, LDB, 371 );
        PLASMA_slacpy( PlasmaUpperLower, N, NRHS, B1, LDB, B2, LDB );

        /* Plasma routines */
        PLASMA_spotrf(uplo[u], N, A2, LDA);
        PLASMA_spotrs(uplo[u], N, NRHS, A2, LDA, B2, LDB);

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

        /* Check the factorization and the solution */
        info_factorization = check_factorization( N, A1, A2, LDA, uplo[u], eps);
        info_solution = check_solution(N, NRHS, A1, LDA, B1, B2, LDB, eps);

        if ((info_solution == 0)&(info_factorization == 0)){
            printf("***************************************************\n");
            printf(" ---- TESTING SPOTRF + SPOTRS (%s)............ PASSED !\n", uplostr[u]);
            printf("***************************************************\n");
        }
        else{
            printf("****************************************************\n");
            printf(" - TESTING SPOTRF + SPOTRS (%s)... FAILED !\n", uplostr[u]);
            printf("****************************************************\n");
        }

        /*-------------------------------------------------------------
         *  TESTING SPOTRF + ZPTRSM + STRSM
         */

        /* Initialize A1 and A2 for Symmetric Positif Matrix */
        PLASMA_splgsy( (float)N, N, A1, LDA, 51 );
        PLASMA_slacpy( PlasmaUpperLower, N, N, A1, LDA, A2, LDA );

        /* Initialize B1 and B2 */
        PLASMA_splrnt( N, NRHS, B1, LDB, 371 );
        PLASMA_slacpy( PlasmaUpperLower, N, NRHS, B1, LDB, B2, LDB );

        /* PLASMA routines */
        PLASMA_spotrf(uplo[u], N, A2, LDA);
        PLASMA_strsm(PlasmaLeft, uplo[u], trans1, PlasmaNonUnit,
                     N, NRHS, 1.0, A2, LDA, B2, LDB);
        PLASMA_strsm(PlasmaLeft, uplo[u], trans2, PlasmaNonUnit,
                     N, NRHS, 1.0, A2, LDA, B2, LDB);

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

        /* Check the factorization and the solution */
        info_factorization = check_factorization( N, A1, A2, LDA, uplo[u], eps);
        info_solution = check_solution(N, NRHS, A1, LDA, B1, B2, LDB, eps);

        if ((info_solution == 0)&(info_factorization == 0)){
            printf("***************************************************\n");
            printf(" ---- TESTING SPOTRF + STRSM + STRSM (%s)..... PASSED !\n", uplostr[u]);
            printf("***************************************************\n");
        }
        else{
            printf("***************************************************\n");
            printf(" - TESTING SPOTRF + STRSM + STRSM (%s)... FAILED !\n", uplostr[u]);
            printf("***************************************************\n");
        }

        /*-------------------------------------------------------------
         *  TESTING ZPOCON on the last call
         */
        {
            float Anorm = PLASMA_slansy( PlasmaOneNorm, uplo[u], N, A1, LDA );
            float Acond;

            info_solution = PLASMA_spocon(uplo[u], N, A2, LDA, Anorm, &Acond);
            if ( info_solution == 0 ) {
                info_solution = check_estimator(uplo[u], N, A1, LDA, A2, Anorm, Acond, eps);
            } else {
                printf(" PLASMA_spocon returned info = %d\n", info_solution );
            }
            if ((info_solution == 0)){
                printf("***************************************************\n");
                printf(" ---- TESTING SPOTRF + ZPOCON (%s) ........... PASSED !\n", uplostr[u]);
                printf("***************************************************\n");
            }
            else{
                printf("**************************************************\n");
                printf(" - TESTING SPOTRF + ZPOCON (%s) ... FAILED !\n", uplostr[u]);
                printf("**************************************************\n");
            }
        }
    }

    free(A1); free(A2); free(B1); free(B2);

    return 0;
}

/*------------------------------------------------------------------------
 *  Check the factorization of the matrix A2
 */
static int check_factorization(int N, float *A1, float *A2, int LDA, int uplo, float eps)
{
    float Anorm, Rnorm;
    float alpha;
    int info_factorization;
    int i,j;

    float *Residual = (float *)malloc(N*N*sizeof(float));
    float *L1       = (float *)malloc(N*N*sizeof(float));
    float *L2       = (float *)malloc(N*N*sizeof(float));
    float *work              = (float *)malloc(N*sizeof(float));

    memset((void*)L1, 0, N*N*sizeof(float));
    memset((void*)L2, 0, N*N*sizeof(float));

    alpha= 1.0;

    LAPACKE_slacpy_work(LAPACK_COL_MAJOR,' ', N, N, A1, LDA, Residual, N);

    /* Dealing with L'L or U'U  */
    if (uplo == PlasmaUpper){
        LAPACKE_slacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L1, N);
        LAPACKE_slacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L2, N);
        cblas_strmm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, N, N, (alpha), L1, N, L2, N);
    }
    else{
        LAPACKE_slacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L1, N);
        LAPACKE_slacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L2, N);
        cblas_strmm(CblasColMajor, CblasRight, CblasLower, CblasTrans, CblasNonUnit, N, N, (alpha), L1, N, L2, N);
    }

    /* Compute the Residual || A -L'L|| */
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
           Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i];

    BLAS_sge_norm( blas_colmajor, blas_inf_norm, N, N, Residual, N, &Rnorm );
    BLAS_sge_norm( blas_colmajor, blas_inf_norm, N, N, A1, LDA, &Anorm );

    printf("============\n");
    printf("Checking the Cholesky Factorization \n");
    printf("-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));

    if ( isnan(Rnorm/(Anorm*N*eps)) || isinf(Rnorm/(Anorm*N*eps)) || (Rnorm/(Anorm*N*eps) > 60.0) ){
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else{
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(Residual); free(L1); free(L2); free(work);

    return info_factorization;
}

/*------------------------------------------------------------------------
 *  Check the accuracy of the solution of the linear system
 */

static int check_solution(int N, int NRHS, float *A1, int LDA, float *B1, float *B2, int LDB, float eps )
{
    int info_solution;
    float Rnorm, Anorm, Xnorm, Bnorm, result;
    float alpha, beta;
    float *work = (float *)malloc(N*sizeof(float));

    alpha = 1.0;
    beta  = -1.0;

    BLAS_sge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B2, LDB, &Xnorm );
    BLAS_sge_norm( blas_colmajor, blas_inf_norm, N, N,    A1, LDA, &Anorm );
    BLAS_sge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B1, LDB, &Bnorm );

    cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, N, NRHS, N, (alpha), A1, LDA, B2, LDB, (beta), B1, LDB);
    BLAS_sge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B1, LDB, &Rnorm );

    if (getenv("PLASMA_TESTING_VERBOSE"))
      printf( "||A||_oo=%f\n||X||_oo=%f\n||B||_oo=%f\n||A X - B||_oo=%e\n", Anorm, Xnorm, Bnorm, Rnorm );

    result = Rnorm / ( (Anorm*Xnorm+Bnorm)*N*eps ) ;
    printf("============\n");
    printf("Checking the Residual of the solution \n");
    printf("-- ||Ax-B||_oo/((||A||_oo||x||_oo+||B||_oo).N.eps) = %e \n", result);

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

    free(work);

    return info_solution;
}

/*------------------------------------------------------------------------
 *  Check the accuracy of the condition estimator
 */
static int check_estimator(PLASMA_enum uplo, int N, float *A1, int LDA,
                           float *A2, float Anorm, float Acond,
                           float eps)
{
    int info_solution;
    float result, Acond_lapack;
    float invcond, invcond_lapack;

    info_solution = LAPACKE_spocon(LAPACK_COL_MAJOR, lapack_const(uplo), N, A2, LDA, Anorm, &Acond_lapack);

    if ( info_solution != 0 ) {
        printf(" PLASMA_sgecon returned info = %d\n", info_solution );
        return info_solution;
    }

    invcond_lapack = 1. / ( Acond_lapack );
    invcond        = 1. / ( Acond );

    printf("============\n");
    printf("Checking the condition number \n");
    printf("-- Acond_plasma = %e, Acond_lapack = %e \n"
           "-- Ainvcond_plasma = %e, Ainvcond_lapack = %e \n",
           Acond, Acond_lapack, invcond, invcond_lapack );

    result = fabs( Acond_lapack - Acond ) / eps;
    if ( result > 60. ) {
        printf("-- The solution is suspicious ! \n");
        info_solution = 1;
     }
    else{
        printf("-- The solution is CORRECT ! \n");
        info_solution = 0;
    }
    return info_solution;
}

ext: c


text:
#ifndef PyGSL_TRANSFORM_TYPES_H
#define PyGSL_TRANSFORM_TYPES_H 1
/*
 * All different transforms are handled by only two functions.
 *       PyGSL_transform_
 * and
 *       PyGSL_transform_2d_.
 *
 * This was done as similar functionallity is needed by all the transforms.
 * Copy the input data (transforms are often in place), check or allocate
 * help spaces. Call the function and make clean up. To allow that the
 * functions need to do quite some decisions to prepare all in proper order
 * and to do the clean up.
 * In this file various data types are defined. The first two are needed by the
 * python type which allows to allocate the proper help space and pass it to
 * the function. All the others are needed to store the information about the
 * transform.
 */

/* ------------------------------------------------------------------------- */
/*
 * One Python type is used to encapsulate all the different workspaces. This
 * enum allows the C Code to destinquish between the different types.
 */

enum pygsl_transform_space_type{
	NOSPACE=0,
	COMPLEX_WORKSPACE,
	REAL_WORKSPACE,
	COMPLEX_WAVETABLE,
	REAL_WAVETABLE,
	HALFCOMPLEX_WAVETABLE,
	COMPLEX_WORKSPACE_FLOAT,
	REAL_WORKSPACE_FLOAT,
	COMPLEX_WAVETABLE_FLOAT,
	REAL_WAVETABLE_FLOAT,
	HALFCOMPLEX_WAVETABLE_FLOAT,
	WAVELET_WORKSPACE
};

#include <pygsl/pygsl_features.h>
#include <gsl/gsl_fft.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>
#include <gsl/gsl_fft_complex_float.h>
#include <gsl/gsl_fft_real_float.h>
#include <gsl/gsl_fft_halfcomplex_float.h>
#include <gsl/gsl_blas.h>
#ifdef  _PYGSL_GSL_HAS_WAVELET
#define forward forward_wavelet
#define backward backward_wavelet
#include <gsl/gsl_wavelet.h>
#include <gsl/gsl_wavelet2d.h>
#undef forward
#undef backward
#else /* _PYGSL_GSL_HAS_WAVELET */
#endif /* _PYGSL_GSL_HAS_WAVELET */

/*
 * Here the corresponding union to address the pointer as the proper type.
 */
union pygsl_transform_space_t{
	gsl_fft_complex_workspace           *cws;
	gsl_fft_complex_wavetable           *cwt;
	gsl_fft_real_workspace              *rws;
	gsl_fft_real_wavetable              *rwt;
	gsl_fft_halfcomplex_wavetable       *hcwt;
	gsl_fft_complex_workspace_float     *cwsf;
	gsl_fft_complex_wavetable_float     *cwtf;
	gsl_fft_real_workspace_float        *rwsf;
	gsl_fft_real_wavetable_float        *rwtf;
	gsl_fft_halfcomplex_wavetable_float *hcwtf;
#ifdef _PYGSL_GSL_HAS_WAVELET
	gsl_wavelet_workspace               *wws;
#endif
	void                                *v;
};

/* ------------------------------------------------------------------------- */
/*
 * From here all data types are defined needed to guide the transform
 * functions.
 */
/*
 * Transformation in float or double mode
 */
enum pygsl_transform_mode{
	MODE_DOUBLE = 1,
	MODE_FLOAT
};

/*
 * Input data in ... output data in ...
 */
enum transform_mode{
	ComplexComplex = 1,
	RealReal,
	RealHalfcomplex,
	HalfComplexReal

};

/*
 * Should be renamed. Originally the transform mode would only handle the FFT.
 * There one only needed to destinquish between the radix. Wavelets are also a
 * bit different. Others to come.
 */
enum radix_mode{
	RADIX_TWO = 1,
	RADIX_FREE,
	WAVELET
};

/*
 * Does the basis type consist of only one machine type or two. e.g real
 * is SINGLE_TYPE, complex PACKED_TYPE (a real and an imaginary part).
 */
enum pygsl_packed_type{
	SINGLE_TYPE=1,
	PACKED_TYPE=2
};

#define PyGSL_TRANSFORM_MODE_SWITCH(mode, double_element, float_element) \
        ((mode == MODE_DOUBLE) ? double_element : float_element)

typedef int transform(void * data, size_t stride, size_t N, const void *, void *);
typedef int transform_r2(void * data, size_t stride, size_t N);
#ifdef  _PYGSL_GSL_HAS_WAVELET
typedef int wavelet(const gsl_wavelet * W, double *data, size_t stride, size_t N, void *);
typedef int wavelet2d(const gsl_wavelet * W, gsl_matrix *m, void *);
#endif
typedef void * pygsl_transform_helpn_t(int);
typedef void * pygsl_transform_help_t(void *);

/*
 * The different transform types.
 */
union transforms{
	transform *free;
	transform_r2 *radix2;
#ifdef  _PYGSL_GSL_HAS_WAVELET
	wavelet *wavelet;
	wavelet2d *wavelet2d;
#endif
	void *v;
};

/*
 * Functions to construct and destruct the helpers space and the helpers table.
 * The user can pass them as a proper python type. If not, the proper space
 * will be allocated using these functions.
 */
struct _pygsl_transform_func_rf_s {
	pygsl_transform_helpn_t  * space_alloc;
	pygsl_transform_help_t   * space_free;
	pygsl_transform_helpn_t  * table_alloc;
	pygsl_transform_help_t   * table_free;
	enum pygsl_transform_space_type space_type;
	enum pygsl_transform_space_type table_type;
};

/*
 * Not all transforms need additional workspace. The ones that need it will put
 * an instance somewhere and point the func to the approbriate initalisers.
 * free_space and free_table is used to store if the space and table have to be
 * freeded at the end of the transform.
 */
struct _pygsl_transform_help_rf_s{
	const struct _pygsl_transform_func_rf_s *func;
	void * space;
	void * table;
	int free_space;
	int free_table;
};

/*
 * The info about which transform, what input and output
 * arrays and so forth. Used by _pygsl_transform_help_s to inform the transform
 * function about the various properties.
 */
struct _pygsl_transform_info_s {
	/*
	 *  complex -> complex
	 *  real    -> half_complex
	 *  half_complex -> real
	 *  real -> real
	 */
	enum transform_mode mode;
        /* float or double ? */
	enum pygsl_transform_mode datatype;
	enum PyArray_TYPES input_array_type;
	enum PyArray_TYPES output_array_type;
        /*
	 *  Do I need to add some additional offset to the data
	 *  pointer before
	 *  calling the function. A trick needed for the real to
	 *  halfcomplex transform.
	 */
	int data_offset;
	/* Type of the transform. Should be renamed */
	enum radix_mode radix2;
	/* one or two machine types make on basis type ? */
	enum pygsl_packed_type packed;
};

/*
 * Finally the struct which gathers all the information. The split up was made, as the
 * info struct is always needed, but not the helpers.
 */
struct _pygsl_transform_help_s{
	struct _pygsl_transform_info_s *info;
	union transforms  transform;
	struct _pygsl_transform_help_rf_s *helpers;
};
typedef  struct _pygsl_transform_help_s pygsl_transform_help_s;

/*
 * The main function doing all the work. Its a function for the ffts and a
 *  method for the wavelets!
 */
static PyObject *
PyGSL_transform_(PyObject *self, PyObject *args, pygsl_transform_help_s *helps);

/*
 * Currently only wavelets provide a 2 dimensional transform
 */
#ifdef _PYGSL_GSL_HAS_WAVELET
static PyObject *
PyGSL_transform_2d_(PyObject *self, PyObject *args, pygsl_transform_help_s *helps);
#endif

#endif /* PyGSL_TRANSFORM_TYPES_H */
/*
 * Local Variables:
 * mode: C
 * c-file-style: "python"
 * End:
 */

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
#pragma once

#include <gsl/span>
#include <array>
#include "halley/core/api/audio_api.h"

namespace Halley
{
	using AudioSourceData = std::array<gsl::span<AudioConfig::SampleFormat>, AudioConfig::maxChannels>;

	class AudioSource {
	public:
		virtual ~AudioSource() {}

		virtual size_t getNumberOfChannels() const = 0;
		virtual bool isReady() const { return true; }
		virtual bool getAudioData(size_t numSamples, AudioSourceData& dst) = 0;
	};
}

max_stars_repo_name: lye/halley

max_issues_repo_name: lye/halley

max_forks_repo_name: lye/halley


text:
/*
 * C version of Diffusive Nested Sampling (DNest4) by Brendon J. Brewer
 *
 * Yan-Rong Li, liyanrong@mail.ihep.ac.cn
 * Jun 30, 2016
 *
 */
#ifndef _DNESTVARS_H
#define _DNESTVARS_H

#ifdef __cplusplus
extern "C" {
#endif

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>
#include <gsl/gsl_rng.h>

#include "dnest.h"

enum PRIOR_TYPE {UNIFORM=0, GAUSSIAN=1, LOG=2};

/* output files */
extern FILE *fsample, *fsample_info;

/* random number generator */
extern const gsl_rng_type * dnest_gsl_T;
extern gsl_rng * dnest_gsl_r;

typedef struct
{
  double value;
  double tiebreaker;
}LikelihoodType;

typedef struct
{
  LikelihoodType log_likelihood;
  double log_X;
  unsigned long long int visits, exceeds;
  unsigned long long int accepts, tries;
}Level;

extern DNestOptions options;
extern char options_file[STR_MAX_LENGTH];

typedef struct
{
  int id;
  void *addr;
  char tag[50];
  int isset;
}DNestPARDICT;

extern void *particles;
extern int dnest_size_of_modeltype;
extern int particle_offset_size, particle_offset_double;

// sampler
extern bool save_to_disk;
extern unsigned int num_threads;
extern double compression;
extern unsigned int regularisation;

extern void *particles;
extern LikelihoodType *log_likelihoods;
extern unsigned int *level_assignments;

// number account of unaccepted times
extern unsigned int *account_unaccepts;

extern int size_levels, size_levels_combine;
extern Level *levels;
extern Level *copies_of_levels, *levels_combine;
extern LikelihoodType *all_above;
extern unsigned int count_saves, num_saves, num_saves_restart;
extern unsigned long long int count_mcmc_steps;
extern LikelihoodType *above;
extern unsigned int size_above, size_all_above;

extern int dnest_flag_restart, dnest_flag_postprc, dnest_flag_sample_info, dnest_flag_limits;
extern double dnest_post_temp;
extern char file_restart[STR_MAX_LENGTH], file_save_restart[STR_MAX_LENGTH];

extern double post_logz;
extern int dnest_num_params;
extern char dnest_sample_postfix[STR_MAX_LENGTH], dnest_sample_tag[STR_MAX_LENGTH], dnest_sample_dir[STR_MAX_LENGTH];
extern double *dnest_param_range, *dnest_prior_info;
extern int *dnest_prior_type;
extern void *dnest_args;

//the limits of parameters for each level;
extern double *limits, *copies_of_limits;

extern int dnest_which_particle_update; // which particle to be updated
extern int dnest_which_level_update;    // which level to be updated;
extern int dnest_thistask, dnest_totaltask;
extern int *dnest_perturb_accept;
extern int dnest_root;
//***********************************************
/*                  functions                  */
extern double mod(double y, double x);
extern void wrap_limit(double *x, double min, double max);
extern int mod_int(int y, int x);
extern int dnest_cmp(const void *pa, const void *pb);

extern void options_load(char *optfile, DNestOptions *opts);
extern void setup(int argc, char** argv, DNestFptrSet *fptrset, int num_params,
           double *param_range, int *prior_type, double *prior_info,
           char *sample_dir, char *optfile, DNestOptions *opts, void *args);
extern void finalise();

extern void dnest_run();
extern void dnest_mcmc_run();
extern void update_particle(unsigned int which);
extern void update_level_assignment(unsigned int which);
extern double log_push(unsigned int which_level);
extern bool enough_levels(Level *l, int size_l);
extern void do_bookkeeping();
extern void save_levels();
extern void save_particle();
extern void save_limits();
extern void kill_lagging_particles();
extern void renormalise_visits();
extern void recalculate_log_X();
extern void dnest_postprocess(double temperature,char *optfile, DNestOptions *opts);
extern void postprocess(double temperature);
extern void initialize_output_file();
extern void close_output_file();
extern void dnest_save_restart();
extern void dnest_restart();
extern void dnest_restart_action(int iflag);
extern void dnest_accept_action();
extern void dnest_kill_action(int i, int i_copy);
extern void dnest_from_prior(void *model);
extern double dnest_perturb(void *model);
extern void dnest_print_particle(FILE *fp, const void *model);
extern void dnest_read_particle(FILE *fp, void *model);
extern void dnest_check_directory(char *sample_dir);
/*=====================================================*/
// users responsible for following functions
extern void (*print_particle)(FILE *fp, const void *model);
extern void (*read_particle)(FILE *fp, void *model);
extern void (*from_prior)(void *model);
extern double (*log_likelihoods_cal)(const void *model);
extern double (*log_likelihoods_cal_initial)(const void *model);
extern double (*log_likelihoods_cal_restart)(const void *model);
extern double (*perturb)(void *model);
extern void (*restart_action)(int iflag);
extern void (*accept_action)();
extern void (*kill_action)(int i, int i_copy);
/*=====================================================*/

#ifdef __cplusplus
}
#endif

#endif

max_stars_repo_name: LiyrAstroph/CDNest


text:
#ifndef CWANNIER_PARTIALDOS_VALUES_H
#define CWANNIER_PARTIALDOS_VALUES_H

#include <gsl/gsl_permutation.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include "HTightBinding.h"
#include "ctetra/partial.h"

double** PartialDosValues(HTightBinding *Hrs, gsl_matrix *R, int na, int nb, int nc, double sigma, double **Es, double num_dos);

void sort_evals_evecs(gsl_vector *energies, gsl_matrix_complex *evecs, int num_bands);

#endif // CWANNIER_PARTIALDOS_VALUES_H

max_stars_repo_name: tflovorn/cwannier

max_issues_repo_name: tflovorn/cwannier

max_forks_repo_name: tflovorn/cwannier