1 ex_linalg_8c example
gullradriel edited this page 2026-10-29 20:54:40 +01:00

ex_linalg.c

n_linalg_solve() Solves a 3 by 3 system with a known answer, one whose first pivot is zero and needs the row swap, the 10 by 10 normal equations of a cubic surface fitted to points it passes through exactly, and refuses a singular system, a zero matrix and NULL arguments. ```c

#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;

}


_Filename: ex_linalg.c_