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

src/n_linalg.c

Small dense linear systems. More...

Detailed Description

Small dense linear systems.

Author: Castagnier Mickael

Version: 1.0

Date: 28/10/2026

Source code

/*
 * Nilorea Library
 * Copyright (C) 2005-2026 Castagnier Mickael
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
 * implied. See the License for the specific language governing
 * permissions and limitations under the License.
 *
 * SPDX-License-Identifier: Apache-2.0
 */


#include "nilorea/n_linalg.h"

#include "nilorea/n_common.h"

#include <math.h>

int n_linalg_solve(double* a, double* b, size_t n) {
    double scale = 0.0;
    size_t i = 0;

    if (!a || !b || n == 0) return FALSE;
    /* the singularity threshold is relative to the matrix, so that a
     * system's answer does not depend on the units it was written in */
    for (i = 0; i < n * n; i++)
        if (fabs(a[i]) > scale) scale = fabs(a[i]);
    if (scale == 0.0) return FALSE;
    for (i = 0; i < n; i++) {
        size_t pivot = i;

        /* the largest entry of the column below takes the diagonal, which
         * keeps the multipliers at most 1 */
        for (size_t k = i + 1; k < n; k++)
            if (fabs(a[k * n + i]) > fabs(a[pivot * n + i])) pivot = k;
        if (fabs(a[pivot * n + i]) < 1e-12 * scale) return FALSE;
        if (pivot != i) {
            double t = 0.0;

            for (size_t j = 0; j < n; j++) {
                t = a[i * n + j];
                a[i * n + j] = a[pivot * n + j];
                a[pivot * n + j] = t;
            }
            t = b[i];
            b[i] = b[pivot];
            b[pivot] = t;
        }
        for (size_t k = i + 1; k < n; k++) {
            double f = a[k * n + i] / a[i * n + i];

            for (size_t j = i; j < n; j++) a[k * n + j] -= f * a[i * n + j];
            b[k] -= f * b[i];
        }
    }
    /* back substitution, the last unknown first */
    for (i = n; i-- > 0;) {
        double s = b[i];

        for (size_t j = i + 1; j < n; j++) s -= a[i * n + j] * b[j];
        b[i] = s / a[i * n + i];
    }
    return TRUE;
}