Go to the first, previous, next, last section, table of contents.


Matrices and Tensors

In JACAL, a matrix is just a bunch of equal length bunchs, and this is the structure that the matrix operations currently supported by JACAL (ncmult(), ^^, transpose(), etc.) expect.

Generating Matrices

Operator: bunch elt_1 elt_2 ...
[elt_1, elt_2, ...]

To collect any number of Jacal objects into a bunch, simply enclose them in square brackets. For example, to make the bunch whose elements are 1, 2, 4, type [1, 2, 4]. One can also nest bunches, for example, [1, [[1, 3], [2, 5]], [1, 4]]. Note however that the bunch whose only element is [1, 2, 3] is [1 2 3]. It is importance to notice that one has commas and the other doesn't.

e3 : a:bunch(1, 2, 3);

e3: [1, 2, 3]

e4 : b:[a];

e4: [1  2  3]

e5 : c:[b];

e5: [[1, 2, 3]]

e6 : [[[1, 2, 3]]];

e6: [[1, 2, 3]]

Command: ident n
The command ident takes as argument a positive integer n and returns an nxn identity matrix. This is sometimes more convenient than obtaining this same matrix using the command scalarmatrix.

e6 : ident(4);

    [1  0  0  0]
    [          ]
    [0  1  0  0]
e6: [          ]
    [0  0  1  0]
    [          ]
    [0  0  0  1]

Command: scalarmatrix size entry
The command scalarmatrix takes as inputs a positive integer size and an algebraic expression entry and returns an n * n diagonal matrix whose diagonal entries are all equal to entry, where n = size.

e1 : scalarmatrix(3, 6);

    [6  0  0]
    [       ]
e1: [0  6  0]
    [       ]
    [0  0  6]

Command: diagmatrix list
The Jacal command diagmatrix takes as input a list of objects and returns the diagonal matrix having those objects as diagonal entries. In case one wants all of the diagonal entries to be equal, it is more convenient to use the command scalarmatrix.

e3 : diagmatrix(12,3,a,s^2);

    [12  0  0  0 ]
    [            ]
    [0   3  0  0 ]
e3: [            ]
    [0   0  a  0 ]
    [            ]
    [0   0  0   2]
    [          s ]

e4 : diagmatrix([1,2],2);

    [[1, 2]  0]
e4: [         ]
    [  0     2]

Command: sylvester poly_1 poly_2 var
Here, poly_1 and poly_2 are polynomials and var is a variable. The function sylvester returns the matrix introduced by Sylvester (A Method of Determining By Mere Inspection the Derivatives from Two Equations of Any Degree, Phil.Mag. 16 (1840) pp. 132-135, Mathematical Papers, vol. I, pp. 54-57) for computing the resultant of the two polynomials poly_1 and poly_2 with respect to the variable var. If one wants to compute the resultant itself, one can simply use the command resultant with the same syntax.

e5 : sylvester(a0 + a1*x + a2*x^2 + a3*x^3, b0 + b1*x + b2*x^2, x);

    [a3  a2  a1  a0  0 ]
    [                  ]
    [0   a3  a2  a1  a0]
    [                  ]
e5: [b2  b1  b0  0   0 ]
    [                  ]
    [0   b2  b1  b0  0 ]
    [                  ]
    [0   0   b2  b1  b0]

Command: genmatrix function rows cols
The function genmatrix takes as arguments a function function of two variables and two positive integers, rows and cols. It returns a matrix with the indicated numbers of rows and columns in which the $(i,j)$th entry is obtained by evaluating function at $(i,j)$. The function may be defined in any of the ways available in Jacal, i.e previously by an explicit algebraic definition, by an explicit lambda expression or by an implicit lambda expression.

e4 : @1^2+@2^2;

                       2      2
e4: lambda([@1, @2], @1  + @2 )

e5 : genmatrix(e4,3,5);

    [2   5   10  17  26]
    [                  ]
e5: [5   8   13  20  29]
    [                  ]
    [10  13  18  25  34]

Matrix Parts

Command: row matrix i
The command row returns the ith row of the matrix matrix, where i = int. If int is larger than the number of rows of matrix, then Jacal prints an error message. The corresponding command for columns of a matrix is col.

e3 : u:[[1, 2, 3], [1, 5, 3]];

    [1  2  3]
e3: [       ]
    [1  5  3]

e4 : row(u, 2);

e4: [1, 5, 3]

Command: col matrix integer
The command col is used to extract a column of a matrix. Here, matrix is a matrix and integer is a positive integer. If that integer exceeds the number of columns, an error message such as

ERROR: list-ref: Wrong type in arg1 ()

appears. Here is an example of correct use of the command col:

e19 : a:[[1,2,4],[2,5,6]];

     [1  2  4]
e19: [       ]
     [2  5  6]

e20 : col(a,2);

     [2]
e20: [ ]
     [5]

Command: minor matrix row col
The command minor returns the submatrix of matrix obtained by deleting the ith row and the jth column, where i = row and j = col.

e21 : b:[[1,2,3],[3,1,5],[5,2,7]];

     [1  2  3]
     [       ]
e21: [3  1  5]
     [       ]
     [5  2  7]

e22 : minor(b,3,1);

     [2  3]
e22: [    ]
     [1  5]

Command: rapply bunch int_1 int_2 ...
The function rapply is used to access elements of bunches. It can also access elements nested at lower levels in a bunch. In particular, it can also access matrix elements. In the above syntax, bunch is the bunch whose parts one wishes to access, and n, int_1, int_2, ..., int_n are positive integers. It returns the int_n-th element of the int_{n-1}-th element of ... of the int_2-th element of the int_1-th element of bunch. One can have n = 0. In that case, rapply simply returns the bunch.

e2 : rapply([[1,2,3],[1,4,6],3],2,3);

e2: 6

e6 : rapply([a,b],2);

e6: b

e7 : rapply([a,b]);

e7: [a, b]

Matrix commands

Command: . matrix1 matrix2

Matrix multiplication.

e1 : a:[[1, 2, 3], [5, 2, 7]];

    [1  2  3]
e1: [       ]
    [5  2  7]

e2 : b:[[3, 2], [6, 4]];

    [3  2]
e2: [    ]
    [6  4]

e3 : b . a;

    [13  10  23]
e3: [          ]
    [26  20  46]

Command: ^^ matrix exponent

The infix operator ^^ is used for raising a square matrix to an integral power. It can also be used on nonsquare matrices, but the results are not documented.

e8 : a:[[1, 0], [-1, 1]];

    [1   0]
e8: [     ]
    [-1  1]

e9 : a^^3;

    [1   0]
e9: [     ]
    [-3  1]

Command: dotproduct vector_1 vector_2
The Jacal function dotproduct returns the dot product of two row vectors of the same length. It will also give the dot product of two matrices of the same size by computing the sum of the dot products of the corresponding rows or, what is the same, the trace of one matrix times the transpose of the other one.

e28 : a:[1,2,3]; b:[3,1,5];

e28: [1, 2, 3]

e29 :
e29: [3, 1, 5]

e30 : dotproduct(a,b);

e30: 20

Command: crossproduct vector_1 vector_2
The Jacal command crossproduct computes the cross product of two vectors. By definition, the two vectors must each have three components.

e24: [2 x, y + 2 x y ]

e25 : crossproduct([1,2,3],[4,2,5]);

e25: [4, 7, -6]

Command: determinant matrix
The Jacal command determinant computes the determinant of a square matrix. Attempting to take the determinant of a non-square matrix will produce an error message.

e1 : a:[[1,2],[6,7]];

    [1  2]
e1: [    ]
    [6  7]

e2 : determinant(a);

e2: -5

Command: transpose matrix
Computes the transpose of (matrix).

Tensors

The tensors supported by JACAL are an extension of the matrix structure (i.e., a bunch of bunches of bunches ...) with the added stipulation that all dimensions of the tensor be the same length (e.g., 4x4x4). The number of dimensions (indices) in a tensor is its rank: A scalar is a tensor of rank 0; a vector is a rank 1 tensor; a matrix has rank 2; and so on.

Further, just as matrix binary operations place restrictions on the matrices involved (e.g., the row/column length requirement for matrix multiplication), the tensor binary operations require that the dimensions of each tensor be of the same length. For example, you could not multiply a 3x3 tensor and a 4x4x4 tensor.

JACAL's tensors do not support the construct of contravariant and covariant indices. Users must keep track of this information themselves, and perform the necessary operations with an appropriate metric so that the "index gymnastics" is performed correctly.

Before using any of JACAL's tensor operations, execute the following command from the JACAL prompt:

require("tensor");

This loads the file `tensor.scm' into JACAL, and makes the tensor operations available for use.

JACAL currently supports four tensor operations: tmult, contract, indexshift, and indexswap. Each of these is described in detail below.

Tensor Multiplication

Command: tmult matrix_1 matrix_2 index_1 index_2
tmult takes a minimum of two arguments which are the tensors on which the multiplication operation is to be performed.

With no additional arguments, tmult will produce the outer product of the two input tensors. The rank of the resulting tensor is the sum of the inputs' ranks, and the components of the result are formed from the pair-wise products of components of the inputs. For example, for the input tensors x[a,b] and y[c]

@result{z:tmult(x,y);} z[a,b,c] = x[a,b]*y[c]

With an additional argument, tmult will produce the inner product of the two tensors on the specified index. For example, given x[i,j] and y[k,l,m]

z:tmult(x,y,3);
=>
                        length
                        -----
                        \
             z[a,b,c] =  >   x[a,q] * y[b,c,q]
                        /
                        -----
                        q = 1

Note that in this case x only has 2 indices. All of JACAL's tensor operations modify index inputs to be between 1 and the rank of the tensor. Thus, in this example, the 3 is modified to 2 in the case of x. As another example, with x[i,j,k] and y[l,m,n]

z:tmult(x,y,2);
=>
                          length
                          -----
                          \
             z[a,b,c,d] =  >   x[a,q,b] * y[c,q,d]
                          /
                          -----
                          q = 1

With four arguments, tmult produces an inner product of the two tensors on the specified indices. For example, for x[i,j] and y[k,l,m]

z:tmult(x,y,1,3);
=>
                        length
                        -----
                        \
             z[a,b,c] =  >   x[q,a] * y[b,c,q]
                        /
                        -----
                        q = 1

Note that matrix multiplication is the special case of an inner product (of two "two dimensional matrices") on the second and first indices, respectively: tmult(x,y,2,1) == ncmult(x,y)

Finally, tmult handles the case of a scalar times a tensor, in which case each component of the tensor is multiplied by the scalar.

Tensor contraction

Command: contract matrix index1 ...

The contraction operation produces a tensor of rank two less than a given tensor. It does this by performing a summation over two of the indices of the given tensor, as clarified in the examples below.

contract takes at least one argument which is the tensor on which the contraction operation is to be performed. One or two additional arguments may be provided to specify the indices to be used in the summation. If no additional arguments are provided, the summation is performed over the first and second indices. With one additional argument, the summation is over the specified index and the one following it (e.g., if 3 is specified, the third and fourth indices are used). With two additional arguments, the summation is performed over the indices specified. The actual indices used will be constrained to be between 1 and the rank of the tensor.

Examples:

1) For a square matrix (tensor of rank 2), contract returns a scalar that is the sum of the diagonal elements of the matrix.

2) Given x[i,j,k,l], the command

y:contract(x,2,4);

produces:

                        length
                        -----
                        \
               y[a,b] =  >   x[a,q,b,q]
                        /
                        -----
                        q = 1

Special cases: If contract is given a scalar (rank 0 tensor) as input, it just returns the scalar. For a vector (tensor of rank 1), contract returns a scalar that is the sum of the elements of the vector.

Shifting of Tensor Indices

Command: indexshift matrix index1 ...

indexshift rearranges the indices of a tensor. It is one of two generalizations of the matrix transpose operation (cf. indexswap).

indexshift takes at least one argument which is the tensor on which the index shifting is to be performed. One or two additional arguments may be provided to specify the index and the position to which it is to be shifted. If no additional arguments are provided, the first index of the tensor is shifted to the second position (equivalent the matrix transpose operation). If one additional argument is provided, it specifies the index to be shifted, and that index will be shifted "to the right" one position (e.g., if 3 is specified, the third index will be shifted to the forth position). If two additional arguments are provided, the first specifies the index and the second specifies the position to which it is to be shifted. The actual index shifted and its shifted position will be constrained to be between 1 and the rank of the tensor.

For example, given x[a,b,c,d], the command y:indexshift(x,1,3); produces a tensor y such that y[a,b,c,d] == x[b,c,a,d]. In this example, the element that was in position [a,b,c,d] in x will be in position [b,c,a,d] in y.

Special cases: If indexshift is given a scalar (rank 0 tensor) as input, it just returns the scalar. For a vector (tensor of rank 1), indexshift transposes the 1-by-n matrix (row vector) to an n-by-1 matrix (column vector).

Swapping of Tensor Indices

Command: indexswap tensor ...

indexswap rearranges the indices of a tensor. It is one of two generalizations of the matrix transpose operation (cf. indexshift).

indexswap takes at least one argument which is the tensor on which index swapping is to be performed. One or two additional arguments may be provided to specify the indices to be swapped. If no additional arguments are provided, the first and second indices of the tensor are swapped (equivalent the matrix transpose operation). With one additional argument, the specified index is swapped with the one following it (e.g., if 2 is specified, the second and third indices will be swapped). If two additional arguments are provided, they specify the indices to be swapped. The actual indices used will be constrained to be between 1 and the rank of the tensor.

For example, given x[a,b,c,d], the command y:indexswap(x,2,4); produces a tensor y such that y[a,b,c,d] = x[a,d,c,b]. In this example, the element that was in position [a,b,c,d] in x will be in position [a,d,c,b] in y.

Special cases: If indexswap is given a scalar (rank 0 tensor) as input, it just returns the scalar. For a vector (tensor of rank 1), indexswap transposes the 1-by-n matrix (row vector) to an n-by-1 matrix (column vector).


Go to the first, previous, next, last section, table of contents.