Calcpad - Online help

This site uses cookies to offer you a better browsing experience. Learn More.

Topic: Expressions

Online Help

Matrices

Internal implementation and types of matrices

Calcpad includes different types of matrices: general (rectangular) and special (column, symmetric, diagonal. upper/lower triangular). Internally, each type is implemented in a different way that benefits from the specific structure for better efficiency. Externally, all of them behave in a similar way, except for a few cases.

Each matrix type is implemented as an array of vectors, as displayed on the figure below. Vectors normally represent matrix rows, except for diagonal and column matrices.

Calcpad uses large vectors to contain the values. So, it does not store the extra zero elements for partially filled (banded) matrices. The indexing operator for each type is internally redefined in order to return directly zero when we try to read a value outside the matrix structure or bandwidth.

a) diagonal matrix -	M[i, j] = d[i], if i = j and 0, if i ≠ j;
b) column matrix -	M[i, j] = c[i], if j = 1, otherwise – error;
c) upper triangular matrix -	M[i, j] = r_i[j – i + 1], if j ≥ i, otherwise – 0;
d) lower triangular matrix -	M[i, j] = r_i[j], if j ≤ i, otherwise – 0;
e) row matrix -	M[i, j] = r[j], if i = 1, otherwise – error;
f) symmetric matrix -	M[i, j] = r_i[j – i + 1], if i ≥ j, otherwise = r_i [i – j + 1];
g) rectangular matrix -	M[i, j] = r_i[j];

If we try to write a non-zero value outside the matrix structure, we will get an "Index out of range" error. For example, you cannot assign a non-zero value to an element outside the main diagonal of a diagonal type of matrix.

Definition

Similar to vectors, you can define matrices by using the "square brackets" syntax, but the rows must be separated by vertical bars " | ", as follows:

A_m⨯n = [a_1,1; a_1,2; … ; a_1,n | a_2,1; a_2,2; … ; a_2,n | … | a_m,1; a_m,2; … ; a_m,n]

In this way, you can create only general (rectangular) types of matrices. For special types of matrices, you have to use the respective creational functions as described further in this manual. If you have rows with different lengths, the number of columns n is assumed to be the maximum number of elements in a row. The missing elements in other rows are assumed to be zero. For example:

A = [1|2; 3|4; 5; 6|7; 8]' = 100 230 456 780

You can use expressions for matrix elements that contain variables, operators, functions, vectors, other matrices, etc. For example, the following code will create a matrix of three rows, by applying a different expression for each row on a single vector:

a = [1; 2; 4]
A = [a|2*a + 1|3*a + 2]' = 124 359 5814

Just like vectors, matrices can also be defined as functions to create them dynamically on demand. The following function generates a 4×4 Vandermonde matrix from a vector containing the elements for the first column:

A(x) = transp([x^0|x|x^2|x^3|x^4])
x = [1; 2; 3; 4]
A = A(x)' = 11111 124816 1392781 141664256

Indexing

You can use indexing to access individual matrix elements for reading and writing their values. Similar to vectors, this is performed by the dot operator, but you have to specify two indexes, as follows:

A.(i; j), where:
i - the index of the row where the element is located,
j - the index of the column.

Indexes must be enclosed in brackets and divided by a semicolon. Row and column numbering start from one. For the Vandermonde matrix from the above example:

A.(3; 2)' = 3.

You can have expressions inside the brackets to calculate the indexes in place:

i = 2', 'j = 3
A.(2*i - 1; j + 1)' A_{3, 4} = 27.

In this way, you can iterate through matrix elements in a loop by including the loop counters in the respective indexes. You can use both inline and block loops for that purpose. The code below creates a Vandermonde matrix from vector ⃗x with the specified number of columns (6):

x = [1; 2; 3; 4]
A = matrix(len(x); 7)
#hide
#for i = 1 : n_rows(A)
#for j = 1 : n_cols(A)
A.(i; j) = x.i^(j - 1)
#loop
#loop
#show
A' = 1111111 1248163264 1392781243729 14166425610244096

The inline equivalent of the above loop is the following:

$Repeat{$Repeat{A.(i; j) = x.i^(j - 1) @ j = 1 : n_cols(A)} @ i = 1 : n_rows(A)}

Creational functions

The "square brackets" syntax is very powerful and flexible for creating small matrices with predefined sizes. However, it also has a lot of limitations. For example, it cannot create special types of matrices and cannot specify the matrix dimensions. That is why, Calcpad also includes various functions for creating matrices, as follows:

matrix(m; n)

Parameters: m - (positive integer) number of rows;
n - (positive integer) number of columns.
Return value: an empty rectangular matrix with dimensions m⨯n.
Notes: m and n must be between 1 and 1 000 000. This also applies to all other kinds of matrices below.
Example: matrix(3; 4)' = 0000 0000 0000

identity(n)

Parameters: n - number of rows/columns.
Return value: an identity matrix with dimensions n⨯n.
Notes: Represents is a diagonal matrix, filled with one along the main diagonal.
This function is equivalent to diagonal ( n ; 1).
Example: identity(3)' = 100 010 001

diagonal(n; d)

Parameters: n - number of rows/columns;
d - the value along the main diagonal.
Return value: an n⨯n diagonal matrix, filled with value d along the main diagonal.
Notes: It is internally different and more efficient than an n⨯n symmetric matrix.
Example: diagonal(3; 2)' = 200 020 002

column(m; c)

Parameters: m - number of rows;
c - a value to fill the matrix with.
Return value: an m⨯1 column matrix, filled with value c.
Notes: It is internally different and more efficient than an mx1 rectangular matrix.
Example: column(3; 2)' = 2 2 2

utriang(n)

Parameters: n - number of rows/columns.
Return value: an empty upper triangular matrix with dimensions n⨯n.
Notes: It is internally different and more efficient than a general n⨯n matrix.
Example: U = utriang(3)
mfill(U; 1)' = 111 011 001

ltriang(n)

Parameters: n - number of rows/columns.
Return value: an empty lower triangular matrix with dimensions n⨯n.
Notes: It is internally different and more efficient than general n⨯n matrix.
Example: L = ltriang(3)
mfill(L; 1)' = 100 110 111

symmetric(n)

Parameters: n - number of rows/columns.
Return value: an empty n⨯n symmetric matrix.
Notes: It is internally different and more efficient than general n⨯n matrix. Only the filled side of the upper-right half of the matrix is stored, forming a skyline structure. If you change a value on either of both sides also changes the respective value on the other side, keeping the matrix always symmetric.
Example: A = symmetric(4)
A.(1; 1) = 5', 'A.(1; 2) = 4
A.(2; 2) = 3', 'A.(2; 3) = 2
A.(4; 2) = 1', 'A.(4; 4) = 1
A' = 5400 4321 0200 0101

vec2diag(⃗v)

Parameters: ⃗v - a vector containing the diagonal elements.
Return value: a diagonal matrix from the elements of vector ⃗v.
Notes: The size of the matrix is equal to the number of elements in ⃗v.
Example: vec2diag([1; 2; 3]' = 100 020 003

vec2row(⃗v)

Parameters: ⃗v - a vector containing the elements of the row matrix.
Return value: a row matrix from the elements of vector ⃗v.
Notes: The number of columns of the matrix is equal to the size of ⃗v.
Example: vec2row([1; 2; 3]' = 1 2 3

vec2col(⃗v)

Parameters: ⃗v - a vector containing the elements of the column matrix.
Return value: a column matrix from the elements of vector ⃗v.
Notes: The number of rows of the matrix is equal to the size of ⃗v.
Example: vec2col([1; 2; 3]' = 1 2 3

join_cols(⃗c₁; ⃗c₂; ⃗c₃ …)

Parameters: ⃗c₁; ⃗c₂; ⃗c₃ … - vectors.
Return value: a new rectangular matrix, obtained by joining the column vectors.
Notes: You can specify arbitrary number of input vectors with different lengths. The number of rows is equal to the maximum vector length and the other columns are filled down with zeros to the end. The vectors are joined sequentially from left to right.
Example: join_cols([1]; [4; 5; 6]; [7; 8]; [10; 11; 12]' = 14710 05811 06012

join_rows(⃗r₁; ⃗r₂; ⃗r₃ …)

Parameters: ⃗r₁; ⃗r₂; ⃗r₃ … - vectors.
Return value: a new rectangular matrix, which rows are the specified vector arguments.
Notes: You can have an arbitrary number of input vectors with different lengths. The number of columns is equal to the maximum vector length. The other rows are filled up with zeros to the end. The vectors are collected from left to right and arranged into rows from top to bottom.
Example: join_rows([1; 2; 3; 4]; [6; 7; 8; 9; 10]' = 12340 678910

augment(A; B; C…)

Parameters: A, B, C … - matrices.
Return value: a rectangular matrix obtained by joining the input matrices side to side.
Notes: You can specify arbitrary number of input matrices with different number of rows. The largest number is relevant and smaller matrices are filled down with zeros to the end. Matrices are joined sequentially from left to right. If the arguments contain vectors, they are treated as columns.
Example: A = [1; 2|3; 4]
B = [1; 2; 3|4; 5; 6|7; 8; 9]
c = [10; 11; 12; 13]
D = augment(A; B; c)' = 1212310 3445611 0078912 0000013

stack(A; B; C…)

Parameters: A, B, C … - matrices.
Return value: a rectang. matrix obtained by stacking the input matrices one below the other.
Notes: You can specify an arbitrary number of input matrices with different row lengths. The largest number is relevant and smaller matrices are filled up with zeros to the end. Matrices are joined sequentially from top to bottom. If the arguments contain vectors, they are treated as columns.
Example: A = [1; 2|3; 4]
B = [1; 2; 3|4; 5; 6|7; 8; 9]
c = [10; 11]
D = stack(A; B; c)' = 120 340 123 456 789 1000 1100

Structural functions

Structural functions are related only to the matrix structure. Unlike data and math functions, the result does not depend much on the values of the elements.

n_rows(M)

Parameters: M - matrix.
Return value: the number of rows in matrix M.
Example: n_rows([1; 2|3; 4|5; 6])' = 3

n_cols(M)

Parameters: M - matrix.
Return value: the number of columns in matrix M.
Example: n_cols([1|2; 3|4; 5; 6])' = 3

mresize(M; m; n)

Parameters: M - the matrix to be resized;
m - the new number of rows;
n - the new number of columns.
Return value: the matrix M resized to m rows and n columns.
Notes: If the new dimensions are larger, the added elements are initialized with zeroes. If smaller, the matrix is truncated. If the new dimensions are compatible to the matrix type, then the original matrix is resized in place and a reference to it is returned. Otherwise, a new rectangular type of matrix is returned with the specified dimensions and the original matrix remains unchanged.
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789 - general type.
B = mresize(A; 2; 4)' = 1230 4560 - also a general type.
A' = 1230 4560 - the original matrix was changed.
A = diagonal(3; 5)' = 500 050 005 - diagonal type.
B = mresize(A; 2; 4)' = 5000 0500 - general type.
A' = 500 050 005 - the original matrix remains unchanged.

mfill(M; x)

Parameters: M - the matrix to be filled;
x - the value to fill the matrix with.
Return value: the matrix M filled with value x.
Notes: The matrix M is modified in place and a reference to it is returned. This function is structurally compliant. Special matrices are filled only in the designated area, preserving the type.
Example: A = matrix(2; 3)' = 000 000
B = mfill(A; 1)' = 111 111
A' = 111 111
L = mfill(ltriang(4); 2)' = 2000 2200 2220 2222

fill_row(M; i; x)

Parameters: M - the matrix which row is to be filled;
i - the index of the row to be filled;
x - the value to fill the row with.
Return value: the matrix M with the i-th row filled with value x.
Notes: The matrix M is modified in place and a reference to it is returned. This function is structurally compliant. Special matrices are filled only inside the designated area, preserving the type.
Example: A = matrix(3; 4)' - general matrix
fill_row(A; 2; 1)' = 0000 1111 0000
L = utriang(4)' - upper triangular matrix
$Repeat{fill_row(L; k; k) @ k = 1 : n_rows(L)}
L' = 1111 0222 0033 0004 - the lower triangle remains unchanged.

fill_col(M; j; x)

Parameters: M - the matrix which column is to be filled;
j - the index of the column to be filled;
x - the value to fill the column with.
Return value: the matrix M with the j-th column filled with value x.
Notes: The matrix M is modified in place and a reference to it is returned. This function is structurally compliant. Special matrices are filled only inside the designated area, preserving the type.
Example: A = matrix(3; 4)
fill_col(A; 2; 1)' = 0100 0100 0100
B = symmetric(4)' - symmetric matrix.
fill_col(B; 2; 1)' = 0000 0111 0100 0100 - the symmetry was preserved.

copy(A; B; i; j)

Parameters: A - source matrix;
B - destination matrix where the elements will be copied;
i - starting row index in matrix B;
j - starting column index in matrix B.
Return value: the matrix B, after the elements from A are copied, starting at indexes i and j in B.
Notes: Copying the elements from A to B modifies the matrix B in place. The existing elements in B are replaced by the respective elements from A. A reference to B is returned as result.
Example: A = [1; 2; 3|4; 5; 6]
B = mfill(matrix(3; 4); -1)' = -1-1-1-1 -1-1-1-1 -1-1-1-1
copy(A; B; 1; 1)
copy(A; B; 2; 2)
B' = 123-1 4123 -1456

add(A; B; i; j)

Parameters: A - source matrix;
B - destination matrix where the elements will be added;
i - starting row index in matrix B;
j - starting column index in matrix B.
Return value: the matrix B with the added elements from A to the respective elements in B, starting at indexes i and j in B.
Notes: Adding the elements from A to B modifies the matrix B in place. A reference to B is returned as result.
Example: A = [1; 2; 3|4; 5; 6]
B = mfill(matrix(3; 4); -1)' = -1-1-1-1 -1-1-1-1 -1-1-1-1
add(A; B; 1; 1)
add(A; B; 2; 2)
B' = 012-1 3572 -1345

row(M; i)

Parameters: M - matrix;
i - the index of the row to be extracted.
Return value: the i-th row of matrix M as a vector.
Notes: The full width row is returned, even for special matrices.
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
row(A; 3)' = [7 8 9]

col(M; j)

Parameters: M - matrix;
j - the index of the column to be extracted.
Return value: the j-th column of matrix M as a vector.
Notes: The full height column is returned, even for special matrices.
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
col(A; 2)' = [2 5 8]

extract_rows(M; ⃗i)

Parameters: M - matrix;
⃗i - vector containing the indexes of the rows to be extracted.
Return value: a new matrix containing the rows of matrix M whose indexes are in vector ⃗i.
Notes: The rows are extracted in the order, specified in vector ⃗i. It can be a result from the order_rows function.
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
extract_rows(A; [1; 2; 1])' = 123 456 123

extract_cols(M; ⃗j)

Parameters: M - matrix;
⃗j - vector containing the indexes of the columns to be extracted.
Return value: a new matrix containing the columns of M whose indexes are in vector ⃗j.
Notes: The columns are extracted in the order, specified in vector ⃗j. It can be a result from the order_cols function.
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
extract_cols(A; [3; 2; 1])' = 321 654 987

diag2vec(M)

Parameters: M - matrix.
Return value: a vector containing the diagonal elements of M.
Notes: The matrix M is not required to be square.
Example: A = [1; 2|3; 4|5; 6; 7|8; 9; 10]' = 120 340 567 8910
diag2vec(A)' = [1 4 7]

submatrix(M; i₁; i₂; j₁; j₂)

Parameters: M - matrix;
i₁ - starting row index;
i₂ - ending row index;
j₁ - starting column index;
j₂ - ending column index.
Return value: a new matrix containing the submatrix of M bounded by rows i₁ to i₂ and columns j₁ to j₂, inclusive.
Notes: Starting indexes can be greater than ending indexes. However, the same part of the matrix is returned. The result is always a general (rectangular) type of matrix, even if the source matrix is special.
Example: A = [ 1; 2; 3; 4| _
5; 6; 7; 8| _
9; 10; 11; 12| _
13; 14; 15; 16]
submatrix(A; 2; 3; 2; 4)' = 678 101112

Data functions

Data functions are not related to the structure of matrices but to the values of the elements. The following functions are available for use:

sort_cols(M; i)
rsort_cols(M; i)

Parameters: M - the matrix to be sorted;
i - the index of the row on which the sorting will be based.
Return value: a new matrix with the columns of M, sorted by the values in row i.
Notes: sort_cols sorts in ascending order, and rsort_cols - in descending order. The original matrix is not modified.
Example: A = [5; 2; 3|4; 9; 1|6; 8; 7]' = 523 491 687
B = sort_cols(A; 2)' = 352 149 768
C = rsort_cols(A; 2)' = 253 941 867
A' = 523 491 687 - the original matrix is unchanged.

sort_rows(M; j)
rsort_rows(M; j)

Parameters: M - the matrix to be sorted;
j - the index of the column on which the sorting will be based.
Return value: A new matrix with the rows of M, sorted by values in column j.
Notes: sort_rows sorts in ascending order, and rsort_rows sorts in descending order. The original matrix is not modified.
Example: A = [5; 2; 3|4; 9; 1|6; 8; 7]' = 523 491 687
B = sort_rows(A; 2)' = 523 687 491
C = rsort_rows(A; 2)' = 491 687 523
A' = 523 491 687 - the original matrix is unchanged.

order_cols(M; i)
revorder_cols(M; i)

Parameters: M - the matrix to be ordered;
i - the index of the row on which the ordering will be based.
Return value: a vector containing the indexes of the columns of matrix M ordered by the values in row i.
Notes: order_cols returns the indexes in ascending order, and revorder_cols - in descending order. The matrix is not modified. Each index in the output vector ⃗j shows which column in M should be placed at the current position to obtain a sorted sequence in row i. You can get the sorted matrix by calling extract_cols (M ; ⃗j ).
Example: A = [5; 2; 3|4; 9; 1|6; 8; 7]' = 523 491 687
b = order_cols(A; 2)' = [3 1 2]
B = extract_cols(A; b)' = 352 149 768
c = revorder_cols(A; 2)' = [2 1 3]
C = extract_cols(A; c)' = 253 941 867

order_rows(M; j)
revorder_rows(M; j)

Parameters: M - the matrix to be ordered;
j - the index of the column on which the ordering will be based.
Return value: a vector containing the indexes of the rows of matrix M ordered by the values in column j.
Notes: order_rows returns the indexes in ascending order, and revorder_rows - in descending order. The matrix is not modified. Each index in the output vector ⃗i shows which row in M should be placed at the current position to obtain a sorted sequence in column j. You can get the sorted matrix by calling extract_rows (M ; ⃗i ).
Example: A = [5; 2; 3|4; 9; 1|6; 8; 7]' = 523 491 687
b = order_rows(A; 2)' = [1 3 2]
B = extract_rows(A; b)' = 523 687 491
c = revorder_rows(A; 2)' = [2 3 1]
B = extract_rows(A; c)' = 491 687 523

mcount(M; x)

Parameters: M - matrix;
x - the value to count the occurrences of.
Return value: the number of occurrences of value x in matrix M.
Example: A = [1; 0; 1|2; 1; 2|1; 3; 1]' = 101 212 131
n = mcount(A; 1)' = 5

msearch(M; x; i; j)

Parameters: M - matrix;
x - the value to search for;
i - starting row index;
j - starting column index.
Return value: a vector containing the row and column indexes of the first occurrence of x in matrix M after indexes i and j, inclusive.
Notes: The searching is performed row by row, consecutively from left to right. If nothing is found, [0 0] is returned as result.
Example: A = [1; 2; 3|1; 5; 6|1; 8; 9]' = 123 156 189
b = msearch(A; 1; 1; 1)' = [1 1]
c = msearch(A; 1; 2; 2)' = [3 1]
d = msearch(A; 4; 1; 1)' = [0 0]

mfind(M; x)

Parameters: M - matrix;
x - the value to search for.
Return value: a two-row matrix containing the indexes of all elements in matrix M that are equal to x.
Notes: The top row in the result contains the row indexes, and the bottom row - the respective column indexes of the elements in M equal to x. If nothing is found, a 2×1 matrix with zeros is returned.
Example: A = [1; 2; 3|4; 1; 6|1; 8; 9]' = 123 416 189
B = mfind(A; 1)' = 123 121
C = mfind(A; 5)' = 0 0

hlookup(M; x; i₁; i₂)

Parameters: M - the matrix in which to perform the lookup;
x - the value to look for;
i₁ - the index of the row where to search for value x;
i₂ - the index of the row from which to return the corresponding values.
Return value: a vector containing the values from row i₂ of matrix M for which the respective elements in row i₁ are equal to x.
Notes: The values are collected consequently from left to right. If nothing is found, an empty vector [] is returned.
Example: A = [0; 1; 0; 1|1; 2; 3; 4; 5|6; 7; 8; 9; 10] ' = 01010 12345 678910
b = hlookup(A; 0; 1; 3)' = [6 8 10]
c = hlookup(A; 2; 1; 3)' = []

vlookup(M; x; j₁; j₂)

Parameters: M - the matrix in which to perform the lookup;
x - the value to look for;
j₁ - the index of the column where to search for value x;
j₂ - the index of the column from which to return the values.
Return value: a vector containing the values from column j₂ of matrix M for which the respective elements in column j₁ are equal to x.
Notes: The values are collected consequently from top to bottom. If nothing is found, an empty vector [] is returned.
Example: A = [1; 2|3; 4; 1|5; 6|7; 8; 1|9; 10]' = 120 341 560 781 9100
b = vlookup(A; 0; 3; 1)' = [1 5 9]
c = vlookup(A; 1; 3; 2)' = [4 8]
d = vlookup(A; 2; 3; 1)' = []

The find, hlookup and vlookup functions have variations with suffixes. Different suffixes refer to different comparison operators. They replace the equality in the original functions while the rest of the behavior remains unchanged. The possible suffixes are given in the table below:

suffix	mfind	hlookup	vlookup	comparison operator
`_eq`	`mfind_eq(M; x)`	`hlookup_eq(M; x; i₁; i₂)`	`vlookup_eq(M; x; j₁; j₂)`	`=` - equal
`_ne`	`mfind_ne(M; x)`	`hlookup_ne(M; x; i₁; i₂)`	`vlookup_ne(M; x; j₁; j₂)`	`≠` - not equal
`_lt`	`mfind_lt(M; x)`	`hlookup_lt(M; x; i₁; i₂)`	`vlookup_lt(M; x; j₁; j₂)`	`<` - less than
`_le`	`mfind_le(M; x)`	`hlookup_le(M; x; i₁; i₂)`	`vlookup_le(M; x; j₁; j₂)`	`≤` - less than or equal
`_gt`	`mfind_gt(M; x)`	`hlookup_gt(M; x; i₁; i₂)`	`vlookup_gt(M; x; j₁; j₂)`	`>` - greater than
`_ge`	`mfind_ge(M; x)`	`hlookup_ge(M; x; i₁; i₂)`	`vlookup_ge(M; x; j₁; j₂)`	`≥` - greater than or equal

Math functions

All standard scalar math functions accept matrix arguments as well. The function is applied separately to all the elements in the input matrix, even if it is of a special type. The results are returned in a corresponding output matrix. It is always a general (rectangular) type, so the structure is not preserved. For example:

#rad
M = diagonal(3; π/2)' = 1.57100 01.5710 001.571
cos(M)' = 011 101 110

Calcpad also includes a lot of math functions that are specific for matrices, as follows:

hprod(A; B)

Parameters: A - first matrix;
B - second matrix.
Return value: (matrix) the Hadamard product of matrices A and B.
Notes: Both matrices must be of the same size. The Hadamard (a.k.a. Schur) product C = A⊙B is an element-wise product of two matrices. The elements of the resulting matrix are obtained by the following equation:
C_ij = A_ij B_ij. If both matrices are of equal type, the type is preserved in the output, otherwise, the returned matrix if of general type.
Example: A = [1; 2|3; 4|5; 6]' = 12 34 56
B = [9; 8|7; 6|5; 4]' = 98 76 54
C = hprod(A; B)' = 916 2124 2524

fprod(A; B)

Parameters: A - first matrix;
B - second matrix.
Return value: (scalar) the Frobenius product of matrices A and B.
Notes: Both matrices must be of the same size. The result is obtained by summing the products of the corresponding element pairs of the input matrices:
p = m∑i= 1n∑j= 1A_ij B_ij
Example: A = [1; 2|3; 4|5; 6]' = 12 34 56
B = [9; 8|7; 6|5; 4]' = 98 76 54
C = fprod(A; B)' = 119

kprod(A; B)

Parameters: A - first matrix;
B - second matrix.
Return value: (matrix) the Kronecker product of matrices A and B.
Notes: If A is m×n and B is p×q matrix, the result is a mp×nq block matrix C, where each block is obtained by the equation: [C ]_ij = A_ij B.
Example: A = [1; 2|3; 4]' = 12 34
B = [5; 6; 7|8; 9; 10]' = 567 8910
C = kprod(A; B)' = 567101214 8910161820 151821202428 242730323640

mnorm_1(M)

Parameters: M - matrix.
Return value: scalar representing the L₁ (Manhattan, a.k.a. taxicab) norm of matrix M.
Notes: It is obtained as the maximum of all L₁ column vector norms by the equation: ||M ||₁ = max1 ≤ j ≤ n m∑i = 1| M_ij |
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
mnorm_1(A)' = 18

mnorm(M) or mnorm_2(M)

Parameters: M - an m×n matrix where m ≥ n.
Return value: scalar representing the L₂ (spectral) norm of matrix M.
Notes: It is obtained as the maximum singular value of matrix M :
||M ||₂ = σ_max (M ).
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
mnorm_2(A)' = 16.8481

mnorm_e(M)

Parameters: M - matrix.
Return value: scalar representing the Frobenius (Euclidean) norm of matrix M.
Notes: It is similar to vector Euclidian norm, assuming that the matrix is linearized by row concatenation. It calculated as the square root of sum of the squares of all elements, as follows:
||M ||_e = m∑i = 1n∑j = 1M_ij²
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
mnorm_e(A)' = 16.8819

mnorm_i(M)

Parameters: M - matrix.
Return value: scalar representing the L_∞ (infinity) norm of matrix M.
Notes: It is evaluated as the maximum of all L₁ row vector norms by the equation:
||M ||_∞ = max1 ≤ i ≤ m n∑j = 1| M_ij |
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
mnorm_i(A)' = 24

cond_1(M)

Parameters: M - square matrix.
Return value: (scalar) the condition number of M based on the L₁ norm.
Notes: It is calculated by the equation:
κ ₁(M ) = ||M ||₁· || M ^-1 ||₁
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
cond_1(A)' = 6.4852×10¹⁷

cond(M) or cond_2(M)

Parameters: M - an m×n matrix where m ≥ n.
Return value: (scalar) the condition number of M based on the L₂ norm.
Notes: The condition number shows how sensitive is the matrix A for solving the equation A⃗x = ⃗b or inverting the matrix. The higher is the number, the lower is the accuracy of the obtained solution. In theory, singular matrices have infinite condition numbers. However, in floating point environment, one can obtain very large, but finite number, due to floating point errors. The condition number is calculated by the following expression:
κ ₂(M ) = σ_max(M ) / σ_min(M )
Since this is computationally expensive, the other functions can be used instead, providing similar values but at lower computational cost.
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
cond_2(A)' = 1.7159×10¹⁷

cond_e(M)

Parameters: M - square matrix.
Return value: (scalar) the condition number of M based on the Frobenius norm.
Notes: It is calculated by the expression:
κ _e(M ) = ||M ||_e· || M ^-1 ||_e
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
cond_e(A)' = 4.5618×10¹⁷

cond_i(M)

Parameters: M - square matrix.
Return value: (scalar) the condition number of M based on the L_∞ norm.
Notes: It is obtained by the equation:
κ _∞(M ) = ||M ||_∞· || M ^-1 ||_∞
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
cond_i(A)' = 8.6469×10¹⁷

det(M)

Parameters: M - square matrix.
Return value: (scalar) the determinant of matrix M.
Notes: To evaluate the result, the matrix is first decomposed to lower triangular (L) and upper triangular (U) matrices by LU factorization. Then, the determinant is obtained as the product of the elements along the main diagonal of the U matrix. Theoretically, for singular matrices, the determinant should be exactly zero. However, it is not recommended to use this criterion in practice. Due to floating point round-off errors, it can return a "small" but non-zero value. It is better to use rank or condition number instead.
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
det(A)' = 6.6613×10^-16

rank(M)

Parameters: M - an m×n matrix where m ≥ n.
Return value: (scalar) the rank of matrix M.
Notes: Rank represents the number of linearly independent rows in a matrix. It is evaluated by performing an SVD decomposition of the matrix and counting the non-zero singular values.
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
rank(A)' = 2

trace(M)

Parameters: M - square matrix.
Return value: (scalar) the trace of matrix M.
Notes: Trace is defined as the sum of the elements along the main diagonal of a square matrix.
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
trace(A)' = 15

transp(M)

Parameters: M - matrix.
Return value: (matrix) the transpose of M.
Notes: The transpose is obtained by swapping the rows and the columns of the matrix. If M is symmetric, the transpose is equal to the matrix itself, so just a copy of M is returned.
Example: A = [1; 2; 3|4; 5; 6|7; 8; 9]' = 123 456 789
transp(A)' = 147 258 369

adj(M)

Parameters: M - square matrix.
Return value: the adjugate of matrix M.
Notes: The adjugate of a square matrix is the transpose of the cofactor matrix. It is obtained by multiplying the inverse of M by the determinant of M : adj(M ) = M ^-1 · |M |. However, this is not applicable when the matrix is singular. In this case, if the size of the matrix is n ≤ 10, Laplace expansion is used: A_ij = C_ji = M_ji · (-1)^i + j. For n > 10, it is reported that the matrix is singular instead, because Laplace expansion has a complexity of O(n!) and would take unreasonably long.
Example: A = [1; 2|3; 4]' = 12 34
adj(A)' = -42 3-1
det(A)' = 2
inverse(A)' = -21 1.5-0.5

cofactor(A)

Parameters: A - square matrix.
Return value: the cofactor matrix of A.
Notes: The cofactor C_ij is defined as the respective minor M_ij is multiplied by (-1)^i + j. The minor M_ij represents the determinant of the submatrix, obtained by removing the i-th row and the j-th column. In Calcpad, the cofactor matrix is calculated by transposing the adjugate: C = adj(A )^T.
Example: cofactor([1; 2|3; 4])' = -43 2-1

eigenvals(M; n_e)

Parameters: M - symmetric matrix;
n_e - the number of eigenvalues to extract (optional).
Return value: a vector containing the ne lowest (or largest) eigenvalues of matrix M.
Notes: If n_e > 0 the first n_e lowest eigenvalues are returned in ascending order. If n_e < 0 the first n_e largest eigenvalues are returned in descending order. If n_e = 0 or omitted, all eigenvalues are returned in ascending order. If M is a high-performance symmetric matrix and the number of rows is > 200 calculations are performed by using iterative symmetric Lanczos solver. Otherwise, direct symmetric QL algorithm with implicit shifts is used.
Example: A = copy( _ [4; 12; -16| _ 12; 37; -43| _ -16; -43; 98]; _ symmetric(3); 1; 1)' = 412-16 1237-43 -16-4398
eigenvals(A)' = [0.0188 15.5 123.48]

eigenvecs(M; n_e)

Parameters: M - symmetric matrix;
n_e - the number of eigenvalues to extract (optional).
Return value: a n_e×n matrix which rows represent the eigenvectors of M.
Notes: The same rules as for eigenvals(M ; n_e ) apply.
Example: A = copy( _ [4; 12; -16| _ 12; 37; -43| _ -16; -43; 98]; _ symmetric(3); 1; 1)' = 412-16 1237-43 -16-4398
eigenvecs(A)' = 0.9634‑0.26480.0411 ‑0.2127-0.849‑0.4838 ‑0.163‑0.45730.8742

eigen(M; n_e)

Parameters: M - symmetric matrix;
n_e - the number of eigenvalues to extract (optional).
Return value: a n_e×(n + 1) matrix were each row contains one eigenvalue, followed by the elements of the respective eigenvector of matrix M.
Notes: In case both are needed, this function is more efficient than calling eigenvals(M ; n_e ) and eigenvecs(M ; n_e ) separately. The same rules as for eigenvals(M ; n_e ) apply.
Example: A = copy( _ [4; 12; -16| _ 12; 37; -43| _ -16; -43; 98]; _ symmetric(3); 1; 1)' = 412-16 1237-43 -16-4398
eigen(A)' = 0.01880.9634‑0.26480.0411 15.504‑0.2127-0.849‑0.4838 123.48‑0.163‑0.45730.8742

cholesky(M)

Parameters: M - symmetric, positive-definite matrix.
Return value: (upper triangular matrix) the Cholesky decomposition of M.
Notes: This kind of decomposition turns the matrix into a product of two triangular matrices: M = L·L ^T. The current function returns only the upper part - L ^T. Cholesky decomposition is faster and more stable than the respective LU decomposition of the same matrix.
Example: A = copy( _ [4; 12; -16| _ 12; 37; -43| _ -16; -43; 98]; _ symmetric(3); 1; 1)' = 412-16 1237-43 -16-4398
LT = cholesky(A)' = 26-8 015 003 - the upper triangular matrix L^T
L = transp(LT)' = 200 610 -803 - the lower triangular matrix L
L*LT' = 412-16 1237-43 -16-4398 - check

lu(M)

Parameters: M - square matrix.
Return value: (matrix) the LU decomposition of M.
Notes: The LU decomposition factors the matrix as a product of two matrices: lower triangular L and upper triangular U (M = L·U ). It does not require the matrix to be symmetric. Since LU decomposition is not unique, there is an infinite number of solutions. Here is assumed that all elements along the main diagonal of L are equal to one. The function returns both matrices L and U, packed in a single square matrix. The elements along the main diagonal belong to U, since it is known that those of L are ones. The solution is performed by using the Crout's algorithm with partial pivoting. Instead of building the permutation matrix P, Calcpad internally creates a vector ⃗ind containing the indexes of the rows after reordering.
If you need both matrices L and U separately, you can extract them by a Hadamard product of the combined matrix by the respective lower/upper triangular matrix. After that, you have to reset the diagonal elements of L to ones.
If the type of M is symmetric matrix, the LDL^T decomposition is returned instead of LU. It is similar to Cholesky decomposition, but avoids taking square roots of diagonal elements. For that reason, the matrix is not required to be positive-definite. However, it makes it necessary to store the diagonal elements in a separate diagonal matrix D. Therefore, matrix M is represented as a product of three matrices:
M = L·D·L ^T. They are also packed in a single square matrix.
Example: A = [4; 12; -16|12; 37; -43|-16; -43; 98]' = 412-16 1237-43 -16-4398
LU = lu(A)' = 1237-43 -1.336.3340.67 0.333-0.05260.474 - the combined matrix
ind' = [2 3 1]
D = not(identity(3))' = 011 101 110 - the index permutation vector
L = hprod(mfill(ltriang(3); 1); LU)^D' =
100 -1.33310 0.3333-0.052631 - extracts the lower triangular matrix
U = hprod(mfill(utriang(3); 1); LU)' =
1237-43 06.33340.667 000.4737 - extracts the upper triangular matrix
extract_rows(L*U; order(ind))' = 412-16 1237-43 -16-4398 - check

qr(M)

Parameters: M - square matrix.
Return value: (matrix) the QR decomposition of matrix M.
Notes: As the name implies, the matrix is factored as a product (M = Q·R ) of an orthonormal matrix Q and upper triangular matrix R. The Householder's method is used for that purpose. The algorithm is stable and does not need pivoting. Both matrices are packed in a single n×2n block rectangular matrix [Q,R] and returned as a result. You can each of the two matrices using the submatrix function.
Example: A = [4; 12; -16|12; 37; -43|-16; -43; 98]' = 412-16 1237-43 -16-4398
QR = qr(A)' =
-0.196-0.1690.966-20.4-57.85105.31 -0.588-0.768-0.2540-3.86-24.85 0.784-0.6180.0508000.457 - the combined QR matrix
Q = submatrix(QR; 1; 3; 1; 3)' =
-0.1961-0.16950.9658 -0.5883-0.7676-0.2542 0.7845-0.61810.05083 - extracts the Q matrix
R = submatrix(QR; 1; 3; 4; 6)' =
-20.396-57.854105.314 0-3.858-24.853 000.4575 - extracts the R matrix
Q*R' = 412-16 1237-43 -16-4398 - check

svd(M)

Parameters: M - an m×n matrix, where m ≥ n.
Return value: the singular value decomposition (SVD) of matrix M.
Notes: The matrix is factored as a product of three matrices: M = U·Σ·V ^T, where Σ is a diagonal matrix containing the singular values of M, and U and V ^T are orthonormal matrices which columns represent the left and right singular vectors, respectively. The result is returned as a single m×(2n + 1) matrix [U, Σ, V], where Σ is a single column containing all singular values. They are sorted in descending order and the singular vectors are reordered respectively to match the corresponding singular values. Note that the returned V ^T matrix is already transposed, so you do not need to transpose it again.
Occasionally, some singular vectors may flip signs so that U·Σ·V ^T will not give M after multiplying the obtained matrices. Sign ambiguity is a well-known and common problem of most SVD algorithms. For symmetric matrices the singular values are equal to the absolute eigenvalues: σ_i = | λ_i|
Example: A = [4; 12; -16|12; 37; -43|-16; -43; 98]' = 412-16 1237-43 -16-4398
SVD = svd(A)' =
-0.1630.2127-0.9634123.477-0.163-0.45730.8742 -0.45730.8490.264815.5040.21270.8490.4838 0.87420.4838-0.04110.0188-0.96340.2648-0.0411 - the combined matrix
U = submatrix(SVD; 1; 3; 1; 3)' =
-0.1630.2127-0.9634 -0.45730.8490.2648 0.87420.4838-0.0411 - extracts the U matrix
V = submatrix(SVD; 1; 3; 5; 7)' =
-0.163-0.45730.8742 0.21270.8490.4838 -0.96340.2648-0.0411 - extracts the V matrix
σ = col(SVD; 4)' = [123.477 15.504 0.0188] - extracts the singular values
Σ = vec2diag(σ)' =
123.47700 015.5040 000.0188 - composes the singular value matrix

inverse(M)

Parameters: M - square matrix.
Return value: the inverse of matrix M.
Notes: The inverse is obtained by LU decomposition for non-symmetric matrices and LDL^T decomposition for symmetric ones. If the matrix is singular, the inverse does not exist. If it is ill conditioned, the result will be distorted by large errors. This is detected during the LU decomposition by watching for a zero or tiny pivot element. If one is detected, an appropriate error message is returned, instead of erroneous solution.
Example: A = [4; 12; -16|12; 37; -43|-16; -43; 98]' = 412-16 1237-43 -16-4398
B = inverse(A)' = 49.361-13.5562.111 -13.5563.778-0.5556 2.111-0.55560.1111
A*B' = 100 010 001 - check

lsolve(A; ⃗b)

Parameters: A - a square matrix with the equation coefficients;
⃗b - the right-hand side vector.
Return value: the solution vector ⃗x of the system of linear equations A⃗x = ⃗b.
Notes: Calculations are performed by using LU decomposition for non-symmetric matrices and LDL^T for symmetric. That is why the matrix is not required to be positive-definite. If A is singular or ill-conditioned, an error message is returned.
Example: A = [8; 6; -4|6; 12; -3|-4; -3; 9]' = 86-4 612-3 -4-39
b = [10; 20; 30]' = [10 20 30]
x = lsolve(A; b)' = [2.5 1.667 5] - the solution vector
A*x' = [10 20 30] - check

clsolve(A; ⃗b)

Parameters: A - symmetric, positive-definite coefficient matrix;
⃗b - the right-hand side vector.
Return value: the solution vector ⃗x of the system of linear equations A⃗x = ⃗b using Cholesky decomposition.
Notes: Cholesky decomposition is faster than LU and LDL^T, so this function should be preferred over lsolve whenever the matrix is symmetric and positive-definite.
Example: A = copy([8; 6; -4|6; 12; -3|-4; -3; 9]; _ symmetric(3); 1; 1)' = 86-4 612-3 -4-39
b = [10; 20; 30]' = [10 20 30]
x = clsolve(A; b)' = [2.5 1.667 5]
A*x' = [10 20 30] - check

slsolve(A; ⃗b)

Parameters: A - high-performance symmetric, positive-definite coefficient matrix;
⃗b - high-performance right-hand side vector.
Return value: the solution vector ⃗x of the system of linear equations A⃗x = ⃗b using preconditioned conjugate gradient (PCG) method.
Notes: The PCG method is iterative and is faster than the direct Cholesky decomposition method. Whenever you have larger systems with tens to hundreds of thousands equations it is recommended to use the slsolve function with high performance matrices and vectors.
Example: A = copy([8; 6; -4|6; 12; -3|-4; -3; 9]; _ symmetric_hp(3); 1; 1)' = 86-4 612-3 -4-39
b = hp([10; 20; 30])' = [10 20 30]
x = slsolve(A; b)' = [2.5 1.667 5]
A*x' = [10 20 30] - check

msolve(A; B)

Parameters: A - a square matrix with the equation coefficients;
B - the right-hand side matrix.
Return value: the solution matrix X of the generalized matrix equation AX = B.
Notes: Similar to lsolve, having that B matrix columns contain multiple right-hand side vectors and X matrix columns represent the respective solution vectors. In this way, the function can solve multiple systems of linear equations in parallel. The LU/LDLT decomposition of A is performed only once in the beginning and the result is reused for backsubstitution multiple times.
Example: A = [8; 6; -4|6; 12; -3|-4; -3; 9]' = 86-4 612-3 -4-39
B = join_cols([10; 20; 30]; [40; 50; 60])' = 1040 2050 3060
X = msolve(A; B)' = 2.58.71 1.672.67 511.43
A*X' = 1040 2050 3060 - check

cmsolve(A; B)

Parameters: A - symmetric, positive-definite coefficient matrix;
B - the right-hand side matrix.
Return value: the solution matrix X of the generalized matrix equation AX = B using Cholesky decomposition.
Notes: Similar to clsolve, having that B matrix columns contain multiple right-hand side vectors and X matrix columns represent the respective solution vectors. In this way, the function can solve multiple systems of linear equations in parallel. The Cholesky decomposition of A is performed only once in the beginning and the result is reused for backsubstitution multiple times.
Example: A = copy([8; 6; -4|6; 12; -3|-4; -3; 9]; _ symmetric(3); 1; 1)' = 86-4 612-3 -4-39
B = join_cols([10; 20; 30]; [40; 50; 60])' = 1040 2050 3060
X = cmsolve(A; B)' = 2.58.71 1.672.67 511.43
A*X' = 1040 2050 3060 - check

smsolve(A; B)

Parameters: A - high-performance symmetric, positive-definite coefficient matrix;
B - high-performance right-hand side matrix.
Return value: the solution matrix X of the generalized matrix equation AX = B using preconditioned conjugate gradient (PCG) method.
Notes: Similar to slsolve, having that B matrix columns contain multiple right-hand side vectors and X matrix columns represent the respective solution vectors.
Example: A = copy([8; 6; -4|6; 12; -3|-4; -3; 9]; _ symmetric_hp(3); 1; 1)' = 86-4 612-3 -4-39
B = hp(join_cols([10; 20; 30]; [40; 50; 60]))' = 1040 2050 3060
X = smsolve(A; B)' = 2.58.71 1.672.67 511.43
A*X' = 1040 2050 3060 - check

fft(M)

Parameters: M - an hp matrix containing the input signal in time domain. It must have one row for real only values or two rows for complex values. The first row must contain the real part and the second row - the corresponding imaginary part.
Return value: the fast Fourier transform of the input signal in M. Similar to the input, the real part of the output is stored in the first row and the imaginary part - in the second row.
Notes: The fast Fourier transform is performed by using the classic Cooley-Turkey algorithm.
Example: A = hp([1; 2; 3; 4|0])' = 1234 0000
B = fft(A)' = 10-2-2-2 0-202

ift(M)

Parameters: M - an hp matrix containing the input signal in frequency domain. It must have one row for real only values or two rows for complex values. The first row must contain the real part and the second row - the corresponding imaginary part.
Return value: the inverse Fourier transform of the input signal in M from frequency domain to time domain. Similar to the input, the real part of the output is stored in the first row and the imaginary part - in the second row.
Notes: The inverse Fourier transform is performed by using the classic Cooley-Turkey algorithm.
Example: A = hp([10; -2; -2; -2|0; -2; 0; 2])' = 10-2-2-2 0-202
B = ift(A)' = 1234 0000

Aggregate and interpolation functions

All aggregate functions can work with matrices. Since they are multivariate, each of them can accept a single matrix, but also a list of scalars, vectors and matrices, mixed in arbitrary order. For example:

A = [0; 2| 4; 8]
b = [5; 3; 1]
sum(10; A; b; 11)' = 44

Interpolation functions behave similarly, when invoked with a mixed list of arguments. In this case matrices are linearized by rows into the common array of scalars at the respective places. Then, the interpolation is executed over the entire array. However, there is also a "matrix" version of these functions, that can perform double interpolation over a single matrix. In this case you must specify exactly three arguments: the first two must be scalars (the interpolation variables) and the third one must be a matrix, as follows:

take(x; y; M)

Parameters: x - row index;
y - column index;
M - matrix.
Return value: (scalar) the element of matrix M at indexes x and y.
Notes: If x and y are not integer, they are rounded to the nearest integer.
Example: A = [1; 2; 5|3; 6; 15|5; 10; 25]' = 125 3615 51025
take(2; 3; A)' = 15

line(x; y; M)

Parameters: x - interpolation variable across rows;
y - interpolation variable across columns;
M - matrix.
Return value: (scalar) the value obtained by double linear interpolation from the elements of matrix M based on the values x and y.
Example: A = [1; 2; 5|3; 6; 15|5; 10; 25]' = 125 3615 51025
line(1.5; 2.5; A)' = 7

spline(x; y; M)

Parameters: x - interpolation variable across rows;
y - interpolation variable across columns;
M - matrix.
Return value: (scalar) the value obtained by double Hermite spline interpolation from the elements of matrix M based on the values x and y.
Example: A = [1; 2; 5|3; 6; 15|5; 10; 25]' = 125 3615 51025
spline(1.5; 2.5; A)' = 6.625

You can use interpolation functions to plot matrix data, as in the example below:

$Map{spline(x; y; A) @ x = 1 : n_rows(A) & y = 1 : n_cols(A)}

A full list of the available aggregate and interpolation functions is provided earlier in this manual (see "Expressions/Functions" above).

Operators

All operators can be used with matrix operands. Both matrices must be of the same sizes. Operations are performed element-by-element and the results are returned in an output matrix. For example:

A = [0; 1; 2|3; 4; 5]' = 012 345
B = [11; 10; 9|8; 7; 6]' = 11109 876
A + B' = 111111 111111

The only exception is multiplication operator "*" which represents the standard matrix multiplication. The element-wise (Hadamard or Schur) product is implemented in Calcpad as a function: hprod. Matrix multiplication C_m×p = A_m×n B_n×p is defined so that each element c_i j is obtained as the dot product of the i-th row of A and the j-th column of B:

c_ij = n∑k = 1a_ik b_kj

For example:

A = [1; 2; 3|4; 5; 6]' = 123 456
B = [6; 5|4; 3|2; 1]' = 65 43 21
C = A*B' = 2014 5641

All binary operators are supported for mixed matrix-vector and vector-matrix operands. In this case, the vector is treated as column matrix. For example:

A = [1; 2; 3|4; 5; 6]' = 123 456
b = [6; 5; 4]
A*b' = [28 73]

Matrix-scalar and scalar-matrix operators are also implemented in a element-wise manner. The operation with the provided scalar is performed for each element and the result is returned as matrix. For example:

[1; 2|3; 4]*5' = 510 1520

Matrices

Table of contents