AlgebraLib

Linear algebra domain providing typed dense, structured, sparse, and matrix-free operations, direct/iterative solvers, partial eigensystems, and vector operations for Free Pascal.

Depends on: MathBase

Learning routes

Beginner route

New code should copy and run the typed double-real solve, which prints solution = 2.0000, 3.0000. Solve allocates a result and private factor storage without mutating its inputs. The lower quick start documents the retained compatibility API, not the primary beginner path.

Common tasks and algorithm choice

Advanced route

Run example 16 for QR/SVD/eigen selection and example 22 for sparse storage, preconditioning, and diagnostic statuses. Both keep double-real typed containers; reusable factors, views, destinations, and workspaces are deeper forms of the same documented APIs and require no private conversion.

Units

The typed dense solver-selection and complete factor contracts are in Typed dense decompositions and solvers. The large-problem selection, storage, residual, and reuse contracts are in Sparse, structured, and matrix-free linear algebra. The IMatrix records below remain the source-compatible legacy surface.

---

Core Types

Exception

EMatrixError = class(Exception);

Raised on dimension mismatches, singular matrices, invalid indices, and other matrix errors.

Storage

TMatrixArray = array of array of Double;

Underlying 2-D array type used by the concrete implementation.

Decomposition Records

All records expose a ToString: string method for easy debugging.

Iterative Solver Enum

TIterativeMethod = (imConjugateGradient, imGaussSeidel, imJacobi);

Vector APIs (AlgebraLib.Vectors)

IVector = IMatrix;     // A vector is a 1-row or 1-column IMatrix

TVector = TMatrixKit;  // Concrete type for construction

The matrix-vector API remains the compatibility-oriented API. For contiguous real or complex arrays, AlgebraLib.Vectors also re-exports:

TRealVector = TDoubleArray;

TComplexVector = TComplexArray;

TVectorKit = class;

TVectorKit provides allocation-returning Add, Subtract, ElementWiseMultiply, ElementWiseDivide, Scale, and Axpy operations, plus compensated Sum/Dot, Mean, Min, Max, stable Norm2, and Normalize for real vectors. Its complex overloads add non-conjugating Dot and Hermitian DotConjugate.

Every allocating transform also has an ...Into variant, for example:

var Destination: TRealVector;

begin

  SetLength(Destination, Length(A));

  TVectorKit.AxpyInto(0.5, A, B, Destination); // Destination := 0.5*A + B

end;

Destination is a var parameter so a correctly sized existing dynamic array is reused rather than replaced. Elementwise ...Into transforms may use one of their inputs as the destination. Inputs must be finite and paired vectors must have equal length; empty-vector sums, dots, and norms are zero, while Mean, Min, Max, and zero-vector normalization raise EVectorError.

Decomposition Entry Point (AlgebraLib.Determinants)

TIterSolverMethod = TIterativeMethod;

The decomposition records retain their original names from AlgebraLib.Matrices: TLUDecomposition, TQRDecomposition, TEigenDecomposition, TSVD, TCholeskyDecomposition, and TEigenpair. AlgebraLib.Determinants provides the focused import path and the one explicit enum alias shown above; it does not declare shortened record aliases such as TLUDecomp or TSVDecomp.

---

IMatrix Interface

Arithmetic, transformations, decompositions, and matrix functions return new matrices. SetValue and SetSubMatrix are explicit in-place mutators; the sparse implementation also exposes mutating storage methods.

Dimensions

Basic Arithmetic

Transformations

Matrix Functions

Matrix Properties (scalar results)

Type Checks

function IsSquare: Boolean;

function IsSymmetric: Boolean;

function IsDiagonal: Boolean;

function IsTriangular(Upper: Boolean = True): Boolean;

function IsPositiveDefinite: Boolean;

function IsPositiveSemidefinite: Boolean;

function IsOrthogonal: Boolean;

function IsVector: Boolean;

function IsRowVector: Boolean;

function IsColumnVector: Boolean;

Decompositions

Linear System Solvers

The compatibility IMatrix interface has no direct SolveLinear method. New code can use the typed Solve(A, B) path, which factors the coefficient matrix instead of forming its inverse. Existing code can opt in through TDenseDoubleMatrix.FromIMatrix. SolveIterative raises EMatrixError if the selected method exhausts MaxIterations; it never silently returns an unconverged last iterate.

Vector Operations (single-row or single-column matrices)

Statistical Methods

Submatrices and Dominant Eigenpair

function GetSubMatrix(StartRow, StartCol, NumRows, NumCols: Integer): IMatrix;

procedure SetSubMatrix(StartRow, StartCol: Integer; const SubMatrix: IMatrix);

function PowerMethod(MaxIterations: Integer = 100;

  Tolerance: Double = 1e-10): TEigenpair;

All indices are zero-based. SetSubMatrix mutates the receiving matrix. PowerMethod returns the dominant real eigenpair and requires a square matrix. It validates positive controls and raises EMatrixError if normalization or convergence fails.

---

TMatrixKit — Factory Methods

TMatrixKit is the concrete class that implements IMatrix. Its constructor creates a zero-filled dense matrix; class functions provide specialised forms.

constructor Create(Rows, Cols: Integer);

class function Identity(N: Integer): IMatrix;

class function Zeros(Rows, Cols: Integer): IMatrix;

class function Ones(Rows, Cols: Integer): IMatrix;

class function CreateFromArray(const Data: TMatrixArray): IMatrix;

class function CreateSparse(Rows, Cols: Integer): IMatrix;

class function CreateDiagonal(const Values: array of Double): IMatrix;

class function CreateBandMatrix(Size, LowerBand, UpperBand: Integer): IMatrix;

class function CreateSymmetric(const Data: TMatrixArray): IMatrix;

class function CreateHilbert(Size: Integer): IMatrix;

class function CreateToeplitz(const FirstRow, FirstCol: TDoubleArray): IMatrix;

class function CreateVandermonde(const Vector: TDoubleArray): IMatrix;

class function CreateRandom(Rows, Cols: Integer; Min, Max: Double): IMatrix; overload;

class function CreateRandom(Rows, Cols: Integer; Min, Max: Double;

  Seed: LongWord): IMatrix; overload;

CreateBandMatrix fills the requested band with ones. CreateSymmetric reads the lower triangle of Data and mirrors it. CreateToeplitz returns a Length(FirstCol) by Length(FirstRow) matrix and requires matching first elements.

Sparse Matrices

TMatrixKitSparse stores non-zero entries in row-major order. GetValue and SetValue are available through IMatrix; zero (within 1e-15) removes an entry. The concrete class additionally exposes:

procedure AddElement(Row, Col: Integer; Value: Double); // set/replace, not accumulate

procedure CompactStorage;

Sparse lookup and insertion are linear in the number of stored entries, and Add currently returns a dense matrix.

---

Quick Start

For new dense code, start with the typed dense five-minute guide. The following example documents the retained compatibility API.

uses AlgebraLib.Matrices;



var

  A, B, C: IMatrix;

  LU: TLUDecomposition;

begin

  A := TMatrixKit.CreateFromArray([[1,2],[3,4]]);

  B := TMatrixKit.Identity(2);



  C  := A.Multiply(B);                // C = A × I = A

  LU := A.LU;                         // PA = LU



  Writeln('Det = ', A.Determinant:0:4);  // -2.0

  Writeln('Rank = ', A.Rank);            // 2

end.

Expected output:

Det = -2.0000

Rank = 2

Design Notes

• Mostly value-oriented — calculations return new IMatrix instances; the explicitly named setters mutate their receiver.
• Interface-based — depend on IMatrix, not on TMatrixKit, for flexibility.
• Cache-aware blocking — matrix multiplication uses cache-aware blocking and bounded parallel workers when the operation count justifies thread startup. Unix callers without an installed thread manager automatically use the serial path.
• Numerically stable — LU uses partial pivoting and preserves prior L multipliers across row swaps; singularity and rank tolerances are relative to matrix scale.
• Explicit convergence — iterative solves and matrix exponential evaluation raise EMatrixError when they cannot converge or produce a finite representable result.
• Real eigensystem contract — the API raises EMatrixError for complex spectra, defective matrices, and unsupported nonsymmetric matrices larger than 2×2 instead of returning misleading real approximations.
• Reproducible random matrices — the seeded CreateRandom overload uses local state and does not change the process-wide RandSeed; the compatibility overload uses caller-managed global state and never calls Randomize.