Creating external types: second example

Topics covered in this chapter:

Mosel defines the types integer, real and boolean on which arithmetic operations may be used. By creating modules it is possible to add other types, such as complex numbers, to this list. In the previous chapters we have already seen an example of how to define a new type in a module, but this new type task was not suited to be used with arithmetic operations. In this chapter we shall therefore give another example of the definition of a type, this time of a type to which such operations may sensibly be applied.

Example

In this chapter we are going to define the type complex to represent complex numbers. The following example demonstrates the typical uses that one may wish to make of a mathematical type like complex numbers in a model:

• use of data structures
• various types of initializations and assignments
• products, sums and other arithmetic operations
• comparison
• printed output on screen and to a file.

The following model shows how one might work with a new type complex in Mosel:

model "Test complex"
 uses "complex"

 declarations
  c:complex
  t:array(1..10) of complex
 end-declarations

 forall(j in 1..10) t(j):=complex(j,10-j)
 t(5):=complex("5+5i")

 c:=prod(i in 1..5) t(i)
 if c<>0 then
  writeln("product: ",c)
 end-if

 writeln("sum: ", sum(i in 1..10) t(i))
 c:= t(1)*t(3)/t(4) + if(t(2)=0,t(10),t(8)) + t(5) - t(9)
 writeln("result: ", c)

 initializations to "test.dat"
  c t
 end-initializations

end-model

Structures for passing information

Complex numbers are usually represented as a+bi where a and b are real numbers. a is called the real part and bi the imaginary part. We implement the following C structure to store a complex number:

typedef struct
    {
     unsigned int refcnt;     /* For reference count and shared flag */
     double re, im;           /* Real and imaginary parts */
    } s_complex;

List of subroutines

The main interest of this example lies in the definition of its list of subroutines which actually is a list of operators:

static XPRMdsofct tabfct[]=
    {
     {"@&", 1000, XPRM_TYP_EXTN, 1, "complex:|complex|", cx_new0},
     {"@&", 1001, XPRM_TYP_EXTN, 1, "complex:r", cx_new1},
     {"@&", 1002, XPRM_TYP_EXTN, 2, "complex:rr", cx_new2},
     {"@0", 1003, XPRM_TYP_EXTN, 0, "complex:", cx_zero},
     {"@1", 1004, XPRM_TYP_EXTN, 0, "complex:", cx_one},
     {"@:", 1005, XPRM_TYP_NOT,  2, "|complex||complex|", cx_asgn},
     {"@:", 1006, XPRM_TYP_NOT,  2, "|complex|r", cx_asgn_r},
     {"@+", 1007, XPRM_TYP_EXTN, 2, "complex:|complex||complex|", cx_pls},
     {"@+", 1008, XPRM_TYP_EXTN, 2, "complex:|complex|r", cx_pls_r},
     {"@*", 1009, XPRM_TYP_EXTN, 2, "complex:|complex||complex|", cx_mul},
     {"@*", 1010, XPRM_TYP_EXTN, 2, "complex:|complex|r", cx_mul_r},
     {"@-", 1011, XPRM_TYP_EXTN, 1, "complex:|complex|", cx_neg},
     {"@/", 1012, XPRM_TYP_EXTN, 2, "complex:|complex||complex|", cx_div},
     {"@/", 1013, XPRM_TYP_EXTN, 2, "complex:|complex|r", cx_div_r1},
     {"@/", 1014, XPRM_TYP_EXTN, 2, "complex:r|complex|", cx_div_r2},
     {"@=", 1015, XPRM_TYP_BOOL, 2, "|complex||complex|", cx_eql},
     {"@=", 1016, XPRM_TYP_BOOL, 2, "|complex|r", cx_eql_r}
    };

In the order of their appearance this list defines the following operators:

@&

creation (construction)

@0

zero element for sums

@1

one element for products

@:

assignment

@+

addition

@*

multiplication

@-

negation

@/

division

@=

comparison (test of equality)

For most operators in the list above several versions are defined, with different types or combinations of types. The only type conversion that is carried out automatically by Mosel is from integer to real (but not the other way round), and no conversions involving external types. It is therefore necessary to define all the operations between two numbers for two complex numbers and also for a complex and a real number. For commutative operations (addition, multiplication, comparison) it is only required to define one version combining the two types, the other sense is deduced by Mosel: for example, if complex + real is defined, Mosel `knows' how to calculate real + complex. For division (not commutative) we need to define every case separately.

List of types

The definition of the new type in the list of types that is passed to Mosel looks as follows:

static XPRMdsotyp tabtyp[]=
    {
     {"complex", 1, XPRM_DTYP_PNCTX|XPRM_DTYP_RFCNT|XPRM_DTYP_APPND,
      cx_create, cx_delete, cx_tostr, cx_fromstr, cx_copy, cx_compare}
    };

The type-related functions (cx_create: creation, cx_delete: deletion, cx_tostr: transformation to a string, cx_fromstr: initialization from a string, cx_copy: copying, cx_compare: comparison) could be implemented in a similar way to what has been shown for the task module in the previous chapters. But, for practical purposes, this rudimentary memory management may not be efficient enough. In this chapter we therefore give an example of improved memory management for external types (see Improved memory management for external types) and as a second implementation option, how to use Mosel's memory management functionality (see Mosel memory management). This includes new versions of the type instance creation and deletion functions, an adaptation of the reset service, and the definition of additional list structures for storing information in the module context.

The functions for converting types to or from strings and also the copy and compare functions described for the task module only require minor modifications to adapt them to this example. Their definition will not be repeated in this chapter.
The list of services (merely consisting of the reset service) and the main interface structure are also very similar to those of the task module, and the module initialization function remains the same except for its name. We therefore refrain from printing them here.

The complete source code of the complex module is among the module examples provided with the Mosel distribution and on the Xpress website.

Definition of operators

In this section we show several examples of the implementation of operators. A comprehensive list of all operators that may be defined in Mosel is given in the appendix.

Constructors

In the chapter about the task module we have already seen examples of functions for cloning a new type and constructing it in different ways. Here the cloning operation is implemented as follows:

static int cx_new0(XPRMcontext ctx, void *libctx)
{
 s_complex *complex,*new_complex;

 complex=XPRM_POP_REF(ctx);
 if(complex!=NULL)
 {
  new_complex=cx_create(ctx, libctx, NULL, 0);
  *new_complex=*complex;
  XPRM_PUSH_REF(ctx, new_complex);
 }
 else
  XPRM_PUSH_REF(ctx, NULL);
 return XPRM_RT_OK;
}

A new complex number is constructed from two given real numbers thus:

static int cx_new2(XPRMcontext ctx, void *libctx)
{
 s_complex *complex;

 complex=cx_create(ctx, libctx, NULL, 0);
 complex->re=XPRM_POP_REAL(ctx);
 complex->im=XPRM_POP_REAL(ctx);
 XPRM_PUSH_REF(ctx, complex);
 return XPRM_RT_OK;
}

Comparison operators

Another operation that we have already seen in the task module is the comparison between new types. This can be done in a very similar way for module complex and is not repeated here. In addition, it makes sense to define a comparison between a complex and a real number:

static int cx_eql_r(XPRMcontext ctx,void *libctx)
{
 s_complex *c1;
 double r;
 int b;

 c1=XPRM_POP_REF(ctx);
 r=XPRM_POP_REAL(ctx);
 if(c1!=NULL)
  b=(c1->im==0)&&(c1->re==r);
 else
  b=(r==0);
 XPRM_PUSH_INT(ctx,b);
 return XPRM_RT_OK;
}

Arithmetic operators

The arithmetic operations must implement the rules to perform these operations on complex numbers.

Multiplication

Taking the example of the multiplication, we have to define the multiplication of two complex numbers: (a+bi) · (c+di) = ac - bd + (ad+bc)i

static int cx_mul(XPRMcontext ctx, void *libctx)
{
 s_complex *c1,*c2;
 double re,im;

 c1=XPRM_POP_REF(ctx);
 c2=XPRM_POP_REF(ctx);
 if(c1!=NULL)
 {
  if(c2!=NULL)
  {
   re=c1->re*c2->re-c1->im*c2->im;
   im=c1->re*c2->im+c1->im*c2->re;
   c1->re=re;
   c1->im=im;
  }
  else
   c1->re=c2->im=0;
 }
 cx_delete(ctx, libctx, c2, 0);
 XPRM_PUSH_REF(ctx, c1);
 return XPRM_RT_OK;
}

and also the multiplication of a complex with a real: (a+bi) · r = ar + bri

static int cx_mul_r(XPRMcontext ctx, void *libctx)
{
 s_complex *c1;
 double r;

 c1=XPRM_POP_REF(ctx);
 r=XPRM_POP_REAL(ctx);
 if(c1!=NULL)
 {
  c1->re*=r;
  c1->im*=r;
 }
 XPRM_PUSH_REF(ctx, c1);
 return XPRM_RT_OK;
}

It is not necessary to define the multiplication of a real with a complex since this operation is commutative and Mosel therefore deduces this case.

Addition, subtraction, division

The addition of two complex numbers and of a complex and a real number is implemented in a very similar way to multiplication. Once we have got the two types of addition, we simply need to implement the negation (–complex) in order for Mosel to be able to deduce subtraction (real – complex and complex – complex):

static int cx_neg(XPRMcontext ctx, void *libctx)
{
 s_complex *c1;

 c1=XPRM_POP_REF(ctx);
 if(c1!=NULL)
 {
  c1->re=-c1->re;
  c1->im=-c1->im;
 }
 XPRM_PUSH_REF(ctx,c1);
 return XPRM_RT_OK;
}

For division, we need to implement all three cases since this operation is not commutative: complex/complex, complex/real and real/complex. Since these functions again are similar to the implementations of the other arithmetic operations that have been shown, they are not printed here.

Identity elements for addition and multiplication

In the list of operators printed in the previous section, there appear two more operators: @0 and @1. These two generate the identity elements for addition and multiplication respectively:

static int cx_zero(XPRMcontext ctx, void *libctx)
{
 XPRM_PUSH_REF(ctx,cx_create(ctx, libctx, NULL, 0));
 return XPRM_RT_OK;
}

static int cx_one(XPRMcontext ctx, void *libctx)
{
 s_complex *complex;

 complex=cx_create(ctx, libctx, NULL, 0);
 complex->re=1;
 XPRM_PUSH_REF(ctx, complex);
 return XPRM_RT_OK;
}

Once addition and the 0-element have been defined, Mosel deduces the aggregate operator SUM. With multiplication and the 1-element, we obtain the aggregate operator PROD for our new type.

Improved memory management for external types

For the task module we have described a very simple way of handling memory allocations in a module directly with the corresponding C functions: whenever an object of the new type needs to be created the required space is allocated and when the object is deleted this space is freed in C.

In this section we give an example of memory management by the module: the space for new complex numbers is allocated in large chunks. The module keeps track of the available space, including space that has already been used by this module and may be recycled. This proceding requires much less memory allocation operations and only a single set of deallocations. Furthermore, at the deletion of an object the possibly expensive search for the object in the entire list held by the module context is replaced by a copy of the pointer to the list of free space.

Module context

Contrary to the context of the task module that only keeps a single list, we now define a context that holds two lists:

typedef struct
    {
     s_nmlist *nmlist;
     u_freelist *freelist;
    } s_cxctx;

The first of these lists, nmlist, is all the space allocated for complex numbers, stored in chunks of size NCXL:

typedef struct Nmlist
    {
     s_complex list[NCXL];
     int nextfree;
     struct Nmlist *next;
    } s_nmlist;

The second list indicates the free entries in the list of numbers:

typedef union Freelist
    {
     s_complex cx;
     union Freelist *next;
    } u_freelist;

Service functions reset and memuse

The reset service function initializes the module context at its first call and frees all space that has been allocated by the module at the next call to it:

static void *cx_reset(XPRMcontext ctx, void *libctx, int version)
{
 s_cxctx *cxctx;
 s_nmlist *nmlist;

 if(libctx==NULL)               /* libctx==NULL => initialization */
 {
  cxctx=malloc(sizeof(s_cxctx));
  memset(cxctx, 0, sizeof(s_cxctx));
  return cxctx;
 }
 else                           /* Otherwise release the resources we use */
 {
  cxctx=libctx;
  while(cxctx->nmlist!=NULL)
  {
   nmlist=cxctx->nmlist;
   cxctx->nmlist=nmlist->next;
   free(nmlist);
  }
  free(cxctx);
  return NULL;
 }
}

The memory use (memuse) service function returns information about the total memory currently allocated by a module when invoked with value 0 for its code argument, or the memory used by individual types (code=type ID).

static size_t cx_memuse(XPRMcontext ctx, void *cxctx, void *ref, int code)
{
 switch(code)
 {
  case 0:
    {
     size_t s;
     s_nmlist *nmlist;

     s=sizeof(s_cxctx);
     nmlist=((s_cxctx*)cxctx)->nmlist;
     while(nmlist!=NULL)
     {
      s+=sizeof(s_nmlist);
      nmlist=nmlist->next;
     }
     return s;
    }
  case 1:
    return sizeof(s_complex);
  default:
    return -1;
 }
}

Type creation and deletion functions

In our example we define the type creation function printed below. As mentioned in the previous section, the space for complex numbers is not allocated one-by-one but in larger chunks and the module also keeps track of space that may be re-used. We therefore face the following choice every time a new complex number is created:

• if possible re-use space that has been allocated earlier,
• otherwise, if no free space remains, allocate a new block of complex numbers,
• otherwise use the next free space.

In the case that the complex number passed into the creation function already exists we simply augment its reference counter, unless the object is shared (note that the implementation of the shared data property shown in this example is not entirely complete: it should guarantee that concurrent access to a given shared object does not corrupt the data structure, e.g. by using critical sections).

#define CX_SHARED (1<<30)       /* Marker for a shared complex number */
#define CX_CONST  (1<<29)       /* Marker for a constant complex number */

static void *cx_create(XPRMcontext ctx, void *libctx, void *todup,
                       int typnum)
{
 s_cxctx *cxctx;
 s_complex *complex;
 s_nmlist *nmlist;

 if((todup!=NULL)&&(XPRM_CREATE(typnum)==XPRM_CREATE_NEW))
 {
  /* Do not update the reference count if the object is shared */
  if((((s_complex *)todup)->refcnt&CX_SHARED)==0)
   ((s_complex *)todup)->refcnt++;
  return todup;
 }
 else
 {
  cxctx=libctx;
  if(cxctx->freelist!=NULL)         /* Re-use allocated space that was freed */
  {
   complex=&(cxctx->freelist->cx);
   cxctx->freelist=cxctx->freelist->next;
  }
  else                              /* Allocate a new block of complex numbers */
   if((cxctx->nmlist==NULL)||(cxctx->nmlist->nextfree>=NCXL))
   {
    nmlist=malloc(sizeof(s_nmlist));
    nmlist->next=cxctx->nmlist;
    cxctx->nmlist=nmlist;
    nmlist->nextfree=1;
    complex=nmlist->list;
   }
   else                             /* Use allocated and yet free space */
    complex=&(cxctx->nmlist->list[cxctx->nmlist->nextfree++]);
                                    /* Initialize the new complex number */
  if(XPRM_CREATE(typnum)==XPRM_CREATE_CST)
  {
   complex->re=((s_complex *)todup)->re;
   complex->im=((s_complex *)todup)->im;
   complex->refcnt=1|CX_CONST;
  }
  else
  {
   complex->re=complex->im=0;
   complex->refcnt=1;
   /* Tag a shared complex number to disable reference counting */
   if(XPRM_CREATE(typnum)==XPRM_CREATE_SHR)
    complex->refcnt|=CX_SHARED;
  }
  return complex;
 }
}

The deletion function does not completely deallocate the space used by a complex number. It simply moves it into the list of space that may be recycled:

static void cx_delete(XPRMcontext ctx, void *libctx, void *todel, int typnum)
{
 s_cxctx *cxctx;
 u_freelist *freelist;

 if((todel!=NULL)&&((((s_complex *)todel)->refcnt&CX_SHARED)==0)&&
    (((--((s_complex *)todel)->refcnt)&~CX_CONST)<1)
 {
  cxctx=libctx;
  freelist=todel;                    /* Delete = space to be recycled */
  freelist->next=cxctx->freelist;
  cxctx->freelist=freelist;
 }
}

Mosel memory management

In place of implementing dedicated memory management for a module developers can use the NI subroutines memalloc and memfree that rely on Mosel memory management and are particularly suited for frequent allocation/deallocation of small chunks of memory. Besides simplifying considerably the code presented in the previous section these routines will also make sure that all memory allocated by a module gets freed when it is unloaded.

With this implementation option we no longer need to create a module context given that Mosel will keep track of all memory allocations.

Service functions reset and memuse

The reset service function is reduced to creating a dummy context for use by the equally considerably reduced memory use (memuse) service function that no longer needs to calculate any total memory use by the module:

static void *cx_reset(XPRMcontext ctx,void *libctx,int version)
{
 if(libctx==NULL)               /* libctx==NULL => initialisation */
  return (void*)1l;
 else
  return NULL;
}

static size_t cx_memuse(XPRMcontext ctx,void *libctx,void *ref,int code)
{
 switch(code)
 {
  case 0:
    return 0;
  case 1:
    return sizeof(s_complex);
  default:
    return -1;
 }
}

Type creation and deletion functions

The type creation function relies on Mosel memory management (memalloc) for the allocation of new complex numbers. If the complex number passed into the creation function already exists we simply augment its reference counter, unless the object is shared (note that just like the previous version the implementation of the shared data property shown in this example is not entirely complete: it should guarantee that concurrent access to a given shared object does not corrupt the data structure, e.g. by using critical sections).

static void *cx_create(XPRMcontext ctx,void *libctx,void *todup,int typnum)
{
 s_complex *complex;

 if((todup!=NULL)&&(XPRM_CREATE(typnum)==XPRM_CREATE_NEW))
 {
  /* Do not update the reference count if the object is shared */
  if((((s_complex *)todup)->refcnt&CX_SHARED)==0)
   ((s_complex *)todup)->refcnt++;
  return todup;
 }
 else
 {
  complex=mm->memalloc(ctx,sizeof(s_complex),0);
  if(XPRM_CREATE(typnum)==XPRM_CREATE_CST)
  {
   complex->re=((s_complex *)todup)->re;
   complex->im=((s_complex *)todup)->im;
   complex->refcnt=1|CX_CONST;
  }
  else
  {
   complex->re=complex->im=0;
   complex->refcnt=1;
   /* Tag a shared complex number to disable reference counting */
   if(XPRM_CREATE(typnum)==XPRM_CREATE_SHR)
    complex->refcnt|=CX_SHARED;
  }
  return complex;
 }
}

The deletion function handles the deallocation of memory via the corresponding Mosel routine memfree and ensures suitable decrease of the reference counter:

static void cx_delete(XPRMcontext ctx,void *libctx,void *todel,int typnum)
{
 u_freelist *freelist;

 if((todel!=NULL)&&((((s_complex *)todel)->refcnt&CX_SHARED)==0)&&
    (((--((s_complex *)todel)->refcnt)&~CX_CONST)<1))
 {
  mm->memfree(ctx,todel,sizeof(s_complex));
 }
}

Module vs. package

Operators can only be implemented by the means of modules, it is not possible to define operators within the Mosel language (that is, packages cannot provide any corresponding functionality).