examples/ex_linalg.c

Functions

	Name
int	main(void )
void	expect_true(const char * what, int ok) report a check

Attributes

	Name
int	failures record one check

Functions Documentation

function main

int main(
    void 
)

function expect_true

static void expect_true(
    const char * what,
    int ok
)

report a check

Attributes Documentation

variable failures

static int failures = 0;

record one check

Source code

#include <math.h>
#include <stdio.h>
#include <string.h>

#include "nilorea/n_common.h"
#include "nilorea/n_linalg.h"
#include "nilorea/n_log.h"

static int failures = 0;

static void expect_true(const char* what, int ok) {
    if (!ok) {
        n_log(LOG_ERR, "FAIL: %s", what);
        failures++;
    }
}

int main(void) {
    set_log_level(LOG_NOTICE);
    {
        /* 2x + y - z = 8, -3x - y + 2z = -11, -2x + y + 2z = -3: x 2, y 3, z -1 */
        double a[9] = {2, 1, -1, -3, -1, 2, -2, 1, 2};
        double b[3] = {8, -11, -3};

        expect_true("a 3 by 3 system", n_linalg_solve(a, b, 3) == TRUE && fabs(b[0] - 2) < 1e-12 &&
                                           fabs(b[1] - 3) < 1e-12 && fabs(b[2] + 1) < 1e-12);
    }
    {
        /* a zero where the first pivot would be: only the swap solves it */
        double a[4] = {0, 1, 1, 0};
        double b[2] = {5, 7};

        expect_true("a zero first pivot, swapped", n_linalg_solve(a, b, 2) == TRUE && b[0] == 7.0 && b[1] == 5.0);
    }
    {
        /* the normal equations of a cubic in u, v fitted to 200 points of
         * a known cubic: the ten coefficients come back */
        const double want[10] = {0.5, -1.25, 2.0, 0.3, -0.7, 0.11, 0.05, -0.02, 0.09, -0.4};
        double a[100], b[10], t[10], worst = 0.0;
        int i = 0, j = 0, k = 0;

        memset(a, 0, sizeof(a));
        memset(b, 0, sizeof(b));
        for (i = 0; i < 200; i++) {
            double u = -1.0 + 2.0 * ((i * 37) % 200) / 199.0, v = -1.0 + 2.0 * ((i * 91) % 200) / 199.0, z = 0.0;

            t[0] = 1, t[1] = u, t[2] = v, t[3] = u * u, t[4] = u * v, t[5] = v * v;
            t[6] = u * u * u, t[7] = u * u * v, t[8] = u * v * v, t[9] = v * v * v;
            for (j = 0; j < 10; j++) z += want[j] * t[j];
            for (j = 0; j < 10; j++) {
                for (k = 0; k < 10; k++) a[j * 10 + k] += t[j] * t[k];
                b[j] += t[j] * z;
            }
        }
        expect_true("a cubic surface's normal equations", n_linalg_solve(a, b, 10) == TRUE);
        for (j = 0; j < 10; j++)
            if (fabs(b[j] - want[j]) > worst) worst = fabs(b[j] - want[j]);
        n_log(LOG_NOTICE, "cubic fit: worst coefficient off by %g", worst);
        expect_true("its coefficients come back", worst < 1e-9);
    }
    {
        double a[4] = {1, 2, 2, 4};
        double b[2] = {1, 2};
        double z[4] = {0, 0, 0, 0};
        double y[2] = {1, 1};

        expect_true("a singular system is refused", n_linalg_solve(a, b, 2) == FALSE);
        expect_true("a zero matrix is refused", n_linalg_solve(z, y, 2) == FALSE);
        expect_true("NULL and n 0 are refused",
                    n_linalg_solve(NULL, y, 2) == FALSE && n_linalg_solve(z, NULL, 2) == FALSE && n_linalg_solve(z, y, 0) == FALSE);
    }
    n_log(LOG_NOTICE, "ex_linalg: %s", failures ? "FAILED" : "ok");
    return failures ? 1 : 0;
}