Skip to content

cisstNumerical tutorial

Jump to bottom Edit New page
Anton Deguet edited this page May 6, 2014 · 5 revisions

Table of Contents generated with DocToc

Introduction

These examples provide a quick introduction to the features of cisstNumerical. The code can be found in the git repository under cisstNumerical/examples/tutorial. To compile your own code, remember to include cisstNumerical.h.

cisstNumerical contains native functions as well as wrappers for existing numerical routines. As many algorithms can be found in FORTRAN, a significant part of cisstNumerical interfaces with FORTRAN routines. This has numerous consequences listed in section 5.

An illustrative example: SVD

The Singular Value Decomposition is a common algorithm and the cisst implementation illustrates many features of the cisstNumerical FORTRAN wrappers. The goal of SVD is to find the decomposition of a matrix A such as A = U * Sigma * V where both U and V are orthonormal and S is a diagonal matrix composed of singular values.

Most of the FORTRAN routines we are using will not allocate any memory nor check that the parameters are valid (see also section 5). Our wrappers not only check the size of the parameters to verify that enough memory has been allocated but they also provide some flexible mechanisms to allocate the required memory.

For the SVD, the underlying FORTRAN routine requires an input matrix A, two matrices to store U and Vt, a vector for the singular values S and a workspace for temporary variables (a.k.a. scratch space). From now on, we will call the two matrices U and Vt and the vector S the output.

Finally, since cisstVector supports both fixed size and dynamic vectors and matrices, cisstNumerical provides different wrappers and classes for each type of memory allocation. The examples we are providing illustrate different possible configurations, i.e. who allocates memory for what and how.

Example Type Output allocation Workspace allocation
2.1 Dynamic User, manually
2.2 Dynamic User, manually nmrSVD function
2.3 Dynamic User, manually User with method WorkspaceSize
2.4 Dynamic User with nmrSVDDynamicData
2.6 Fixed size User, manually
2.7 Fixed size User with nmrSVDFixedSizeData

2.1 Using nmrSVD with user allocated containers

This example shows how to use the function nmrSVD with a dynamic matrix.

void ExampleSVDUserOutputWorkspace(void) {
    const unsigned int size = 6;
    // fill a matrix with random numbers
    vctDynamicMatrix<double> A(size, size);
    vctRandom(A, -10.0, 10.0);
    // create matrices (U, Vt) and vector (S) for the result
    vctDynamicMatrix<double> U(size, size);
    vctDynamicMatrix<double> Vt(size, size);
    vctDynamicVector<double> S(size);
    // now, create a workspace of the right size
    // we will explain later why this is needed
    vctDynamicVector<double> workspace(size * 10);
    // and we can finally call the nmrSVD function
    // using a copy of A because nmrSVD modifies the input
    vctDynamicMatrix<double> Acopy = A;
    try {
        nmrSVD(Acopy, U, S, Vt, workspace);
    } catch(...) {
        std::cout << "An exception occured, check cisstLog.txt." << std::endl;
    }
    // display the result
    std::cout << "U:\n" << U << "\nS:\n" << S << "\nV:\n"
              << Vt.TransposeRef() << std::endl;
}

In this example, we have used a workspace 10 times bigger than the initial matrix which is large enough. Since the size of the workspace can be determined automatically, cisstNumerical also provides an overloaded version of nmrSVD which doesn't require a workspace.

Using nmrSVD without specifying a workspace

void ExampleSVDImplicitWorkspace(void) {
    const unsigned int size = 6;
    // fill a matrix with random numbers
    vctDynamicMatrix<double> A(size, size);
    vctRandom(A, -10.0, 10.0);
    // create matrices (U, Vt) and vector (S) for the result
    vctDynamicMatrix<double> U(size, size);
    vctDynamicMatrix<double> Vt(size, size);
    vctDynamicVector<double> S(size);
    // and we can finally call the nmrSVD function
    // using a copy of A because nmrSVD modifies the input
    vctDynamicMatrix<double> Acopy = A;
    try {
        nmrSVD(Acopy, U, S, Vt);
    } catch(...) {
        std::cout << "An exception occured, check cisstLog.txt." << std::endl;
    }
}

This is easier to use, but one has to remember that a workspace is created dynamically by nmrSVD, i.e. every time the function is called some memory is allocated and released.

This behavior might not suit everyone, therefore cisstNumerical provides a couple of methods to ease the allocation of the workspace and the output matrices and vectors. All these methods are declared within the scope of a class called "data". For SVD, we have two different classes available, nmrSVDDynamicData and nmrSVDFixedSizeData.

Using nmrSVDDynamicData::WorkspaceSize 

void ExampleSVDWorkspaceSize(void) {
    const unsigned int size = 6;
    // create the input matrix with the correct size
    vctDynamicMatrix<double> A(size, size);
    // now, create a workspace of the right size
    vctDynamicVector<double> workspace;
    workspace.SetSize(nmrSVDDynamicData::WorkspaceSize(A));
    // Allocate U, Vt, S and use the workspace for nmrSVD ...
}

This method simplifies the allocation of the workspace but doesn't solve another problem that we have ignored so far: If the input matrix is not square, the size of the different output containers are a bit trickier to determine, i.e. if the input matrix is m x n then U should be m x m, Vt n x n and the vector S should be min(m, n) long. This is not awfully difficult but still requires some extra attention from the caller. To facilitate the user's work, it is possible the use the class nmrSVDDynamicData to allocate not only the workspace but the output as well.

Using nmrSVDDynamicData to allocate everything

void ExampleSVDDynamicData(void) {
    // fill a matrix with random numbers
    vctDynamicMatrix<double> A(10, 3, VCT_COL_MAJOR);
    vctRandom(A, -10.0, 10.0);
    // create a data object
    nmrSVDDynamicData svdData(A);
    // and we can finally call the nmrSVD function
    vctDynamicMatrix<double> Acopy = A;
    nmrSVD(Acopy, svdData);
    // display the result
    std::cout << "A:\n" << A
              << "\nU:\n" << svdData.U()
              << "\nS:\n" << svdData.S()
              << "\nV:\n" << svdData.Vt().TransposeRef() << std::endl;
}

In this example we have declared a data object based on the input matrix A. The constructor of nmrSVDDynamicData allocates the required memory for all the output containers (U, Vt and S) as well as the workspace. Another overloaded version nmrSVD takes the matrix A and the object svdData to perform the singular value decomposition.

This data object performs a one time memory allocation and can be used multiple times (in a loop for example) without performing any new memory allocation.

Using nmrSVDDynamicData::UpdateMatrixS

Once the decomposition has been performed, nmrSVD stores all the singular values in decreasing order in a vector. This might be convenient for some but one might need a diagonal matrix instead. To update this matrix, the class nmrSVDDynamicData provides the method UpdateMatrixS.

void ExampleSVDUpdateMatrixS(void) {
    // fill a matrix with random numbers
    vctDynamicMatrix<double> A(5, 7, VCT_COL_MAJOR);
    vctRandom(A, -10.0, 10.0);
    // create a data object
    nmrSVDDynamicData svdData(A);
    // and we can finally call the nmrSVD function
    vctDynamicMatrix<double> Acopy = A;
    nmrSVD(Acopy, svdData);
    // compute the matrix S
    vctDynamicMatrix<double> S(5, 7);
    nmrSVDDynamicData::UpdateMatrixS(A, svdData.S(), S);
    // display the initial matrix as well as U * S * V
    std::cout << "A:\n" << A
              << "\nU * S * Vt:\n"
              << svdData.U() * S * svdData.Vt() << std::endl;
}

Note that the method UpdateMatrixS is a static method which can be called even if no data object has been created. Since the method is static, it needs the input matrix A to determine the correct size of the matrix S.

Using fixed size matrices without a data object

void ExampleSVDFixedSize(void) {
    // fill a matrix with random numbers
    vctFixedSizeMatrix<double, 5, 5> A, Acopy;
    vctRandom(A, -10.0, 10.0);
    Acopy = A;
    // create U, S, Vt and a workspace
    vctFixedSizeMatrix<double, 5, 5> U, Vt;
    vctFixedSizeVector<double, 5> S;
    vctFixedSizeVector<double, 50> workspace;
    // and we can finally call the nmrSVD function
    nmrSVD(Acopy, U, S, Vt, workspace);
    // display the result
    std::cout << "U:\n" << U << "\nS:\n" << S << "\nV:\n"
              << Vt.TransposeRef() << std::endl;
}

The example is very similar to the first one, i.e. the user has to create all the containers with the correct sizes before calling nmrSVD.

As for the dynamic matrices, nmrSVD has an overloaded version which doesn't require a user allocated workspace. For the fixed size matrices, the function nmrSVD will create a workspace on the stack if none was provided (not a dynamic memory allocation as seen for dynamic matrices!).

Using fixed size matrices with nmrSVDFixedSizeData 

void ExampleSVDFixedSizeData(void) {
    // fill a matrix with random numbers
    vctFixedSizeMatrix<double, 5, 7, VCT_COL_MAJOR> A, Acopy;
    vctRandom(A, -10.0, 10.0);
    Acopy = A;
    // create a data object
    typedef nmrSVDFixedSizeData<5, 7, VCT_COL_MAJOR> SVDDataType;
    SVDDataType svdData;
    // and we can finally call the nmrSVD function
    nmrSVD(Acopy, svdData);
    // compute the matrix S
    SVDDataType::MatrixTypeS S;
    SVDDataType::UpdateMatrixS(svdData.S(), S);
    // display the initial matrix as well as U * S * V
    std::cout << "A:\n" << A
              << "\nU * S * Vt:\n"
              << svdData.U() * S * svdData.Vt() << std::endl;
}

The interface of the nmrSVDFixedSizeData is pretty much the same as nmrSVDDynamicData except that the size and storage order are now specified using template parameters. To simplify our example, we introduced a type for DataType. This approach is strongly recommended whenever one uses the cisst fixed size vectors and matrices.

Functions with data object

Besides the function nmrSVD, cisstNumerical includes more numerical functions which can be used with either a data object or some vectors and matrices provided by the caller.

FORTRAN based functions

The cisstNumerical FORTRAN wrappers are all written using the approach used for nmrSVD and share the different properties listed for nmrSVD.

nmrInverse

This function computes the inverse of a matrix using an LU decomposition. It can be used for dynamic and fixed size matrices of any storage order. Nevertheless, for fixed size matrices of size 2, 3 or 4, we recommend to use nmrGaussJordanInverse.

void ExampleInverse(void) {
    // Start with a fixed size matrix
    vctFixedSizeMatrix<double, 6, 6> A, AInverse;
    // Fill with random values
    vctRandom(A, -10.0, 10.0);
    AInverse = A;
    // Compute inverse and check result
    nmrInverse(AInverse);
    std::cout << A * AInverse << std::endl;

    // Continue with a dynamic matrix
    vctDynamicMatrix<double> B, BInverse;
    // Fill with random values
    B.SetSize(8, 8, VCT_COL_MAJOR);
    vctRandom(B, -10.0, 10.0);
    BInverse = B;
    // Compute inverse and check result
    nmrInverse(BInverse);
    std::cout << B * BInverse << std::endl;
}

In this example, we used the overloaded version of nmrInverse which doesn't require a data object. This is possible since the data object doesn't provide any useful information or result.

As for most wrappers, using the function nmrInverse without a data object is not optimal if the function is going to be called multiple times. To optimize the memory allocation, one should use nmrInverseFixedSizeData or nmrInverseDynamicData.

nmrLU

The goal of LU is to find the factorization of a matrix A such as A = L x U where U is an upper matrix and V is a lower matrix.

void ExampleLUDynamicData(void) {
    // fill a matrix with random numbers
    vctDynamicMatrix<double> A(5, 7, VCT_COL_MAJOR);
    vctRandom(A, -10.0, 10.0);
    // create a data object
    nmrLUDynamicData luData(A);
    // and we can finally call the nmrLU function
    vctDynamicMatrix<double> Acopy = A;
    nmrLU(Acopy, luData);
    // LAPACK routine store the LU in input A and use
    // a vector to store the permutations P
    vctDynamicMatrix<double> P, L, U;
    P.SetSize(nmrLUDynamicData::MatrixPSize(A));
    L.SetSize(nmrLUDynamicData::MatrixLSize(A));
    U.SetSize(nmrLUDynamicData::MatrixUSize(A));
    nmrLUDynamicData::UpdateMatrixP(Acopy, luData.PivotIndices(), P);
    nmrLUDynamicData::UpdateMatrixLU(Acopy, L, U);
    std::cout << "A:\n" << A
              << "\nP * L * U:\n" << (P * L * U) << std::endl;
}

It is important to notice that in this example we explicitly created the input using VCT_COL_MAJOR. As it is, nmrLU doesn't support the row first storage order.

Besides this constraint, the LU decomposition routine provided by LAPACK stores the result in one single matrix, replacing the input. This is perfectly good for most applications but one can also use the helper methods of nmrLUDynamicData to determine the size and compute the matrices P, L and U.

nmrPInverse

This function actually relies on nmrSVD. The corresponding data object nmrPInverseDynamicData and nmrPInverseFixedSizeData allocate a workspace large enough for nmrSVD.

Native cisst functions

The cisst native functions are more flexible than the FORTRAN wrappers mostly because the restrictions regarding the storage order and the compactness are lifted. The elements might also be different, i.e. one can use single precision floating point numbers if this makes any sense for his/her application.

nmrIsOrthonormal

In this example, we are using nmrSVD to create a couple of orthonormal matrices.

void ExampleIsOrthonormal(void) {
    // fill a matrix with random numbers
    vctDynamicMatrix<double> A(5, 7);
    vctRandom(A, -10.0, 10.0);
    // create a workspace and use it for the SVD data
    vctDynamicVector<double>
	workspace(nmrSVDDynamicData::WorkspaceSize(A));
    nmrSVDDynamicData svdData(A, workspace);
    // we can call the nmrSVD function
    vctDynamicMatrix<double> Acopy = A;
    nmrSVD(Acopy, svdData);
    // check that the output is correct using our workspace
    if (nmrIsOrthonormal(svdData.U(), workspace)) {
        std::cout << "U is orthonormal" << std::endl;
    }
    // same with dynamic creation of a workspace
    if (nmrIsOrthonormal(svdData.Vt())) {
        std::cout << "Vt is orthonormal" << std::endl;
    }
}

This example demonstrates two different ways to use the function nmrIsOrthonormal, one with a user defined workspace and one with no workspace at all (i.e. the function will allocate and free memory on the fly).

Please note in this example how we created a single workspace used by different routines. This is very convenient to avoid any unnecessary memory allocation but one must make sure that this workspace is not being used by two different threads.

It is also possible to create a data object for this problem (see nmrIsOrthonormalDynamicData and nmrIsOrthonormalFixedSizeData).

Others

nmrGaussJordanInverse

The Gauss Jordan inverse methods are implemented for fixed size matrices 2 x 2 , 3 x 3 and 4 x 4. These functions, compiled in release mode, are faster than their FORTRAN counterparts.

void ExampleGaussJordanInverse(void) {
    vctFixedSizeMatrix<double, 4, 4> A, AInverse;
    vctRandom(A, -10.0, 10.0);
    bool nonSingular;
    double tolerance = 10E-6;
    // call nmrGaussJordanInverse4x4
    nmrGaussJordanInverse4x4(A, nonSingular, AInverse, tolerance);
    if (nonSingular) {
	std::cout << "A * AInverse:\n" << A * AInverse << std::endl;
    } else {
	std::cout << "A is a singular matrix" << std::endl;
    }
}

Note that since these functions are fully implemented using the cisst package, any storage order or stride can be used (i.e. the matrices don't need to be compact).

FORTRAN specifics

Compilation

Most of the FORTRAN routines we are using come from the on-line code repository netlib.org. There is no standard binary distribution of these routines therefore we decided to provide a binary version of these routines (library and header files). We have two different versions for historical reasons:

  • CNetlib: This is the oldest version, now deprecated. The main default of this version was that it is not thread safe.

  • cisstNetlib: This version is based on LAPACK routines and is thread safe. We strongly recommend to use this version.

The cisstNumerical API is the same (i.e. your code will be the same) for both binary distributions. Remember that you need to configure your build with the CMake option CISST_HAS_CISSTNETLIB turned ON.

For more details and to download these libraries, see http://github.com/jhu-cisst/cisstNetlib.

Common properties

All our wrappers for FORTRAN routines share the following properties:

The default storage order for matrices is column first in FORTRAN while it is row first in C/C++. Since cisstVector supports both formats, the user has to remember to create his matrices column first (using VCT_COL_MAJOR). This is the default but there are some exceptions. For example, nmrSVD can be used with any storage order. In the case, cisstNumerical uses the fact that changing the storage order is similar to a transpose. For SVD, the problem $A^t = (U * \Sigma * V)^t = V^t * \Sigma^t * U^t$ is basically the same as $A = U * \Sigma * T$. Whenever this is possible, classes and functions of cisstNumerical are implemented to support both storage orders (either row or column major).

Most FORTRAN routines were not written with the concept of stride in mind. This means that all matrices and vectors which are finally used by a FORTRAN routine must be compact (i.e. use a contiguous block of memory).

Most LAPACK routines, will modify the input to avoid unnecessary memory allocation. Since cisstNumerical has been designed to avoid implicit memory allocation and copies as well, it is the caller's responsibility to create a copy of the input for future use.

These functions can only operate on matrices and vectors of doubles. This is because this function is actually a wrapper to a LAPACK routine which requires double precision floating point numbers. For the integers (e.g. vector of pivot indices), FORTRAN uses the equivalent of a C/C++ long int. To enforce this and remind the caller of this subtlety, the cisstNumerical interface defines and uses F_INTEGER. If the matrices or vectors provided by the user are not correct (size, storage order, compact), an exception will occur (std::runtime_error). Since these exceptions are logged, the user might want to look at the file cisstLog.txt if his/her application quits unexpectedly.

JHU - ERC-CISST - LCSR

Clone this wiki locally