std.linalg — dichte lineare Algebra

Matrizen in f64: LU, QR, Cholesky, Gleichungssysteme lösen, Determinante, Inverse, Rang, Ausgleichsrechnung und symmetrische Eigenwerte nach Jacobi.

→ Mathematik-Units · std.decomp · std.sparse · std.vector

import std.linalg;

Schicht 0, hängt an std.math und std.alloc. 35 Funktionen. Alle Beispiele mit lyxc 1.0.21A übersetzt und ausgeführt und gegen eine unabhängige Referenzrechnung geprüft.

Eine m×n-Matrix ist ein Speicherblock aus m·n f64-Werten, zeilenweise abgelegt. Der Zugriff läuft über LinGet/LinSet, die die Spaltenzahl als Argument bekommen:

var n: int64 := 3;
var a: int64 := alloc(n * n * 8);
LinSet(a, n, 0, 1, -2.0);        // Zeile 0, Spalte 1
var v: f64 := LinGet(a, n, 0, 1);

unit main;
import std.io;
import std.string;
import std.alloc;
import std.linalg;

fn ZF(t: pchar, v: f64): void { Print(t); PrintF64(v); }
fn ZI(t: pchar, v: int64): void { Print(t); PrintLn(IntToStr(v)); }

fn main(): int64 {
    var n: int64 := 3;
    var a: int64 := alloc(n * n * 8);
    LinSet(a, n, 0,0,  4.0); LinSet(a, n, 0,1, -2.0); LinSet(a, n, 0,2,  1.0);
    LinSet(a, n, 1,0, -2.0); LinSet(a, n, 1,1,  4.0); LinSet(a, n, 1,2, -2.0);
    LinSet(a, n, 2,0,  1.0); LinSet(a, n, 2,1, -2.0); LinSet(a, n, 2,2,  4.0);

    ZF("Det         : ", LinDet(a, n));
    ZF("Trace       : ", LinTrace(a, n));
    ZI("Rank        : ", LinRank(a, n, n));
    ZI("IsSymmetric : ", LinIsSymmetric(a, n, LinEpsilon()) as int64);
    ZI("IsPosDef    : ", LinIsPosDef(a, n) as int64);
    ZF("NormFro     : ", LinNormFro(a, n, n));
    ZF("Cond2       : ", LinCond2(a, n, n));

    var b: int64 := alloc(n * 8);
    pokef64(b, 1.0); pokef64(b + 8, 2.0); pokef64(b + 16, 3.0);
    var x: int64 := alloc(n * 8);
    var ok: bool := LinSolve(a, n, b, x);
    Print("Solve ok="); Print(IntToStr(ok as int64)); PrintLn("");
    ZF("x0          : ", peekf64(x));
    ZF("x1          : ", peekf64(x + 8));
    ZF("x2          : ", peekf64(x + 16));

    var ev: int64 := alloc(n * 8);
    var evec: int64 := alloc(n * n * 8);
    var ok2: bool := LinEigenSym(a, n, ev, evec);
    LinEigenSortDesc(ev, evec, n);
    Print("Eigen ok="); Print(IntToStr(ok2 as int64)); PrintLn("");
    ZF("lambda0     : ", peekf64(ev));
    ZF("lambda1     : ", peekf64(ev + 8));
    ZF("lambda2     : ", peekf64(ev + 16));
    return 0;
}

Det         : 36.000000
Trace       : 12.000000
Rank        : 3
IsSymmetric : 1
IsPosDef    : 1
NormFro     : 8.124038
Cond2       : 4.529210
Solve ok=1
x0          : 0.666666
x1          : 1.500000
x2          : 1.333333
Eigen ok=1
lambda0     : 7.372281
lambda1     : 2.999999
lambda2     : 1.627718

Sämtliche Werte stimmen mit der Referenzrechnung überein. lambda1 erscheint als 2.999999 statt 3: das Jacobi-Verfahren ist iterativ und endet an einer Schwelle, nicht exakt. Für Eigenwerte gilt daher immer ein Vergleich mit Toleranz (→ std.feq), nie ==.

Die Zerlegung einmal, das Lösen mehrfach. LinSolve zerlegt bei jedem Aufruf neu; wer dasselbe a mit zehn verschiedenen b löst, ruft einmal LinLuDecompose und danach zehnmal LinLuSolve. LinLuDecompose schreibt in die Matrix hinein und braucht zusätzlich ein Permutationsfeld (n Einträge) sowie eine Adresse für das Vorzeichen — wer a danach noch braucht, kopiert vorher mit LinCopy.

LinIsPosDef beantwortet vorab, ob der Cholesky-Weg überhaupt gangbar ist.

LinEquals nimmt ein eps, weil ein Vergleich zweier gerechneter Matrizen mit == praktisch immer scheitert. LinEpsilon() liefert den Vorgabewert 1e−12.

LinEigenSym setzt eine symmetrische Matrix voraus — das Jacobi-Verfahren beruht darauf. Für allgemeine Matrizen und für die vollständige SVD ist std.decomp zuständig.

LinEigenSortDesc sortiert Werte und Vektoren gemeinsam; wer nur die Werte sortiert, ordnet die Vektoren falsch zu.

LinAdd, LinSub, LinScale, LinMul, LinMulVec, LinTranspose — dazu LinFill, LinZero, LinIdentity und LinCopy zum Anlegen.

LinMul(a, b, m, k, n, out) multipliziert eine m×k- mit einer k×n-Matrix. Die drei Maße stehen einzeln in der Signatur, weil sich das aus den Zeigern nicht ableiten lässt. out darf nicht auf a oder b zeigen.

Dünnbesetzte Matrizen und iterative Löser (CG, GMRES) gibt es derzeit nicht. std.matrix aus dem Bestand existiert, hat aber keine eigene Wiki-Seite; für Rechnungen in f64 ist std.linalg der Weg.

Letzte Aktualisierung: 2026-08-16 · alle Beispiele mit lyxc 1.0.21A übersetzt und ausgeführt; Determinante, Norm, Kondition, Lösung und Eigenwerte gegen eine unabhängige Referenzrechnung geprüft.