The mcx matrix calculator
Purpose
mcx is a matrix calculator that integrates with other components of the Questaal package. It can read data files written in several formats, include the standard Questaal format which features programming language capabilities.
Table of Contents
Preliminaries
mcx is required and is assumed to be in your path.
This manual is written for version 1.057. To see what version you are using, do:
$ mcx version
1. Introduction
mcx is a matrix extension of an ordinary calculator. It is mainly designed to work with 2D arrays; it operates on numerical arrays. Usually the arrays are stored in ASCII format, though mcx has the ability to read binary files.
mcx is commandline driven and can efficiently manipulate arrays for convenient analysis. Some Questaal testing scripts make use of this calculator to analyze whether a test passes a certain tolerance.
mcx is stackbased: each time you read an array from a file it is pushed onto the stack. You can assign the toplevel array to a name, in which case an independent copy is made.
The ASCII format of the file is quite flexible; you can specify the number of rows nr and number of columns nc in a number of ways:

You can specify the number of columns on the command line, with −nc=#; similarly with the number of rows (−nr=#). −nc=#; this switch is used in Example 2.2.

You can include a directive at the start of the file, before any data is read. mcx passes the input through a file preprocessor first; in addition a directive specifying some combination nr and nc, e.g.
% rows #1 cols #2
supplies the needed information to mcx.

In the absence of an explicit specification, mcx will infer nc by counting how many numbers (more generally expressions) are on the first line.
If nc is obtained somehow but not nr, mcx determines it by counting the number of expressions in the entire file and dividing by nc.
Note: : if the −nc=# or −nr=# switches are not used, mcx determines them from the standard Questaal algorithm for reading 2D arrays, described in more detail here.
2. Examples
This section develops a few examples to provide an intuitive feel for how mcx works and to illustrate some features.
A systematic description of mcx’s features and arguments is given in Section 3.
Cut and paste the following data into array into file mat1
1.1 2.2
3.3 4.4
and this data into file mat2.
5 6
7 8
Example 2.1. Add mat1 to mat2
If you do the following:
$ mcx mat1 mat2 +
you should see
% rows 2 cols 2 real
6.100000 8.200000
10.300000 12.400000
Any argument that begins with “−“ is a switch or an operator (as distinct from data). The string following “−“ names the operator.
Thus:

mat1 : reads file mat1 from disk and pushes it onto the stack; call it s_{0}.

mat2 : reads file mat2 and pushes it onto the stack so s_{0} → s_{1} and the contents of mat2 → s_{0}.

−+ : is the binary operation “add.” If s_{0} and s_{1} exist, they are summed and popped off the stack. Their sum is pushed onto toplevel array s_{0}.
Note: if operations occur before arrays are given, they are push on an operations stack and execute when arrays become available. The following commands accomplish the same thing
$ mcx + mat1 mat2
$ mcx mat1 + mat2
$ mcx mat1 mat2 +
Use −show to see the stack. Try this
$ mcx f2f8.1 mat1 mat2 show +
You should get
# 0 named arrays, 2 on stack; pending 0 unops 0 bops (vsn 1.057)
# stack nr nc cast
# 2 2 2 real
# 1 2 2 real
% rows 2 cols 2 real
6.1 8.2
10.3 12.4
−show prints out what arrays are in the calculator; the addition is subsequently performed and the result printed.
“−f…“ is a formatting statement; 2f8.1 follows the fortran convention for formatting real numbers.
Example 2.2 Standard input and standard output
After all the commandline arguments are parsed, usually mcx prints the toplevel array s_{0}, and exits silently.
However,

If s_{0} is not present, mcx waits for you to enter an array from standard input.
For example try$ mcx f2f8.1 $ 1 2 $ 3 $ 4 <ctrlD>
mcx should print out
% rows 2 cols 2 real 1.0 2.0 3.0 4.0

If the last command is −show, mcx just prints out information about named arrays and the stack and exits (even if the stack has no arrays). Thus the following command
$ mcx f2f8.1 mat1 mat2 show + show
does the same as it did in Example 2.1 but prints out stack contents a second time instead of s_{0}.
You can also tell mcx to read the next array from standard input by using a full stop (“ . “) in lieu of file name. Thus
$ echo 1 2 3 4  mcx .
pushes a 1×4 array onto s_{0}, while
$ echo 1 2 3 4  mcx nc=2 . mat2 +
pushes a 2×2 array onto s_{0}, then pushes mat2, and adds the two so that a single array remains on the stack. −nc=2 tells mcx to treat the line from stdin as an array with 2 columns.
Example 2.3 Manipulations of eigenvalues and eigenvectors of an array
This example finds the eigenvalues e and eigenvectors z of mat2, and shows that mat2 z = z e.
First try
$ mcx mat2 evl
You should see a 2×1 complex array
% rows 2 cols 1 complex
0.152067
13.152067
0.000000
0.000000
The eigenvalues are complex because the matrix is not hermitian. The imaginary part follows the real part; this is how mcx displays and reads complex arrays in ASCII format.
You can force mat2 to be hermitian (symmetric since mat2 is real) with −herm. Now do
$ mcx mat2 herm a h h evl h evc ap z tog v2dia x h z x 
These instructions do the following:

mat2 : reads file mat2 and pushes it onto s_{0}.

−herm : symmetrizes s_{0}.

−a h : copies s_{0} to a named array h and pops it from the stack. The stack is now empty.

h −evl : pushes h onto s_{0} and replaces s_{0} with its eigenvalues. Because h is hermitian, s_{0} is real.

h −evc : pushes h onto the stack and replaces s_{0} with its eigenvectors. Now the stack has two arrays, s_{0} = z and s_{1} = e.

−ap z : copies s_{0} to a named array z.

−tog : toggles s_{0} and s_{1}, making s_{1} = z and s_{0} = e.

−v2dia : Turns s_{0} (a vector or 2×1 array of eigenvalues) into a diagonal 2×2 array.

−x : multiplies s_{1} × s_{0}. One array remains on the stack, s_{0} = z × e.

h z −x : pushes h and z onto the stack and multiplies them. Now s_{0} = h × z, while s_{1} = z × e

−− : Adds s_{1} to −s_{0}. s_{1} and s_{1} are mathematically identical so the difference should be zero.
Note: this formula should still work even if h is not hermitian.
Example 2.4 Numerical integration, differentiation, and interpolation of a function
Integrate and differentiate the function $y = x^n e^{\lambda x}$, by tabulating it on a mesh and evaluating integrals and derivatives numerically. This example also uses the tabulated data to interpolate it to another mesh.
Copy the contents of the box below into file dat.
% const n=100 p=1 lam=2
% save
% macro iy(z) exp(lam*z)*(lam*z1)/lam^2
% repeat i= 0:n
% var x=10*i/n
{x} {x^p*exp(lam*x)}
% end
This creates 101 rows of xy pairs with x ranging between 0 and 10.
Note: : The number of points n+1, and also p and λ are declared in the file with the const preprocessor directive. You can override the values assigned there with commandline arguments, e.g. −vp=#. The save directive retains the variables declared in this file after it is read. The macro directive will be used for the indefinite integral, below.
Derivative
The derivative is readily found to be
$y' = {x^{p  1}}\,\left( {p  \lambda x} \right)\,{e^{  \lambda x}}$Try some of the following commands
$ mcx vp=2 dat diff
$ mcx vp=2 dat diff e3 x1 x2 'x1^(p1)*(plam*x1)*exp(lam*x1)'
$ mcx vp=2 f3f15.10 dat diff:nord=5 e3 x1 x2 'x1^(p1)*(plam*x1)*exp(lam*x1)' e3 x1 x2 x2x3
All of them differentiate the second column with respect to the first, using p=2.
 The first returns two columns with x and y′.
 The second returns three columns with x, y′, and the analytic derivative of y′.
 The third returns three columns (more decimals) with x, y′, and the error in the numerical estimate for y′.
Unpacking the third command:
Argument  Function 
−vp=2  Declares variable p and assigns the value 2. This overrides the assignment in dat. 
−f3f15.10  Formats output (fortran format 3f15.10) 
dat  read s_{0} from dat. 
−diff:nord=5  Replace column 2 with a numerical estimate for y′ Use a 5point polynomial to interpolate the data; estimate is derivative of polynomial interpolation 
−e3 x1 x2 ‘x1^(p−1)*(p−lam*x1)*exp(−lam*x1)’  replace s_{0} with a three column array consisting of x, y′, and the analytic derivative of y′ 
−e3 x1 x2 x2−x3  replace s_{0} with a three column array consisting of x, y′, and difference between the numerical and analytical y′ 
Some observations:
 The largest error appears for x→0. The interpolation is less accurate when all the data lies on one side of the interpolating point.
 The error improves when higher order polynomials (nord=#) are used, or when the mesh is made finer (−vn=#).
 If you use p<1, y′ diverges at the origin. The numerical derivative cannot reproduce this.
Integral
$I = \int_0^\infty x^p e^{\lambda{x}} \mathrm{d}x = \lambda^{p} \Gamma(p+1) \rightarrow p!/\lambda^p \mathrm{\quad if\ p\ is\ an\ integer}$y is small at x=10 if λ=−2, provided p is 3 or less; so we will use 10 in place of ∞.
Try some of the following commands
$ mcx vp=1 dat int 0 10
$ mcx vp=2 dat int 0 10
$ mcx vp=3 dat int 0 10
These commands calculate I between 0 and 10 for three values of p.
You should find that I is close to $p!/\lambda^p$, i.e. 1/4, 1/4, and 3/8 for p=1,2,3. Some numerical errors appear in the 6^{th} digit. The integral is carried out by fitting the data to a polynomial of order nord−1, and integrating the polynomial. You can reduce the error by using a higher order than the default value of 4 for nord, viz
$ mcx vp=3 dat int:nord=6 0 10
Also for p>3 you should increase the upper bound beyond 10 since the integral from 10 to ∞ is on the order of 10^{−6}.
Indefinite Integral
As noted I must be evaluated between definite limits. However, you can make mcx simulate an indefinite integral by integrating over a range of upper bounds.
To compare with exact results, note that when p=1 the indefinite integral is
$I = \int x e^{\lambda{x}} \mathrm{d}x =  {e^{  \lambda{x}}} \left( \frac{\lambda{x} + 1}{\lambda^2} \right)$The macro in dat evaluates this integral.
Try the following:
$ mcx dat int 0 0:2:.2 e3 x1 x2 'iy(x1)iy(0)'
This evaluates I for a lower bound of 0 and a uniform mesh of points for the upper bound between 0 and 2. In the third column the analytic integral is evaluated from the macro iy at the upper and lower bounds.
You should that the numerical and analytic integrals agree to about 6 decimal places:
% rows 11 cols 3 real
0.000000 0.000000 0.000000
0.200000 0.015387 0.015388
0.400000 0.047801 0.047802
0.600000 0.084343 0.084343
0.800000 0.118767 0.118767
1.000000 0.148498 0.148499
1.200000 0.172890 0.172890
1.400000 0.192230 0.192230
1.600000 0.207200 0.207200
1.800000 0.218578 0.218578
2.000000 0.227105 0.227105
Interpolation
Interpolation proceeds much in the same was as integration; only interpolation has lower bound. Try
$ mcx vp=2 dat intrp .5:1:.05 e3 x1 x2 'x1^p*exp(lam*x1)' e3 x1 x2 x2x3
This returns the abscissa on a mesh twice finer than the original mesh. Every odd point is perfectly interpolated (they lie on the original mesh); the even points reflect the error of the interpolation.
3. mcx manual
mcx is a stackbased, commandline driven calculator for matrices. Matrices reside on the stack, ordered as s_{0}, s_{1}, … . There are unary operators that operate on the toplevel element s_{0}, replacing it with some transformation, and binary operators that operate on s_{1} and s_{0} replacing both of them with the result of some operation, e.g. s_{1} × s_{0}.
Usage:
mcx [−switches] datafile [ops] …
Arguments that do not begin with “−“ must be files, stored in the form of a 2D array. For ASCII files, data is read using the standard Questaal format which features programming language capabilities.
When a file is read, its contents (together with the number of rows nr and columns nc) is pushed onto the stack and becomes the toplevel stack element s_{0}. Elements already existing on the stack get pushed down one level. If there are n such elements, s_{i1} → s_{i} for i = n, n−1, …, and the new element becomes s_{0}.
Data is normally read from a file; however if datafile is a full stop (“ . ”), data is read from standard input in lieu of a file. (It can occur only once). See Example 2.2.
Any argument that begins with “−“ is a switch, a unary operator, or a binary operator.
Switches
 −nc=# (nr=#)
Stipulate that next matrix read has # columns (rows). See Example 2.2 for an illustration.  −vvar=#
define variable −vvar assign value to #. This is the standard way Questaal programs assign variables from the command line. Since data files are parsed by the preprocessor, such variables may enter into preprocessor directives or as part of expressions in the data itself.  −show
show data stack and any operations pending. See Example 2.1.  −w[l=string] fname  −bw[l=string] fname
write s_{0} to file fname. −bw writes to a binary file.  −wap
(complex arrays only) write s_{0} as (amplitude,phase) rather than (Re, Im).  −a[nr=#:nc=#] nam  −ap nam
assign s_{0} to nam, and pop s_{0} off the stack. −ap nam performs the assignment but does not pop s_{0} off the stack.  −av[ir,ic] var
assigns scalar variable var to element from s_{0}(ir,ic). If ir and ic are specified, the (1,1) element is used. 
r~switches
switches are separated by ’~’; you can use a different character instead. Switches are:~qr read with fortran read (fast, no algebra) ~s=# skips # records before reading ~open leaves file open after reading; thus if you read the file again it will read the next array ~br read from binary file, using first record to read nr and nc ~spc load in sparse format ~br,nr,nc binary file has no first record; user supplies nr,nc  −px[:nprec]  −pxc
write in row (column) compressed storage format, displaying only elements larger than a tolerance By default the tolerance is 10^{−8}. If nprec is supplied, the tolerance is By default the tolerance is 10^{nprec}.
Unary operators
These operators act on the top level matrix, s_{0}, and replace it with the result of the unop.
Where expr appears in the switches below, algebraic variables declared from directives may be used; also you can use xn for the contents of the n^{th} column.
 −p:
Pushes top array s_{0} onto stack, duplicating it.  −p+n (pn):
Pushes nth array (from bottom) onto stack  −pop:
Pops s_{0} from stack, shifting the other elements up one level.  −csum[:list]
Sums the columns of s_{0}, replacing it with a matrix of one column.  −rsum[:list]
Sums the rows of s_{0}, replacing it with a matrix of one row.  −s#:
Scale s_{0} by # (may use s#real,#imag)  −shft=#:
Adds a constant # to s_{0}  −sort expr
Sorts rows s_{0}, ordering them by the result of expr  −i  − iq:
Inverts s_{0}  −1:n
Pushes unit matrix, dimension n, onto stack  −tp [nc~]list
Generates matrix from list. Elements in list are real numbers; see here for syntax.  −evl  −evc
Replaces s_{0} by its eigenvalues (eigenvectors)  −t:
Transposes s_{0}  −cc:
Take complex conjugate of s_{0}  −herm: Replaces s_{0} with its hermitian part
 −real:
Replaces s_{0} with its real part  −v2dia:
If s_{0} is a vector, it is expanded into diagonal matrix.
If s_{0} is a square matrix, its diagonal elements are used to form a vector.  −split nam x1,x2…,xn y1,y2,…,yn:
Splits s_{0} into subblocks; assign them names namij.  −rep:n1,n2
Concatenates replicas of s_{0} to create (nr×n1,nc×n2) array.  −roll:#1[,#2]
Cyclically shifts rows (and columns, respectively) by #1 (and #2) elements.  −pwr=#:
Raises (col1) matrix to a power
Row and column Unary Operators
These unops treat s_{0} as an array of nr rows and nc columns.
 −rowl list  −coll list
Creates a new array from a list of rows (columns) of s_{0}  −rowl:mirr  −coll:mirr
Rearranges rows (columns) in reverse order  −rowl:pf=fnam  −coll:pf=fnam  −rowl:ipf=fnam  −coll:ipf=fnam
Same as −rowl (−coll) but list is read from permutation file fnam. ipf reverses the sense of the permutation.
Example: file perm contains 2 3 1 4. Then −rowl:pf=perm returns s_{0} with rows 2,3,1,4 permuted into 1,2,3,4
 −rowl:ipf=perm returns s_{0} w/ rows 3,1,2,4 permuted into 1,2,3,4
 −inc expr
Retains rows from s_{0} for which expr is nonzero.  −sub t,b,l,r  sub t,l
Extracts a subblock of s_{0}. In second form, bottom right corner = (nr,nc).  −subs:# t,b,l,r
Scales a subblock of s_{0} by #.  −subv:# t,b,l,r
Copies # to subblock of s_{0}.  −e# expr1 expr2 … expr#
Create new matrix of # columns with values expr1 expr2 … expr# .
Variables x1, x2, … can be used in the expressions. These variables refer to the elements of s_{0}in columns 1, 2, …, and change with each row.
This switch is used in Example 2.4.
Note: : at present this switch does not work for complex elements.
Unary Operators that treat data as discretized continuous functions of the first column
In these unops, column 1 (x1) is treated as an the independent variable.
Three of these unops fit data in other columns with a polynomial in x1. The polynomial is used to differentiate (−diff), integrate (−int), or interpolate (−intrp) data in the remaining columns. There are some modifiers (called [:opts] in the documentation below):
 :nord=# will change the number of points in fitting polynomial to # (polynomial order is #−1). The default value is #=4.
 :rat will use a rational polynomial, rather than an ordinary one. Good for data with singularities.
 :mesh (applies to −int only) replaces s_{0} with its integral on the same mesh as the original set of values x1. A trapezoidal rule is used at present.
You can string the options together. These three unops are demonstrated in Example 2.4.
 −diff[:opts]
differentiates columns 2…nc wrt column 1, evaluated at each x1.
s_{0} is returned as a table of points with list in the first column and the result of the derivative in the remaining columns.  −int[:opts] xlo list
integrates columns 2…nc wrt column 1 from the lower bound xlo to a set of upper bounds, given by list.
s_{0} is returned as a table of points with list in the first column and the result of the integral in the remaining columns. 
−intrp[:opts] list
interpolates column 2 to points in list.
s_{0} is returned as a table of points with list in the first column and the result of the interpolation in the remaining columns.  −smo width,xmin,xmax,dx:
smooths vector of deltafunctions with gaussians  −abs
takes the absolute value of each element  −max[:ig]
puts the largest element into the first column.
Optional g returns max value of the entire array
Optional i returns index to max value.  −unx[:i1] #1,#2
(uncross) exchanges points in columns #1 and #2 after their point of closest approach  −unx2[:i1] #1,#2
(cross) exchanges points in columns #1 and #2 at their point of closest approach, if they do not cross.  −at[:i1] val expr
find adjacent rows that bracket expr=val. Array contains linearly interpolated x1 and expr.  −nint[#1]
Replaces each element with nearest integer.
Optional #1: scale array by #1 before operation, and by 1/#1 after. Thus nint:1000 rounds to nearest thousandth.
Binary Operators
These operate on the top two matrices in the stack, s_{1} and s_{0}. s_{1} and s_{0} are popped from the stack, and the result of the binary operation is put into s_{0}.
Dimensions for s_{0} and s_{1} are not independent; for example if they are added they must have the same number of rows and columns.
 −tog
Toggle s_{0} and s_{1}.  −+
Add s_{1} + s_{0}  −−
Add s_{1} − s_{0}  −x
Multiply s_{1} × s_{0}  −xe
Multiply s_{1} and s_{0} element by element  −de
Divide s_{1} / s_{0} element by element  −x3
Multiply s_{1} s_{0} are thou they are 3D arrays: s_{1} = s_{1}(n11,n21,n31) s_{0} = s_{0}(n10,n20,n30) where n10=nr(0)/n20,n20=nc(1),n30=nc(0); n11=n10,n21=nr(1)/n11,n31=n20 Result(i,j,k) = sum_m s1(i,j,m) s0(i,m,k) is condensed into 2D (nr(1),nc(0))  −gevl  −gevc
Same as −evl  −evc, but for the generalized eigenvalue problem. s_{1} is the nonorthogonal matrix.  −orthos:
Replace s_{0} with s_{1}^{1/2} s_{0} s_{1}^{1/2}  −ccat:
Concatenate columns of s_{1} and s_{0} into a single array  −rcat:
Concatenate rows of s_{1} and s_{0} into a single array  −cross:
Cross product s_{1}(1,1..3) x s_{0}(:,1..3)  −suba[#] t,b,l,r  suba[#] t,l
Copy s_{1} to subblock of s_{0}. Conventions for subblock are the same as for −suba t,b,l,r  sub t,l.
Optional # copies #×s_{1} into s_{0}.  −index:
Use s_{0} as an row index to s_{1}. s_{0}(i) is overwritten by s_{1}(s_{0}(i)).
s_{1} is preserved. New s_{0} has row dimensions of the original s_{0} and column dimensions of s_{1}.
Other resources
The source code to mcx can be found in the bitbucket repository.
Edit This Page