Go to the previous, next chapter.

PROFIL/BIAS Extension Package

After the first official version of PROFIL/BIAS has been released, several extensions including a set of test matrices, general linear lists, automatic differentiation, and more have been added as a separate package to the PROFIL/BIAS distribution. Some examples and test programs are also contained in this package.

This is a brief description of these extensions and their usage. It is not intended to explain the specific implementation details of the extensions.

In an extra section we will give some hints for the adaptation of the BIAS IEEE version to non-supported architectures. For this section some experience with C and good knowledge about the architecture (that is how to switch the rounding mode of the floating point unit) is required.

Extension Package Installation

The installation of the PROFIL/BIAS extensions is similar to the installation of PROFIL/BIAS. As a prerequisite it is assumed that PROFIL and BIAS are properly installed. See section Installation, for details.

To install the extensions, change to the PROFIL directory and copy the package file into this directory. The package file is called profext.tar.Z under Unix systems and profext.tgz under DOS or OS/2. To unpack the compressed package archive, type

zcat profext.tar.Z | tar xf -

under Unix or

gzip -cd profext.tgz | tar xf -

under DOS or OS/2. After that the compressed package (profext.tar.Z or profext.tgz) may be deleted.

Type cd Tools to change to the directory containing the extensions. To build all libraries, type

dmake architecture

where architecture is the same value used for building PROFIL. After the dmake finishes, the following libraries have been built:

ProfilLIN: Linear subroutines.
ProfilNLI: Non-linear subroutines.
ProfilMISC: miscellaneous subroutines.

If _test is appended to the above build command, additionally a set of test routines is also built. These test routines perform a simple check of some of the new functions. Under Unix, they are executed automatically, but under DOS or OS/2 you have to start them explicitly. If the test routines terminate silently, the test succeeded. Otherwise, some information is given.

Library and Include File Installation

Installing the Libraries and Include Files

Until now, the libraries and the include files with the PROFIL/BIAS extensions are contained in the Tools directory. To install them in a different place, please refer to the installation section of PROFIL/BIAS, where the appropriate steps are described for the PROFIL/BIAS distribution. See section Installation, for more information. Do not forget that files ending with .hgen and .Cgen (.hge and .cge under DOS) must be treated like include files.

Due to limitations in file name length on some systems, some of the file names mentioned below may be different on your system. If any differences occur, those can be identified in the file names where the first name in each line corresponds to the name mentioned in this documentation and the second name is the file name on your system. On such systems for example the library ProfilMISC is called prmis instead (i.e. the file name is libprmis.a).

Some Miscellaneous Functions --- `MiscFunctions.h`

This file contains the definitions of some functions that do not fit into any other header file but are useful enough to provide them for end users.

Function: REAL Rand01 (): Returns a pseudo-random number between 0 and 1.

Function: REAL Random (): Returns a pseudo-random number between -1 and 1.

Function: UNSIGNED Rand (): Returns an integer pseudo-random number between 0 and PROFIL_RAND_MAX.

Function: void Randomize (): Initializes the pseudo-random number generator with a new time dependent random sequence.

Function: INT Binom (INT k, INT p): Returns the binomial coefficient k over p.

Function: INT Legendre (INT k, INT p)

Returns the Legendre coefficient k over p. This coefficient is defined as follows:

0: if p divides k.
1: if k is congruent to a square mod p.
-1: otherwise.

Test Matrices --- `TestMatrices.h`

The extension package provides a set of test matrices, which mainly will be used for testing and comparing linear algorithms. Many of them have been taken from the Gregory and Karney book.

To get a quick overview concerning the matrix properties, each of the test matrices has a four character type field in the first line of the header comments. The values used for each of the four characters are:

S: Denotes a symmetric matrix.
PD: Denotes a positive definite matrix.
M: Denotes a M-matrix.

The symbol * in any position denotes the absence of a certain property.

Except symmetry, any mentioned properties of the matrices below (e.g. being positive definite or a given range for the eigenvalues) are not guaranteed due to the finite precision of the floating point arithmetic used to calculate the matrices.

The following test matrices are provided:

Function: MATRIX Lietzke (INT n): Returns the square Lietze matrix of dimension n.

Function: MATRIX Pascal (INT n [, REAL k]): Returns the square Pascal matrix of dimension n. If k is present, the elements in the first column and row are set to k, otherwise to 1. If k is the reciprocal of an integer, the elements of the inverse matrix are integers.

Function: MATRIX Hilbert (INT n): Returns the square finite segment of the Hilbert matrix of dimension n.

Function: MATRIX ExactHilbert (INT n): Same as above, but the matrix is scaled to an integer matrix.
Only matrices which are exactly representable on the computer are returned. Beyond a machine-dependent matrix dimension, the subroutine aborts with an error message.

Function: MATRIX InverseHilbert (INT n): Returns the square finite segment of the inverse Hilbert matrix of dimension n.

Function: MATRIX Lotkin (INT n): Returns the square Lotkin matrix of dimension n.

Function: MATRIX Westlake (INT n): Returns the square Westlake matrix of dimension n.

Function: MATRIX Newman (INT n): Returns the square Newman matrix of dimension n. This matrix is orthogonal an the inverse of the matrix is the matrix itself.

Function: MATRIX Frank (INT n): Returns the square Frank matrix of dimension n.

Function: MATRIX BoothroydMax (INT n): Returns a square matrix of dimension n where the matrix element in the i-th row and j-th column contains max(i,j).

Function: MATRIX Hadamard (INT n): Returns a square matrix of dimension n where the matrix element in the i-th row and j-th column contains the Legendre coefficient of (i+j) and (n+1).

Function: MATRIX Binomial (INT n): Returns the square Binomial matrix of dimension n.

Function: MATRIX Aergerter (INT n): Returns the square Aergerter matrix of dimension n.

Function: MATRIX Todd (INT n): Returns a square matrix of dimension n where the matrix element in the i-th row and j-th column contains the absolute value of (i-j).

Function: MATRIX Milnes (INT n [, VECTOR u]): Returns the square Milnes matrix of dimension n. The optional vector u of dimension (n-1) denotes an additional parameter used for the matrix construction. The matrix is singular, if at least one component of u is equal to 1.

Function: MATRIX Combinatorial (INT n [, REAL x, REAL y]): Returns the square combinatorial matrix of dimension n. If the optional parameters x and y are missing, random numbers are chosen. The diagonal contains (x+y) and all other elements of the matrix are set to y.

Function: MATRIX Cauchy (INT n [, REAL x, REAL y]): Returns the square Cauchy matrix of dimension n. If the optional parameters x and y are missing, random numbers are chosen.

Function: MATRIX Vandermonde (INT n, VECTOR x): Returns the square Vandermonde matrix of dimension n.

Function: MATRIX Boothroyd (INT n): Returns the square Boothroyd matrix of dimension n.

Function: MATRIX InverseBoothroyd (INT n): Returns the square inverse of the Boothroyd matrix of dimension n.

Function: MATRIX RandomM_Matrix (INT n): Returns a square random M-matrix of dimension n.

Function: MATRIX RandomPD (INT n): Returns a square random positive definite matrix of dimension n.

Function: MATRIX RandomSPD (INT n [, REAL lower bound, REAL upper bound]): Returns a square random symmetric positive definite matrix of dimension n with optional specified bounds for the eigenvalues of the generated matrix.

Function: MATRIX RandomTeoplitz (INT n): Returns a square random Toeplitz matrix of dimension n.

Function: MATRIX RandomMatrix (INT n [, INT m): Returns a random matrix with n rows and m columns. If m is not specified, a square matrix is generated.

Function: MATRIX RandomSymmetric (INT n): Returns a square random symmetric matrix of dimension n.

Function: MATRIX RandomPerSymmetric (INT n): Returns a square random persymmetric matrix of dimension n.

Subroutines for Local Minimization

Two methods for finding a local minimum of a function with n variables and one REAL result without using derivatives are contained in the package. The first method is the downhill simplex method due to Nelder and Mead. It is a robust but rather slow method in terms of required function evaluations. The implementation has been taken from the Numerical Receipes in C with some slight modifications. The file NelderMead.h contains the definitions needed to use this method.

The second method is due to Brent and defined in Brent.h. Both methods are available as easy-to-use and as expert version. The latter having more parameters allowing more freedom in customization.

Function: REAL NelderMead (VECTOR & x, VECTOR & diam, REAL tol, REAL abstol, INT & ok, INT niter, REAL (*pfunc)(VECTOR &))

This is the easy-to-use version of the Nelder and Mead minimization method. It returns an approximation of the minimum value found. The parameters are:

x: On input, contains the start point. On output, contains the approximation of the minimum point.
diam: On input, contains the approx. diameter of the box with center x, which contains a minimum point. This parameter is used to construct the start simplex. Unchanged on output.
tol: On input, relative termination tolerance. Unchanged on output.
abstol: On input, absolute termination tolerance. Unchanged on output.
ok: On output, success flag. A value of 1 indicates that the algorithm terminates normally. Otherwise, a value of 0 is used.
niter: On input, denotes the maximal number of iterations to be performed. Unchanged on output.
pfunc: On input, pointer to the function to be minimized. Unchanged on output.

Function: VOID NelderMead_Expert (MATRIX & p, VECTOR & y, REAL tol, REAL abstol, INT & ok, INT niter, PREAL parms, REAL (*pfunc)(VECTOR &))

This is the expert version of the Nelder and Mead minimization method. The parameters are:

p

On input, the rows of the matrix contain the start simplex. The number of rows must be (1+n) where n is the dimension of the problem. On output, contains the last calculated simplex.

y

On input/output, contains the function values of each simplex point given by p.

tol

On input, relative termination tolerance. Unchanged on output.

abstol

On input, absolute termination tolerance. Unchanged on output.

ok

On output, success flag. A value of 1 indicates that the algorithm terminates normally. Otherwise, a value of 0 is used.

niter

On input, denotes the maximal number of iterations to be performed. Unchanged on output.

parms

On input, points to an array with 3 elements containing the parameters of the minimization method:

parms[0]: expansion factor for new trial point
parms[1]: shrink factor
parms[2]: expansion factor

Unchanged on output.

pfunc On input, pointer to the function to be minimized. Unchanged on output.

Function: REAL BrentMinimize (VECTOR & p, VECTOR & diam, REAL rtol, REAL atol, REAL (*pfunc)(VECTOR &))

This is the easy-to-use version of Brent's algorithm. It returns an approximation of the minimum value found. The parameters are:

p: On input, contains the start point. On output, contains the approximation of the minimum point.
diam: On input, contains the approx. diameter of the box with center x, which contains a minimum point. This parameter is used to set up some initial parameters. Unchanged on output.
rtol: On input, relative termination tolerance. Unchanged on output.
atol: On input, absolute termination tolerance. Unchanged on output.
pfunc: On input, pointer to the function to be minimized. Unchanged on output.

Function: REAL BrentMinimize_Expert (VECTOR & p, REAL h, REAL rtol, REAL atol, REAL scbd, INT illc, INT ktm, REAL (*pfunc)(VECTOR &))

This is the expert version of Brent's minimization method. It returns an approximation of the minimum value found. The parameters are:

p: On input, contains the start point. On output, contains the approximation of the minimum point.
h: On input, contains an estimation of the distance from p to the minimum point. Unchanged on output.
rtol: On input, relative termination tolerance. Unchanged on output.
atol: On input, absolute termination tolerance. Unchanged on output.
scbd: On input, contains the scale factor. Unchanged on output.
illc: On input, contains the ill-conditioned flag. The value 1 denotes an ill-conditioned problem. Otherwise the value must be 0. Unchanged on output.
ktm: On input, the number of iterations minus one which are allowed to be taken until an improvement must take place. In most cases 1 should do. Unchanged on output.
pfunc: On input, pointer to the function to be minimized. Unchanged on output.

General Linear Singly Linked Lists

The extension package of PROFIL/BIAS supports the creation and usage of linear singly linked lists of arbitrary data types. The general list definitions and implementations are contained in the two generic files LinearList.hgen and LinearList.Cgen.

The creation of the new list data types is done by defining some macros with the names of the data types to be used and including the generic files mentioned above.

The following macros control the generation of the linear lists:

LIST_ELEMENT: Defines the type name of the list elements (e.g. a C++ class name).
LIST: Defines the type name of the list to be generated.
LISTOBJECT: Defines the type name of one list node (only used internally).

After defining these macros, the file LinearList.hgen can be included to generate all definitions of the new list data type. By including the file LinearList.Cgen (the file LinearList.hgen must have been included before), all necessary routines for the list implementation are generated.

Example

As an example, we consider a list of real vectors (data type VECTOR). The new list data type should be named VECTORLIST. With VECTORLISTOBJECT we want to denote the list node data type.

To create the definitions (headers) and implementations (subroutines) of the new list type, use the following code fragment:

#define LIST_ELEMENT VECTOR #define LIST VECTORLIST #define LISTOBJECT VECTORLISTOBJECT

#include "LinearList.hgen"

#include "LinearList.Cgen"

#undef LISTOBJECT #undef LIST #undef LIST_ELEMENT

The #undef commands may be necessary, if more than one list type is defined. After the generation step, you can define a new list li of vectors by

VECTORLIST li;

Sorted Lists

The general linear list package also supports sorted lists where the sorting criterion is given by a comparison function provided by the user. The sorting is performed when inserting a new list element, i.e. the new element is inserted into an already sorted list.

The comparison function used for the sorted insert must have the following prototype

INT SortCompare (LIST_ELEMENT & e1, LIST_ELEMENT & e2);

where the name of the function can be chosen arbitrarily. The return value of the comparison function must be 1, if the element given by e1 should precede the element given by e2 in the list. Otherwise, a value of 0 must be returned. To create the list considered in the above example as sorted list with the comparison function called SortComp, the list is defined as:

VECTORLIST li(SortComp);

The elements are inserted into this list by using the <<= operator (see below).

Supported List Functions

Function: LIST_ELEMENT First (LIST & li): Returns the first element of the list and sets the first list element as current.

Function: LIST_ELEMENT Next (LIST & li): Returns the next element of the list and sets this element as current list element.

Function: LIST_ELEMENT Last (LIST & li): Returns the last element of the list.

Function: LIST_ELEMENT Current (LIST & li): Returns the current element of the list.

Function: VOID RemoveCurrent (LIST & li): Removes the current list element from the list and sets the new current element to the following element.

Function: BOOL Finished (LIST & li): Returns TRUE if the current element is undefined, i.e. if the end of the list is reached. Otherwise, FALSE is returned.

Function: BOOL IsEmpty (LIST & li): Returns TRUE if the list is empty. Otherwise, FALSE is returned.

Function: INT Length (LIST & li): Returns the current number of list elements.

Function: INT MaxLength (LIST & li): Returns the maximal number of list elements since the last call of ResetLength or the definition of the list variable.

Function: VOID ResetLength (LIST & li): Resets the maximal list length to the current length.

To append the element e at the end of a list li, use

li += e

To prepend the element e at the beginning of a list li, use

li *= e

The element e is inserted into the sorted list li by

li <<= e

and the first element of a list li can be removed by li--.

A complete list can be written to output by using the << operator, e.g.

cout << li;

Automatic Differentiation

For functions with n variables and one REAL or INTERVAL result, an automatic computation of the gradient and optionally the Hessian value can be performed during the evaluation of the functional expression in parallel to the computation of the function value.

If only the first derivative (the gradient) has to be computed, the necessary header files are (depending whether you have a real or interval valued function) Gradient.h or IntervalGradient.h.

If the Hessian is also needed, the header files are AutoDiff.h or IntervalAutoDiff.h.

With automatic differentiation, new data types which contain the current function value as well as the current derivative values, do exist. Their names depend on the header file used

GRADIENT: Is defined by Gradient.h
INTERVAL_GRADIENT: Is defined by IntervalGradient.h
AUTODIFF: Is defined by AutoDiff.h
INTERVAL_AUTODIFF: Is defined by IntervalAutoDiff.h

In the following, we will use always AUTODIFF as new type name.

To use automatic differentiation, you first have to define your independent variable by creating a new variable of the type AUTODIFF containing the values of your independent variable. After that, you write your functional expression as you would do in the real case, using variables of the type AUTODIFF instead of REAL or INTERVAL for your intermediate results. Last, you can assign your result (either the function value, the gradient value, or the Hessian value).

Components of a variable of the type AUTODIFF can be accessed by the following functions:

Function: REAL FunctionValue (AUTODIFF a): Returns a reference to the function value. For interval functions, the return type is INTERVAL, instead.
As an alternative, a.fkt() can be used.

Function: VECTOR GradientValue (AUTODIFF a): Returns a reference to the gradient value. For interval functions, the return type is INTERVAL_VECTOR, instead.
As an alternative, a.grd() can be used.

Function: MATRIX HessianValue (AUTODIFF a): Returns a reference to the Hessian value. For interval functions, the return type is INTERVAL_MATRIX, instead.
As an alternative, a.hessian() can be used.

In most cases, the above functions are not needed, because the conversion is performed automatically. For example:

VECTOR x(2), y(2);

x(1) = 1.0; x(2) = 2.0;

GRADIENT xg(x); // x resp. xg is the independent variable
  
y = 2.0 * xg(1) * xg(2); // automatic call of GradientValue(...)

Caution: This automatic conversion will fail if the expression does not depend on the independent variable. In this case, the compiler will generate an error message about incompatible types. Using an intermediate result of an automatic differentiation type or the explicit assignment of zero to the derivatives will be a work-around. This is not a bug of the automatic differentiation package but due to the limited capabilities of C++.

The file AutoDiffExample.C contains a simple example of how to use the automatic differentiation package. If you want to use other than the predefined set of functions with automatic differentiation, you have to modify the files AutoDiff.hgen and AutoDiff.Cgen. These files contain the generic parts for all automatic differentiation files.

The optional computation of the Hessian is controlled by two variables:

Variable: BOOL AUTODIFF::ComputeHessian: If this variable contains TRUE, the evaluation of the Hessian is performed for AUTODIFF objects. Otherwise, the Hessian is not computed.

Variable: BOOL INTERVAL_AUTODIFF::ComputeHessian: If this variable contains TRUE, the evaluation of the Hessian is performed for INTERVAL_AUTODIFF objects. Otherwise, the Hessian is not computed.

Sample Programs

The subdirectory Examples contains sample programs for some of the new features discussed in this documentation. The examples are kept quite simple to make the main ideas of the new features clear.

All sample programs can simply be build by appending _example to the architecture of the dmake command discussed in the installation section. For example,

dmake xlC_examples

builds all libraries (if not already done) and all sample programs.

Adaptation of BIAS to non-supported Architectures

In this section, we will give some information of how to adapt the distributed BIAS version to other non-supported architectures which support the IEEE-standard with directed roundings. For other architectures, a different version of BIAS is available on request.

The adaptation is performed in 4 steps:

The file Portab.h in the BIAS subdirectory must be modified to define some general information about your system.
Subroutines for switching the rounding mode must be written or taken from the system.
If the routines for switching the rounding mode are given in an extra file, the Makefile has to be modified to include this file into the BIAS library.
If the rounding mode routines are taken from the system or you can use inline assembler statements, the files BiasInt.h and Bias0.c have to be modified.

For the correct automatic compilation of BIAS, your system architecture must be identified at compilation time. For this purpose, most compilers support pre-defined macros which can be used for that purpose. For example, the xlc on RS/6000 architectures running under AIX defines the macro _AIX.

To identify your architecture, you have to modify the file Portab.h. The necessary insertions should be behind the comment starting with

Insert settings for new compilers / machines here!

The most important setting is whether you have a big or little endian architecture on your machine. Big endian architectures are e.g. Sparc, RS/6000. Little endian architectures are e.g. PCs with 80x86. As a help, the following C program may be used to identify your machine:

#include 

union {
  struct { char a; char b; } bytes;
  short word;
} un;

main()
{
  un.bytes.a = 1; un.bytes.b = 2;
  if (un.word == 258)
    printf ("big endian machine\n");
  else if (un.word == 513)
    printf ("little endian machine\n");
  else
    printf ("unknown machine\n");
}

Define the macro __HIGHBYTEFIRST__ for big endian and __LOWBYTEFIRST__ for little endian machines. Additionally, you have to define __32BIT__ or __16BIT__ if your system word length is either 32 or 16 bit. Define __ANSI__ if your compiler supports the ANSI keyword volatile. If your compiler supports the library function memmove with overlapping memory fields and does not define the macro __STDC__, you should define it yourself.

As an example, consider the xlc compiler on RS/6000 architectures running AIX. For this case, the settings to be inserted into Portab.h should be:

#if defined (_AIX)
#ifndef __ANSI__
#define __ANSI__
#endif
/* __STDC__ predefined */
#define __32BIT__
#define __HIGHBYTEFIRST__
#endif

After that, you have to define the subroutines for switching the rounding mode. These subroutines should be called

BiasRoundUp()
BiasRoundDown()
BiasRoundNear()

for the different rounding mode settings. Then the assemblation of the rounding mode switching routines must be inserted into the Makefile of BIAS. For this purpose, you need the name of your compiler, the necessary command line flags and the name of the assembler file, where the rounding mode setting routines are contained. If you use C routines instead of assembler subroutines, you have to change in the Makefile in the line

$(ROUNDOBJ)    : $(ROUND).s

the Extension .s into .c.

The BIAS Makefile is prepared for one new architecture/compiler to be simply added. Search for the line

.ELIF $(YOURCOMP)

and modify the following lines for your purposes, if you consider that the compiler name is defined by the macro CC, the command line flags by the macro CFLAGS and the file with the new rounding mode subroutines is defined by ROUND.

The compilation then can be started by

dmake yourcomp

dmake yourcomp_test

If you can use assembler inline sequences to switch the rounding mode, more modifications which are beyond the scope of this report are necessary. Take a look at the files BiasInt.h and Bias0.c. Look for the keyword __I386__ together with __GNUC__ for a machine which only uses inline assembler sequences. If you provide subroutines as well as inline sequences, look for the keyword _AIX or sparc together with __GNUC__ as example.

If you have adapted the BIAS IEEE version to any non-supported architecture and you think that your adaptation is of general interest, you may send it with all modified and/or new files via email with the subject field PROFIL/BIAS to the author. Please include a declaration that your adaptation may be distributed for non-commercial use free of charge.