la #

fn den_solve(mut x []f64, a &Matrix[f64], b []f64, preserve_a bool)

den_solve solves dense linear system using LAPACK (OpenBLaS)

Given: a ⋅ x = b find x such that x = a⁻¹ ⋅ b

fn jacobi(mut q Matrix[f64], mut v []f64, mut a Matrix[f64]) !

jacobi performs the Jacobi transformation of a symmetric matrix to find its eigenvectors and eigenvalues.

The Jacobi method consists of a sequence of orthogonal similarity transformations. Each transformation (a Jacobi rotation) is just a plane rotation designed to annihilate one of the off-diagonal matrix elements. Successive transformations undo previously set zeros, but the off-diagonal elements nevertheless get smaller and smaller. Accumulating the product of the transformations as you go gives the matrix of eigenvectors (Q), while the elements of the final diagonal matrix (A) are the eigenvalues.

The Jacobi method is absolutely foolproof for all real symmetric matrices.

A = Q ⋅ L ⋅ Qᵀ

Input: A -- matrix to compute eigenvalues (SYMMETRIC and SQUARE) Output: A -- modified Q -- matrix which columns are the eigenvectors v -- vector with the eigenvalues

Note: for matrices of order greater than about 10, say, the algorithm is slower,by a significant constant factor, than the QR method.

fn matrix_add(mut res Matrix[f64], alpha f64, a &Matrix[f64], beta f64, b &Matrix[f64])

matrix_add adds the scaled components of two matrices res := alpha⋅a + beta⋅b ⇒ result[i][j] := alpha⋅a[i][j] + beta⋅b[i][j]

fn matrix_cond_num(mut a Matrix[f64], normtype string) f64

matrix_cond_num returns the condition number of a square matrix using the inverse of this matrix; thus it is not as efficient as it could be, e.g. by using the SV decomposition. normtype -- Type of norm to use: "I" => Infinite "F" or "" (default) => Frobenius

fn matrix_det(o &Matrix[f64]) f64

det computes the determinant of matrix using the LU factorization

Note: this method may fail due to overflow...

fn matrix_inv(mut ai Matrix[f64], mut a Matrix[f64], calc_det bool) f64

matrix_inv computes the inverse of a general matrix (square or not). It also computes the pseudo-inverse if the matrix is not square. Input: a -- input matrix (M x N) Output: ai -- inverse matrix (N x M) det -- determinant of matrix (ONLY if calc_det == true and the matrix is square)

Note: the dimension of the ai matrix must be N x M for the pseudo-inverse

fn matrix_inv_small(mut ai Matrix[f64], a Matrix[f64], tol f64) f64

matrix_inv_small computes the inverse of small matrices of size 1x1, 2x2, or 3x3. It also returns the determinant. Input: a -- the matrix tol -- tolerance to assume zero determinant Output: ai -- the inverse matrix det -- determinant of a

fn matrix_matrix_mul(mut c Matrix[f64], alpha f64, a &Matrix[f64], b &Matrix[f64])

matrix_matrix_mul returns the matrix multiplication (scaled)

c := alpha⋅a⋅b ⇒ cij := alpha * aik * bkj

fn matrix_matrix_muladd(mut c Matrix[f64], alpha f64, a &Matrix[f64], b &Matrix[f64])

matrix_matrix_muladd returns the matrix multiplication (scaled)

c += alpha⋅a⋅b ⇒ cij += alpha * aik * bkj

fn matrix_matrix_tr_mul(mut c Matrix[f64], alpha f64, a &Matrix[f64], b &Matrix[f64])

matrix_matrix_tr_mul returns the matrix multiplication (scaled) with transposed(b)

c := alpha⋅a⋅bᵀ ⇒ cij := alpha * aik * bjk

fn matrix_matrix_tr_muladd(mut c Matrix[f64], alpha f64, a &Matrix[f64], b &Matrix[f64])

matrix_matrix_tr_muladd returns the matrix multiplication (scaled) with transposed(b)

c += alpha⋅a⋅bᵀ ⇒ cij += alpha * aik * bjk

fn matrix_svd(mut s []f64, mut u Matrix[f64], mut vt Matrix[f64], mut a Matrix[f64], copy_a bool)

matrix_svd performs the SVD decomposition Input: a -- matrix a copy_a -- creates a copy of a; otherwise 'a' is modified Output: s -- diagonal terms [must be pre-allocated] s.len = imin(a.m, a.n) u -- left matrix [must be pre-allocated] u is (a.m x a.m) vt -- transposed right matrix [must be pre-allocated] vt is (a.n x a.n)

fn matrix_tr_matrix_mul(mut c Matrix[f64], alpha f64, a &Matrix[f64], b &Matrix[f64])

matrix_tr_matrix_mul returns the matrix multiplication (scaled) with transposed(a)

c := alpha⋅aᵀ⋅b ⇒ cij := alpha * aki * bkj

fn matrix_tr_matrix_muladd(mut c Matrix[f64], alpha f64, a &Matrix[f64], b &Matrix[f64])

matrix_tr_matrix_muladd returns the matrix multiplication (scaled) with transposed(a)

c += alpha⋅aᵀ⋅b ⇒ cij += alpha * aki * bkj

fn matrix_tr_matrix_tr_mul(mut c Matrix[f64], alpha f64, a &Matrix[f64], b &Matrix[f64])

matrix_tr_matrix_tr_mul returns the matrix multiplication (scaled) with transposed(a) and transposed(b)

c := alpha⋅aᵀ⋅bᵀ ⇒ cij := alpha * aki * bjk

fn matrix_tr_matrix_tr_mul_add(mut c Matrix[f64], alpha f64, a &Matrix[f64], b &Matrix[f64])

matrix_tr_matrix_tr_mul_add returns the matrix multiplication (scaled) with transposed(a) and transposed(b)

c += alpha⋅aᵀ⋅bᵀ ⇒ cij += alpha * aki * bjk

fn matrix_tr_vector_mul[T](alpha T, a &Matrix[T], u []T) []T

matrix_tr_vector_mul returns the transpose(matrix)-vector multiplication

v = alpha⋅aᵀ⋅u ⇒ vi = alpha * aji * uj = alpha * uj * aji

fn matrix_vector_mul[T](alpha T, a &Matrix[T], u []T) []T

matrix_vector_mul returns the matrix-vector multiplication

v = alpha⋅a⋅u ⇒ vi = alpha * aij * uj

fn matrix_vector_mul_add(alpha f64, a &Matrix[f64], u []f64) []f64

matrix_vector_mul_add returns the matrix-vector multiplication with addition

v += alpha⋅a⋅u ⇒ vi += alpha * aij * uj

fn safe_print[T](format string, message T) string

fn vector_accum[T](o []T) T

vector_accum sum/accumulates all components in a []T sum := Σ_i v[i]

fn vector_add[T](alpha T, u []T, beta T, v []T) []T

vector_add adds the scaled components of two vectors res := alpha⋅u + beta⋅v ⇒ result[i] := alpha⋅u[i] + beta⋅v[i]

fn vector_apply[T](mut o []T, a T, another []T)

vector_apply sets this []T with the scaled components of another []T this := a * another ⇒ this[i] := a * another[i]

Note: "another" may be "this"

fn vector_apply_func[T](mut o []T, f fn (i int, x T) T)

vector_apply_func runs a function over all components of a []T vi = f(i,vi)

fn vector_dot[T](u []T, v []T) T

vector_dot returns the dot product between two vectors: s := u・v

fn vector_largest[T](o []T, den T) T

vector_largest returns the largest component |u[i]| of this []T, normalised by den largest := |u[i]| / den

fn vector_max_diff[T](u []T, v []T) T

vector_max_diff returns the maximum absolute difference between two vectors maxdiff = max(|u - v|)

fn vector_norm(o []f64) f64

vector_norm returns the Euclidean norm of a []T: nrm := ‖v‖

fn vector_norm_diff[T](o []T, v []T) T

vector_norm_diff returns the Euclidean norm of the difference: nrm := ||u - v||

fn vector_rms[T](o []T) T

vector_rms returns the root-mean-square of this []T

fn vector_rms_error[T](u []T, v []T, a T, m T, s []T) T

fn vector_scale_abs[T](a T, m T, x []T) []T

vector_scale_abs creates a "scale" vector using the absolute value of another vector scale := a + m ⋅ |x| ⇒ scale[i] := a + m ⋅ |x[i]|

fn vector_unit(mut o []f64) []f64

vector_unit returns the unit []f64 parallel to this []f64 b := a / norm(a)

fn vector_vector_tr_mul[T](alpha T, u []T, v []T) &Matrix[T]

vector_vector_tr_mul returns the matrix = vector-transpose(vector) multiplication (e.g. dyadic product)

a = alpha⋅u⋅vᵀ ⇒ aij = alpha * ui * vj

fn Matrix.deep2[T](a [][]T) &Matrix[T]

Matrix.deep2 allocates a new Matrix from given (Deep2) nested slice.

Note: make sure to have at least 1x1 item

fn Matrix.new[T](m int, n int) &Matrix[T]

Matrix.new allocates a new (empty) Matrix with given (m,n) (row/col sizes)

fn Matrix.raw[T](m int, n int, rawdata []T) &Matrix[T]

Matrix.raw creates a new Matrix using given raw data Input: rawdata -- data organized as column-major; e.g. Fortran format

Note: (1) rawdata is not copied! (2) the external slice rawdata should not be changed or deleted

fn SparseConfig.new() SparseConfig

SparseConfig.new returns a new SparseConfig Input: comm -- may be nil

fn SparseConfig.with_comm(comm &mpi.Communicator) SparseConfig

SparseConfig.with_comm returns a new SparseConfig

fn Triplet.new[T](m int, n int, max int) &Triplet[T]

Triplet.new returns a new Triplet. This is a wrapper to new(Triplet) followed by init()

fn (mut o Matrix[T]) set_from_deep2(a [][]T)

set_from_deep2 sets matrix with data from a nested slice (Deep2) structure

fn (mut o Matrix[T]) set_diag(val T)

set_diag sets diagonal matrix with diagonal components equal to val

fn (mut o Matrix[T]) set(i int, j int, val T)

set sets value

fn (o &Matrix[T]) get(i int, j int) T

get gets value

fn (o &Matrix[T]) get_deep2() [][]T

get_deep2 returns nested slice representation

fn (o &Matrix[T]) clone() &Matrix[T]

clone returns a copy of this matrix

fn (o &Matrix[T]) transpose() &Matrix[T]

transpose returns the transpose matrix

fn (o &Matrix[T]) copy_into(mut result Matrix[T], alpha T)

copy_into copies the scaled components of this matrix into another one (result) result := alpha * this ⇒ result[ij] := alpha * this[ij]

fn (mut o Matrix[T]) add(i int, j int, val T)

add adds value to (i,j) location

fn (mut o Matrix[T]) fill(val T)

fill fills this matrix with a single number val aij = val

fn (mut o Matrix[T]) clear_rc(rows []int, cols []int, diag T)

fn (mut o Matrix[T]) clear_bry(diag T)

fn (o &Matrix[T]) max_diff(another &Matrix[T]) T

max_diff returns the maximum difference between the components of this and another matrix

fn (o &Matrix[T]) largest(den T) T

largest returns the largest component |a[ij]| of this matrix, normalised by den largest := |a[ij]| / den

fn (o &Matrix[T]) col(j int) []T

col access column j of this matrix. No copies are made since the internal data are in row-major format already.

Note: this method can be used to modify the columns; e.g. with o.col(0)[0] = 123

fn (o &Matrix[T]) get_row(i int) []T

get_row returns row i of this matrix

fn (o &Matrix[T]) get_col(j int) []T

get_col returns column j of this matrix

fn (o &Matrix[T]) extract_cols(start int, endp1 int) !&Matrix[T]

extract_cols returns columns from j=start to j=endp1-1 start -- first column endp1 -- "end-plus-one", the number of the last requested column + 1

fn (o &Matrix[T]) extract_rows(start int, endp1 int) !&Matrix[T]

extract_rows returns rows from i=start to i=endp1-1 start -- first column endp1 -- "end-plus-one", the number of the last requested column + 1

fn (mut o Matrix[T]) set_row(i int, value T)

set_row sets the values of a row i with a single value

fn (mut o Matrix[T]) set_col(j int, value T)

set_col sets the values of a column j with a single value

fn (o &Matrix[T]) split_by_col(j int) !(&Matrix[T], &Matrix[T])

split_by_col splits this matrix into two matrices at column j j -- column index

fn (o &Matrix[T]) split_by_row(i int) !(&Matrix[T], &Matrix[T])

split_by_row splits this matrix into two matrices at row i i -- row index

fn (o &Matrix[T]) norm_frob() T

norm_frob returns the Frobenius norm of this matrix nrm := ‖a‖_F = sqrt(Σ_i Σ_j a[ij]⋅a[ij]) = ‖a‖_2

fn (o &Matrix[T]) norm_inf() T

norm_inf returns the infinite norm of this matrix nrm := ‖a‖_∞ = max_i ( Σ_j a[ij] )

fn (mut o Matrix[T]) apply(alpha T, another &Matrix[T])

apply sets this matrix with the scaled components of another matrix this := alpha * another ⇒ this[i] := alpha * another[i]

Note: "another" may be "this"

fn (o &Matrix[T]) equals(another &Matrix[T]) bool

equals returns true if this matrix is equal to another matrix this == another ⇒ this[i] == another[i]

Note: "another" may be "this"

fn (o &Matrix[T]) str() string

fn (o &Matrix[T]) print(nfmt_ string) string

print prints matrix (without commas or brackets)

fn (o &Matrix[T]) print_v(nfmt_ string) string

print_v prints matrix in V format

fn (o &Matrix[T]) print_py(nfmt_ string) string

print_py prints matrix in Python format

fn (mut o Triplet[T]) init(m int, n int, max int)

init allocates all memory required to hold a sparse matrix in triplet form

fn (mut o Triplet[T]) put(i int, j int, x T) !

put inserts an element to a pre-allocated (with init) triplet matrix

fn (mut o Triplet[T]) put_matrix_and_matrix_t(a &Triplet[T]) !

fn (mut o Triplet[T]) put_cc_matrix_and_matrix_t(a &CCMatrix[T]) !

fn (mut o Triplet[T]) start()

start (re)starts index for inserting items using the put command

fn (mut o Triplet[T]) len() int

Len returns the number of items just inserted in the triplet

fn (o Triplet[T]) max() int

max returns the maximum number of items that can be inserted in the triplet

fn (o Triplet[T]) size() (int, int)

size returns the row/column size of the matrix represented by the Triplet

fn (o Triplet[T]) to_dense() &Matrix[T]

to_dense returns the dense matrix corresponding to this Triplet

struct CCMatrix[T] {
mut:
	// matrix dimension (rows, columns)
	m int
	n int
	// number of non-zeros
	nnz int
	// pointers and row indices (len(p)=n+1, len(i)=nnz)
	p []int
	i []int
	// values (len(x)=nnz)
	x []T
}

CCMatrix represents a sparse matrix using the so-called "column-compressed format".

struct Matrix[T] {
pub mut:
	m    int
	n    int
	data []T
}

struct SparseConfig {
mut:
	communicator   ?&mpi.Communicator // MPI communicator for parallel solvers [may be none]
	mumps_ordering int  // ICNTL(7) default = "" == "auto"
	mumps_scaling  int  // Scaling type (check MUMPS solver) [may be empty]
	symmetric      bool // indicates symmetric system. NOTE: when using MUMPS, only the upper or lower part of the matrix must be provided
	sym_pos_def    bool // indicates symmetric-positive-defined system. NOTE: when using MUMPS, only the upper or lower part of the matrix must be provided
pub mut:
	verbose                             bool // run on verbose mode
	mumps_increase_of_working_space_pct int  // MUMPS control parameters (check MUMPS solver manual) ICNTL(14) default = 100%
	mumps_max_memory_per_processor      int  // MUMPS control parameters (check MUMPS solver manual) ICNTL(23) default = 2000Mb
}

The SparseConfig structure holds configuration arguments for sparse solvers

fn (mut o SparseConfig) set_mumps_symmetry(only_upper_or_lower_given bool, positive_defined bool)

set_mumps_symmetry sets symmetry options for MUMPS solver

fn (mut o SparseConfig) set_umfpack_symmetry()

set_umfpack_symmetry sets symmetry options for UMFPACK solver

fn (mut o SparseConfig) set_mumps_ordering(ordering string)

set_mumps_ordering sets ordering for MUMPS solver ordering -- "" or "amf" [default] "amf", "scotch", "pord", "metis", "qamd", "auto" ICNTL(7) 0: "amd" Approximate Minimum Degree (AMD) 2: "amf" Approximate Minimum Fill (AMF) 3: "scotch" SCOTCH5 package is used if previously installed by the user otherwise treated as 7. 4: "pord" PORD6 is used if previously installed by the user otherwise treated as 7. 5: "metis" Metis7 package is used if previously installed by the user otherwise treated as 7. 6: "qamd" Approximate Minimum Degree with automatic quasi-dense row detection (QAMD) is used. 7: "auto" automatic choice by the software during analysis phase. This choice will depend on the ordering packages made available, on the matrix (type and size), and on the number of processors.

fn (mut o SparseConfig) set_mumps_scaling(scaling string)

set_mumps_scaling sets scaling for MUMPS solver scaling -- "" or "rcit" [default] "no", "diag", "col", "rcinf", "rcit", "rrcit", "auto" ICNTL(8) 0: "no" No scaling applied/computed. 1: "diag" Diagonal scaling computed during the numerical factorization phase, 3: "col" Column scaling computed during the numerical factorization phase, 4: "rcinf" Row and column scaling based on infinite row/column norms, computed during the numerical factorization phase, 7: "rcit" Simultaneous row and column iterative scaling based on [41] and [15] computed during the numerical factorization phase. 8: "rrcit" Similar to 7 but more rigorous and expensive to compute; computed during the numerical factorization phase. 77: "auto" Automatic choice of the value of ICNTL(8) done during analy

struct Triplet[T] {
mut:
	// matrix dimension (rows, columns)
	m int
	n int
	// current position and max number of entries allowed (non-zeros, including repetitions)
	pos int
	max int
	// indices for each x values (size=max)
	i []int
	j []int
	// values for each i, j (size=max)
	x []T
}

Triplet is a simple representation of a sparse matrix, where the indices and values of this matrix are stored directly.