By Vadim Engelson and Tommy Persson. Some material is adapted from Anders Gertz' master thesis.
LightMat is a collection of C++ classes that provide storage in memory, access and mathematical operations on arrays or tensors. LightMat currently supports arrays of rank 1,2,3 and 4 and can be extended to higher dimensions if needed. Each array is stored as a C++ object, which can be defined and manipulated by corresponding C++ methods. To use LightMat you need just to write the following lines:
#define LM_ALL
#include "lightmat.h"
at the start of your code and add some object files when linking the application.
If you are interested in more details, read the Implementation notes and Optimization (in a separate page).
Some options are not implemented yet, but will be available in the next versions.
Compilers | Preprocessor definitions set by the compiler and used in code |
WorkShop Compilers 4.2 C++ 4.2 (Sun / Solaris 5.6) | __SUNPRO_CC |
SC4.0 / C++ 4.1 (Sun / Solaris 5.5.1) | __SUNPRO_CC |
GNU C++ V.2.6.3(Sun / Solaris 5.5.1, 5.6) | __GNUG__ |
GNU C++ 2.7.2.1 (Sun / Solaris 5.5.1, 5.6, Linux 2.0.30) | __GNUG__ |
GNU C++ 2.7.2 (Digital UNIX V3.2) | __GNUG__ |
DEC C++ V1.3B-0 DEC OSF/1 (Alpha) | __alpha |
MicroSoft Visual C++ 4.0, Win95/NT4.0 | _WIN32 |
MicroSoft Visual C++ 5.0, Win95/NT4.0 | _WIN32 |
If you use one of these compilers there should be no problems in compilation. Otherwise there is a chance that some small modifications should be introduced in the code. Please, let us know about your modifications.
You download the library from the LightMat home page at http://www.ida.liu.se/~pelab/lightmat.
After un-zipping you can go to the "src" subdirectory. Compilation can take some 4-5 minutes on a 133MHz Pentium processor. You must have 70MB free on the Windows disk where your swap file is located (usually in C:\WINDOWS).
make -f make.unx
nmake -f make.win
The code will not work on 8.3 file systems like Windows 3.1 (long file names are used). As result of compilation all necessary object files in the LightMat library are created.
It is important to test the library before using it on your platform. A number of predefined test cases are included with the library.
make -f make.unx all
nmake -f make.win all
Several executable files are created. It can take 4-5 minutes on a 133MHz
Pentium laptop. The make
utility automatically runs the test. If
there were no assertion violations and no core dumps, then, presumably you can
use the code further, e.g. writing new applications from scratch or modifying
the supplied test cases.
The test cases in subdirectory tst
of the distribution show two
different ways of using the library: using inlined functions or non-inlined
functions.
The execution is identical, but compilation in different.
This variant gives the fastest code, but slow compilation. In your files you specify:
#define LM_ALL
#include "lightmat.h"
You compile your files with the usual compilation flags. Compilation usually
takes one or more minutes. You link your object files with cstring.o
and extwin.o
. For example the file foo.cpp
(the
default extension .cpp is used for MicroSoft Visual C++ files) contains:
#define LM_ALL
#include <lightmat.h>
int main(){
doubleNNN f(5,6,7,3.11);
cout << f;
return 0;
}
Compile with GNU C++ (assuming that LightMat is installed in the directory ???
)
g++ foo.cpp -I../include ../obj/cstring.o ../obj/extwin.o -lm
or Visual C++:
foo.cpp -I..\include ..\obj\cstring.o ..\obj\extwin.o
Note: Both forward slash and backslash can be used as directory separators with
cl
.
This variant gives slower code (although still reasonably fast), but significantly faster compilation. In your files you specify:
#define LM_ALL
#include <lightmat.h>
You compile your files with the usual flags and -DNOT_INLINED
. You
link your files with cstring.o
, extwin.o
and not_inlined.o
.
All the library functions are non inlined. This gives reasonable, but not the
highest performance. Each time you compile your file, it will take less than 1
minute. For example the file foo.cpp
contains:
#define LM_ALL
#include <lightmat.h>
int main(){
doubleNNN f(5,6,7,3.11);
cout <<f;
return 0;
}
Compile with
g++ foo.cpp -I../include -DNOT_INLINED ../obj/cstring.o
../obj/extwin.o ../obj/not_inlined.o -lm
or
cl foo.cpp -I..\include -DNOT_INLINED ..\obj\cstring.o
..\obj\extwin.o ..\obj\not_inlined.o -lm
If you use only a subset of the LightMat classes you can reduce compilation time
for your code. Use a preprocessor definition #define LM_name
for each class you use. Names of these definitions are given in the
table of classes. Important: This option can be used only with
the flag -DNOT_INLINED
. #define LM_ALL
means that you
use all available classes. For example the file foo1.cpp
contains:
#define LM_N
#include <lightmat.h>
int main(){
doubleN x(5,2.0);
cout << x;
//doubleNN y(3,4);- Forbidden
// This is not allowed because this class requires LM_NN.
// Compiler reports an error.
return 0;
}
Compile with
g++ foo1.cpp -I../include -DNOT_INLINED ../obj/cstring.o
../obj/extwin.o ../obj/not_inlined.o -lm
or
cl foo1.cpp -I..\include -DNOT_INLINED ..\obj\cstring.o
..\obj\extwin.o ..\obj\not_inlined.o
If you compile several files into a single executable, the following holds:
-DNOT_INLINED
flag
-DNOT_INLINED
flagThe library includes implementation of complex numbers. The class for a stand-alone complex number used is called lm_complex. This should not be mixed with <complex> template from Standard C++ library.
The library includes implementation of strings. The class for a stand-alone string used is called string. This is the same as std::string from standard C++ library.
There is no need in a special boolean type - they are implemented via integer types. However all boolean operations (not, and, or) should be defined for all intege-based lightmat types. The value 1 stands for true, 0 for false in boolean operation results.
In the file lightmat.h
(or lightmat_id.h
) there are a
number of preprocessor definitions that can be changed. Those definitions are:
|
|
|
Meaning if undefined |
LIGHTMAT_LIMITS_CHECKING | defined | LightMat does limits checking during calculations. I.e. LightMat will abort computations and dump core if an index is out of bounds. Computation speed is reduced by these checks. | An operation with bad index can cause unpredictable results. |
LIGHTMAT_INIT_ZERO | undefined | LightMat will initialize all array elements to zero when constructing a new object. Computation speed is slightly reduced. | Array elements initially may contain any value. Care should be taken if non-initialized values are used in some other interfaces, like MathLink. |
LIGHTMAT_DONT_USE_BLAS | defined | LightMat uses its own function for double array multiplication | LightMat calls BLAS libraryroutines. The
library should be linked with the application.
(N/A) |
LIGHTMAT_OUTPUT_FUNCS | defined |
The ToStr (conversion to Tools.h++ class RWCString) functions are defined and operator<< functions are defined for all arrays. In version 0.80 std::string is used instead of RWCString. |
(N/A) |
LIGHTMAT_TEMPLATE | undefined | LightMat classes can be used as templates. This option is currently not available due to difficulties with various compilers.(N/A) | LightMat classes are used as traditional C++ classes. |
If preprocessor definitions change, all application files referring to LightMat should be recompiled.
We will need some definitions first:
double
or int
lightN
.
lightNN
.
lightNNN
.
lightNNNN
.There exist 8 different classes for vectors, matrices, tensors of rank 3 and tensors of rank 4. There are 4 universal classes and 4 specially optimized classes. The specially optimized classes can be used in particular application areas: for instance matrix of dimension 3x3 is used in geometry and 4x4 in computer graphics. If you are not sure which classes to use then use universal classes.
Dimension
suffix |
Preprocessor
symbol for inclusion |
Use | Default stack allocation size | Internal class names | Class names | Explanation |
3 | LM_3 | special | 3 |
light3double
|
double3
int3 |
Vector of 3 elements |
4 | LM_4 | special | 4 |
light4double
|
double4
int4 |
Vector of 4 elements |
N | LM_N | universal | 10 |
lightNdouble
|
doubleN
intN lm_complexN stringN |
Vector of arbitrary size |
33 | LM_33 | special | 3 x 3 | light33double
light33int |
double33
int33 |
Matrix 3 by 3 elements |
44 | LM_44 | special | 4 x 4 | light44double
light44int |
double44
int44 |
Matrix 4 by 4 elements |
NN | LM_NN | universal | 10 x 10 | lightNNdouble
lightNNint lightNNlm_complex lightNNstring |
doubleNN
intNN lm_complexNN stringNN |
Matrix of arbitrary size |
NNN | LM_NNN | universal | 5 x 5 x 5 | lightNNNdouble
lightNNNint lightNNNlm_complex lightNNNstring |
doubleNNN
intNNN lm_complexNNN stringNNN |
Tensor of arbitrary size |
NNNN | LM_NNNN | universal | 4 x 4 x 4 x 4 | lightNNNNdouble
lightNNNNint lightNNNNlm_complex lightNNNNstring |
doubleNNNN
intNNNN lm_complexNNNN stringNNNN |
Tensor??? of arbitrary size |
Note on higher ranks: There are no classes for tensors of higher ranks larger than 4. It is possible to add support for that in the future though. Having different classes for vectors, matrices etc is quite natural, and that's what most other packages also have chosen. Other packages sometimes have one class for tensors of high ranks, but that does not exist in LightMat. Having different classes for the different ranks also avoids run-time choices of what algorithm, for the different ranks, that should be used.
Note on internal class names: The class names in the ...???
Array classes for elements of type int
and double
are
implemented. The user declares a variable as follows:
double4 foo
This means that foo
is a vector with 4 elements of type double
.
All available class names are given in the class table
above. The properties of classes with elements of type int
and double
are very similar. For instance, by classes with suffix 3
is
usually meant both classes int3
and double3
.
Extra optimized classes (suffixes 3, 4, 33, 44
) have been
implemented for some specific sizes of vectors and matrices for double
and int. Those are vectors of length 3 and 4, and matrices of
size 3x3 and 4x4. They have, of course, a stack allocated memory area for the
elements. That means that the compiler can generate code which need fewer
references through pointers. Code for calculations is written especially suited
to those sizes. Calculation loops are completely unrolled. This makes these
classes extra fast computationally.
These classes, though, can also interface with the other existing classes in lightmat.
Since the sizes of objects of these classes are partially or totally fixed they cannot dynamically change their sizes as freely as the universal LightMat classes.
All LightMat classes have some basic member functions like constructors, destructors, assignment operators and a few more.
There exist a number of constructors, with different arguments, with which a new object can be created. The reasons for having several different constructors are efficiency and ease of use.
intNN a;
doubleNN b;
stringNNN c;
lm_complexNNNN p;
a
is of type doubleNN
and d
is of type intNNN
:
doubleNN b = a;
intNNN c = d;
doubleNN a(3,5);
intNNNN c(3,5,2,4);
double33 a;
doubleNN c = a;
Such constructors are also called converters. There are converters for the following cases:
intN
to doubleN
)
intN
to complexN
)
doubleN
to complexN
)
int3 |
intN |
int4 |
intN |
int33 |
intNN |
int44 |
intNN |
double3 |
doubleN |
double4 |
doubleN |
double33 |
doubleNN |
double44 |
doubleNN |
double arr[] = { 1.3, 2.6, 3, 4, 5.7, 6 };
doubleNN a(2,3,arr);
The same result can be achieved by writing a loop:
doubleNN a(2,3);
idx=0;
for(inti=1;i<=2;i++)< br> for(intj=1;j<=3;j++)
a(i,j)=arr[idx++];
doubleNN
a(3, 4, 17.7);
lm_complexNN b(3, 5, lm_complex(7.0,6.0));
light3
,
light4
, light33
and light44
). The values
for the elements should be given in row major order. double33 a(1.1,
2, 3.5, 4, 5, 6.9, 7, 8.3, 9);
If values for the elements are not supplied in the constructor, the
initialization behaviour depends on preprocessor definition LIGHTMAT_INIT_ZERO
(see preprocessor definition table above.) If it is
defined, the elements are initialized to zero. Otherwise, the elements are not
initialized. This saves some time. It can be a good idea not todefine LIGHTMAT_INIT_ZERO
if the user knows that objects are never used without explicit initialization.
Functions within the classes that do calculations and need a local object use a
special internal protected constructor that never initializes the elements.
The LightMat constructors automatically allocate dynamic memory on the heap if the pre-allocated memory area (area allocated by default on the stack) is not sufficiently large. The class table above describes the size of memory allocated by default on the stack.
The destructors free any memory which might have been allocated on the heap.
The assignment operators copy values from the elements in one array object to another. The array objects can be of different size, but the rank must be the same. The right-hand operand can however be scalar. The size of the left-hand operand is changed to the size of the right-hand operand if needed. More memory is also allocated if needed.
Note: Memory is not deallocated even if the objects need of memory decreases. The reason is that deallocation and the sometimes following allocation of a smaller memory area takes time. There may also, later in the calculations, once again be a need of a larger memory area in which case another reallocation would be needed once again.
doubleNN a(2,3), b(3,2);
a = b; // a becomes 3x2
a = 5.7; // all elements of a become 5.7
Assignment from a scalar will set all the elements in the object to that value. Read the section about indexing on how to assign a value to an element of array.
All classes have member functions Get
and Set
which
takes a pointer to a memory area as argument. Set
will set the
elements in the object from the values which are stored in the memory area. Get
does the opposite. The elements are stored in row major order (C++ style) in
the memory area. No bounds check is performed. For example:
doubleNN a(5,6,17.7), b(5,6);
double arr[30];
a.Get(arr); // a --> arr
b.Set(arr); // arr --> b
The size of a dimension of an array object can be found with the function dimension()
.
It takes an integer argument which should specify the dimension. The argument
should, for example, be 1 for the number of rows in a matrix, and 2 for the
number of columns. Thus, a a.dimensions(1)
returns the size of the
list dimension of a
, a.dimensions(2)
returns the size
of the second dimension, etc.
The usual operators for equality and inequality exist for all classes. Two objects are considered equal if they have the same size and if the values of all corresponding elements are equal. Only arrays of the same rank can be compared. For example:
if(a == b ?? c != d) { /* something */ }
A function named reshape
is used to promote an object of rank n to
another object of rank n+1. E.g. a vector can be promoted to a matrix and all
the columns in the matrix will then be given the values of the vector. For
example:
doubleN v(3);
doubleNN a;
a.reshape(3,5,v); // a becomes 3x5
The function make_lightN()
is not a method of a class, but a friend
function. You can use it for construction of an array from its components.
Given several arrays of rank n it constructs an array of rank n+1. The number
of components should be specified as the first argument. Obviously, the
dimensions and element types of all arrays should be the same. Not more than 10
parameters can be given in the argument list.
lightNN a;
a=make_lightN(2, make_lightN( 3, 11,12,13), make_lightN( 3, 14,15,16));
This code creates array 2x3 of integers from 11 to 16. Specification is always
given in row-major order. Note: This function is very time and space consuming.
Therefore use Set()
if you are interested in higher performance.
The indexes in LightMat arrays starts at 1, like in Fortran, Matlab and Mathematica, not at 0 like in C. For matrices the first coordinate is usually named "row" and the second "column". Bound check in indexing operations can be turned on or off by preprocessor definitions. There is indexing for individual elements, indexing for array extraction and special indexing for Set().
The index argument to operator ()
and Set
can be of
the following types:
R
.The class R
contains a start index and an end index. Both are
integers. This specifies a range (which includes the indexes). There is a
constructor taking two arguments, that is start index and end index. For
example:
lightNN a(7, 9);
lightN b(3);
b = c(3, R(4, 6));
// Equivalent to b = make_lightN (c(3,4), c(3,5), c(3,6))
Suppose that the range operator is used to index a range in an object where the
dimension of the indexed position is size
. Then the following
trasformations are done:
All -> R(1, -1)
R(imin, imax) -> R(imin+size+1, imax)
if imin
<0
R(imin, imax) -> R(imax+size+1, imax)
if imax
<0After these transformations a limit check is done (if required). If imax == 0
or imin == size+1
then the range is assumed to contain zero
elements and the limit check does not fail. Otherwise it is required that imax-imin+1
> size
.
If the size of the range is zero then the Set
method does not do
anything and the operator()
method return an object with size
zero.
Suppose that a
is intNNNN
and mat
is intNN
and we have the following method call: Set (2, range1, 3, range2, mat)
.
Then if limit checks are done it is checked that the size of range1
is less than mat.dimension (1)
and that the size of range2
is less than mat.dimension (2)
.
Individual elements can be indexed with the parentheses operator. The indexed element can be either an l-value or an r-value. The number of indexes should be equal to the rank of the array.
doubleNN a(5,3);
a(1,1) = a(2,3);
Arrays can be extracted from arrays of higher dimensions using the range class
as index argument for operator()
. The result can only be used as
an r-value, and the indexed elements are copied in the process. The elements
need to be copied since the concept of different views of the same data does
not exist in this design. Note that this is relatively expensive operation. If
the result is used as a l-value, the assignment result is simply discarded.
doubleN v;
doubleNN a(5, 5);
doubleNNN c(7, 8, 9);
v = a(4, All);
v = c(4, 5, All);
a = c(All, 3, All);
a(1, All) = v; // The result is discarded, and matrix a does not change.
double3 w;
m(3) = w + m(2); // This class is exception, m changes.
The method Set()
is used to set a slice of an object to the value
of another object. For example:
doubleN v(3);
doubleNN a(5, 5);
doubleNNN c(7, 8, 9);
a.Set(1, All, v); // Error because a is of size 5x5 and5>3
a.Set(1, R(1, 2), v); // OK, because the range size is 2
SetShape()
This function changes the shape of an array. (Available for universal
classes (doubleN
, ..., doubleNNNN
,) only). If
necessary, new memory is allocated. The elements are not initialized. They can
be accessed by the usual indexing operation.
doubleNN a; // Allocates 10x10, has
shape 0x0
doubleNN b(5,6); // Allocates 10x10, has shape 5x6
a.SetShape(20,40); // Allocates 20x40, has shape 20x40
b(70,30)=1; // Causes error
b.SetShape(70,30); // Allocates 70x30, has shape 70x43
b(70,30)=1; // OK
data()
This function returns the address of internal presentation of the array. (Available for universal classes only) Internally the arrays are stored in column major order. This order is used for interfacing with Fortran routines.
All arrays can be printed out to an output stream by the usual << operator: If the NOPAR preprocessor definition is defined, no parentheses are used in printouts.
doubleNN a;
cout << a;
May produce:
( 5.0 7.0
2.6 7.2 )
Complex numbers are written as 5+7 I, 5-7 I
Strings are written without quotes.
-,+
+, -, pow, Mod
*
+=, -=, *=, /=
/
ElemProduct, ElemQuotient
Cross, OuterProduct, Apply
sign, abs, LightMax, LightMin,
ifloor, iceil, irint, IntegerPart, FractionalPart, sqrt, exp, log, sin, cos,
tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh
A collection of operators and functions is defined for every class in the LightMat library. The operators and function names are overloaded and perform relevant operations on various arguments. All the operators and functions are available as friend functions. Therefore the notation is similar for scalar arguments and for arrays of different ranks.
Operators available for all classes are:
-x
, unary minus
+x
, unary plusOne of the arguments can be scalar. If both are arrays they must be of the same size and shape. This applies to int, double and complex.
x+y
, Elementwise addition.
x-y
, Elementwise subtraction.
pow(x,y)
, Elementwise power. Extended for complex
numbers.
Mod(x,y)
, Elementwise modulo. The definition is extended for real and
complex numbers. If x and y are integers of different sign (x%y+y) will
be the result.This applies for int, double and lm_complex.
xarr*yscalar
or xscalar*yarr
: If one of the
arguments is scalar, all elements of other argument are multiplied
by the scalar. The result has the same rank and shpe as the array.
xarr*yarr
: If both arguments are arrays, the inner
products is computed. (For vector also called dot product.). If the arguments
have rank m and n, the result will be of rank m+n-2. The last dimension of the
first argument should be equal to the first dimension of the second argument.
This is checked at runtime.All the possible variants are given in the table (s is a scalar; other symbols are suffixes of class names):
|
|
|
|
If y is an array, it must have the same dimensions as x. Result has the same dimension as array x. This applies to int, double and complex.
x+=y
, add array or scalar y
to x
.
x-=y
, subtract array or scalar y
from x
.
x*=y
, multiply elements of x
by scalar y
x/=y
, divide elements of x
by scalar y
The result array has the same dimension as the array argument. This applies to int, double and complex.
x/y
, divide elements of array x
by scalar y
x/y
, divide scalar x
by array y
. For
each element a of y
, form an array of elements x/a
.
The x and y are arrays. Result has the same dimension. This applies to int, double and complex.
ElemProduct(x,y)
- element-wise multiplication
ElemQuotient(x,y)
- element-wise divisionCross(x,y)
- defined for int3
and double3 and
lm_complexN
of length 3 only. Returns the same type. More
???.
OuterProduct(x,y)
- defined for classes with suffixes 3, 4, N
.
Returns value of corresponding class with suffix NN
. For example
double3, doubleN
gives a result value of doubleNN
.
Apply(arr,f)
- applies function f to every element of the array.
This applies to int, double, string and complex. Example: double
add(double x) { return x+1; }
double3 a(1.5,1.6,1.7);
double3 b=Apply(a,add);
cout << b; // prints 2.5, 2.6, 2.7
Apply(arr1,arr2,f)
- applies function f(x,y)
to every
element of array arr1
and arr2
element-wise. This
applies to int, double, string and complex. The arrays arr1
and arr2
have the same dimension sizes. Example: int
sum(int x1, int x2){ return x1+x2;}
int4 s(5,6,7,8);
int4 t(10,11,12,13);
int4 u=Apply(s,t,sum); // 15, 17, 19, 21
sign
(returns integer, -1, 0 or 1; defined for integer and
double). Returns complex for complex as z/abs(z)
.
abs
(returns the same type as argument). For complex returns a double
value |z|.
LightMax
returns the maximal element of an array (for integer and double).
LightMin
returns the minimal element of an array (for integer and double).
ifloor, iceil, irint, IntegerPart
(returns integer).
FractionalPart, sqrt, exp, log, sin, cos, tan, asin, acos, atan, sinh, cosh,
tanh, atan2
asinh, acosh, atanh
(returns double, not defined for Win32 implementation)
Note: (Version 0.76) FractionalPart, IntegerPart, Mod, LightMax, LightMin - implemented for universal types (N,NN,...) only.
All complex operations are needed foe universal types (N,NN,...) only.
These operations can be applied to universal array types (N,NN,...) of int, double, lm_complex, string
Joins two arrays of the same rank. All dimensions except the first must be identical.
Joins two arrays of the rank n and n-1. All dimensions except the first dimension of the first array must be identical.
Joins two arrays of the rank n-1 and n. All dimensions except the first dimension of the first array must be identical.
Accepts an array and a range (object of class R). Removes all elements whose first index is in the range.