Advanced Introduction to C++, Scientific Computing and Machine Learning




Claudius Gros, WS 2019/20

Institut für theoretische Physik
Goethe-University Frankfurt a.M.

Functions & Templates


functions: global & static variables

Copy Copy to clipboad
Downlaod Download
#include <iostream>
#include <cstdio>       // equivalent to <stdio.h>  
using namespace std;
 
int globalVariable = 0;        // for everbody
 
int myFunction(int a)          // call by value
{
 static int storeValue = 0;    // executed only (!) at first function call
 storeValue++;                 // counting number of times called
 printf("myFunction : %3d %5d\n",storeValue, globalVariable);
 globalVariable += 100;
 a = a + 1;
 return a + 1;
}
 
int main()
{
 int mainInt = 10;
 int returnInt = myFunction(mainInt);
//
 printf("mainInt    : %3d\n",mainInt);
 printf("returnInt  : %3d\n",returnInt);
//
 globalVariable += 10;
 returnInt = myFunction(mainInt);
 globalVariable += 10;
 returnInt = myFunction(mainInt);
 globalVariable += 10;
 returnInt = myFunction(mainInt);
 return 1;
}

functions & references

Copy Copy to clipboad
Downlaod Download
#include <iostream>
#include <stdio.h>             // c-type output with printf
using namespace std;
 
int myReferenceCall(int &a)     // call by reference
{
 a = a + 1;                     // changing a
 return a + 1;                 
}
 
int myPointerCall(int *b)       // calling with a pointer
{
 *b = (*b) * 3;                 // (*b) is the value
 return *b;                    
}
 
int main()
{
 int beforeReference = 10;
 int beforePointer   = 10;
 int referenceInt = myReferenceCall(beforeReference);
 int pointerInt   = myPointerCall(&beforePointer);
 
 
 printf("beforeReference : %5d\n",beforeReference);
 printf("referenceInt    : %5d\n",referenceInt);
//
 printf("\n");
//
 printf("beforePointer   : %5d\n",beforePointer);
 printf("pointerInt      : %5d\n",pointerInt);
 
 return 1;
}

passing arrays to functions

Copy Copy to clipboad
Downlaod Download
#include <iostream>
#include <stdio.h>                // c-type output with printf
#include <cstdlib>                // for atof
 
using namespace std;
 
void ordering(int sizeA, int a[])  // note: a is a pointer and hence
 {                                 // sizeof(a) / sizeof(a[0]);
                                   // does not work anymore 
 printf("# in ordering(): sizeof(a)       = %d\n",(int)sizeof(&a));
 printf("# in ordering(): sizeof(a[0])    = %d\n",(int)sizeof(a[0]));
 int xx;
 for (int i=0;i<sizeA;i++)
   for (int j=i+1;j<sizeA;j++)
     if (a[i]>a[j])
       { 
       xx = a[i];
       a[i] = a[j];
       a[j] = xx;
       }
 }                                // end of ordering()
 
int main(int argLength, char* argValues[]) 
{
 size_t arraySize = argLength-1;
 int myArray[arraySize];          
//
 for (int k=0; k<arraySize; k++)
   myArray[k] = atof(argValues[k+1]);
//
 printf("\n");
 printf(" input array: ");
 for (int k=0; k<arraySize; k++)
   printf("%5d",myArray[k]);
 printf("\n");
//
  printf("# in     main(): sizeof(myArray) = %d\n",(int)sizeof(myArray));
  ordering(arraySize,myArray);
//
 printf("\n");
 printf("orderd array: ");
 for (int k=0; k<arraySize; k++)
   printf("%5d",myArray[k]);
 printf("\n");
//
 return 1;
}

constants and function arguments

a plethora of function arguments

C++ namespaces

Copy Copy to clipboad
Downlaod Download
#include <iostream>
#include <stdio.h>             // c-type output with printf
using namespace std;                 // for everthing to come
 
namespace Tom 
{
 void aFunction(int a)        // Tom has a function
   {printf("Hi %d, welcome to the namespace of Tom!\n",a);}
}
 
namespace Lara
{
 void aFunction(int a)        // Lara has a function too
   {printf("Hi %d, welcome to the namespace of Lara!\n",a);}
 inline int plus4(int a)     // code replacement when called
   {return a += 4;}
}
 
int main()
{
 int mainInt = 10;
// 
 Tom::aFunction(mainInt);
 mainInt = Lara::plus4(mainInt);
 Lara::aFunction(mainInt);
 
 return 1;
}

function naming / doxygen documentation

Copy Copy to clipboad
Downlaod Download
/** @brief illustrates function naming and doxygen documentation
  * @file functionNaming.cpp
  */
#include <iostream>
#include <stdio.h>             // c-type output with printf
using namespace std;


/** @return absolute value of an double
  */
double myAbs(double a)
{
 printf("         : using   double myAbs(double a)\n");
 return (a>0.0)?a:(-a);
}
 
/** @return  absolute value of an int
  * @author  Superior Programmer
  * @date    most exciting day of year
  * @version last version before final version
  * @brief   comment block illustration
  */
int myAbs(int a)                 // alowed
{
 printf("         : using   int myAbs(int a)\n");
 return (a>0.0)?a:(-a);
}
 
/** @brief entry point
  */
int main()
{
 double inA;
 int intInA;
 printf("please enter a number : ");
 cin >> inA;
 intInA = (int)inA;
//
 printf("\n");
 printf(" inA     : %8.2f\n",inA);
 printf("|inA|    : %8.2f\n",myAbs(inA));
 printf("\n");
 printf(" intInA  : %8d\n",intInA);
 printf("|intInA| : %8d\n",myAbs(intInA));
 
 return 1;
}

function templates

Copy Copy to clipboad
Downlaod Download
#include <iostream>
#include <stdio.h>             // c-type output with printf
using namespace std;
 
// -- ---------------------------
// -- max for all even data types
// -- ---------------------------
 
template <typename T>          // T can be: int, double, ...
T maxT(T a, T b)             //      for: return and argument types
{
 return (a>b)?a:b;
}
 
// -- -----------------------
// -- swapping two variables
// -- ----------------------
 
template <typename T>       
void swapVariables(T &var1, T &var2)  // call by value
{
  T xx = var1;
  var1 = var2;
  var2 = xx;
}
 
// -- ----
// -- main
// -- ----
 
int main()
{
 double inA, inB;
 printf("please enter two real numbers : ");
 cin >> inA >> inB;
//
 printf("\n");
 printf("your input : %7.3f %7.3f\n",inA,inB);
 swapVariables(inA,inB);
 printf("your input : %7.3f %7.3f\n",inA,inB);
//
 printf("\n");
 printf("the larger real number is : %7.3f\n",maxT(inA,inB));
 int intInA = (int)inA;
 int intInB = (int)inB;
 printf("the larger int  number is : %7d\n",maxT(intInA,intInB));
 printf("\n");
//
 return 1;
}

recursive functions calls

Copy Copy to clipboad
Downlaod Download
#include <iostream>   /* standard IO  */
#include <stdio.h>    /* for printf  */
#include <cstdlib>    /* for atof     */
#include <stdlib.h>   /* srand, rand  */
#include <time.h>     /* for the time */
 
using namespace std;
 
// -- ----------------------------
// -- ordering for all array types
// -- ----------------------------
 
template <typename T>       
void ordering(int sizeA,T a[])    // note: a is a pointer and hence
{                                 
 T xx;
 for (int i=0;i<sizeA;i++)
   for (int j=i+1;j<sizeA;j++)
     if (a[i]>a[j])
       { 
       xx = a[i];
       a[i] = a[j];
       a[j] = xx;
       }
}                                 // end of ordering()
 
// -- ---------------------------
// -- generate a random int array
// -- ---------------------------
 
void randArray(int sizeA, int a[], int range)
{                                 
 srand( (unsigned)time( NULL ) );                  // random seed 
 for (int i=0;i<sizeA;i++)
   a[i] = (int)(range*(rand()/(double)RAND_MAX));  // (..) important!
}
 
// -- ------------------------------------
// -- recursive search of an ordered array
// -- ------------------------------------
 
template <typename T>             // any type: int, double, ..
int recursiveSearch(T sortedArray[], int first, int last, T key) 
{
 if (first <= last) 
   {
   int mid = (first + last) / 2;  // compute mid point.
   if (key == sortedArray[mid]) 
     return mid;                  // found it.
   else 
      if (key < sortedArray[mid]) 
        return recursiveSearch(sortedArray, first, mid-1, key);
                                  // call youself for the lower part of the array
      else
        return recursiveSearch(sortedArray, mid+1, last, key);
                                  // call ourself for the upper part of the array
   }
 return -(first + 1);             // failed to find key
}
 
// -- ----
// -- main
// -- ----
 
int main(int argLength, char* argValues[]) 
{
 const int sizeA = 10;
 int range = 20;
 int a[sizeA] = {0};       // initializing to zero
 randArray(sizeA,a,range); // generating random entries
//
 printf("\n");
 printf(" input array: ");
 for (int k=0; k<sizeA; k++)
   printf("%4d",a[k]);
 printf("\n");
//
  ordering(sizeA,a);       // ordering the array
//
 printf("\n");
 printf("orderd array: ");
 for (int k=0; k<sizeA; k++)
   printf("%4d",a[k]);
 printf("\n");
//
 int guess;
 printf("\n");
 printf("please guess a number: ");
 cin >> guess;           // guessing a number
//
 int searchResult = recursiveSearch(a,0,sizeA-1,guess);
 printf("\n");
 printf("your search result: %4d\n ",searchResult);
//
 return 1;
 }

optional function arguments

$$ \sum_{n=0}^\infty x^n \ =\ \frac{1}{1-x} \qquad\quad |x|<1 $$
Copy Copy to clipboad
Downlaod Download
#include <iostream>       // standard IO
#include <stdio.h>        // for printf 
using namespace std;

double geometricSeries(double x, bool starting=true) // the second argument is optional
 {
 static const double accuracy = 0.000001;
 static double factor;
 static double result;
 if (starting)
   {
   factor = 1.0;             // x*x* ... *x
   result = 0.0;             // the sum over powers
   }
 result += factor;
 factor *= x;
//
// printf("factor, result %12.8f %12.8f\n",factor, result);
//
 if (factor*factor>accuracy*accuracy)
   geometricSeries(x, false);
//
 return result;
 }

int main()
{
 double x = 0.4; 
 double result = geometricSeries(x);              // called only with a single argument
 printf("\n");
 printf("x, result, exact %6.2f %16.8f %12.8f\n",x , result,1.0/(1.0-x));
 return 1;
}

pointers to functions

Copy Copy to clipboad
Downlaod Download
/** @brief function pointers may be passed as arguments to other funcions
  * @file functionPointers.cpp
  */
#include <iostream>  /* standard IO */
#include <stdio.h>   /* for printf  */
#include <stdlib.h>  /* srand, rand */
#include <time.h>    /* for the time */
using namespace std;
 
/** @brief defining a pointer to a void functions with int argument
  */
typedef void (*FF_void_int)(int);  // new types should be global
 
 
/** @brief a dummy first function
  */
void firstFunction(int a)  { printf(" firstFunction():   a : %d\n",a); }

/** @brief a dummy second function
  */
void secondFunction(int a) { printf("secondFunction(): a*a : %d\n",a*a); }
 
/** @brief entry point
  */
int main()
{
 FF_void_int pointerToFunction;
//
 srand( (unsigned)time( NULL ) );  // initializing the random number generator
 for (int i=0;i<10;i++)
   {
   pointerToFunction =  firstFunction;   // functionsa are anyhow pointers
// pointerToFunction = &firstFunction;   // also ok
// pointerToFunction = *firstFunction;   // even this!
//
   if ((rand()%2)==1)
     pointerToFunction = *secondFunction;
//
   pointerToFunction(rand()%10);         // calling like this
//   (*pointerToFunction)(11);           // or like this
   }
//
 return 1;
}

functions as arguments

$$ \int_0^{2\pi} \frac{dx}{a+b\sin(x)}=\frac{2\pi}{\sqrt{a^2-b^2}}, \qquad\quad \int_0^{1} dx\frac{\log(1+x)}{x}=\frac{\pi^2}{12} $$
Copy Copy to clipboad
Downlaod Download
#include <iostream>          // standard IO
#include <stdio.h>           // for printf
#include <math.h>            // math
#include <assert.h>          // (`assert()' calls abort if arg not true)
 
#define VARIABLE_NAME(x) #x  // a C++ macro definition
using namespace std;
 
typedef double (*FF_double_double)(double);  // pointer to function
 
// --- -------------------
// --- defining integrands
// --- -------------------
 
inline double myLog_1(double x)
{
 assert(x>0.0);                           // should always be used
 double xx = log(1.0+x)/x;
 return xx*12.0/(M_PI*M_PI);
}
 
inline double mySin_1(double x)
{
 const double a = 5.0*0.5*M_PI;
 const double b = 3.0*0.5*M_PI;
 return 1.0/(a+b*sin(x));
}
 
// --- --------------------------
// --- defining simple integrator
// --- --------------------------
 
double simpleIntegrate(FF_double_double integrand,  // integrating the integrand
                       double lowerBound,           // from the lower to the
                       double upperBound)           // upper bound
{
 assert(upperBound>lowerBound); 
 int nSteps = 10;
 double x, deltaX = (upperBound-lowerBound)/nSteps;
 double result = 0.0;
//
 for (int iSteps=0;iSteps<nSteps;iSteps++)
   {
   x = lowerBound + (iSteps+0.5)*deltaX;
   result += integrand(x)*deltaX;                   // simplest 
   }
 return result;
} // end of simpleIntegrate()
 
// --- ----
// --- main
// --- ----
 
int main()
{
 double lowerBound = 0.0;
 double upperBound = 1.0;
 double result;
//
 result = simpleIntegrate(myLog_1,lowerBound,upperBound);
 printf("integrating %8s from %5.2f to %5.2f yields: %19.16f\n",
        VARIABLE_NAME(myLog_1), lowerBound, upperBound,result);
 printf("\n");
//
 lowerBound = 0.0;
 upperBound = 2.0*M_PI;
 result = simpleIntegrate(mySin_1,lowerBound,upperBound);
 printf("integrating %8s from %5.2f to %5.2f yields: %19.16f\n",
        VARIABLE_NAME(mySin_1), lowerBound, upperBound,result);
//
 return 1;
}

lambda expressions



  1. capture clause (or lambda-introducer)
    captures variables from the surrounding scope
  2. parameter list (or lambda declarator; optional)
  3. mutable specification (optional)
    allows to change captured variables
  4. exception-specification (optional)
  5. trailing-return-type (optional)
    normally automatically induced
  6. lambda body
Copy Copy to clipboad
Downlaod Download
#include <iostream>          // standard IO
#include <stdio.h>           // for printf
#include <math.h>            // math
#include <assert.h>          // (`assert()' calls abort if arg not true)
 
using namespace std;

/** @brief function type definition
  */
typedef double (*FF_double_double)(double);  
 
/** @brief Euler integrator
  */ 
double simpleIntegrate(FF_double_double integrand, 
                       double lowerBound,         
                       double upperBound)        
{
 assert(upperBound>lowerBound); 
 int nSteps = 10;
 double x, deltaX = (upperBound-lowerBound)/nSteps;
 double result = 0.0;
//
 for (int iSteps=0; iSteps < nSteps; iSteps++)
   {
   x = lowerBound + (iSteps+0.5)*deltaX;
   result += integrand(x)*deltaX;               
   }
 return result;
} 
 
/** @brief entry point
  */
int main()
{
 double lowerBound = 0.0;
 double upperBound = 1.0;
 double result;

 FF_double_double myFunction =         // binding lambda expression to function pointer
   [](double a){return 3.0*a*a;};      // lambda expresion

 result = simpleIntegrate(myFunction,lowerBound,upperBound);

 printf("integrating %18s from %5.2f to %5.2f yields: %12.8f\n\n",
        "f(x)=3*x*x    ", lowerBound, upperBound,result);

 result = simpleIntegrate([](double a){return 6.0*a*(1.0-a);}  // locally defined function
                            ,lowerBound,upperBound);

 printf("integrating %18s from %5.2f to %5.2f yields: %12.8f\n\n",
        "f(x)=6*x*(1-x)", lowerBound, upperBound,result);

// constants, as needd for array delimiters, can be evaluated
// on the fly by lambda expressions
 const int fibo5 = []() {                       // assigning const. variable to lambda expresssion
                        int a=1,b=1,c=0; 
                        for (int i=0;i<5;i++) {c=a+b;a=b;b=c;}
                        return c;
                        }();                    // () call lambda expression

 printf("constant variable value: %4d : \n\n", fibo5);

// capture the counting variable by reference and do stuff with it
 double rr=1.0;
 for (int i=0;                             // for loop start statement
      i<100;                               // for loop condition
      [&i](){printf("%3d :",i);i++;}()     // for loop loop-statement
     )
   { rr = rr*0.999;}                       // for loop body

 cout << endl; 
 return 1;
}

returning values by reference

Copy Copy to clipboad
Downlaod Download
#include <iostream>
 
using namespace std;

/** a global array
  */
double vals[] = {10.1, 20.2, 30.3, 40.4, 50.5};
 
/** a function returning the pointer (reference) 
  * to the ith element of the global array
  */
double& setValue(int i) {return vals[i]; }

/** a standard getter
  */
double  getValue(int i) {return vals[i]; }

/** just printing the array content
  */
template<int size>
void printArray(double inArray[size])
 {
 for (int i=0; i<size; i++) 
  {
  cout << "vals[" << i << "] = ";
  cout << vals[i] << endl;
  }
 }
 
/** entry point
  */
int main () 
{
cout << "value before change" << endl;
printArray<5>(vals);
 
setValue(1) = 3.0*getValue(1);  // change 2nd element
setValue(3) = 44.7;             // change 4th element
 
cout << "value after change" << endl;
printArray<5>(vals);
return 0;
}

outsourcing codes to multiple files

Copy Copy to clipboad
Downlaod Download
// --- -----------------
// --- the file main.cpp
// --- -----------------
#include <iostream>
 
using namespace std;
 
int add(int x, int y);  // forward declaration
 
int main()
{
 cout <<  endl;
 cout << "the sum of 3 and 4 is: " << add(3, 4) << endl;
 cout <<  endl;
 return 1;
}

Copy Copy to clipboad
Downlaod Download
// --- ----------------
// --- the file add.cpp
// --- ----------------
int add(int x, int y)
{
    return x + y;
}

C++ header files



Copy Copy to clipboad
Downlaod Download
// --- -----------------
// --- the file main.cpp
// --- -----------------
#include <iostream>
#include "add.h"         // my own header file
 
using namespace std;
 
int main()
{
 cout <<  endl;
 cout << "the sum of 3 and 4 is: " << add(3, 4) << endl;
 cout <<  endl;
 return 1;
}

Copy Copy to clipboad
Downlaod Download
// --- --------------
// --- the file add.h
// --- --------------
int add(int x, int y);    // forward declaration (function prototype)

Copy Copy to clipboad
Downlaod Download
// --- ----------------
// --- the file add.cpp
// --- ----------------
int add(int x, int y)
{
    return x + y;
}

Copy Copy to clipboad
Downlaod Download
// --- ----------------------------------
// --- the file add.h with a header guard
// --- ----------------------------------
#ifndef add_INCLUDED // is the identifier "add_INCLUDED" not yet defined (if not defined)?
#define add_INCLUDED // define "add_INCLUDED" 
 
int add(int x, int y);  // forward declaration (function prototype)
 
#endif                  // end of if-macro

single file projects

/* July 2001, Claudius Gros
 * stochastic series expansion with worm updates, here 
 * lattice with periodic boundary conditions or infinite lattices  
 */

/* g++     -O3 -DNDEBUG -o stoch.e stoch.c  (on a Linux machine) */
/*             -DNDEBUG    turns the assert() off */

#include <iostream>     // (input/ouput)
#include <iomanip>      // (manipulating input/ouput with `setw()' etc)
#include <fstream>      // (in/output on files)
#include <math.h>       // (math's commands)
#include <assert.h>     // (`assert()' calls abort if arg not true), 
                        // turnd off with -DNDEBUG
#include <time.h>       // ('clock()' measures time in nanoseconds)
using namespace std;

/* contains the classes
 * ====================
 * A_SITE        :  all static information relevant for given site
 * A_BOND        :  all static information relevant for a specfic bond 
 * A_LINK        :  a link in the stochastic series expansion 
 * SYSTEM        :  contains a variable array with all links and other data
 *               :  characterizing the system via stochastic series expansion.
 *               :  SYSTEM.diagonal_update() does the full diagonal update,
 *               :  i.e. adding and removing diagonal links
 *               :  SYSTEM.loop_update() constructs a full loop
 */

/* contains the soubroutines
 * =========================
 * spin_output()   : diagnostic output: spin-configuration and bond-types
 * link_output()   : diagnostic output: links in between the links
 * bond_output()   : diagnostic output: bonds
 * initial_hyper() : intializes the the lattice (for hypercubic lattices)
 * initial_fromFile() : reads lattice from file 'stoch.input'
 * expectationValues() : evaluates for given data mean, varianz, auto-tau
 * measure_chi()   : measures the magnetization (q=0 and q=pi) and chi's
 * measure_cV()    : measures the internal energy and spezific heat
 * timeInSeconds() : returns running time
 * ggubs()         : random-number generator
 * min()           : returns the minimum of to numbers (double,double)--> double
 * abs()           : returns the absolute value (double)--> double
 * main()          : here the program starts
 */ 

extern double ggubs();              // forward declaration, neccessary in old-fashion c++
extern double min(double,double);   // forward declaration, old-fashion c++

/* the following global constants are defined
 * ==========================================
 *  xL, yL, zL  : # unit-cells in x-/y-/z-direction  
 *  dimension   : dimension of the lattice 
 *  nSites      : total number if sites 
 *  nBonds      : total number of bonds 
 *  nHamil      : Hamiltonian = sum_{i=0}%{N_hamil-1} H_i 
 *                = # of spin configuration  on a cluster 
 *  nParHamil   : number of parameters in Hamiltonian
 *  Jxx,Jzz,Bz  : global coupling constants, longitudinal magnetic field
 *  dimerization: dimerization of interaction in between bonds
 *  nBondTypes  : number of different bonds (==spin configurations)
 *              : == 2 (in general) == 1 (for B==0)
 *  nBranchings : number of possible worm branchings at a given bond [1,2]
 *  maxMCsteps  : maximum number of MC-steps 
 *  smallNumber : a small number == zero
 */
/*  The following global variables are defined just behind the 
 *  respective class definition  : 
 *
 *  A_SITE allSites[nSites];     : array with all sites 
 *  A_BOND allBonds[nBonds];     : array with all bonds 
 *  SYSTEM system;               : class which all links and other system information
 */

const int xL = 10;        // if (xL*yL*zL==1) --> read from file 'stoch.input'
const int yL = 10;     
const int zL = 1;    
const int dimension = (xL*yL*zL==1)?3:( (yL*zL==1)?1:( (zL==1)?2:3 ) );
//
const int nSites    = (xL*yL*zL==1)?8:( xL*yL*zL);         
const int nBonds    = (xL*yL*zL==1)?12:( dimension*nSites ); 
//
const int nParHamil = 6;
const double Jxx          =  1.00;            
const double Jzz          =  1.00;            
const double Bz           =  0.0;            
const double dimerization =  0.0;
//
const int nBondTypes = 6;                          // : 'iBondTypes'
const int nBranchings = 1;                         // : 'ibranchings'
const int maxMCsteps =  20000;                     // : 'iMC' 
const double smallNumber = 0.0000000001;  
//
// nBondTypes = 6  corresponds to the following six spin configurations
// ====================================================================
//   
//  ferro(1) = + +   anti_ferro(1) = + -   spin-flip(1) = - +
//       [0]   + +             [2]   + -            [4]   + -    
//               
//  ferro(2) = - -   anti_ferro(2) = - +   spin_flip(2) = + - 
//       [1]   - -             [3]   - +            [5]   - +
//
//
/* ******************************* */
/* *****   class A_SITE    ******* */
/* ******************************* */
//
/*  contains for following member functions
 *  =======================================
 *  set/get_coordinate_X/Y/Z() : set/get the X/Y/Z coordinate of the lattice site
 */
//
class A_SITE
{
public:
//
inline void set_coordinate_X(int XX)
    {assert(dimension>=1); coordinates[0]=XX; }
inline  int get_coordinate_X()
    {assert(dimension>=1); return coordinates[0]; }
//
inline void set_coordinate_Y(int YY)
    {assert(dimension>=2); coordinates[1]=YY; }
inline  int get_coordinate_Y()
    {assert(dimension>=2); return coordinates[1]; }
//
inline void set_coordinate_Z(int ZZ)
    {assert(dimension>=3); coordinates[2]=ZZ; }
inline  int get_coordinate_Z()
    {assert(dimension>=3); return coordinates[2]; }
//
inline void set_sublattice(int AorB)
    {assert(AorB*AorB==1);sublattice = AorB; }
inline  int get_sublattice()
    {return sublattice; }
//
//  =======================================
//  the following data characterizes a site
//  =======================================
//
private:
  int coordinates[dimension];    // vector with the real-space postion 
  int sublattice;                // +1/-1  for  A-/B-sublattice
};   // end of class A_SITE
//
/* ***************************************** */
/* ***   definitions of global variable  *** */
/* ***        A_SITE allSites[]          *** */
/* ***************************************** */
//
   A_SITE allSites[nSites]; 

...