OpenNL

Solving Linear Systems with OpenNL

OpenNL is a library for solving large sparse linear systems. It includes an easy-to-use API for assembling matrices, and various iterative solvers for symmetric and non-symmetric systems. OpenNL API is declared in geogram/NL/nl.h.

Solving a linear system.

Let us start with a simple example. Suppose you want to solve the following linear system:

   [ 1 2 ] [x]   [5]
   [ 3 4 ] [y] = [6]

First, you need to declare the variables that will store the solution of the system. Then, we create an OpenNL context, and declare the number of variables:

  double x,y;
  nlNewContext();
  nlSolverParameteri(NL_NB_VARIABLES, 2);

Now we can built the linear system. The linear system is built row by row. OpenNL has a state machine, controlled by nlBegin(),nlEnd() calls (that are similar to OpenGL 2.x, except that they build sparse matrices instead of graphic primitives !).

  nlBegin(NL_SYSTEM);
  nlBegin(NL_MATRIX);
  nlBegin(NL_ROW);
  nlCoefficient(0, 1.0);
  nlCoefficient(1, 2.0);
  nlRightHandSide(5.0);
  nlEnd(NL_ROW);
  nlBegin(NL_ROW);
  nlCoefficient(0, 3.0);
  nlCoefficient(1, 4.0);
  nlRightHandSide(6.0);
  nlEnd(NL_ROW);
  nlEnd(NL_MATRIX);
  nlEnd(NL_SYSTEM);

Next, we solve the system and get the solution:

  nlSolve();
  x = nlGetVariable(0);
  y = nlGetVariable(1);

And finally, we delete the OpenNL context.

  nlDeleteContext(nlGetCurrent());

Least squares regression

OpenNL has a special least squares mode, that assembles the normal equation (least squares). Suppose that you have a file that contains X,Y point coordinates on each line, and you want to find the line equation y=ax+b that best fits the data points. Like before, we create an OpenNL context and declare the number of variables. In addition, we activate least squares mode:

   nlNewContext();
   nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
   nlSolverParameteri(NL_NB_VARIABLES, 2);

Then, we read the file and assemble the matrix from the file data:

  FILE* input = fopen("datapoints.dat", "r");
  double X,Y; // current datapoint
  nlBegin(NL_SYSTEM);
  nlBegin(NL_MATRIX);
  while(!feof(input)) {
     fread(input, "%f %f", &X, &Y);
     nlBegin(NL_ROW);
     nlCoefficient(0, X);
     nlCoefficient(1, 1.0);
     nlRightHandSide(Y);
     nlEnd(NL_ROW);
  }
  nlEnd(NL_MATRIX);
  nlEnd(NL_SYSTEM);

Just like in the previous example, we can now solve the system and get the solution:

  nlSolve();
  a = nlGetVariable(0);
  b = nlGetVariable(1);

And do not forget to close the input file and to delete the OpenNL context:

  fclose(input);
  nlDeleteContext(nlGetCurrent());

Least squares regression with locked variables

As in the previous example, we Suppose that we have a file that contains X,Y point coordinates on each line, and we still want to find the line equation y=ax+b that best fits the data points, but this time we have an additional constraint: the slope a of the line should be equal to 1. The OpenNL context is initialized just like before:

  FILE* input = fopen("datapoints.dat", "r");
  double X,Y; 
  nlNewContext();
  nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
  nlSolverParameteri(NL_NB_VARIABLES, 2);

When we build the system, we lock variable 0 (that corresponds to a) and set its value (1.0) as follows:

  nlBegin(NL_SYSTEM);
  nlLockVariable(0);
  nlSetVariable(0, 1.0);
  nlBegin(NL_MATRIX);
  while(!feof(input)) {
     fread(input, "%f %f", &X, &Y);
     nlBegin(NL_ROW);
     nlCoefficient(0, X);
     nlCoefficient(1, 1.0);
     nlRightHandSide(Y);
     nlEnd(NL_ROW);
  }
  nlEnd(NL_MATRIX);
  nlEnd(NL_SYSTEM);

Then the rest is just the same as before:

  nlSolve();
  a = nlGetVariable(0);
  b = nlGetVariable(1);
  fclose(input);
  nlDeleteContext(nlGetCurrent());

Fine tuning

WIP

OpenNL extensions

• SuperLU
• CUDA

WIP

Eigen solver

WIP