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/libStGermain/src/StGermain.h>
#include <StgDomain/libStgDomain/src/StgDomain.h>
#include "StgFEM/Discretisation/src/types.h"

#include "ElementType.h"
#include "ElementType_Register.h"
#include "Element.h"
#include "FeMesh.h"
#include "FeEquationNumber.h"
#include "FeVariable.h"
#include "LinkedDofInfo.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <mpi.h>
#include <petsc.h>
#include <petscvec.h>

// Private prototype
int GenerateEquationNumbering(
		int NX, int NY, int NZ,
		int nlocal, int g_node_id[],
		int dof, int nglobal,
		PetscTruth periodic_x, PetscTruth periodic_y, PetscTruth periodic_z,
		int npx, int npy, int npz,
		int periodic_x_gnode_id[], int periodic_y_gnode_id[], int periodic_z_gnode_id[],
		int eqnums[], int *neqnums );

int stgCmpInt( const void *l, const void *r ) {
   return *(int*)l - *(int*)r;
}

/*###### Typedefs and Structs ######*/

/** Textual name of this class */
const Type FeEquationNumber_Type = "FeEquationNumber";

/** struct to store sub-totals: what is the equation number up to at a given
    node? These can then be exchanged between processors */
typedef struct CritPointInfo {
      Node_GlobalIndex	index;
      Dof_EquationNumber	eqNum;
} CritPointInfo;

/** An element of linked list of critical point info. Several of the functions
    use this to keep track of key points */
typedef struct AddListEntry {
      CritPointInfo*		critPointInfo;
      struct AddListEntry*	next;
} AddListEntry;

/** Enum to say whetehr values at crit. nodes should be printed */
typedef enum PrintValuesFlag {
   DONT_PRINT_VALUES,
   PRINT_VALUES
} PrintValuesFlag;

/** MPI datatype handle for efficiently exchanging CritPointInfo structures.
	see FeEquationNumber_Create_CritPointInfo_MPI_Datatype() for where this
	handle is defined. */
MPI_Datatype MPI_critPointInfoType;

/*###### Private Function Declarations ######*/

/*###### Function Definitions ######*/

/** Public constructor */

FeEquationNumber* FeEquationNumber_New(
	Name						name,
	DomainContext*			context,
	void*						mesh,
	DofLayout*				dofLayout,
	VariableCondition*	bcs,
	LinkedDofInfo*			linkedDofInfo )
{
   FeEquationNumber* self = FeEquationNumber_DefaultNew( name );

	self->isConstructed = True;
	_FeEquationNumber_Init( self, context, mesh, dofLayout, bcs, linkedDofInfo );

	return self;
}

void* FeEquationNumber_DefaultNew( Name name ) {
	/* Variables set in this function */
	SizeT                                              _sizeOfSelf = sizeof(FeEquationNumber);
	Type                                                      type = FeEquationNumber_Type;
	Stg_Class_DeleteFunction*                              _delete = _FeEquationNumber_Delete;
	Stg_Class_PrintFunction*                                _print = _FeEquationNumber_Print;
	Stg_Class_CopyFunction*                                  _copy = _FeEquationNumber_Copy;
	Stg_Component_DefaultConstructorFunction*  _defaultConstructor = (Stg_Component_DefaultConstructorFunction*)FeEquationNumber_DefaultNew;
	Stg_Component_ConstructFunction*                    _construct = _FeEquationNumber_AssignFromXML;
	Stg_Component_BuildFunction*                            _build = _FeEquationNumber_Build;
	Stg_Component_InitialiseFunction*                  _initialise = _FeEquationNumber_Initialise;
	Stg_Component_ExecuteFunction*                        _execute = _FeEquationNumber_Execute;
	Stg_Component_DestroyFunction*                        _destroy = _FeEquationNumber_Destroy;
	AllocationType                              nameAllocationType = NON_GLOBAL;

   return _FeEquationNumber_New(  FEEQUATIONNUMBER_PASSARGS  );
}
/** Constructor implementation. */
FeEquationNumber* _FeEquationNumber_New(  FEEQUATIONNUMBER_DEFARGS  ){
   FeEquationNumber* self;

   /* Allocate memory */
   assert( _sizeOfSelf >= sizeof(FeEquationNumber) );
   self = (FeEquationNumber*)_Stg_Component_New(  STG_COMPONENT_PASSARGS  );

   /* General info */

   /* Virtual info */
   self->_build = _build;
   self->_initialise = _initialise;

   /* Mesh info */

   return self;
}

/** Constructor variables initialisation. Doesn't allocate any
    memory, just saves arguments and sets counters to 0. */
void _FeEquationNumber_Init(
   FeEquationNumber*		self,
	DomainContext*			context,
   void*						feMesh,
   DofLayout*				dofLayout,
   VariableCondition*	bcs,
   LinkedDofInfo*			linkedDofInfo )
{
   /* General and Virtual info should already be set */

   /* FinteElementMesh info */
	self->context = context;
   self->feMesh = (FeMesh*)feMesh;
   self->globalSumUnconstrainedDofs = 0;
   self->isBuilt = False;
   self->locationMatrixBuilt = False;
   self->_lowestLocalEqNum = -1;
   self->_lowestGlobalEqNums = NULL;
   self->_highestLocalEqNum = -1;
   self->dofLayout = dofLayout;
   self->bcs = bcs;
   self->linkedDofInfo = linkedDofInfo;
   self->remappingActivated = False;
   /* register streams */
   self->debug = Stream_RegisterChild( StgFEM_Discretisation_Debug, FeEquationNumber_Type );
   self->debugLM = Stream_RegisterChild( self->debug, "LM" );
   self->warning = Stream_RegisterChild( StgFEM_Warning, FeEquationNumber_Type );
   self->removeBCs = True;
   self->bcEqNums = STree_New();
   STree_SetIntCallbacks( self->bcEqNums );
   STree_SetItemSize( self->bcEqNums, sizeof(int) );
   self->ownedMap = STreeMap_New();
   STreeMap_SetItemSize( self->ownedMap, sizeof(int), sizeof(int) );
   STree_SetIntCallbacks( self->ownedMap );

   Stream_SetPrintingRank( self->debug, 0 );
}

void _FeEquationNumber_AssignFromXML( void* feEquationNumber, Stg_ComponentFactory *cf, void* data ) {
}

void _FeEquationNumber_Execute( void* feEquationNumber, void *data ){

}

void _FeEquationNumber_Destroy( void* feEquationNumber, void *data ){
   FeEquationNumber* self = (FeEquationNumber*) feEquationNumber;
   Index ii;

   /* free destination array memory */
   Journal_DPrintfL( self->debug, 2, "Freeing I.D. Array\n" );
   FreeArray( self->mapNodeDof2Eq );

   if (self->locationMatrix) {
      Journal_DPrintfL( self->debug, 2, "Freeing Full L.M. Array\n" );
      for( ii = 0; ii < self->nDomainEls; ii++ )
         FreeArray( self->locationMatrix[ii] );

      FreeArray( self->locationMatrix );
   }

   if( self->bcEqNums )
      Stg_Class_Delete( self->bcEqNums );

   if( self->ownedMap )
      Stg_Class_Delete( self->ownedMap );
}

/* Copy */

/** Stg_Class_Delete implementation. */
void _FeEquationNumber_Delete( void* feEquationNumber ) {
   FeEquationNumber* self = (FeEquationNumber*) feEquationNumber;

   // This calls the Destroy() function
   _Stg_Component_Delete( self );
}

/** Print implementation */
void _FeEquationNumber_Print( void* mesh, Stream* stream ) {
   FeEquationNumber* self = (FeEquationNumber*)mesh;

   /* General info */
   Journal_Printf( stream, "FeEquationNumber (ptr): %p\n", self );

   /* Print parent */
   _Stg_Class_Print( self, stream );

   /* Virtual info */
   Journal_Printf( stream,  "\t_build (func ptr): %p\n", self->_build );
   Journal_Printf( stream,  "\t_intialise (func ptr): %p\n", self->_initialise );

   /* FeEquationNumber info */
   /* Don't print dofs or bcs as these will be printed by FeVariable */

   if ( self->mapNodeDof2Eq ) {
      FeEquationNumber_PrintmapNodeDof2Eq( self, stream );
      FeEquationNumber_PrintLocationMatrix( self, stream );
   }
   else {
      Journal_Printf( stream, "\tmapNodeDof2Eq: (null)... not built yet\n" );
      Journal_Printf( stream, "\tlocationMatrix: (null)... not built yet\n" );
   }
}

void* _FeEquationNumber_Copy( void* feEquationNumber, void* dest, Bool deep, Name nameExt, PtrMap* ptrMap ) {

   abort();
}

void _FeEquationNumber_Build( void* feEquationNumber, void* data ) {
   FeEquationNumber* self = (FeEquationNumber*) feEquationNumber;

   assert(self);

   Journal_DPrintf( self->debug, "In %s:\n",  __func__ );
   Stream_IndentBranch( StgFEM_Debug );

   Stg_Component_Build( self->feMesh   , data, False );
   Stg_Component_Build( self->dofLayout, data, False );
   if ( self->linkedDofInfo ) Stg_Component_Build( self->dofLayout, data, False );
   if ( self->bcs )           Stg_Component_Build( self->bcs      , data, False );
	/* If we have new mesh topology information, do this differently. */
   /* if( self->feMesh->topo->domains && self->feMesh->topo->domains[MT_VERTEX] ) { */

   if( Mesh_HasExtension( self->feMesh, "vertexGrid" ) )
      FeEquationNumber_BuildWithDave( self );
   else
      FeEquationNumber_BuildWithTopology( self );

   /* If not removing BCs, construct a table of which equation numbers are actually BCs. */
   if( !self->removeBCs ) {
      FeMesh* mesh = self->feMesh;
      DofLayout* dofLayout = self->dofLayout;
      VariableCondition* bcs = self->bcs;
      int nDofs, varInd;
      int ii, jj;

      for( ii = 0; ii < FeMesh_GetNodeLocalSize( mesh ); ii++ ) {
         nDofs = dofLayout->dofCounts[ii];
         for( jj = 0; jj < nDofs; jj++ ) {
            varInd = dofLayout->varIndices[ii][jj];
            if( bcs && VariableCondition_IsCondition( bcs, ii, varInd ) ) {
               if( !STree_Has( self->bcEqNums, self->mapNodeDof2Eq[ii] + jj ) )
                  STree_Insert( self->bcEqNums, self->mapNodeDof2Eq[ii] + jj );
            }
         }
      }
   }

   if ( Stream_IsPrintableLevel( self->debug, 3 ) ) {
      FeEquationNumber_PrintmapNodeDof2Eq( self, self->debug );
   }

   Stream_UnIndentBranch( StgFEM_Debug );
}

/** Initialise implementation. Currently does nothing. */
void _FeEquationNumber_Initialise( void* feEquationNumber, void* data ) {
   FeEquationNumber* self = (FeEquationNumber*) feEquationNumber;

   Stg_Component_Initialise( self->feMesh   , data, False );
   Stg_Component_Initialise( self->dofLayout, data, False );
   if ( self->linkedDofInfo ) Stg_Component_Initialise( self->dofLayout, data, False );
   if ( self->bcs )           Stg_Component_Initialise( self->bcs      , data, False );
}

Index FeEquationNumber_CalculateActiveEqCountAtNode(
   void*			feEquationNumber,
   Node_DomainIndex	dNode_I,
   Dof_EquationNumber*	lowestActiveEqNumAtNodePtr )
{
   FeEquationNumber*	self = (FeEquationNumber*) feEquationNumber;
   Dof_Index		nodalDof_I = 0;
   Index			activeEqsAtCurrRowNode = 0;
   Dof_EquationNumber	currEqNum;
   Bool			foundLowest = False;

   for ( nodalDof_I = 0; nodalDof_I < self->dofLayout->dofCounts[dNode_I]; nodalDof_I++ ) {
      currEqNum = self->mapNodeDof2Eq[dNode_I][nodalDof_I];
      if ( currEqNum != -1 ) {
         activeEqsAtCurrRowNode++;
         if ( False == foundLowest ) {
            (*lowestActiveEqNumAtNodePtr) = currEqNum;
            foundLowest = True;
         }
      }
   }

   return activeEqsAtCurrRowNode;
}

void FeEquationNumber_BuildLocationMatrix( void* feEquationNumber ) {
   FeEquationNumber*	self = (FeEquationNumber*)feEquationNumber;
   FeMesh*			feMesh;
   unsigned		nDomainEls;
   unsigned*		nNodalDofs;
   unsigned		nElNodes;
   int*       elNodes;
   int**			dstArray;
   int***			locMat;
   IArray*			inc;
   unsigned		e_i, n_i, dof_i;

   assert( self );

   /* Don't build if already done. */
   if( self->locationMatrixBuilt ) {
      Journal_DPrintf( self->debugLM, "In %s: LM already built, so just returning.\n",  __func__ );
      Stream_UnIndentBranch( StgFEM_Debug );
      return;
   }

   inc = IArray_New();

   /* Shortcuts. */
   feMesh = self->feMesh;
   nDomainEls = FeMesh_GetElementDomainSize( feMesh );
   nNodalDofs = self->dofLayout->dofCounts;
   dstArray = self->mapNodeDof2Eq;

   /* Allocate for the location matrix. */
   locMat = AllocArray( int**, nDomainEls );
   for( e_i = 0; e_i < nDomainEls; e_i++ ) {
      FeMesh_GetElementNodes( feMesh, e_i, inc );
      nElNodes = IArray_GetSize( inc );
      elNodes = IArray_GetPtr( inc );
      locMat[e_i] = AllocArray( int*, nElNodes );
      for( n_i = 0; n_i < nElNodes; n_i++ )
         locMat[e_i][n_i] = AllocArray( int, nNodalDofs[elNodes[n_i]] );
   }

   /* Build location matrix. */
   for( e_i = 0; e_i < nDomainEls; e_i++ ) {
      FeMesh_GetElementNodes( feMesh, e_i, inc );
      nElNodes = IArray_GetSize( inc );
      elNodes = IArray_GetPtr( inc );
      for( n_i = 0; n_i < nElNodes; n_i++ ) {
         for( dof_i = 0; dof_i < nNodalDofs[elNodes[n_i]]; dof_i++ )
            locMat[e_i][n_i][dof_i] = dstArray[elNodes[n_i]][dof_i];
      }
   }

   Stg_Class_Delete( inc );

   /* Store result. */
   self->locationMatrix = locMat;
}

FeEquationNumber* _FeEquationNumber_Create( void* _self, Bool removeBCs ) {
  /* Build a n equation number for a given feVariable */
  FeVariable* feVar = Stg_CheckType( (FeVariable*)_self, FeVariable) ;
  FeEquationNumber* eqNum=NULL;
  Stg_Component_Build( feVar, NULL, False );
  Stg_Component_Initialise( feVar, NULL, False );
  eqNum = FeEquationNumber_New( "",
                                  feVar->context,
                                  feVar->feMesh,
                                  feVar->dofLayout,
                                  feVar->bcs,
                                  NULL );

  eqNum->removeBCs = removeBCs;
  memcpy( eqNum->periodic, feVar->periodic, 3*sizeof(Bool) );
  Stg_Component_Build( eqNum, NULL, False );
  Stg_Component_Initialise( eqNum, NULL, False );
  // feVar->eqNum = self->eqNum; // We do it for now
  return eqNum;
}

void FeEquationNumber_PrintmapNodeDof2Eq( void* feFeEquationNumber, Stream* stream ) {
   FeEquationNumber* self = (FeEquationNumber*) feFeEquationNumber;
   FeMesh*		feMesh = self->feMesh;
   MPI_Comm comm = Comm_GetMPIComm( Mesh_GetCommTopology( feMesh, MT_VERTEX ) );
   unsigned rank;
   Node_GlobalIndex gNode_I;
   Node_GlobalIndex nodeGlobalCount = FeMesh_GetNodeGlobalSize( feMesh );

   MPI_Comm_rank( comm, (int*)&rank );
   Journal_Printf( stream, "%d: *** Printing destination array ***\n", rank );

   for (gNode_I =0; gNode_I < nodeGlobalCount; gNode_I++) {
      Node_DomainIndex dNode_I;

      if ( !Mesh_GlobalToDomain( feMesh, MT_VERTEX, gNode_I, &dNode_I ) ) {
         Journal_Printf( stream, "\tmapNodeDof2Eq[(gnode)%2d]: on another proc\n", gNode_I);
      }
      else {
         Dof_Index currNodeNumDofs = self->dofLayout->dofCounts[ dNode_I ];
         Dof_Index nodeLocalDof_I;

         Journal_Printf( stream, "\tmapNodeDof2Eq[(gnode)%2d][(dof)0-%d]:",gNode_I, currNodeNumDofs );
         for( nodeLocalDof_I = 0; nodeLocalDof_I < currNodeNumDofs; nodeLocalDof_I++ ) {
            Journal_Printf( stream, "%3d, ", self->mapNodeDof2Eq[dNode_I][nodeLocalDof_I] );
         }
         Journal_Printf( stream, "\n" );
      }
   }
}

void FeEquationNumber_PrintmapNodeDof2EqBox( void* feFeEquationNumber, Stream* stream ) {
	FeEquationNumber*   self = (FeEquationNumber*) feFeEquationNumber;
	FeMesh*             feMesh = self->feMesh;
	MPI_Comm            comm = Comm_GetMPIComm( Mesh_GetCommTopology( feMesh, MT_VERTEX ) );
	unsigned            rank;
	Grid*               vGrid;
	int                 ijk[3];
	Element_LocalIndex  lNode_I;
	Node_GlobalIndex    gNode_I;
	Node_Index          sizeI, sizeJ, sizeK;
	Dof_Index           nDofs;
	Dof_Index           dof_I;

	MPI_Comm_rank( comm, (int*)&rank );
	Journal_Printf( stream, "%d: *** Printing destination array ***\n", rank );

	vGrid = *Mesh_GetExtension( feMesh, Grid**,  feMesh->vertGridId );
	nDofs = self->dofLayout->dofCounts[0];
	sizeI = vGrid->sizes[ I_AXIS ];
	sizeJ = vGrid->sizes[ J_AXIS ];
	sizeK = vGrid->sizes[ K_AXIS ] ? vGrid->sizes[ K_AXIS ] : 1;

	for ( ijk[2] = 0 ; ijk[2] < sizeK ; ijk[2]++ ) {
		if ( sizeK != 1 )
			Journal_Printf( stream, "\nk = %d\n", ijk[2] );
		for ( ijk[1] = sizeJ - 1 ; ijk[1] >= 0 ; ijk[1]-- ) {
			Journal_Printf( stream, "%2d - ", ijk[1] );
			for ( ijk[0] = 0 ; ijk[0] < sizeI ; ijk[0]++ ) {
				gNode_I = Grid_Project( vGrid, ijk );
				Journal_Printf( stream, "{ " );
				if ( Mesh_GlobalToDomain( feMesh, MT_VERTEX, gNode_I, &lNode_I ) ) {
					for ( dof_I = 0 ; dof_I < nDofs ; dof_I++ )
						Journal_Printf( stream, "%3d ", self->mapNodeDof2Eq[lNode_I][dof_I] );
				}
				else {
					for ( dof_I = 0 ; dof_I < nDofs ; dof_I++ )
						Journal_Printf( stream, " XX " );
				}
				Journal_Printf( stream, "}" );
			}
			Journal_Printf( stream, "\n" );
		}
		/* Bottom row */
		Journal_Printf( stream, "    " );
		for ( ijk[0] = 0 ; ijk[0] < sizeI ; ijk[0]++ ) {
			Journal_Printf( stream, "    %3d    ", ijk[0] );
			if ( nDofs == 3 )
				Journal_Printf( stream, "    " );
		}
		Journal_Printf( stream, "\n" );
	}
}

void FeEquationNumber_PrintLocationMatrix( void* feFeEquationNumber, Stream* stream ) {
   FeEquationNumber* self = (FeEquationNumber*) feFeEquationNumber;
   FeMesh*		feMesh = self->feMesh;
   MPI_Comm comm = Comm_GetMPIComm( Mesh_GetCommTopology( feMesh, MT_VERTEX ) );
   unsigned rank;
   Element_GlobalIndex gEl_I;
   unsigned nDims = Mesh_GetDimSize( feMesh );
   Element_GlobalIndex elementGlobalCount = FeMesh_GetElementGlobalSize( feMesh );
   unsigned nLocalEls = FeMesh_GetElementLocalSize( feMesh );

   Journal_Printf( stream, "%d: *** Printing location matrix ***\n", rank  );

   MPI_Comm_rank( comm, (int*)&rank );

   for (gEl_I =0; gEl_I < elementGlobalCount; gEl_I++ ) {
      Element_LocalIndex lEl_I;

      if ( !Mesh_GlobalToDomain( feMesh, nDims, gEl_I, &lEl_I ) || lEl_I >= nLocalEls ) {
         Journal_Printf( stream, "\tLM[(g/l el)%2d/XXX]: on another proc\n", gEl_I);
      }
      else {
         Node_LocalIndex numNodesAtElement;
         Node_LocalIndex elLocalNode_I;
         unsigned*	incNodes;
         IArray*		inc;

         inc = IArray_New();
         FeMesh_GetElementNodes( self->feMesh, lEl_I, inc );
         numNodesAtElement = IArray_GetSize( inc );
         incNodes = IArray_GetPtr( inc );

         Journal_Printf( stream, "\tLM[(g/l el)%2d/%2d][(enodes)0-%d]", gEl_I, lEl_I, numNodesAtElement);
         /* print the nodes and dofs */
         for ( elLocalNode_I = 0; elLocalNode_I < numNodesAtElement; elLocalNode_I++ ) {
            /* look up processor local node number. */
            Element_LocalIndex currNode = incNodes[elLocalNode_I == 2 ? 3 :
                                                   elLocalNode_I == 3 ? 2 :
                                                   elLocalNode_I == 6 ? 7 :
                                                   elLocalNode_I == 7 ? 6 :
                                                   elLocalNode_I];
            /* get the number of dofs at current node */
            Dof_Index currNodeNumDofs = self->dofLayout->dofCounts[ currNode ];
            Dof_Index nodeLocalDof_I;

            Journal_Printf( stream, "({%2d}", currNode );
            for( nodeLocalDof_I = 0; nodeLocalDof_I < currNodeNumDofs; nodeLocalDof_I++ ) {
               Journal_Printf( stream, "%3d,", self->mapNodeDof2Eq[currNode][nodeLocalDof_I] );
            }
            Journal_Printf( stream, "), " );
         }

         Journal_Printf( stream, "\n" );

         Stg_Class_Delete( inc );
      }

   }
}

Partition_Index FeEquationNumber_CalculateOwningProcessorOfEqNum( void* feEquationNumber, Dof_EquationNumber eqNum ) {
   FeEquationNumber* self = (FeEquationNumber*)feEquationNumber;
   /* Partition_Index ownerProc = (unsigned int)-1; */
   Comm*	comm = Mesh_GetCommTopology( self->feMesh, MT_VERTEX );
   MPI_Comm	mpiComm = Comm_GetMPIComm( comm );
   unsigned	nProcs;
   unsigned	p_i;

   MPI_Comm_size( mpiComm, (int*)&nProcs );
   for( p_i = 1; p_i < nProcs; p_i++ ) {
      if( eqNum < self->_lowestGlobalEqNums[p_i] )
         break;
   }

   return p_i - 1;

}

void FeEquationNumber_BuildWithTopology( FeEquationNumber* self ) {
   Stream*		stream;
   double		startTime, endTime, time, tmin, tmax;
   FeMesh*		feMesh;
   Sync*		sync;
   Comm*		comm;
   MPI_Comm		mpiComm;
   unsigned		rank, nProcs;
   unsigned		nDims;
   unsigned		nDomainNodes;
   unsigned		nLocalNodes;
   unsigned*		nNodalDofs;
   unsigned		nElNodes, *elNodes;
   int**		dstArray;
   int			*nLocMatDofs, ***locMat;
   unsigned		varInd;
   unsigned		curEqNum;
   unsigned		base;
   unsigned		subTotal;
   MPI_Status		status;
   unsigned		maxDofs;
   unsigned*		tuples;
   LinkedDofInfo*	links;
   unsigned		highest;
   IArray*		inc;
   unsigned             e_i, n_i, dof_i, s_i;
   int			ii;

   assert( self );

   inc = IArray_New();

//   stream = Journal_Register( Info_Type, (Name)self->type  );
//   Stream_SetPrintingRank( stream, 0 );
//
//   Journal_RPrintf( stream, "FeEquationNumber: '%s'\n", self->name );
//   Stream_Indent( stream );
//   Journal_RPrintf( stream, "Generating equation numbers...\n" );
//   Stream_Indent( stream );
//   if( self->removeBCs )
//      Journal_RPrintf( stream, "BCs set to be removed.\n" );
//   else
//      Journal_RPrintf( stream, "BCs will not be removed.\n" );

   startTime = MPI_Wtime();

   /* Shortcuts. */
   feMesh = self->feMesh;
   comm = Mesh_GetCommTopology( feMesh, MT_VERTEX );
   mpiComm = Comm_GetMPIComm( comm );
   MPI_Comm_size( mpiComm, (int*)&nProcs );
   MPI_Comm_rank( mpiComm, (int*)&rank );
   nDims = Mesh_GetDimSize( feMesh );
   nDomainNodes = FeMesh_GetNodeDomainSize( feMesh );
   self->nDomainEls = FeMesh_GetElementDomainSize( feMesh );
   nLocalNodes = FeMesh_GetNodeLocalSize( feMesh );
   nNodalDofs = self->dofLayout->dofCounts;
   links = self->linkedDofInfo;

   /* Allocate for destination array. */
   dstArray = Memory_Alloc_2DComplex( int, nDomainNodes, nNodalDofs,
                                      "FeEquationNumber::mapNodeDof2Eq" );

   /* If needed, allocate for linked equation numbers. */
   if( links ) {
      unsigned	s_i;

      links->eqNumsOfLinkedDofs = ReallocArray( links->eqNumsOfLinkedDofs, int, links->linkedDofSetsCount );
      for( s_i = 0; s_i < links->linkedDofSetsCount; s_i++ )
         links->eqNumsOfLinkedDofs[s_i] = -1;
   }

   /* Allocate for the location matrix. */
   nLocMatDofs = NULL;
   locMat = AllocArray( int**, self->nDomainEls );
   for( e_i = 0; e_i < self->nDomainEls; e_i++ ) {
      FeMesh_GetElementNodes( feMesh, e_i, inc );
      nElNodes = IArray_GetSize( inc );
      elNodes = IArray_GetPtr( inc );
      nLocMatDofs = ReallocArray( nLocMatDofs, int, nElNodes );
      for( n_i = 0; n_i < nElNodes; n_i++ )
         nLocMatDofs[n_i] = nNodalDofs[elNodes[n_i]];
      locMat[e_i] = AllocComplex2D( int, nElNodes, nLocMatDofs );
   }
   FreeArray( nLocMatDofs );

   /* Build initial destination array and store max dofs. */
   curEqNum = 0;
   maxDofs = 0;
   for( n_i = 0; n_i < nLocalNodes; n_i++ ) {
      if( nNodalDofs[n_i] > maxDofs )
         maxDofs = nNodalDofs[n_i];

      for( dof_i = 0; dof_i < nNodalDofs[n_i]; dof_i++ ) {
         varInd = self->dofLayout->varIndices[n_i][dof_i];
         if( !self->bcs || !VariableCondition_IsCondition( self->bcs, n_i, varInd ) ||
             !self->removeBCs )
         {
            if( links && links->linkedDofTbl[n_i][dof_i] != -1 ) {
               if( rank > 0 ) {
                  dstArray[n_i][dof_i] = -2;
                  continue;
               }
               if( links->eqNumsOfLinkedDofs[links->linkedDofTbl[n_i][dof_i]] == -1 )
                  links->eqNumsOfLinkedDofs[links->linkedDofTbl[n_i][dof_i]] = curEqNum++;
               dstArray[n_i][dof_i] = links->eqNumsOfLinkedDofs[links->linkedDofTbl[n_i][dof_i]];
            }
            else
               dstArray[n_i][dof_i] = curEqNum++;
         }
         else
            dstArray[n_i][dof_i] = -1;
      }
   }

   /* Order the equation numbers based on processor rank; cascade counts forward. */
   base = 0;
   if( rank > 0 )
      (void)MPI_Recv( &base, 1, MPI_UNSIGNED, rank - 1, 6669, mpiComm, &status );
   subTotal = base + curEqNum;
   if( rank < nProcs - 1 )
      (void)MPI_Send( &subTotal, 1, MPI_UNSIGNED, rank + 1, 6669, mpiComm );

   if( links ) {
      /* Reduce to find lowest linked DOFs. */
      for( s_i = 0; s_i < links->linkedDofSetsCount; s_i++ ) {
         if( links->eqNumsOfLinkedDofs[s_i] != -1 )
            links->eqNumsOfLinkedDofs[s_i] += base;
/*
  MPI_Allreduce( links->eqNumsOfLinkedDofs + s_i, &lowest, 1, MPI_UNSIGNED, MPI_MAX, mpiComm );
*/
         MPI_Allreduce( links->eqNumsOfLinkedDofs + s_i, &highest, 1, MPI_INT, MPI_MAX, mpiComm );
/*
  assert( (lowest == (unsigned)-1) ? lowest == highest : 1 );
*/
         links->eqNumsOfLinkedDofs[s_i] = highest;
      }
   }

   /* Modify existing destination array and dump to a tuple array. */
   tuples = AllocArray( unsigned, nDomainNodes * maxDofs );
   for( n_i = 0; n_i < nLocalNodes; n_i++ ) {
      for( dof_i = 0; dof_i < nNodalDofs[n_i]; dof_i++ ) {
         varInd = self->dofLayout->varIndices[n_i][dof_i];
         if( !self->bcs || !VariableCondition_IsCondition( self->bcs, n_i, varInd ) ||
             !self->removeBCs )
         {
            if( links && links->linkedDofTbl[n_i][dof_i] != -1 ) {
               highest = links->eqNumsOfLinkedDofs[links->linkedDofTbl[n_i][dof_i]];
               dstArray[n_i][dof_i] = highest;
            }
            else
               dstArray[n_i][dof_i] += base;
         }
         tuples[n_i * maxDofs + dof_i] = dstArray[n_i][dof_i];
      }
   }

   /* Update all other procs. */
   sync = Mesh_GetSync( feMesh, MT_VERTEX );
   Sync_SyncArray( sync, tuples, maxDofs * sizeof(unsigned),
                   tuples + nLocalNodes * maxDofs, maxDofs * sizeof(unsigned),
                   maxDofs * sizeof(unsigned) );

   /* Update destination array's domain indices. */
   for( n_i = nLocalNodes; n_i < nDomainNodes; n_i++ ) {
      for( dof_i = 0; dof_i < nNodalDofs[n_i]; dof_i++ ) {
         varInd = self->dofLayout->varIndices[n_i][dof_i];
         if( !self->bcs || !VariableCondition_IsCondition( self->bcs, n_i, varInd ) ||
             !self->removeBCs )
         {
            dstArray[n_i][dof_i] = tuples[n_i * maxDofs + dof_i];
         }
         else
            dstArray[n_i][dof_i] = -1;
      }
   }

   /* Destroy tuple array. */
   FreeArray( tuples );

   /* Build location matrix. */
   for( e_i = 0; e_i < self->nDomainEls; e_i++ ) {
      FeMesh_GetElementNodes( feMesh, e_i, inc );
      nElNodes = IArray_GetSize( inc );
      elNodes = IArray_GetPtr( inc );
      for( n_i = 0; n_i < nElNodes; n_i++ ) {
         for( dof_i = 0; dof_i < nNodalDofs[elNodes[n_i]]; dof_i++ )
            locMat[e_i][n_i][dof_i] = dstArray[elNodes[n_i]][dof_i];
      }
   }

   /* Store stuff on class. */
   self->mapNodeDof2Eq = dstArray;
   self->locationMatrix = locMat;
   self->locationMatrixBuilt = True;
   self->remappingActivated = False;
   self->localEqNumsOwnedCount = curEqNum;
   self->firstOwnedEqNum = base;
   self->lastOwnedEqNum = subTotal - 1;
   self->_lowestLocalEqNum = self->firstOwnedEqNum;

   /* Setup owned mapping. */
   STree_Clear( self->ownedMap );
   for( ii = self->firstOwnedEqNum; ii <= self->lastOwnedEqNum; ii++ ) {
      int val = ii - self->firstOwnedEqNum;
      STreeMap_Insert( self->ownedMap, &ii, &val );
   }

   /* Bcast global sum from highest rank. */
   if( rank == nProcs - 1 )
      self->globalSumUnconstrainedDofs = self->lastOwnedEqNum + 1;
   (void)MPI_Bcast( &self->globalSumUnconstrainedDofs, 1, MPI_UNSIGNED, nProcs - 1, mpiComm );

   /* Construct lowest global equation number list. */
   self->_lowestGlobalEqNums = AllocArray( int, nProcs );
   (void)MPI_Allgather( &self->firstOwnedEqNum, 1, MPI_UNSIGNED, self->_lowestGlobalEqNums, 1, MPI_UNSIGNED, mpiComm );

//   endTime = MPI_Wtime();

//   Journal_RPrintf( stream, "Assigned %d global equation numbers.\n", self->globalSumUnconstrainedDofs );
//   Journal_Printf( stream, "[%u] Assigned %d local equation numbers, within range %d to %d.\n",
//                   rank, self->lastOwnedEqNum - self->firstOwnedEqNum + 1, self->firstOwnedEqNum, self->lastOwnedEqNum + 1 );
//   Stream_UnIndent( stream );

//   time = endTime - startTime;
//   (void)MPI_Reduce( &time, &tmin, 1, MPI_DOUBLE, MPI_MIN, 0, mpiComm );
//   (void)MPI_Reduce( &time, &tmax, 1, MPI_DOUBLE, MPI_MAX, 0, mpiComm );
//   Journal_RPrintf( stream, "... Completed in %g [min] / %g [max] seconds.\n", tmin, tmax );
//   Stream_UnIndent( stream );

   Stg_Class_Delete( inc );
}

void FeEquationNumber_BuildWithDave( FeEquationNumber* self ) {
   int nLocals, *locals;
   Grid *vGrid;
   int varInd;
   int nEqNums, **dstArray;
   IArray *inc;
   int nDofs;
   int *periodic;
   int ***locMat;
   int nDims;
   int *elNodes;
   Comm *comm;
   MPI_Comm mpiComm;
   int nRanks, rank;
   Sync *sync;
   Bool isCond;
   int nPeriodicInds[3];
   int *periodicInds[3];
   int inds[3];
   Bool usePeriodic;
   int *tmpArray, nLocalEqNums;
   int lastOwnedEqNum, ind;
   STree *doneSet;
   int ii, jj, kk;

   comm = Mesh_GetCommTopology( self->feMesh, 0 );
   mpiComm = Comm_GetMPIComm( comm );
   MPI_Comm_size( mpiComm, &nRanks );
   MPI_Comm_rank( mpiComm, &rank );

   /* Setup an array containing global indices of all locally owned nodes. */
   nLocals = Mesh_GetLocalSize( self->feMesh, 0 );
   locals = AllocArray( int, nLocals );
   for( ii = 0; ii < nLocals; ii++ )
      locals[ii] = Mesh_DomainToGlobal( self->feMesh, 0, ii );

   /* Allocate for destination array. */
   nDofs = self->dofLayout->dofCounts[0];
   dstArray = AllocArray2D( int, Mesh_GetDomainSize( self->feMesh, 0 ), nDofs );

   /* Get the vertex grid extension and any periodicity. */
   nDims = Mesh_GetDimSize( self->feMesh );
   vGrid = *Mesh_GetExtension( self->feMesh, Grid**,  self->feMesh->vertGridId );
   periodic = Mesh_GetExtension( self->feMesh, int*, self->feMesh->periodicId );

   /* Fill destination array with initial values, setting dirichlet BCs as we go. */
   for( ii = 0; ii < nLocals; ii++ )
   {
/*
      Grid_Lift( vGrid, locals[ii], inds );
      usePeriodic = False;
      for( jj = 0; jj < nDims; jj++ )
      {
         if( (periodic[jj] || self->periodic[jj]) && (inds[jj] == 0 || inds[jj] == vGrid->sizes[jj] - 1) )
         {
            usePeriodic = True;
            break;
         }
      }
*/
      for( jj = 0; jj < nDofs; jj++ )
      {
         varInd = self->dofLayout->varIndices[ii][jj];
         if( self->bcs )
            isCond = VariableCondition_IsCondition( self->bcs, ii, varInd );
         else
            isCond = False;

         if( isCond && self->removeBCs )
            dstArray[ii][jj] = -1;
         else
            dstArray[ii][jj] = 0;
      }
   }

   /* Generate opposing indices for periodicity. */
   for( ii = 0; ii < nDims; ii++ )
   {
      nPeriodicInds[ii] = 0;
      periodicInds[ii] = NULL;
      if( periodic[ii] || self->periodic[ii] )
      {
         periodicInds[ii] = AllocArray( int, nLocals );
         for( jj = 0; jj < nLocals; jj++ )
         {
            Grid_Lift( vGrid, locals[jj], inds );
            if( inds[ii] != vGrid->sizes[ii] - 1 )
               continue;
            /*
            for( kk = 0; kk < nDofs; kk++ )
               if( dstArray[jj][kk] == -1 )
                  break;
               if( kk < nDofs )
                  continue;
            */
            periodicInds[ii][nPeriodicInds[ii]++] = locals[jj];
         }
      }
   }

   /* Call Dave's equation number generation routine. */
   if( nDims == 2 ) {
      GenerateEquationNumbering( vGrid->sizes[0], vGrid->sizes[1], 1,
				 nLocals, locals,
				 nDofs, Mesh_GetGlobalSize( self->feMesh, 0 ),
				 periodic[0] || self->periodic[0], periodic[1] || self->periodic[1], False,
				 nPeriodicInds[0], nPeriodicInds[1], 0,
				 periodicInds[0], periodicInds[1], NULL,
				 dstArray[0], &nEqNums );
   }
   else if( nDims == 3 ) {
      GenerateEquationNumbering( vGrid->sizes[0], vGrid->sizes[1], vGrid->sizes[2],
				 nLocals, locals,
				 nDofs, Mesh_GetGlobalSize( self->feMesh, 0 ),
				 periodic[0] || self->periodic[0], periodic[1] || self->periodic[1], periodic[2] || self->periodic[2],
				 nPeriodicInds[0], nPeriodicInds[1], nPeriodicInds[2],
				 periodicInds[0], periodicInds[1], periodicInds[2],
				 dstArray[0], &nEqNums );
   }
   else if( nDims == 1 ) {
      GenerateEquationNumbering( vGrid->sizes[0], 1, 1,
				 nLocals, locals,
				 nDofs, Mesh_GetGlobalSize( self->feMesh, 0 ),
				 periodic[0], False, False,
				 nPeriodicInds[0], 0, 0,
				 periodicInds[0], NULL, NULL,
				 dstArray[0], &nEqNums );
   }
   else
      abort();

   /* Free periodic arrays. */
   for( ii = 0; ii < nDims; ii++ )
   {
      if( periodicInds[ii] )
         FreeArray( periodicInds[ii] );
   }

   /* Setup owned mapping part 1. */
   STree_Clear( self->ownedMap );
   for( ii = 0; ii < nLocals; ii++ )
   {
      Grid_Lift( vGrid, locals[ii], inds );
      for( jj = 0; jj < nDims; jj++ )
      {
         if( (periodic[jj] || self->periodic[jj]) && inds[jj] == vGrid->sizes[jj] - 1 )
         {
            inds[jj] = 0;
            ind = Grid_Project( vGrid, inds );
            if( !FeMesh_NodeGlobalToDomain( self->feMesh, ind, &ind ) )
               break;
         }
      }
      if( jj < nDims )
         continue;
      for( jj = 0; jj < nDofs; jj++ )
      {
         if( dstArray[ii][jj] == -1 || STreeMap_HasKey( self->ownedMap, dstArray[ii] + jj ) )
            continue;
         STreeMap_Insert( self->ownedMap, dstArray[ii] + jj, &ii );
      }
   }

   /* Setup owned mapping. */
   tmpArray = AllocArray( int, nLocals * nDofs );
   memcpy( tmpArray, dstArray[0], nLocals * nDofs * sizeof(int) );
   qsort( tmpArray, nLocals * nDofs, sizeof(int), stgCmpInt );
   doneSet = STree_New();
   STree_SetItemSize( doneSet, sizeof(int) );
   STree_SetIntCallbacks( doneSet );
   for( nLocalEqNums = 0, ii = 0; ii < nLocals * nDofs; ii++ )
   {
      if( tmpArray[ii] != -1 && STreeMap_HasKey( self->ownedMap, tmpArray + ii ) && !STree_Has( doneSet, tmpArray + ii ) )
      {
         if( !nLocalEqNums )
            self->_lowestLocalEqNum = tmpArray[ii];
         *(int*)STreeMap_Map( self->ownedMap, tmpArray + ii ) = nLocalEqNums;
         STree_Insert( doneSet, tmpArray + ii );
         nLocalEqNums++;
      }
   }
   lastOwnedEqNum = -1; /* Don't need this anymore. */
   FreeArray( tmpArray );
   Stg_Class_Delete( doneSet );

   /* Transfer remote equation numbers. */
   sync = Mesh_GetSync( self->feMesh, 0 );
   Sync_SyncArray( sync, dstArray[0], nDofs * sizeof(int),
                   dstArray[0] + nLocals * nDofs, nDofs * sizeof(int),
                   nDofs * sizeof(int) );

   /* Allocate for location matrix. */
   /* first store nDomainEls for usage during destroy */
   self->nDomainEls = Mesh_GetDomainSize( self->feMesh, nDims );
   locMat = AllocArray( int**, self->nDomainEls );
   for( ii = 0; ii < self->nDomainEls; ii++ )
      locMat[ii] = AllocArray2D( int, FeMesh_GetElementNodeSize( self->feMesh, 0 ), nDofs );

   /* Fill in location matrix. */
   inc = IArray_New();
   for( ii = 0; ii < Mesh_GetDomainSize( self->feMesh, nDims ); ii++ )
   {
      FeMesh_GetElementNodes( self->feMesh, ii, inc );
      elNodes = IArray_GetPtr( inc );
      for( jj = 0; jj < FeMesh_GetElementNodeSize( self->feMesh, 0 ); jj++ )
      {
         for( kk = 0; kk < nDofs; kk++ )
         locMat[ii][jj][kk] = dstArray[elNodes[jj]][kk];
      }
   }
   Stg_Class_Delete( inc );

   /* Fill in our other weird values. */
   self->mapNodeDof2Eq = dstArray;
   self->locationMatrix = locMat;
   self->locationMatrixBuilt = True;
   self->remappingActivated = False;
   self->localEqNumsOwnedCount = nLocalEqNums;

   /* Bcast global sum from highest rank. */
   self->globalSumUnconstrainedDofs = nEqNums;

   /* Construct lowest global equation number list. */
   self->_lowestGlobalEqNums = AllocArray( int, nRanks );
   (void)MPI_Allgather( &self->_lowestLocalEqNum, 1, MPI_UNSIGNED,
                        self->_lowestGlobalEqNums, 1, MPI_UNSIGNED,
                        mpiComm );

   FreeArray( locals );

   /*
   printf( "%d: localEqNumsOwned = %d\n", rank, self->localEqNumsOwnedCount );
   printf( "%d: globalSumUnconstrainedDofs = %d\n", rank, self->globalSumUnconstrainedDofs );
   */
}

/*

Input:
  nlocal - number of locally ownded nodes
  g_node_id - global indices of nodes owned locally. Size nlocal
  dof - degrees of freedom per node
  nglobal - number of global nodes
  npx - number of consider to be periodic in x (local to this proc)
  npy - number of consider to be periodic in y (local to this proc)
  periodic_x_gnode_id - global indices of nodes (on this proc) which are on right hand side boundary
  periodic_y_gnode_id - global indices of nodes (on this proc) whipch are on top boundary
  eqnums - contains any dirichlet boundary conditions. Size nlocal*dof

Output;
  eqnums - contains full list of eqnums

Assumptions:
- Ordering eqnums[] = { (node_0,[dof_0,dof_1,..,dof_x]), (node_1,[dof_0,dof_1,..,dof_x]), ... }
- Any dirichlet set along a boundary deemed to be periodic will be clobbered.
- Processors may have duplicate nodes in the g_node_id[] list.
- A number in the corner is considered part of both boundaries (horiz and vert)
- If npx is not 0, then periodicity is assumed in x
- If npy is not 0, then periodicity is assumed in y
- Dofs constrained to be dirichlet must be marked with a negative number.
- Dofs are NOT split across processors.
- We can define a logical i,j,k ordering to uniquely identify nodes.

*/

PetscErrorCode _VecScatterBeginEnd( VecScatter vscat, Vec FROM, Vec TO, InsertMode addv,ScatterMode mode )
{
#if( (PETSC_VERSION_MAJOR==2) && (PETSC_VERSION_MINOR==3) && (PETSC_VERSION_SUBMINOR==2) )
	// 2.3.2 ordering of args
	VecScatterBegin( FROM, TO, addv, mode, vscat );
	VecScatterEnd( FROM, TO, addv, mode, vscat );
#else
	// 2.3.3 or 3.0.0
	VecScatterBegin( vscat, FROM, TO, addv, mode );
	VecScatterEnd( vscat, FROM, TO, addv, mode );
#endif

	PetscFunctionReturn(0);
}

int GenerateEquationNumbering(
		int NX, int NY, int NZ,
		int nlocal, int g_node_id[],
		int dof, int nglobal,
		PetscTruth periodic_x, PetscTruth periodic_y, PetscTruth periodic_z,
		int npx, int npy, int npz,
		int periodic_x_gnode_id[], int periodic_y_gnode_id[], int periodic_z_gnode_id[],
		int eqnums[], int *neqnums )
{
	PetscErrorCode ierr;
	PetscInt periodic_mask;
	Vec global_eqnum, g_ownership;
	PetscInt i;
	PetscMPIInt rank;        /* processor rank */
	PetscMPIInt size;        /* size of communicator */
	Vec local_ownership, local_eqnum;
	PetscInt *_g_node_id;
	IS is_gnode, is_eqnum;
	VecScatter vscat_ownership, vscat_eqnum;
	PetscScalar *_local_ownership, *_local_eqnum;
	PetscInt local_eqnum_count,global_eqnum_count;
	PetscScalar val[10];
	PetscInt d,idx[10];
	PetscInt *to_fetch,cnt,number_to_fetch;
	PetscInt eq_cnt;
	Vec offset_list;
	VecScatter vscat_offset;
	Vec seq_offset_list;
	PetscInt offset, inc;

	PetscInt spanx,spany,spanz,total;
	PetscInt loc;
	PetscReal max;
	PetscInt n_inserts;

	ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
	ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);

	if( dof>=10 ) {
		Stg_SETERRQ(PETSC_ERR_SUP, "Max allowable degrees of freedom per node = 10. Change static size" );
	}

	/*
	Claim locally owned nodes. Duplicate nodes on the interior will be resolved by the processor
	which inserts last.
	*/
	VecCreate( PETSC_COMM_WORLD, &g_ownership );
	VecSetSizes( g_ownership, PETSC_DECIDE, nglobal );
	VecSetFromOptions( g_ownership );

	for( i=0; i<nlocal; i++ ) {
		VecSetValue( g_ownership, g_node_id[i], rank, INSERT_VALUES );
	}
	VecAssemblyBegin(g_ownership);
	VecAssemblyEnd(g_ownership);

	/* Mask out the periodic boundaries. */
	periodic_mask = -6699.0;
	if (periodic_x_gnode_id!=NULL) {
		for( i=0; i<npx; i++ ) {
			VecSetValue( g_ownership, periodic_x_gnode_id[i], periodic_mask, INSERT_VALUES );
		}
	}
	if (periodic_y_gnode_id!=NULL) {
		for( i=0; i<npy; i++ ) {
			VecSetValue( g_ownership, periodic_y_gnode_id[i], periodic_mask, INSERT_VALUES );
		}
	}
	if (periodic_z_gnode_id!=NULL) {
		for( i=0; i<npz; i++ ) {
			VecSetValue( g_ownership, periodic_z_gnode_id[i], periodic_mask, INSERT_VALUES );
		}
	}
	VecAssemblyBegin(g_ownership);
	VecAssemblyEnd(g_ownership);

	/*
	PetscPrintf(PETSC_COMM_WORLD, "g_ownership \n");
	VecView( g_ownership, PETSC_VIEWER_STDOUT_WORLD );
	*/

	/* Get all locally owned nodes */
	VecCreate( PETSC_COMM_SELF, &local_ownership );
	VecSetSizes( local_ownership, PETSC_DECIDE, nlocal );
	VecSetFromOptions( local_ownership );

	PetscMalloc( sizeof(PetscInt)*nlocal, &_g_node_id);
	for( i=0; i<nlocal; i++ ) {
		_g_node_id[i] = g_node_id[i];
	}
	ISCreateGeneralWithArray( PETSC_COMM_WORLD, nlocal, _g_node_id, &is_gnode );
	VecScatterCreate( g_ownership, is_gnode, local_ownership, PETSC_NULL, &vscat_ownership );

	/* assign unique equation numbers */
	VecSet( local_ownership, -6699 );
	_VecScatterBeginEnd( vscat_ownership, g_ownership, local_ownership, INSERT_VALUES, SCATTER_FORWARD );

	/* Count instances of rank in the local_ownership vector */
	VecGetArray( local_ownership, &_local_ownership );
	local_eqnum_count = 0;
	for( i=0; i<nlocal; i++ ) {
		if( ((PetscInt)_local_ownership[i]) == rank ) {
			local_eqnum_count++;
		}
	}

	(void)MPI_Allreduce( &local_eqnum_count, &global_eqnum_count, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD );
	/* PetscPrintf( PETSC_COMM_SELF,
	    "[%d] number of local,global equations (without dofs) %d,%d \n", rank, local_eqnum_count, global_eqnum_count ); */
	/* check */
	spanx = NX;
	spany = NY;
	spanz = NZ;
	if( periodic_x==PETSC_TRUE ) {
		spanx--;
	}
	if( periodic_y==PETSC_TRUE ) {
		spany--;
	}
	if( periodic_z==PETSC_TRUE ) {
		spanz--;
	}
	total = spanx*spany*spanz;
	if( total!=global_eqnum_count ) {
		Stg_SETERRQ(PETSC_ERR_SUP, "Something stinks. Computed global size for nodes does not match expected" );
	}

	VecCreate( PETSC_COMM_WORLD, &global_eqnum );
	VecSetSizes( global_eqnum, PETSC_DECIDE, nglobal*dof );
	VecSetFromOptions( global_eqnum );
	VecSet( global_eqnum, 0.0 );

	/* Load existing eqnums in */
	for( i=0; i<nlocal; i++ ) {
		n_inserts = 0;
		for( d=0; d<dof; d++ ) {
/*
			idx[d] = -(g_node_id[i]*dof + d);
			val[d] = 0.0;
*/
			if( eqnums[ i*dof + d ] < 0.0 ) {
				idx[n_inserts] = g_node_id[i]*dof + d;
				val[n_inserts] = eqnums[ i*dof + d ];
				n_inserts++;
			}
		}
		/* only insert dirichlet bc's */
		VecSetValues( global_eqnum, n_inserts, idx, val, INSERT_VALUES );
	}
	VecAssemblyBegin(global_eqnum);
	VecAssemblyEnd(global_eqnum);

	/* Generate list of eqnums to get */
	PetscMalloc( sizeof(PetscInt)*nlocal*dof, &to_fetch );
	cnt = 0;
	for( i=0; i<nlocal; i++ ) {
		if( _local_ownership[i]==rank ) {
			for( d=0; d<dof; d++ ) {
				to_fetch[cnt] = g_node_id[i]*dof + d;
				cnt++;
			}
		}
	}
	number_to_fetch = cnt;

	VecCreate( PETSC_COMM_SELF, &local_eqnum );
	VecSetSizes( local_eqnum, PETSC_DECIDE, number_to_fetch);
	VecSetFromOptions( local_eqnum );

	ISCreateGeneralWithArray( PETSC_COMM_SELF, number_to_fetch, to_fetch, &is_eqnum );
	VecScatterCreate( global_eqnum, is_eqnum, local_eqnum, PETSC_NULL, &vscat_eqnum );

	VecSet( local_eqnum, -6699 );
	_VecScatterBeginEnd( vscat_eqnum, global_eqnum, local_eqnum, INSERT_VALUES, SCATTER_FORWARD );

	/* compute offset */
	/* count how many entries there are */
	VecGetArray( local_eqnum, &_local_eqnum );
	eq_cnt = 0;
	for( i=0; i<number_to_fetch; i++ ) {
		if( (PetscInt)_local_eqnum[i]==0 ) {
			eq_cnt++;
		}
	}
	VecRestoreArray( local_eqnum, &_local_eqnum );

    VecCreate( PETSC_COMM_WORLD, &offset_list );
    VecSetSizes( offset_list, PETSC_DECIDE, size );
    VecSetFromOptions( offset_list );
    VecSetValue( offset_list, rank, eq_cnt, INSERT_VALUES );
    VecAssemblyBegin(offset_list);
    VecAssemblyEnd(offset_list);
    /*
    PetscPrintf(PETSC_COMM_WORLD, "offset_list \n");
    VecView( offset_list, PETSC_VIEWER_STDOUT_WORLD );
    */

    VecScatterCreateToAll(offset_list,&vscat_offset,&seq_offset_list);
    _VecScatterBeginEnd( vscat_offset, offset_list, seq_offset_list, INSERT_VALUES, SCATTER_FORWARD );

    {
        PetscScalar *_seq_offset_list;

        VecGetArray( seq_offset_list, &_seq_offset_list );
        offset = 0;
        for (i=0; i<rank; i++) {
            offset+=_seq_offset_list[ i ];
        }
        VecRestoreArray( seq_offset_list, &_seq_offset_list );
    }
    Stg_VecScatterDestroy(&vscat_offset);
    Stg_VecDestroy(&offset_list);
    Stg_VecDestroy(&seq_offset_list);

//    PetscPrintf( PETSC_COMM_SELF, "[%d]: offset = %d \n", rank, offset );

	VecGetArray( local_eqnum, &_local_eqnum );
	inc = 0;
	for( i=0; i<number_to_fetch; i++ ) {
		if( (PetscInt)_local_eqnum[i]==0 ) {
			_local_eqnum[i] = offset+inc;
			inc++;
		}
	}
	VecRestoreArray( local_eqnum, &_local_eqnum );

	_VecScatterBeginEnd( vscat_eqnum, local_eqnum, global_eqnum, INSERT_VALUES, SCATTER_REVERSE );
	/*
	PetscPrintf(PETSC_COMM_WORLD, "global_eqnum \n");
	VecView( global_eqnum, PETSC_VIEWER_STDOUT_WORLD );
	*/

	/* For each periodic boundary, get the mapped nodes */

	if( periodic_x==PETSC_TRUE ) {
		VecScatter vscat_p;
		IS is_from;
		PetscInt *from, *to;
		Vec mapped;
		PetscScalar *_mapped;
		PetscInt c;

		PetscMalloc( sizeof(PetscInt)*npx*dof, &from );
		PetscMalloc( sizeof(PetscInt)*npx*dof, &to );

		c = 0;
		for( i=0; i<npx; i++ ) {
			PetscInt I,J,K,gid,from_gid;

			gid = periodic_x_gnode_id[i];
			K = gid/(NX*NY);
			J = (gid - K*(NX*NY))/NX;
			I = gid - K*(NX*NY) - J*NX;
			from_gid = (I-(NX-1)) + J*NX + K*(NX*NY);

			for( d=0; d<dof; d++ ) {
				to[c] = gid * dof + d;
				from[c] = from_gid * dof + d;
				c++;
			}
		}

		VecCreate( PETSC_COMM_SELF, &mapped );
		VecSetSizes( mapped, PETSC_DECIDE, npx*dof );
		VecSetFromOptions( mapped );

		ISCreateGeneralWithArray( PETSC_COMM_SELF, npx*dof, from, &is_from );
		VecScatterCreate( global_eqnum, is_from, mapped, PETSC_NULL, &vscat_p );

		_VecScatterBeginEnd( vscat_p, global_eqnum, mapped, INSERT_VALUES, SCATTER_FORWARD );
		if( npx>0 ) {
			VecGetArray( mapped, &_mapped );
			VecSetValues( global_eqnum, npx*dof, to, _mapped,  INSERT_VALUES );
			VecRestoreArray( mapped, &_mapped );
		}

		VecAssemblyBegin(global_eqnum);
		VecAssemblyEnd(global_eqnum);

		Stg_VecScatterDestroy(&vscat_p );
		Stg_ISDestroy(&is_from );
		Stg_VecDestroy(&mapped );
		PetscFree( from );
		PetscFree( to );
	}

	if( periodic_y==PETSC_TRUE ) {
		VecScatter vscat_p;
		IS is_from;
		PetscInt *from, *to;
		Vec mapped;
		PetscScalar *_mapped;
		PetscInt c;

		PetscMalloc( sizeof(PetscInt)*npy*dof, &from );
		PetscMalloc( sizeof(PetscInt)*npy*dof, &to );

		c = 0;
		for( i=0; i<npy; i++ ) {
			PetscInt I,J,K,gid,from_gid;

			gid = periodic_y_gnode_id[i];
			K = gid/(NX*NY);
			J = (gid - K*(NX*NY))/NX;
			I = gid - K*(NX*NY) - J*NX;
			from_gid = I + (J - (NY - 1))*NX + K*(NX*NY);

			for( d=0; d<dof; d++ ) {
				to[c] = gid * dof + d;
				from[c] = from_gid * dof + d;
				c++;
			}
		}

		VecCreate( PETSC_COMM_SELF, &mapped );
		VecSetSizes( mapped, PETSC_DECIDE, npy*dof );
		VecSetFromOptions( mapped );

		ISCreateGeneralWithArray( PETSC_COMM_SELF, npy*dof, from, &is_from );
		VecScatterCreate( global_eqnum, is_from, mapped, PETSC_NULL, &vscat_p );

		_VecScatterBeginEnd( vscat_p, global_eqnum, mapped, INSERT_VALUES, SCATTER_FORWARD );
		if( npy>0 ) {
			VecGetArray( mapped, &_mapped );
			VecSetValues( global_eqnum, npy*dof, to, _mapped,  INSERT_VALUES );
			VecRestoreArray( mapped, &_mapped );
		}

		VecAssemblyBegin(global_eqnum);
		VecAssemblyEnd(global_eqnum);

		Stg_VecScatterDestroy(&vscat_p );
		Stg_ISDestroy(&is_from );
		Stg_VecDestroy(&mapped );
		PetscFree( from );
		PetscFree( to );
	}

	if( periodic_z==PETSC_TRUE ) {
		VecScatter vscat_p;
		IS is_from;
		PetscInt *from, *to;
		Vec mapped;
		PetscScalar *_mapped;
		PetscInt c;

		PetscMalloc( sizeof(PetscInt)*npz*dof, &from );
		PetscMalloc( sizeof(PetscInt)*npz*dof, &to );

		c = 0;
		for( i=0; i<npz; i++ ) {
			PetscInt I,J,K,gid,from_gid;

			gid = periodic_z_gnode_id[i];
			K = gid/(NX*NY);
			J = (gid - K*(NX*NY))/NX;
			I = gid - K*(NX*NY) - J*NX;
			from_gid = I + J*NX + (K - (NZ-1))*(NX*NY);

			for( d=0; d<dof; d++ ) {
				to[c] = gid * dof + d;
				from[c] = from_gid * dof + d;
				c++;
			}
		}

		VecCreate( PETSC_COMM_SELF, &mapped );
		VecSetSizes( mapped, PETSC_DECIDE, npz*dof );
		VecSetFromOptions( mapped );

		ISCreateGeneralWithArray( PETSC_COMM_SELF, npz*dof, from, &is_from );
		VecScatterCreate( global_eqnum, is_from, mapped, PETSC_NULL, &vscat_p );

		_VecScatterBeginEnd( vscat_p, global_eqnum, mapped, INSERT_VALUES, SCATTER_FORWARD );
		if( npz>0 ) {
			VecGetArray( mapped, &_mapped );
			VecSetValues( global_eqnum, npz*dof, to, _mapped,  INSERT_VALUES );
			VecRestoreArray( mapped, &_mapped );
		}

		VecAssemblyBegin(global_eqnum);
		VecAssemblyEnd(global_eqnum);

		Stg_VecScatterDestroy(&vscat_p );
		Stg_ISDestroy(&is_from );
		Stg_VecDestroy(&mapped );
		PetscFree( from );
		PetscFree( to );
	}

	/*
	PetscPrintf(PETSC_COMM_WORLD, "global_eqnum following periodic \n");
	VecView( global_eqnum, PETSC_VIEWER_STDOUT_WORLD );
	*/

	/* Finally, scatter stuff from global_eqnums into MY array */
	{
		IS is_all_my_eqnums;
		PetscInt *all_my_eqnum_index;
		PetscInt _I,_D,CNT;
		Vec all_my_eqnums;
		PetscScalar *_all_my_eqnums;
		VecScatter vscat_p;

		PetscMalloc( sizeof(PetscInt)*dof*nlocal, &all_my_eqnum_index );

		CNT = 0;
		for( _I=0; _I<nlocal; _I++ ) {
			for( _D=0; _D<dof; _D++ ) {
				all_my_eqnum_index[CNT] = g_node_id[_I]*dof + _D;
				CNT++;
			}
		}

		ISCreateGeneralWithArray( PETSC_COMM_SELF, nlocal*dof, all_my_eqnum_index, &is_all_my_eqnums );
		VecCreate( PETSC_COMM_SELF, &all_my_eqnums );
		VecSetSizes( all_my_eqnums, PETSC_DECIDE, nlocal*dof );
		VecSetFromOptions( all_my_eqnums );
		VecScatterCreate( global_eqnum, is_all_my_eqnums, all_my_eqnums, PETSC_NULL, &vscat_p );
		_VecScatterBeginEnd( vscat_p, global_eqnum, all_my_eqnums, INSERT_VALUES, SCATTER_FORWARD );
		VecGetArray( all_my_eqnums, &_all_my_eqnums );

		for( i=0; i<nlocal*dof; i++ ) {
			eqnums[i] = (int)_all_my_eqnums[i];
		}
		VecRestoreArray( all_my_eqnums, &_all_my_eqnums );

		Stg_VecScatterDestroy(&vscat_p );
		Stg_VecDestroy(&all_my_eqnums );
		Stg_ISDestroy(&is_all_my_eqnums );
		PetscFree( all_my_eqnum_index );
	}

	VecMax( global_eqnum, &loc, &max );
	*neqnums = (int)max;
	(*neqnums)++;

	/* tidy up */
	VecRestoreArray( local_ownership, &_local_ownership );

	Stg_VecDestroy(&g_ownership );
	Stg_VecDestroy(&local_ownership );
	PetscFree( _g_node_id );
	Stg_ISDestroy(&is_gnode );
	Stg_VecScatterDestroy(&vscat_ownership );

	Stg_VecDestroy(&global_eqnum );
	Stg_VecDestroy(&local_eqnum );
	PetscFree( to_fetch );
	Stg_ISDestroy(&is_eqnum );
	Stg_VecScatterDestroy(&vscat_eqnum );

	return 0;
}


text:
/**
 * A metropolis hasting sampler for bremsstrahlung and pair production cross sections
 *
 * copyright  (C) 2007
 * the icecube collaboration
 *
 * $Id: I3MetropolisHastings.h 143040 2016-03-11 17:35:25Z nega $
 *
 * @version $Revision: 143040 $
 * @date $LastChangedDate: 2016-03-11 10:35:25 -0700 (Fri, 11 Mar 2016) $
 * @author Bernhard Voigt <bernhard.voigt@desy.de>   $LastChangedBy$
 */

#ifndef I3METROPOLISHASTINGS_h_
#define I3METROPOLISHASTINGS_h_

// icetray includes
#include "icetray/I3Logging.h"
#include "cmc/I3CascadeMCCommon.h"

// gsl includes
#include <gsl/gsl_rng.h>

// c++ includes
#include <string>

/**
 * @brief This class implements a Metropolis Hastings sampler
 *
 * It uses a pdf and a proposal distribution to efficiently
 * produce random samples from the pdf
 *
 * The proposal distribution used is a beta distribution, which
 * matches almost the shape of the bremsstrahlung and pair production
 * cross sections, using the right parameters. See I3CascadeSimulation
 * for the usage.
 *
 * For details on the algorithm, take a look at :
 * http://en.wikipedia.org/wiki/Metropolis-Hastings_algorithm
 */
class I3MetropolisHastings : public I3CascadeMCCommon::Sampler {

  // set logging name for icetray
  SET_LOGGER("I3MetropolisHastings");

 public:

  /**
   * Default Constructor
   */
  I3MetropolisHastings() {}

  /**
   * Constructor
   *
   * @param rng      a gsl_random instance
   * @param pdf      a function with two parameters, ie. the energy and
   *                 fractional energy of the cross section
   * @param proposal a gsl random instance sampling from the proposal
   *                 distribution which is a beta distribution
   * @param proposalDistribution a beta function in x
   * @param name     name of this instance for logging purposes
   *
   * The proposal gsl random instance as well as the proposalDistribution are implemented
   * in a struct beta_dist in the header file I3CascadeSimulation.h
   */
  I3MetropolisHastings(I3RandomServicePtr rng,
                       double (*pdf)(const double&, const double&),
                       double (*proposal)(const gsl_rng*),
                       double (*proposalDistribution)(const double&),
                       std::string name);

  /**
   * Get the total sample loop count
   */
  unsigned long GetCounter() {
    return counter_;
  }

  /**
   * Reset the total sample loop count
   */
  void ResetCounter() {
    counter_ = 0;
  }

  /**
   * Get the number of returned samples
   */
  unsigned long GetYield() {
    return yield_;
  }

  /**
   * Reset the number of returned samples
   */
  void ResetYield() {
    yield_ = 0;
  }

  /**
   * Samples n times starting at start position
   *
   * Used to forget the initial state
   *
   * @param energy parameter of underlying pdf
   * @param start starting position
   * @param n number of samples
   */
  void BurnIn(double energy, double start, unsigned int n=100);

  /**
   * Samples the underlying pdf in the given range
   *
   * @param energy parameter to the underlying pdf
   * @param lower edge of sample range
   * @param higher edge of sample range
   * @return random number
   */
  double Sample(double energy, double lower=0, double higher=1);

 private:

  /**
   * @brief pdf, ie a cross section with energy and fractional energy parameters
   */
  double (*pdf_)(const double&, const double&);

  /**
   * @brief a gls random instance drawing samples from the proposal distribution
   */
  double (*proposal_)(const gsl_rng*);

  /**
   * @brief a proposal distribution, in this case a beta distribution in x
   */
  double (*proposalDistribution_)(const double&);

  /**
   * @brief loop counter
   */
  unsigned long counter_;

  /**
   * @brief counts how many successful attemps have been made to draw a new random number
   */
  unsigned long yield_;

  /**
   * @brief name of this instance, for logging purpose
   */
  std::string name_;

  /**
   * @brief last drawn sample
   */
  double x_;

  /**
   * @brief last value of pdf(x)
   */
  double y_;

  /**
   * @brief last value of proposalDist(x)
   */
  double propDistY_;

};

#endif

ext: h


text:
// Copyright (c) 2021 Stig Rune Sellevag
//
// This file is distributed under the MIT License. See the accompanying file
// LICENSE.txt or http://www.opensource.org/licenses/mit-license.php for terms
// and conditions.

#ifndef SCILIB_LINALG_BLAS1_VECTOR_NORM2_H
#define SCILIB_LINALG_BLAS1_VECTOR_NORM2_H

#ifdef USE_MKL
#include <mkl.h>
#else
#include <cblas.h>
#endif

#include <experimental/mdspan>

namespace Sci {
namespace Linalg {

namespace stdex = std::experimental;

template <stdex::extents<>::size_type ext_x, class Layout_x, class Accessor_x>
inline double
norm2(stdex::mdspan<double, stdex::extents<ext_x>, Layout_x, Accessor_x> x)
{
    const BLAS_INT n = static_cast<BLAS_INT>(x.size());
    const BLAS_INT incx = static_cast<BLAS_INT>(x.stride(0));

    return cblas_dnrm2(n, x.data(), incx);
}

template <stdex::extents<>::size_type ext_x, class Layout_x, class Accessor_x>
inline double norm2(
    stdex::mdspan<const double, stdex::extents<ext_x>, Layout_x, Accessor_x> x)
{
    const BLAS_INT n = static_cast<BLAS_INT>(x.size());
    const BLAS_INT incx = static_cast<BLAS_INT>(x.stride(0));

    return cblas_dnrm2(n, x.data(), incx);
}

template <class Layout, class Allocator>
inline double norm2(const Sci::Vector<double, Layout, Allocator>& x)
{
    return norm2(x.view());
}

} // namespace Linalg
} // namespace Sci

#endif // SCILIB_LINALG_BLAS1_VECTOR_NORM2_H

max_stars_repo_name: stigrs/scilib

max_issues_repo_name: stigrs/scilib

max_forks_repo_name: stigrs/scilib


text:
#ifndef COMMAND_LINE_INCLUDE
#define COMMAND_LINE_INCLUDE

#include <string>
#include <vector>

#include <gsl/string_span>

#include "config/fleetingOptionsInterface.h"
#include "config/variablesMap.h"
#include "core/task.h"

namespace execHelper {
namespace plugins {
using CommandLineArg = std::string;
using CommandLineArgs = std::vector<CommandLineArg>;

static const gsl::czstring<> COMMAND_LINE_KEY = "command-line";

/**
 * \brief Extends the functionality to include the _command-line_ config parameter and processes this parameter
 */
struct CommandLine {
    /*! @copydoc AddEnvironment::getVariables(config::VariablesMap&, const config::FleetingOptionsInterface&)
     */
    static void
    getVariables(config::VariablesMap& variables,
                 const config::FleetingOptionsInterface& options) noexcept;

    /*! @copydoc JobsLong::apply(core::Task&, const config::VariablesMap&)
     */
    inline static void apply(core::Task& task,
                             const config::VariablesMap& variables) noexcept {
        task.append(*(variables.get<CommandLineArgs>(COMMAND_LINE_KEY)));
    }
};

} // namespace plugins
} // namespace execHelper

#endif /* COMMAND_LINE_INCLUDE */

max_stars_repo_name: exec-helper/source

max_issues_repo_name: exec-helper/source

max_forks_repo_name: exec-helper/source


text:
#pragma once

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <cblas.h>
#include <omp.h>
#include "matrix.h"

double l1_distance(double *vec_1, double *vec_2, int N);
void naive_ranking(double *matrix, double *result, int N);
void pagerank(double *matrix, double *result, int N,
              double tolerance, double damping_factor, int max_iter);


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>

#include <omp.h>

#include "util.h"

int main(int argc, char *argv[]) {
    char* input_fileName1 = argv[1];
    char* input_fileName2 = argv[2];
    int N_doc_bg = atoi(argv[3]);
    int N_kw = atoi(argv[4]);
    int N_obs = atoi(argv[5]);
    int N_iter = atoi(argv[6]);
    char* output_fileName = argv[7];
    int N_doc = 480000 / 2;

    int* matrix = (int*) malloc(sizeof(int) * N_obs * N_obs);
    int* matrix_bg = (int*) malloc(sizeof(int) * N_kw * N_kw);
    int* matrix_padded = (int*) malloc(sizeof(int) * N_obs * N_obs);
    int* true_index = (int*) malloc(sizeof(int) * N_kw);
    int* permutation = (int*) malloc(sizeof(int) * N_obs);
    gsl_matrix* matrix_obs;

    for (int round = 0; round < 10; round++)
    {
        char input_fileName1_extend[40];
        char input_fileName2_extend[40];
        sprintf(input_fileName1_extend, "%s%d", input_fileName1, round);
        sprintf(input_fileName2_extend, "%s%d", input_fileName2, round);

        // Setup
        struct timeval tv1,tv2;
        gettimeofday(&tv1, NULL);
        read_matrix(&true_index, &matrix_bg, 1.0*N_doc/N_doc_bg, N_kw, input_fileName2_extend);
        read_matrix(&true_index, &matrix, 1.0, N_obs, input_fileName1_extend);
        gettimeofday(&tv2, NULL);
        printf("Reading done: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));
        fflush(stdout);

        // compute the maximum frequency of a keyword
        int freq_max = 0;
        for (int ii = 0; ii < N_kw; ii++)
            if (matrix_bg[ii*N_kw+ii] > freq_max)
                freq_max = matrix_bg[ii*N_kw+ii];
        for (int ii = 0; ii < N_obs; ii++)
            if (matrix[ii*N_obs+ii] > freq_max)
                freq_max = matrix[ii*N_obs+ii];
        printf("Max frequency: %d\n", freq_max);

        for (int iter = 0; iter < 10; iter++)
        {
            printf("Run %d\n", iter);
            matrix_obs = gsl_matrix_alloc(N_obs, N_obs);

            gettimeofday(&tv1, NULL);
            pad_matrix(&matrix_padded, &matrix, N_obs, N_doc, freq_max);
            gettimeofday(&tv2, NULL);
            printf("Padding done: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));
            fflush(stdout);

            gettimeofday(&tv1, NULL);
            observe_matrix(matrix_obs, &matrix_padded, N_obs);
            gettimeofday(&tv2, NULL);
            printf("Observed matrix generated: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));
            fflush(stdout);

            // Permute observed matrix randomly and attack
            gettimeofday(&tv1, NULL);
            attack(matrix_obs, &matrix_bg, &permutation, N_kw, N_obs, N_doc, freq_max, N_iter);
            gettimeofday(&tv2, NULL);
            printf("Main attack done: %d.\n", (tv2.tv_sec - tv1.tv_sec));
            fflush(stdout);

            char output_fileName_full[40];
            sprintf(output_fileName_full, "%s%d-%d", output_fileName, round, iter);
            print_result(output_fileName_full, &permutation, &true_index, N_obs);
            //sprintf(output_fileName_full, "%s%d-%d-full", output_fileName, round, iter);
            //print_full_result(output_fileName_full, &permutation, &true_index, N_obs);
        }
    }

    free(matrix);
    free(matrix_padded);
    gsl_matrix_free(matrix_obs);
    return(0);
}

double log_score(int idx1, int idx2, gsl_matrix* matrix_obs, int** matrix, int** permutation, int N_kw, int N_doc, int freq_max)
{
    if (idx1 == idx2)
        return(0.0);

    int idx1_m = (*permutation)[idx1];
    int idx2_m = (*permutation)[idx2];

    int n1 = 2*freq_max - (*matrix)[idx1_m*N_kw + idx1_m];
    int n2 = 2*freq_max - (*matrix)[idx2_m*N_kw + idx2_m];

    /*
    int N_upper = 3 * sqrt((*matrix)[idx1_m*N_kw + idx2_m] * (N_doc - (*matrix)[idx1_m*N_kw + idx2_m]) / (double) N_doc);
    int N_lower = (*matrix)[idx1_m*N_kw + idx2_m] - N_upper;
    N_upper = (*matrix)[idx1_m*N_kw + idx2_m] + N_upper;
    int count_total = (int) gsl_matrix_get(matrix_obs, idx1, idx2);

    double score = 0;
    double prob = (*matrix)[idx1_m*N_kw + idx2_m] / (double) N_doc;
    for (int kk = N_lower; kk < N_upper+1; kk+=20)
        score +=  gsl_ran_binomial_pdf(kk, prob, N_doc) * gsl_ran_hypergeometric_pdf(count_total - kk, n1, 2*N_doc - n1, n2);
    */

    double mean = 1.0 * freq_max / N_doc * (n1 + n2);
    double var = 1.0 * (*matrix)[idx1_m*N_kw + idx2_m] / N_doc * (N_doc - (*matrix)[idx1_m*N_kw + idx2_m]);
    var += 1.0 * n1 / N_doc * freq_max * (2 * N_doc - 2*freq_max) / N_doc;
    var += 1.0 * n2 / N_doc * freq_max * (2 * N_doc - 2*freq_max) / N_doc;
    int count_total = (int) gsl_matrix_get(matrix_obs, idx1, idx2);

    int N_upper = 3 * sqrt((*matrix)[idx1_m*N_kw + idx2_m] * (N_doc - (*matrix)[idx1_m*N_kw + idx2_m]) / (double) N_doc);
    int N_lower = (*matrix)[idx1_m*N_kw + idx2_m] - N_upper;
    N_upper = (*matrix)[idx1_m*N_kw + idx2_m] + N_upper;

    mean += (*matrix)[idx1_m*N_kw + idx2_m];
    var += (*matrix)[idx1_m*N_kw + idx2_m] / N_doc * (N_doc - (*matrix)[idx1_m*N_kw + idx2_m]);

    double score = gsl_ran_gaussian_pdf(count_total - mean, sqrt(var));

    if (score == 0)
        return(-500.0);
    return(log(score));
}

void attack(gsl_matrix* matrix_obs, int** matrix, int** permutation, int N_kw, int N_obs, int N_doc, int freq_max, int N_iter)
{
    // Initialise data structures
    double* score_matrix = (double*) malloc(sizeof(double) * N_obs * N_obs);
    double* score_row1 = (double*) malloc(sizeof(double) * N_obs);
    double* score_row2 = (double*) malloc(sizeof(double) * N_obs);

    int* permutation_tmp = (int*) malloc(sizeof(int) * N_obs);
    int* permutation_inv = (int*) malloc(sizeof(int) * N_kw);

    // Initialise permutations
    for (int ii = 0; ii < N_obs; ii++)
        (*permutation)[ii] = ii;
    for (int ii = 0; ii < N_obs; ii++)
        permutation_tmp[ii] = ii;
    for (int ii = 0; ii < N_kw; ii++)
        permutation_inv[ii] = -1;
    for (int ii = 0; ii < N_obs; ii++)
        permutation_inv[permutation_tmp[ii]] = ii;

    // Initialising RNG
    const gsl_rng_type * T;
    gsl_rng * r;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    struct timeval tv1,tv2;
    gettimeofday(&tv1, NULL);

    // Compute initial score
    #pragma omp parallel for shared(score_matrix, matrix_obs, matrix)
    for (int ii = 0; ii < N_obs * N_obs; ii++)
        score_matrix[ii] = log_score((int) (ii / N_obs), ii % N_obs, matrix_obs, matrix, permutation, N_kw, N_doc, freq_max);
    gettimeofday(&tv2, NULL);
    printf("Initial score computed: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));

    // Iterations of simulated annealing
    double temp = (double) N_kw;
    int N_stuck = 0;
    for (int iter = 0; iter < N_iter; iter++)
    {
        /* Status code */
        if (iter % (N_iter / 10) == 0)
        {
            gettimeofday(&tv1, NULL);
            printf("Iteration: %d, %d, %d.\n", iter, N_stuck, (int) (tv1.tv_sec - tv2.tv_sec));
            fflush(stdout);
            gettimeofday(&tv2, NULL);
        }

	    // used to be 20k
        if (N_stuck >= 80000)
            iter = N_iter;

        /* Main code */
        int idx1, idx2;
        permutation_generation(&idx1, &idx2, &permutation_tmp, permutation, &permutation_inv, N_kw, N_obs);

        int ii = 0;
        #pragma omp parallel for shared(score_row1)
        for (ii = 0; ii < N_obs; ii++)
            score_row1[ii] = log_score(idx1, ii, matrix_obs, matrix, &permutation_tmp, N_kw, N_doc, freq_max);

        if (idx2 >= 0)
            #pragma omp parallel for shared(score_row2)
            for (ii = 0; ii < N_obs; ii++)
                score_row2[ii] = log_score(idx2, ii, matrix_obs, matrix, &permutation_tmp, N_kw, N_doc, freq_max);

        double score_diff = 0;
        for (int ii = 0; ii < N_obs; ii++)
            score_diff += score_row1[ii];
        for (int ii = 0; ii < N_obs; ii++)
            score_diff -= score_matrix[idx1*N_obs + ii];
        if (idx2 >= 0)
        {
            for (int ii = 0; ii < N_obs; ii++)
                score_diff += score_row2[ii];
            for (int ii = 0; ii < N_obs; ii++)
                score_diff -= score_matrix[idx2*N_obs + ii];
        }

        // compute difference in score, with exponentiation
        score_diff = score_diff / temp;

        if (score_diff < -40)
            score_diff = 0;
        else if (score_diff > 0)
            score_diff = 1.01;
        else
            score_diff = exp(score_diff);

        if (gsl_ran_flat(r, 0, 1) < score_diff)
        {
            // Update the scores
            for (int ii = 0; ii < N_obs; ii++)
                score_matrix[idx1*N_obs + ii] = score_row1[ii];
            for (int ii = 0; ii < N_obs; ii++)
                score_matrix[ii*N_obs + idx1] = score_row1[ii];
            if (idx2 >= 0)
            {
                for (int ii = 0; ii < N_obs; ii++)
                    score_matrix[idx2*N_obs + ii] = score_row2[ii];
                for (int ii = 0; ii < N_obs; ii++)
                    score_matrix[ii*N_obs + idx2] = score_row2[ii];
            }

            // Update the permutation
            permutation_inv[(*permutation)[idx1]] = -1;
            (*permutation)[idx1] = permutation_tmp[idx1];
            permutation_inv[permutation_tmp[idx1]] = idx1;
            if (idx2 >= 0)
            {
                (*permutation)[idx2] = permutation_tmp[idx2];
                permutation_inv[permutation_tmp[idx2]] = idx2;
            }
            N_stuck = 0;
        }
        else
        {
            // Update the permutation
            permutation_tmp[idx1] = (*permutation)[idx1];
            if (idx2 >= 0)
                permutation_tmp[idx2] = (*permutation)[idx2];
            N_stuck += 1;
        }

        temp *= 0.995;
    }

    free(score_matrix);
    free(score_row1);
    free(score_row2);
    gsl_rng_free(r);
}

void print_result(char* output_fileName, int** permutation, int** true_index, int N_obs)
{
    FILE* fp = fopen(output_fileName, "w");
    int count = 0;
    for (int ii = 0; ii < N_obs; ii++)
        if ((*permutation)[ii] == (*true_index)[ii])
            count++;

    fprintf(fp, "%d\n", count);
    fclose(fp);
    printf("Success: %d/%d.\n", count, N_obs);
}

void print_full_result(char* output_fileName, int** permutation, int** true_index, int N_obs)
{
    FILE* fp = fopen(output_fileName, "w");
    int count = 0;
    for (int ii = 0; ii < N_obs; ii++)
        if ((*permutation)[ii] == (*true_index)[ii])
            count++;

    printf("Success: %d/%d.\n", count, N_obs);
    fprintf(fp, "%d\n", count);
    for (int ii = 0; ii < N_obs; ii++)
        fprintf(fp, "%d,%d\n", (*permutation)[ii], (*true_index)[ii]);

    fclose(fp);
}

ext: c

max_stars_repo_name: RethinkingSSE/Attacks-on-SSE

max_issues_repo_name: RethinkingSSE/Attacks-on-SSE

max_forks_repo_name: RethinkingSSE/Attacks-on-SSE


text:
/*==========================================================
 * Multinomial_gsl.c
 *
 * Return multinomially distributed integers based on input probability vector
 *
 *========================================================*/
/*  */

#include "mex.h"
#include "matrix.h"
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h> /* For EXIT_FAILURE, EXIT_SUCCESS */
#include "mat.h"
// #include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_randist.h>
#include <time.h>
#define pi acos(-1.0)

int rlxd_seed(){
    int seed;
    FILE *file_seeds;
    file_seeds=fopen("/dev/urandom","rb"); // /dev/urandom is designed so that each time it is opened it produces a different random seed
    fread(&seed, sizeof(int), 1, file_seeds);
    fclose(file_seeds);

    return abs(seed);
}

void Multinomial(unsigned int *n, double N, double *x, mwSize K)
{

//     Multinomial(n,N,x,K);

    //    /* calculate eigenvalues and eigenvectors */
    //mexCallMATLAB(2, lhs, 1, &x, "eig");

//    printf("N =%f\n",N);

//    printf("K =%d\n",K);

    //Initialisation for using gls multinomial function

    const gsl_rng_type *Type;         //initialise a constant Type which holds type of random number generator
    gsl_rng *r;                              //Initialise pointer to r which is an "instance" of random number generator

    //Set seed for random numbers from current time
//     srand(time(0));
//     srand(rlxd_seed());
    unsigned long int seed = rlxd_seed();         //initialise integer seed and set to value shown
//     unsigned long int seed = random_seed();         //initialise integer seed and set to value shown

    gsl_rng_env_setup();                        //read in environmental variables set at command line

    Type = gsl_rng_default;                        //set random number generator to default, which is mt19937
    r = gsl_rng_alloc (Type);                       //Allocate memory for type T generator and get pointer
    gsl_rng_set (r, seed);                      //initialise random number generator

    unsigned int NN = N;                    //Sample size
    size_t KK =K;
//     double mx[KK];                    //Initialise pointer to mx real positive array of probabilities

//     unsigned int n[KK];           //Initialise output integer array

//     printf ("generator type: %s\n", gsl_rng_name (r));
//     printf ("seed = %d\n", seed);

//     printf("NN=%d\n",NN);
//     printf("KK=%d\n",KK);

    gsl_ran_multinomial(r, KK, NN, x, n);

    gsl_rng_free(r);

//     free(n);
//     free(mx);

 }

/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{

     /* check for proper number of arguments */
    if(nrhs!=3) {
        mexErrMsgIdAndTxt("MyToolbox:Multinomial_gsl:nrhs","Three inputs required.");
    }
    if(nlhs!=1) {
        mexErrMsgIdAndTxt("MyToolbox:Multinomial_gsl:nlhs","One output required.");
    }

    /* get input variables */
    double N=mxGetScalar(prhs[0]);
    double *x = mxGetPr(prhs[1]);
    double K=mxGetScalar(prhs[2]);

    int *n;

    mwSize Ka=K;

//     plhs[0] = mxCreateDoubleMatrix(1, Ka,mxREAL);
    //mxCreateNumericMatrix is used here as we want to create a vector of integers
    plhs[0] = mxCreateNumericMatrix(1, Ka, mxUINT32_CLASS, mxREAL);

   /* get a pointer to the real data in the output matrix */
   //As arrays is not DOUBLE then we can't use mxGetPr
   n = (int*) mxGetData(plhs[0]);

   Multinomial(n,N,x,Ka);

}

max_stars_repo_name: BhavKhatri/MutliplexDriveResistanceSims

max_issues_repo_name: BhavKhatri/MutliplexDriveResistanceSims

max_forks_repo_name: BhavKhatri/MutliplexDriveResistanceSims


text:
/*
  NAME:
     proj_EM_step
  PURPOSE:
     one proj_EM step
  CALLING SEQUENCE:
     proj_EM_step(struct datapoint * data, int N, struct gaussian * gaussians,
     int K,bool * fixamp, bool * fixmean, bool * fixcovar,
     double * avgloglikedata, bool likeonly, double w, bool noproj,
     bool diagerrs, bool noweight)
  INPUT:
     data         - the data
     N            - number of data points
     gaussians    - model gaussians
     K            - number of model gaussians
     fixamp       - fix the amplitude?
     fixmean      - fix the mean?
     fixcovar     - fix the covar?
     likeonly     - only compute likelihood?
     w            - regularization parameter
     noproj       - don't perform any projections
     diagerrs     - the data->SS errors-squared are diagonal
     noweight     - don't use data-weights
  OUTPUT:
     avgloglikedata - average loglikelihood of the data
  REVISION HISTORY:
     2008-09-21 - Written Bovy
     2010-03-01 Added noproj option - Bovy
     2010-04-01 Added noweight option - Bovy
*/
#ifdef _OPENMP
#include <omp.h>
#endif

#include <math.h>
//#include <time.h>
//#include <sys/time.h>
#include <float.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <proj_gauss_mixtures.h>

#define CHUNKSIZE 1

void proj_EM_step(struct datapoint * data, int N,
		  struct gaussian * gaussians, int K,bool * fixamp,
		  bool * fixmean, bool * fixcovar, double * avgloglikedata,
		  bool likeonly, double w, bool noproj, bool diagerrs,
		  bool noweight){
  *avgloglikedata = 0.0;
  //struct timeval start,time1, time2, time3, time4, time5,end;
  struct datapoint * thisdata;
  struct gaussian * thisgaussian;
  struct gaussian * thisnewgaussian;
  int signum,di;
  double exponent;
  double currqij;
  struct modelbs * thisbs;
  int d = (gaussians->VV)->size1;//dim of mm

  //gettimeofday(&start,NULL);
  //Initialize new parameters
  int kk;
  for (kk=0; kk != K*nthreads; ++kk){
    newgaussians->alpha = 0.0;
    gsl_vector_set_zero(newgaussians->mm);
    gsl_matrix_set_zero(newgaussians->VV);
    ++newgaussians;
  }
  newgaussians= startnewgaussians;

  //gettimeofday(&time1,NULL);
  //check whether for some Gaussians none of the parameters get updated
  double sumfixedamps= 0;
  bool * allfixed = (bool *) calloc(K, sizeof (bool) );
  double ampnorm;
  for (kk=0; kk != K; ++kk){
    if (*fixamp == true){
      sumfixedamps += gaussians->alpha;
    }
    ++gaussians;
    if (*fixamp == true && *fixmean == true && *fixcovar == true)
      *allfixed= true;
    ++allfixed;
    ++fixamp;
    ++fixmean;
    ++fixcovar;
  }
  gaussians -= K;
  allfixed -= K;
  fixamp -= K;
  fixmean -= K;
  fixcovar -= K;

  //gettimeofday(&time2,NULL);

  //now loop over data and gaussians to update the model parameters
  int ii, jj, ll;
  double sumSV;
  int chunk;
  chunk= CHUNKSIZE;
#pragma omp parallel for schedule(static,chunk) \
  private(tid,di,signum,exponent,ii,jj,ll,kk,Tij,Tij_inv,wminusRm,p,VRTTinv,sumSV,VRT,TinvwminusRm,Rtrans,thisgaussian,thisdata,thisbs,thisnewgaussian,currqij) \
  shared(newgaussians,gaussians,bs,allfixed,K,d,data,avgloglikedata)
  for (ii = 0 ; ii < N; ++ii){
    thisdata= data+ii;
#ifdef _OPENMP
    tid= omp_get_thread_num();
#else
    tid = 0;
#endif
    di = (thisdata->SS)->size1;
    //printf("Datapoint has dimension %i\n",di);
    p = gsl_permutation_alloc (di);
    wminusRm = gsl_vector_alloc (di);
    TinvwminusRm = gsl_vector_alloc (di);
    Tij = gsl_matrix_alloc(di,di);
    Tij_inv = gsl_matrix_alloc(di,di);
    if ( ! noproj ) VRT = gsl_matrix_alloc(d,di);
    VRTTinv = gsl_matrix_alloc(d,di);
    if ( ! noproj ) Rtrans = gsl_matrix_alloc(d,di);
    for (jj = 0; jj != K; ++jj){
      //printf("%i,%i\n",(thisdata->ww)->size,wminusRm->size);
      gsl_vector_memcpy(wminusRm,thisdata->ww);
      //fprintf(stdout,"Where is the seg fault?\n");
      thisgaussian= gaussians+jj;
      //prepare...
      if ( ! noproj ) {
	if ( diagerrs ) {
	  gsl_matrix_set_zero(Tij);
	  for (ll = 0; ll != di; ++ll)
	    gsl_matrix_set(Tij,ll,ll,gsl_matrix_get(thisdata->SS,ll,0));}
	else
	  gsl_matrix_memcpy(Tij,thisdata->SS);
      }
      //Calculate Tij
      if ( ! noproj ) {
	gsl_matrix_transpose_memcpy(Rtrans,thisdata->RR);
	gsl_blas_dsymm(CblasLeft,CblasUpper,1.0,thisgaussian->VV,Rtrans,0.0,VRT);//Only the upper right part of VV is calculated --> use only that part
	gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,thisdata->RR,VRT,1.0,Tij);}//This is Tij
      else {
	if ( diagerrs ) {
	  for (kk = 0; kk != d; ++kk){
	    gsl_matrix_set(Tij,kk,kk,gsl_matrix_get(thisdata->SS,kk,0)+gsl_matrix_get(thisgaussian->VV,kk,kk));
	    for (ll = kk+1; ll != d; ++ll){
	      sumSV= gsl_matrix_get(thisgaussian->VV,kk,ll);
	      gsl_matrix_set(Tij,kk,ll,sumSV);
	      gsl_matrix_set(Tij,ll,kk,sumSV);}}}
	else {
	  for (kk = 0; kk != d; ++kk){
	    gsl_matrix_set(Tij,kk,kk,gsl_matrix_get(thisdata->SS,kk,kk)+gsl_matrix_get(thisgaussian->VV,kk,kk));
	    for (ll = kk+1; ll != d; ++ll){
	      sumSV= gsl_matrix_get(thisdata->SS,kk,ll)+gsl_matrix_get(thisgaussian->VV,kk,ll);
	      gsl_matrix_set(Tij,kk,ll,sumSV);
	      gsl_matrix_set(Tij,ll,kk,sumSV);}}}}
      //gsl_matrix_add(Tij,thisgaussian->VV);}
      //Calculate LU decomp of Tij and Tij inverse
      gsl_linalg_LU_decomp(Tij,p,&signum);
      gsl_linalg_LU_invert(Tij,p,Tij_inv);
      //Calculate Tijinv*(w-Rm)
      if ( ! noproj ) gsl_blas_dgemv(CblasNoTrans,-1.0,thisdata->RR,thisgaussian->mm,1.0,wminusRm);
      else gsl_vector_sub(wminusRm,thisgaussian->mm);
      //printf("wminusRm = %f\t%f\n",gsl_vector_get(wminusRm,0),gsl_vector_get(wminusRm,1));
      gsl_blas_dsymv(CblasUpper,1.0,Tij_inv,wminusRm,0.0,TinvwminusRm);
      //printf("TinvwminusRm = %f\t%f\n",gsl_vector_get(TinvwminusRm,0),gsl_vector_get(TinvwminusRm,1));
      gsl_blas_ddot(wminusRm,TinvwminusRm,&exponent);
      //printf("Exponent = %f\nDet = %f\n",exponent,gsl_linalg_LU_det(Tij,signum));
      gsl_matrix_set(qij,ii,jj,log(thisgaussian->alpha) - di * halflogtwopi - 0.5 * gsl_linalg_LU_lndet(Tij) -0.5 * exponent);//This is actually the log of qij
      //printf("Here we have = %f\n",gsl_matrix_get(qij,ii,jj));
      //Now calculate bij and Bij
      thisbs= bs+tid*K+jj;
      gsl_vector_memcpy(thisbs->bbij,thisgaussian->mm);
      //printf("%i,%i,%i\n",tid,ii,jj);
      if ( ! noproj ) gsl_blas_dgemv(CblasNoTrans,1.0,VRT,TinvwminusRm,1.0,thisbs->bbij);
      else gsl_blas_dsymv(CblasUpper,1.0,thisgaussian->VV,TinvwminusRm,1.0,thisbs->bbij);
      //printf("bij = %f\t%f\n",gsl_vector_get(bs->bbij,0),gsl_vector_get(bs->bbij,1));
      gsl_matrix_memcpy(thisbs->BBij,thisgaussian->VV);
      if ( ! noproj ) {
	gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,VRT,Tij_inv,0.0,VRTTinv);
	gsl_blas_dgemm(CblasNoTrans,CblasTrans,-1.0,VRTTinv,VRT,1.0,thisbs->BBij);}
      else {
	gsl_blas_dsymm(CblasLeft,CblasUpper,1.0,thisgaussian->VV,Tij_inv,0.0,VRTTinv);
	gsl_blas_dsymm(CblasRight,CblasUpper,-1.0,thisgaussian->VV,VRTTinv,1.0,thisbs->BBij);}
      gsl_blas_dsyr(CblasUpper,1.0,thisbs->bbij,thisbs->BBij);//This is bijbijT + Bij, which is the relevant quantity
      }
    gsl_permutation_free (p);
    gsl_vector_free(wminusRm);
    gsl_vector_free(TinvwminusRm);
    gsl_matrix_free(Tij);
    gsl_matrix_free(Tij_inv);
    if ( ! noproj ) gsl_matrix_free(VRT);
    gsl_matrix_free(VRTTinv);
    if ( ! noproj ) gsl_matrix_free(Rtrans);
    //Again loop over the gaussians to update the model(can this be more efficient? in any case this is not so bad since generally K << N)
#pragma omp critical
    {
      //Normalize qij properly
      *avgloglikedata += normalize_row(qij,ii,true,noweight,thisdata->logweight);
    }
      //printf("qij = %f\t%f\n",gsl_matrix_get(qij,ii,0),gsl_matrix_get(qij,ii,1));
      //printf("avgloglgge = %f\n",*avgloglikedata);
      for (jj = 0; jj != K; ++jj){
	currqij = exp(gsl_matrix_get(qij,ii,jj));
	//printf("Current qij = %f\n",currqij);
	thisbs= bs+tid*K+jj;
	thisnewgaussian= newgaussians+tid*K+jj;
	gsl_vector_scale(thisbs->bbij,currqij);
	gsl_vector_add(thisnewgaussian->mm,thisbs->bbij);
	gsl_matrix_scale(thisbs->BBij,currqij);
	gsl_matrix_add(thisnewgaussian->VV,thisbs->BBij);
	//printf("bij = %f\t%f\n",gsl_vector_get(bs->bbij,0),gsl_vector_get(bs->bbij,1));
	//printf("Bij = %f\t%f\t%f\n",gsl_matrix_get(bs->BBij,0,0),gsl_matrix_get(bs->BBij,1,1),gsl_matrix_get(bs->BBij,0,1));
      }
  }
  *avgloglikedata /= N;
  if (likeonly)
    return;
  //gettimeofday(&time3,NULL);

  //gather newgaussians
  if ( nthreads != 1 )
#pragma omp parallel for schedule(static,chunk) \
  private(ll,jj)
    for (jj = 0; jj < K; ++jj)
      for (ll = 1; ll != nthreads; ++ll) {
	gsl_vector_add((newgaussians+jj)->mm,(newgaussians+ll*K+jj)->mm);
	gsl_matrix_add((newgaussians+jj)->VV,(newgaussians+ll*K+jj)->VV);
      }

  //gettimeofday(&time4,NULL);

  //Now update the parameters
  //Thus, loop over gaussians again!
  double qj;
#pragma omp parallel for schedule(dynamic,chunk) \
  private(jj,qj)
  for (jj = 0; jj < K; ++jj){
    if (*(allfixed+jj)){
      continue;
    }
    else {
      qj = exp(logsum(qij,jj,false));
      (qj < DBL_MIN) ? qj = 0: 0;
      //printf("qj = %f\n",qj);
      if (*(fixamp+jj) != true) {
	(gaussians+jj)->alpha = qj;
	if (qj == 0) {//rethink this
	  *(fixamp+jj)=1;
	  *(fixmean+jj)=1;
	  *(fixcovar+jj)=1;
	  continue;
	}
      }
      gsl_vector_scale((newgaussians+jj)->mm,1.0/qj);
      if (*(fixmean+jj) != true){
	gsl_vector_memcpy((gaussians+jj)->mm,(newgaussians+jj)->mm);
    }
      if (*(fixcovar+jj) != true){
	//	if (*(fixmean+jj) != true)
	//  gsl_blas_dsyr(CblasUpper,-qj,(gaussians+jj)->mm,(newgaussians+jj)->VV);
	//else {
	gsl_blas_dsyr(CblasUpper,qj,(gaussians+jj)->mm,(newgaussians+jj)->VV);
	gsl_blas_dsyr2(CblasUpper,-qj,(gaussians+jj)->mm,(newgaussians+jj)->mm,(newgaussians+jj)->VV);
	//}
      if (w > 0.){
	gsl_matrix_add((newgaussians+jj)->VV,I);
	gsl_matrix_scale((newgaussians+jj)->VV,1.0/(qj+1.0));
      }
      else gsl_matrix_scale((newgaussians+jj)->VV,1.0/qj);
      gsl_matrix_memcpy((gaussians+jj)->VV,(newgaussians+jj)->VV);
      }
    }
  }
  //gettimeofday(&time5,NULL);

  //normalize the amplitudes
  if ( sumfixedamps == 0. && noweight ){
    for (kk=0; kk != K; ++kk){
      if ( noweight ) (gaussians++)->alpha /= (double) N;
    }
  }
  else {
    ampnorm= 0;
    for (kk=0; kk != K; ++kk){
      if (*(fixamp++) == false) ampnorm += gaussians->alpha;
      ++gaussians;
    }
    fixamp -= K;
    gaussians -= K;
    for (kk=0; kk != K; ++kk){
      if (*(fixamp++) == false){
	gaussians->alpha /= ampnorm;
	gaussians->alpha *= (1. - sumfixedamps);
      }
      ++gaussians;
    }
    fixamp -= K;
    gaussians -= K;
  }
  //gettimeofday(&end,NULL);
  //double diff, diff1, diff2, diff3, diff4, diff5,diff6;
  //diff= difftime (end.tv_sec,start.tv_sec)+difftime (end.tv_usec,start.tv_usec)/1000000;
  //diff1= (difftime(time1.tv_sec,start.tv_sec)+difftime(time1.tv_usec,start.tv_usec)/1000000)/diff;
  //diff2= (difftime(time2.tv_sec,time1.tv_sec)+difftime(time2.tv_usec,time1.tv_usec)/1000000)/diff;
  //diff3= (difftime(time3.tv_sec,time2.tv_sec)+difftime(time3.tv_usec,time2.tv_usec)/1000000)/diff;
  //diff4= (difftime(time4.tv_sec,time3.tv_sec)+difftime(time4.tv_usec,time3.tv_usec)/1000000)/diff;
  //diff5= (difftime(time5.tv_sec,time4.tv_sec)+difftime(time5.tv_usec,time4.tv_usec)/1000000)/diff;
  //diff6= (difftime(end.tv_sec,time5.tv_sec)+difftime(end.tv_usec,time5.tv_usec)/1000000)/diff;
  //printf("%f,%f,%f,%f,%f,%f,%f\n",diff,diff1,diff2,diff3,diff4,diff5,diff6);

  free(allfixed);

  return;
}


text:
/**
 *
 * @file qwrapper_slange.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Julien Langou
 * @author Henricus Bouwmeester
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:44:56 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_slange(Quark *quark, Quark_Task_Flags *task_flags,
                       int norm, int M, int N,
                       const float *A, int LDA, int szeA,
                       int szeW, float *result)
{
    szeW = max(1, szeW);
    DAG_CORE_LANGE;
    QUARK_Insert_Task(quark, CORE_slange_quark, task_flags,
        sizeof(PLASMA_enum),                &norm,  VALUE,
        sizeof(int),                        &M,     VALUE,
        sizeof(int),                        &N,     VALUE,
        sizeof(float)*szeA,     A,             INPUT,
        sizeof(int),                        &LDA,   VALUE,
        sizeof(float)*szeW,                 NULL,          SCRATCH,
        sizeof(float),                      result,        OUTPUT,
        0);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_slange_f1(Quark *quark, Quark_Task_Flags *task_flags,
                          PLASMA_enum norm, int M, int N,
                          const float *A, int LDA, int szeA,
                          int szeW, float *result,
                          float *fake, int szeF)
{
    szeW = max(1, szeW);
    DAG_CORE_LANGE;

    if ( result == fake ) {
        QUARK_Insert_Task(quark, CORE_slange_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(int),                        &M,     VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(float)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(float)*szeW,                 NULL,          SCRATCH,
            sizeof(float),                      result,        OUTPUT | GATHERV,
            0);
    } else {
        QUARK_Insert_Task(quark, CORE_slange_f1_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(int),                        &M,     VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(float)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(float)*szeW,                 NULL,          SCRATCH,
            sizeof(float),                      result,        OUTPUT,
            sizeof(float)*szeF,                 fake,          OUTPUT | GATHERV,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slange_quark = PCORE_slange_quark
#define CORE_slange_quark PCORE_slange_quark
#endif
void CORE_slange_quark(Quark *quark)
{
    float *normA;
    int norm;
    int M;
    int N;
    float *A;
    int LDA;
    float *work;

    quark_unpack_args_7(quark, norm, M, N, A, LDA, work, normA);
    *normA = LAPACKE_slange_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm),
        M, N, A, LDA, work);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slange_f1_quark = PCORE_slange_f1_quark
#define CORE_slange_f1_quark PCORE_slange_f1_quark
#endif
void CORE_slange_f1_quark(Quark *quark)
{
    float *normA;
    int norm;
    int M;
    int N;
    float *A;
    int LDA;
    float *work;
    float *fake;

    quark_unpack_args_8(quark, norm, M, N, A, LDA, work, normA, fake);
    *normA = LAPACKE_slange_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm),
        M, N, A, LDA, work);
}


text:
// SPDX-License-Identifier: Unlicense

/*
 * libinput calibration calculation
 */

#include <stdio.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

/*
 * Compute the (Moore-Penrose) pseudo-inverse of a matrix.
 *
 * PUBLIC DOMAIN by Charl Linssen <charl@itfromb.it>
 *
 * If the singular value decomposition (SVD) of A = UΣVᵀ then the
 * pseudoinverse A⁻¹ = VΣ⁻¹Uᵀ, where ᵀ indicates transpose and
 * Σ⁻¹ is obtained by taking the reciprocal of each nonzero element on
 * the diagonal, leaving zeros in place. Elements on the diagonal smaller
 * than ``rcond`` times the largest singular value are considered zero.
 *
 * @parameter A		Input matrix. **WARNING**: the input matrix ``A`` is
 * 			destroyed. However, it is still the responsibility of
 * 			the caller to free it.
 *
 * @parameter rcond	A real number specifying the singular value threshold
 * 			for inclusion. NumPy default for ``rcond`` is 1E-15.
 *
 * @returns A_pinv	Matrix containing the result. ``A_pinv`` is allocated
 * 			in this function and it is the responsibility of the
 * 			caller to free it.
 *
 *
 */
gsl_matrix* moore_penrose_pinv(gsl_matrix *A, const double rcond) {

	gsl_matrix *V, *Sigma_pinv, *U, *A_pinv;
	gsl_matrix *_tmp_mat = NULL;
	gsl_vector *_tmp_vec;
	gsl_vector *u;
	double x, cutoff;
	size_t i, j;
	unsigned int n = A->size1;
	unsigned int m = A->size2;
	unsigned int was_swapped = 0;

	if (m > n) {
		/* libgsl SVD can only handle the case m <= n - transpose matrix */
		was_swapped = 1;
		_tmp_mat = gsl_matrix_alloc(m, n);
		gsl_matrix_transpose_memcpy(_tmp_mat, A);
		A = _tmp_mat;
		i = m;
		m = n;
		n = i;
	}

	/* do SVD */
	V = gsl_matrix_alloc(m, m);
	u = gsl_vector_alloc(m);
	_tmp_vec = gsl_vector_alloc(m);
	gsl_linalg_SV_decomp(A, V, u, _tmp_vec);
	gsl_vector_free(_tmp_vec);

	/* compute Σ⁻¹ */
	Sigma_pinv = gsl_matrix_alloc(m, n);
	gsl_matrix_set_zero(Sigma_pinv);
	cutoff = rcond * gsl_vector_max(u);

	for (i = 0; i < m; ++i) {
		if (gsl_vector_get(u, i) > cutoff) {
			x = 1. / gsl_vector_get(u, i);
		}
		else {
			x = 0.;
		}
		gsl_matrix_set(Sigma_pinv, i, i, x);
	}

	/* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
	U = gsl_matrix_alloc(n, n);
	gsl_matrix_set_zero(U);

	for (i = 0; i < n; ++i) {
		for (j = 0; j < m; ++j) {
			gsl_matrix_set(U, i, j, gsl_matrix_get(A, i, j));
		}
	}

	if (_tmp_mat != NULL) {
		gsl_matrix_free(_tmp_mat);
	}

	/* two dot products to obtain pseudoinverse */
	_tmp_mat = gsl_matrix_alloc(m, n);
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., _tmp_mat);

	if (was_swapped) {
		A_pinv = gsl_matrix_alloc(n, m);
		gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1., U, _tmp_mat, 0., A_pinv);
	}
	else {
		A_pinv = gsl_matrix_alloc(m, n);
		gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1., _tmp_mat, U, 0., A_pinv);
	}

	gsl_matrix_free(_tmp_mat);
	gsl_matrix_free(U);
	gsl_matrix_free(Sigma_pinv);
	gsl_vector_free(u);
	gsl_matrix_free(V);

	return A_pinv;
}

int main(int argc, char **argv) {
	const unsigned int col = 4;
	const unsigned int row = 3;
	const double rcond = 1E-15;

	const double oneEigth = 1.0/8.0;
	const double sevenEight = 7.0/8.0;

	if (argc < 11) {
		printf("usage: %s x-res y-res x0 y0 x1 y1 x2 y2 x3 y3 [swap-x-y]\n", argv[0]);
		return EXIT_FAILURE;
	}

        int xres = atoi(argv[1]);
        int yres = atoi(argv[2]);
        int x0in = atoi(argv[3]);
        int y0in = atoi(argv[4]);
        int x1in = atoi(argv[5]);
        int y1in = atoi(argv[6]);
        int x2in = atoi(argv[7]);
        int y2in = atoi(argv[8]);
        int x3in = atoi(argv[9]);
        int y3in = atoi(argv[10]);

	int swapxy = 0;
	if (argc > 11)
		swapxy = 1;

	double x0 = (double)x0in / (double)xres;
	double y0 = (double)y0in / (double)yres;
	double x1 = (double)x1in / (double)xres;
	double y1 = (double)y1in / (double)yres;
	double x2 = (double)x2in / (double)xres;
	double y2 = (double)y2in / (double)yres;
	double x3 = (double)x3in / (double)xres;
	double y3 = (double)y3in / (double)yres;

	gsl_matrix *A = gsl_matrix_alloc(row, col);
	gsl_matrix *C = gsl_matrix_alloc(row, col);
	gsl_matrix *R = gsl_matrix_alloc(3, 3);
	gsl_matrix *S;
	gsl_matrix *T;
	gsl_matrix *A_pinv;

	gsl_matrix_set(A, 0, 0, x0);
	gsl_matrix_set(A, 0, 1, x1);
	gsl_matrix_set(A, 0, 2, x2);
	gsl_matrix_set(A, 0, 3, x3);
	gsl_matrix_set(A, 1, 0, y0);
	gsl_matrix_set(A, 1, 1, y1);
	gsl_matrix_set(A, 1, 2, y2);
	gsl_matrix_set(A, 1, 3, y3);
	gsl_matrix_set(A, 2, 0, 1.0);
	gsl_matrix_set(A, 2, 1, 1.0);
	gsl_matrix_set(A, 2, 2, 1.0);
	gsl_matrix_set(A, 2, 3, 1.0);

	gsl_matrix_set(C, 0, 0, oneEigth);
	gsl_matrix_set(C, 0, 1, sevenEight);
	gsl_matrix_set(C, 0, 2, oneEigth);
	gsl_matrix_set(C, 0, 3, sevenEight);
	gsl_matrix_set(C, 1, 0, oneEigth);
	gsl_matrix_set(C, 1, 1, oneEigth);
	gsl_matrix_set(C, 1, 2, sevenEight);
	gsl_matrix_set(C, 1, 3, sevenEight);
	gsl_matrix_set(C, 2, 0, 1.0);
	gsl_matrix_set(C, 2, 1, 1.0);
	gsl_matrix_set(C, 2, 2, 1.0);
	gsl_matrix_set(C, 2, 3, 1.0);

	A_pinv = moore_penrose_pinv(A, rcond);
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, C, A_pinv, 0.0, R);

	if (swapxy) {
		S = gsl_matrix_alloc(3, 3);

		/* extra swapping x-y axes as requested */
		gsl_matrix_set(S, 0, 0,  0.0);
		gsl_matrix_set(S, 0, 1, -1.0);
		gsl_matrix_set(S, 0, 2,  1.0);
		gsl_matrix_set(S, 1, 0, -1.0);
		gsl_matrix_set(S, 1, 1,  0.0);
		gsl_matrix_set(S, 1, 2,  1.0);
		gsl_matrix_set(S, 2, 0,  0.0);
		gsl_matrix_set(S, 2, 1,  0.0);
		gsl_matrix_set(S, 2, 2,  1.0);

		T = R;
		R = gsl_matrix_alloc(3, 3);
		gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, T, S, 0.0, R);
	}

        printf("%f, %f, %f, %f, %f, %f, %f, %f, %f\n",
                gsl_matrix_get((const gsl_matrix *)R, 0, 0), gsl_matrix_get((const gsl_matrix *)R, 0, 1), gsl_matrix_get((const gsl_matrix *)R, 0, 2),
                gsl_matrix_get((const gsl_matrix *)R, 1, 0), gsl_matrix_get((const gsl_matrix *)R, 1, 1), gsl_matrix_get((const gsl_matrix *)R, 1, 2),
                gsl_matrix_get((const gsl_matrix *)R, 2, 0), gsl_matrix_get((const gsl_matrix *)R, 2, 1), gsl_matrix_get((const gsl_matrix *)R, 2, 2));

	gsl_matrix_free(A);
	gsl_matrix_free(C);
	gsl_matrix_free(R);
	gsl_matrix_free(A_pinv);

	if (swapxy) {
		gsl_matrix_free(S);
		gsl_matrix_free(T);
	}
	return EXIT_SUCCESS;
}


text:
#ifndef H_LASER_DATA_JSON
#define H_LASER_DATA_JSON

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

#include <json-c/json.h>
#include <json-c/json_more_utils.h>

#include "laser_data.h"
#include "algos.h"

/* Laserdata to/from json */

JO ld_to_json(LDP);
LDP json_to_ld(JO);

JO corr_to_json(struct correspondence*, int n);
int json_to_corr(JO jo, struct correspondence*, int n);

LDP ld_from_json_stream(FILE*);
LDP ld_from_json_string(const char*s);
void ld_write_as_json(LDP ld, FILE * stream);

/* Other stuff to/from json */

JO matrix_to_json(gsl_matrix*m);
JO vector_to_json(gsl_vector*v);
JO result_to_json(struct sm_params*p, struct sm_result *r);

#endif

max_stars_repo_name: ozaslan/basics

max_issues_repo_name: ozaslan/basics

max_forks_repo_name: ozaslan/basics


text:
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>
#include "timer.c"
#include "timer.h"

double g (double *t, size_t dim, void *params);

void timer_start(void);
double timer_stop(void);

int main (void)
{
  double tvegas, tmonte;
  double res, err;
  double distmin = 1.001;
  double distmax = 4.;
  double dist = distmin;
  int np = 20;
  double nt = (distmax - distmin) / ((double) np - 1.);
  size_t dim = 6;
  double x1[] = { 0., 0., 0., 0., 0., 0. };
  double xu[] = { 1., 1., 1., 1., 1., 1. };

  double vegas[20], dipole[20], distance[20], monte[20];
  int f = 0;
  gsl_rng *r = gsl_rng_alloc (gsl_rng_taus2);
  unsigned long seed = 1UL;

  gsl_rng_set (r, seed);

  size_t calls = 1500000;

  //Vegas Integration

  gsl_monte_function G = { &g, dim, &dist };

  gsl_monte_vegas_state *sv = gsl_monte_vegas_alloc (dim);

  gsl_monte_vegas_init (sv);

  timer_start();

  do
    {

      gsl_monte_vegas_integrate (&G, x1, xu, dim, calls / 5, r, sv, &res,
				 &err);
      do
	{
	  gsl_monte_vegas_integrate (&G, x1, xu, dim, calls, r, sv, &res,
				     &err);
	  fflush (stdout);
	}
      while (fabs (gsl_monte_vegas_chisq (sv) - 1.) > .2);
      vegas[f] = res;
      f++;
      dist = dist + nt;
    }
  while (dist <= distmax);

  tvegas = timer_stop();

  gsl_monte_vegas_free (sv);

  // Montecarlo Integration

  double sum = 0.;
  double x[6];
  long i, j, l;

  l = 1E6;
  dist = distmin;

  timer_start();

  for (j = 0; j < np; j++)
    {
      sum = 0.;
      for (i = 0; i < l; i++)
	{
	  for (int k = 0; k < (int) dim; k++)
	    {
	      x[k] = gsl_rng_uniform (r);
	    }
	  sum += g (x, dim, &dist);
	}

      res = sum / (double) l;
      monte[j] = res;
      dist += nt;
    }

  tmonte =  timer_stop();

  gsl_rng_free (r);
  dist = distmin;

  for (int y = 0; y < np; y++)
    {
      dipole[y] = 2. / pow (dist, 3);
      distance[y] = dist;
      dist = dist + nt;
    }

  printf
    ("#    Dist               Vegas           Monte       Dipolapprox\n");
  for (int q = 0; q < np; q++)
    {
      double dd = distance[q];
      double vv = fabs (vegas[q]);
      double hh = fabs (monte[q]);
      double di = fabs (dipole[q]);
      printf ("   %.6f         %.6f       %.6f      %.6f\n", dd, vv, hh, di);
    }

  printf("\n\n");
  printf("Time for Vegas  = %f\n", tvegas);
  printf("Time for Monte  = %f\n", tmonte);
  printf("Speed up Factor = %f\n", tvegas/tmonte);

  return 0;
}