<< 2-Multi-precision basics Introduction Matrix manipulation >>

mpscilab >> mpscilab > Introduction > 3-MPScilab Basics

3-MPScilab Basics

Introduction to MPScilab

The MPS datatype

The MPScilab functionalities revolves around a new datatype made available in Scilab. The MPS datatype stores a full matrix of arbitrary precision floating point values. Precision is the same for the whole matrix and must be set when first creating the MPS matrix. At this moment resizing or changing the precision of an existing MPS matrix is not possible. Automatic growth of the matrix when inserting/setting values outside the current matrix index is also not possible.

The next paragraph quickly explains the internal structure of an MPS matrix. It is not required to understand that information to use the toolbox but it might help explain a few tricky parts of using the toolbox. Knowledge of the internal structure will also help users who wants to extract the maximum performance out of the toolbox.

Internally, an MPS matrix is composed of three parts. The first part, is a Scilab mlist that resides inside the Scilab stack and contains the basic information about a matrix such as its dimensions and precision. This list also contains a reference to the two other pars of the matrix located outside of Scilab's managed memory. The first of the two external storage regions contains the MPFR headers of all the elements of the matrix. Each of those header stores the precision, exponent and sign of an entry. Additionally the header contains a pointer to its mantissa. All of the mantissas are stored in the third part.

The MPS functions

Functions and operations provided by the toolbox can be classed into three distinct types. The first type, the most abundant are the native MPS functions. All of the MPS functions are prefixed with mps_ and contain all the basic and not so basic functions to manipulate MPS matrices. A lot of MPS functions can be accessed by another overloaded function or operator. This second classed simply called the overloaded functions contains the native Scilab functions and operators which are overloaded by this toolbox to operate on MPS matrices. The third and final class called the MPX functions are prefixed with mpx_ and regroup advanced and generally non-mathematic functions. Most of the MPX functions are undocumented and should not be required by most tasks. They are however useful for debugging and provide insight into storage statistics and the inner working of the garbage collector and memory management.

Most MPS functions use a similar signature and can be used with either an explicitly given result operator or an implicit one. For example let's take the mps_add() function which adds two operands together and returns the result. mps_add() is defined as:

rop = mps_add(op1, op2)

In the example above, the result operand rop is on the left size of the function call. This is the natural way of calling a Scilab function. In MPScilab this usage will create a new MPS matrix to hold the result. The created matrix will have the correct dimension for the operation and the precision will be assumed from the input operands. However MPS matrices can be large and initializing a new matrix at every function call may be costly. To alleviate this most MPS functions can be called with an implicit output operand. For example the same mps_add() function can be called like this :

mps_add(rop, op1, op2)

When called like this the result operand is always the first argument on the left. Here rop must be a preinitialized MPS matrix of the correct type and dimension. The various ways to create an MPS matrix are described a little later in this document.

In the documentation input operands are generally called op, op1, op2... while a strictly output operand is called rop. Every right hand side operands of such a function call must exist, be initialized and of a valid size for the operation in question. For the matrix addition example the computed operation would be rop = op1 + op2.

Except for very specific functions most MPS functions will accept either MPS matrices or a Scilab matrix of double which is the native Scilab datatype. The result operand must always be an MPS matrix however. The documentation page will specify the range of possible types a function can accept. Moreover, most functions can perform an in place operation where the result operand is the same as one or more input operands. Again the various functions documentation will specify cases where this is not possible. The Advanced topics section explains in more details the rational and consequences of in-place operations and other possible optimizations.

The overloaded functions

An alternative way of calling most MPS functions is by using their overloaded counterparts. For example the overloaded addition operator can be used instead of calling mps_add() like this:

rop = op1 + op2

The same way the mps_sin() function that computes the element-wise Sine of its argument is defined as such:

rop = mps_sin(op)
mps_sin(rop, op)

The same result can be done by using the Scilab native sin() function like this:

rop = sin(op)

When using an overloaded function the result operand rop is created on the fly with the correct size. This implicit creation is done for every overloaded call, as such there is a performance penalty in doing so. More informations on this is located in the Advanced topics section.

Creating an MPS matrix

Most of the time, to perform arbitrary precision computations one or more MPS matrix must be explicitly created. This is done using an initialization function. The basic initialization function is mps_init() and is defined as below:

rop = mps_init(m, n, prec)

mps_init() will create and initialize an MPS matrix of size m rows by n columns with precision prec in bits. For example to create a square 10x10 MPS matrix with 100 bits of precision called mpsA one would invoke mps_init() like this:

mpsA = mps_init(10, 10, 100)

Alternatively, one can use mps_init2() which creates an MPS matrix out of a Scilab matrix of double. Only the precision needs to be specified in this case, for example to create a matrix with 100 bits of precision out of the pre-defined matrix A:

A = [1 2; 3 4]
mpsA = mps_init2(A, 100)

A third way of creating an MPS matrix is by either using a named matrix generator or one of the constant generating function.s One of the well-known named matrix is the identity matrix which can be created using the mps_eye() function.

Constant generating functions can be called nearly the same way as mps_init() with the size and precision of the desired matrix. As a shortcut when a scalar constant is required they can be called with a single argument specifying the precision of the wanted constant. For example to generate a scalar MPS variable containing the constant pi accurate to 1000 bits of precision:

pi = mps_pi(1000)

Manipulating MPS data

As described earlier MPS matrix can be manipulated using various functions. Often there is more than one way of performing the same thing, the different choices depends on the situation and a balance between speed and ease of use.

For the first usage example let's view the different ways to perform matrix addition. The canonical example would be to add two MPS matrix using mps_add() like this:

A = [1 2; 3 4]
B = [4 3; 2 1]
mpsA = mps_init2(A, 100)
mpsB = mps_init2(B, 100)

// Perform the addition.
rop = mps_add(mpsA, mpsB)

As described earlier, it is also possible to use an existing matrix to receive the result :

A = [1 2; 3 4]
B = [4 3; 2 1]
mpsA = mps_init2(A, 100)
mpsB = mps_init2(B, 100)
rop = mps_init(2, 2, 100) // Matrix to hold the result.

// Perform the addition.
mps_add(rop, mpsA, mpsB)

Alternatively, in the previous example we could have used the two Scilab matrices A and B directly with the multi-precision result in rop.

A = [1 2; 3 4]
B = [4 3; 2 1]
rop = mps_init(2, 2, 100) // Matrix to hold the result.

// Perform the addition using A and B directly.
mps_add(rop, A, B)

Yet another method to perform the addition is to use the overloaded addition operator.

A = [1 2; 3 4]
B = [4 3; 2 1]
mpsA = mps_init2(A, 100)
mpsB = mps_init2(B, 100)

// Perform the addition using the overloaded operator, no need to define
// rop beforehand.
rop = mpsA + mpsB

Here also we could substitute the MPS operands with Scilab matrices, however when using the overloaded functions at least one of the input must be an MPS matrix if the calculation is to be done with MPScilab. Otherwise a normal Scilab addition would be performed and rop would be a normal matrix.

A = [1 2; 3 4]
B = [4 3; 2 1]
mpsA = mps_init2(A, 100)

// Performs the addition using the overloaded operator, no need to define
// rop beforehand. Here only mpsA is an MPS matrix for the input.
rop = mpsA + B

Right now to plot MPS data it is necessary to convert it back to a Scilab matrix. This can be performed by calling mps_todouble().

A = [1 2; 3 4]
B = [4 3; 2 1]
mpsA = mps_init2(A, 100)
mpsB = mps_init2(B, 100)

// Perform the addition.
rop = mps_add(mpsA, mpsB)
C = mps_todouble(rop) //Convert the result back to a Scilab matrix.

Clearing MPS variables

MPS matrices created either explicitly or implicitly will be automatically cleared when they fall out of use. The exact working of the garbage collector is described in the Advanced topics

.

Since MPS matrices can be very large, only limited by available memory, it is sometime needed to free the memory used by a matrix immediately. This is achieved using the mps_clear() function which removes an MPS from the Scilab stack and immediately reclaims the memory used to store the matrix.

A = [1 2; 3 4]
B = [4 3; 2 1]
mpsA = mps_init2(A, 100)
mpsB = mps_init2(B, 100)

// Perform the addition.
rop = mps_add(mpsA, mpsB)

// Clear mpsA and mpsB since they are no longer needed.
mps_clear(mpsA)
mps_clear(mpsB)

Understanding assignments with MPS matrices

There is sadly one important limitation when using MPScilab. The simple assignment usually considered a copy by Scilab must never be used.

For reference a simple assignment is using the assignment operator = alone in a statement without any other operators like this.

A = [1 2; 3 4]
mpsA = mps_init2(A, 100)

mpsB = mpsA //This is bad...

Doing a simple assignment will result in having more than one named variable in Scilab that refers to the same underlying MPS matrix. For example in the previous example modifying mpsB would have the same effect on mpsA since they are in fact the same matrix. More importantly, clearing one of them would render the other variable invalid and cause unpredictable results. MPS functions will generally warn if an invalid matrix is passed to them, however a more tricky side effect of the simple assignment is that some operation will corrupt the matrix when accessed with its alias. An example of that is the transpose operation. Doing a transpose on mpsB for the last example would corrupt the matrix viewed by mpsA and using mpsA again would likely result in a crash/segfault.

This limitation is both a result of the way Scilab is designed and the internal representation used by MPScilab and at the moment there is no possible workaround.

In a similar fashion, assignment to an existing variable without using the insertion operator () will overwrite the variable. There is no way for a function to query the output arguments as such a new MPS matrix is created every time. For example:

A = [1 2; 3 4]
mpsA = mps_init2(A, 100)

mpsA = mpsA + mpsA

In the previous example the variable mpsA is added to itself and the result is assigned to itself. However the storage area of mpsA won't be reused and a new matrix will be created to store the result. The original variable is overwritten and no longer visible from within Scilab. This is entirely acceptable and is the same behavior seen with Scilab variables. In MPScilab however the memory used by overwritten variable will not be reclaimed immediately but on the next run of the garbage collector.

This is why MPS functions can have the result operand in the input operand list. Since creating and initializing very large matrices is slow and consumes a significant amount of memory, reusing variables is faster and more efficient.

<< 2-Multi-precision basics Introduction Matrix manipulation >>