In Part I , we looked at the most basic operations of the numerical workbenches GNU/Octave 2.1.34, Scilab 2.6, and Tela 1.32. This time we will talk about matrices, have a look at some of the predefined functions, learn how to write our own functions, and introduce flow control statements. The article closes with a brief discussion of the applications' input and output facilities.

Vectors help a lot if data depend on a single parameter. The different parameter values are reflected by different index values. If data depend on two parameters, vectors are a clumsy container and a more general structure, which allows for two independent indices is needed. Such a structure is called a matrix. Matrices are packed like a fresh six-pack: they are rectangular storage containers and no bottle -- oops -- element is missing.

Matrices are, for example, built from scalars as the next transcript of a GNU/Octave session demonstrates.

octave:1> # temperature rain sunshine octave:1> # degF inches hours octave:1> weather_data = [ 73.4, 0.0, 10.8; ... > 70.7, 0.0, 8.5; ... > 65.2, 1.3, 0.7; ... > 68.2, 0.2, 4.1] weather_data =

73.40000 0.00000 10.80000 70.70000 0.00000 8.50000 65.20000 1.30000 0.70000 68.20000 0.20000 4.10000

Three new ideas appear in the example. First, we have introduced some
comments to label the columns of our matrix. A comment starts with a pound
sign ```#`

'' and extends until the end of the line. Second,
the rows of a matrix are separated by semi-colons ```;`

'', and
third, if an expression stretches across two or more lines, the unfinished
lines must end with the line-continuation operator ```...`

''.

Similarly to vectors, matrices can not only be constructed from scalars, but from vectors or other matrices. If we had some variables holding the weather data of each day, like

weather_mon = [73.4, 0.0, 10.8] weather_tue = [70.7, 0.0, 8.5] weather_wed = [65.2, 1.3, 0.7] weather_thu = [68.2, 0.2, 4.1]

we would have defined `weather_data`

with

weather_data = [weather_mon; weather_tue; weather_wed; weather_thu]

or, on the other hand, if we had the data from the various instruments as

temperature = [73.4; 70.7; 65.2; 68.2] rain = [0.0; 0.0; 1.3; 0.2] sunshine = [10.8; 8.5; 0.7; 4.1]

we would have defined `weather_data`

with

weather_data = [temperature, rain, sunshine]

The fundamental rule is: *Commas separate columns, semi-colons separate
rows.*

The scalars living in matrix `m`

are accessed by applying
two indices: `m(row, column)`

, where *row* is the row-index,
and *column* is the column index. Thus, the amount of rain fallen on
Wednesday is fetched with the expression

octave:10> weather_data(3, 2) ans = 1.3000

Entries are changed by assigning to them:

octave:11> weather_data(3, 2) = 1.1 weather_data =

73.40000 0.00000 10.80000 70.70000 0.00000 8.50000 65.20000 1.10000 0.70000 68.20000 0.20000 4.10000

Now that we have defined `weather_data`

we want to work with it.
We can apply all binary operations that we have seen in last month's article
on vectors. However, for this particular example, computing

rain_forest_weather_data = weather_data + 2.1 siberian_summer_weather_data = weather_data / 3.8

does not make much sense, though the computer will not complain at all. In
the first example it would dutifully add `2.1`

to every element of
`weather_data`

, in the second it would -- obedient like a sheepdog
-- divide each element by `3.8`

.

Say we want to do something meaningful to `weather_data`

and
convert all temperatures from degrees Fahrenheit to degrees Celsius. To that
end, we need to access all elements in the first column. The vector of
interest is

octave:16> temp = [weather_data(1, 1); ... > weather_data(2, 1); ... > weather_data(3, 1); ... > weather_data(4, 1)] temp =

73.400 70.700 65.200 68.200

Obviously, the row-indices `[1, 2, 3, 4]`

form a vector themselves. We can use a shortcut and write this vector of
indices as

temp = weather_data([1, 2, 3, 4], 1)

In general, any vector may be used as index vector. Just watch out that no
index is out of range. Ordering of the indices does matter (for example ```
weather_data([2, 1, 4, 3], 1)
```

puts Tuesday's temperature in front) and
repeated indices are permitted (for example ```
weather_data([3, 3, 3, 3, 3,
3, 3], 1)
```

holds Wednesday's temperature seven times).

In our example, the index-vector can be generated by a special built-in,
the range generation operator ```:`

''. To make a vector that
starts at *low* and contains all integers from *low* to *
high*, we say

low:high

For example

octave:1> -5:2 ans =

-5 -4 -3 -2 -1 0 1 2

Our weather data example now simplifies to

temp = weather_data(1:4, 1)

Accessing a complete column or row is so common that further shortcuts
exist. If we drop both, *low* and *high* from the
colon-operator, it will generate all valid column indices for us. Therefore,
we reach at the shortest form to get all elements in the first column.

octave:17> temp = weather_data(:, 1) temp =

73.400 70.700 65.200 68.200

With our new knowledge, we extract the sunshine hours on Tuesday, Wednesday, and Thursday

octave:19> sunnyhours = weather_data(2:4, 3) sunnyhours =

8.50000 0.70000 4.10000

and Tuesday's weather record

octave:20> tue_all = weather_data(2, :) tue_all =

70.70000 0.00000 8.50000

Now it is trivial to convert the data on the rain from inches to
millimeters: Multiply the second column of `weather_data`

by 25.4
(Millimeters per Inch) to get the amount of rain in metric units:

octave:21> rain_in_mm = 25.4 * weather_data(:, 2) rain_in_mm =

0.00000 0.00000 27.94000 5.08000

We have already seen that vectors are compatible with scalars

1.25 + [0.5, 0.75, 1.0]

or

[-4.49, -4.32, 1.76] * 2

Scalars are also compatible with matrices.

octave:1> 1.25 + [ 0.5, 0.75, 1.0; ... > -0.75, 0.5, 1.25; ... > -1.0, -1.25, 0.5] ans =

1.75000 2.00000 2.25000 0.50000 1.75000 2.50000 0.25000 0.00000 1.75000

octave:2> [-4.49, -4.32, 1.76; ... > 9.17, 6.35, 3.27] * 2 ans =

-8.9800 -8.6400 3.5200 18.3400 12.7000 6.5400

In each case the result is the scalar applied to every element in the vector or matrix.

How about vectors and matrices? Obviously, an expressions like

[7, 4, 9] + [3, 2, 7, 6, 6] [2, 4; 1, 6] - [1, 1, 9, 4]

do not make any sense. In the first line the vectors disagree in size (3
vs. 5 elements), in the second line they have different shapes
(2 columns and 2 rows vs. 4 columns and 1 row). To make
sense, vectors or matrices that are used in an addition or subtraction must
have the same shape, which means the same number of rows and the same number
of columns. The technical term for ``shape'' in this context is dimension. We
can query the dimension of anything with the built-in function ```
size()
```

.

octave:22> size(weather_data) ans =

4 3

octave:23> size(sunnyhours) ans =

3 1

The answer is a vector whose first element is the number of rows, and whose second element is the number of columns of the argument.

Multiplication and division of matrices can be defined in two flavors, both of which are implemented in the numerical workbenches.

- Element by element multiplication or division of two vectors or matrices
of same dimensions: The number in the first row and first column of the first
matrix is multiplied by the number in the first row and first column of the
second matrix and so on for every element.
a = [3, 3; ... 6, 4; ... 6, 3] b = [9, 3; ... 8, 2; ... 0, 3]

octave:1> a .* b ans =

27 9 48 8 0 9

The element-by-element operators are preceded by a dot: element-by-element multiplication ``

`.*`

'' and element-by-element division ```./`

''. - Matrix multiplication as known from Linear Algebra:
*c*=*a***b*, where*a*is a*p*-times-*q*matrix and*b*is a*q*-times-*r*matrix. The result*c*is a*p*-times-*r*matrix.Example:

a = [3, 3; ... 6, 4; ... 6, 3]

b = [-4, 0, 1, -4; ... -1, -3, 2, 0]

octave:1> a * b ans =

-15 -9 9 -12 -28 -12 14 -24 -27 -9 12 -24

Although we have not seen

`for`

-loops yet (they will be discussed farther down), I would like to write the code behind the matrix multiplication operator ```*`

'' to give the reader an impression of the operations involved.for i = 1:p for j = 1:r sum = 0 for k = 1:q sum = sum + a(i, k)*b(k, j) end c(i, j) = sum end end

Compare these triply nested

`for`

-loops with the simple expression`c = a * b`

. - Matrix division? You cannot divide by a matrix! However, operator
``
`/`

'' is defined for vectors and matrices. But writing*x*=*b*/*a*, where*a*and*b*are matrices or vectors has nothing to do with division at all! It means: please solve the system of linear equationsx * a = b

for

*x*, given matrix*a*and the right-hand-side(s)*b*. Here ```*`

'' denotes matrix multiplication as defined in the previous item, and the same rules for compatible dimensions of*a*and*b*apply.a = [-2, 3, 1; ... 7, 8, 6; ... 2, 0, -1]

b = [-26, 5, -6; ... 24, 53, 26]

octave:1> x = b / a x =

7.00000 -2.00000 1.00000 7.00000 4.00000 5.00000

Isn't that an easy way to solve a system of linear equations? Imagine you had to write the code which does exactly that.

Finally, let us verify the result by multiplying with

*a*againoctave:2> x*a ans =

-26.0000 5.0000 -6.0000 24.0000 53.0000 26.0000

which, as expected, recovers

*b*.

**Details**

- For convenience GNU/Octave and Scilab define an alternative matrix
division operator ``
`\`

''.*x*=*a*\*b*solves the linear system of equationsa * x = b

for

*x*, given matrix*a*and the right-hand-side(s)*b*. This is the form most users prefer, because here*x*is a column vector, whereas operator ```/`

'' returns*x*as row-vector. - operator ``
`\`

'' has the dotted cousin ```.\`

'' and the relation*a*./*b*==*b*.\*a*holds.

**Differences**

- Scilab and Tela use C++-like comments
// This is a Scilab or a Tela comment

- Tela does not need or understand the line continuation operator
``
`...`

''weather_data = #(73.4, 0.0, 10.8; 70.7, 0.0, 8.5; 65.2, 1.3, 0.7; 68.2, 0.2, 4.1)

In interactive mode, Tela does not handle multi-line expressions as the above. Multi-line expressions must be read from a file (with

`source("filename.t")`

). - In Tela the operators ``
`*`

'' and ```/`

'' work element by element, this is, they work like ```.*`

'' and ```./`

'' do in GNU/Octave and Scilab. Matrix multiplication (*a***b*in GNU/Octave or Scilab) is written asa ** b

or

matmul(a, b)

solving systems of linear equations (

*b*/*a*in Octave or Scilab) aslinsolve(a, b)

Ugh -- far too many to mention! The workbenches supply dozens of predefined functions. Here I can only wet the reader's appetite.

**Generating Special Matrices**

- Several matrices occur often enough in computations that they have been
given their own generating functions. These are for example,
*m*-times-*n*matrices filled with zeros:`zeros(m, n)`

or ones:`ones(m, n)`

, or*n*-times-*n*diagonal matrices, where the diagonal consists entirely of ones:`eye(n)`

or the diagonal is set to numbers supplied in a vector:`diag([a1, a2, ..., I<an>])`

. **Analyzing Matrices**

- Getting the smallest or largest element in matrix
*a*:`min(a)`

,`max(a)`

, or totaling matrix*a*:`sum(a)`

.**Differences:**GNU/Octave's`min(a)`

,`max(a)`

, and`sum(a)`

return the column-wise result as a row vector. To get the minimum, maximum, and sum of all elements in matrix*a*, use`min(min(a))`

,`max(max(a))`

,`sum(sum(a))`

. **Linear Algebra**

- We mentioned that systems of linear equations, like
*x***a*=*b*, are solved for*x*with the slash operator ```/`

''. But many more linear algebra functions exist, for example singular value decomposition:`svd(a)`

, or eigenvalue computation:`eig(a)`

.**Differences:**In Tela uses`SVD(a)`

instead of`svd(a)`

, and instead of`eig(a)`

, Scilab uses`spec(a)`

to compute the eigenvalue spectrum.

One note on performance: basically, all three applications are interpreters. This means that each expression is first parsed, then the interpreter performs desired computations, finally calling the functions inside of the expressions -- all in all a relatively slow process in comparison to a compiled program. However, functions like those shown above are used in their compiled form! They execute almost at top speed. What the interpreter does in these cases is to hand over the complete matrix to a compiled Fortran, C, or C++ function, let it do all the work, and then pick up the result.

Thus we deduce one of the fundamental rules for successful work with numerical workbenches: prefer compiled functions over interpreted code. It makes a tremendous difference in execution speed.

No matter how many functions a program may provide its users, they are never enough. Users always need specialized functions to deal with their problems, or they simply want to group repeated, yet predefined operations. In other words, there always is a need for user-defined functions.

User functions are best defined in files, so that they can be used again in
later sessions. For GNU/Octave, functions files end in *.m*, and are
loaded either automagically or with
`source("`

. Scilab calls its
function files *filename.m*")*.sci*, and requires them to be loaded with ```
getf("
```

. Tela functions are stored
in *filename.sci*")*.t*-files and loaded with ```
source("
```

. As big as the differences
are in loading functions, all workbenches use quite similar syntax for the
definition of functions.*filename.t*")

GNU/Octave and Scilab

function [res1, res2, ..., resM] = foo(arg1, arg2, ..., argN) # function body endfunction

Tela

function [res1, res2, ..., resM] = foo(arg1, arg2, ..., argN) { // function body };

where *arg1* to *argN* are the functions' arguments (also
known as parameters), and *res1* to *resN* are the return
values. Yes, trust your eyes, multiple return values are permitted, what might
come as a surprise to most readers who are acquainted with popular programming
languages. However, this is a necessity, as no function is allowed to change
any of its input arguments.

Enough theory! let us write a function that takes a matrix as input and returns a matrix of the same dimensions, with the entries rescaled to lie in the interval (0, 1).

### Octave

function y = normalize(x) ## Return matrix X rescaled to the interval (0, 1).

minval = min(min(x)) maxval = max(max(x))

y = (x - minval) / (maxval - minval) endfunction

Now define a Scilab function that returns the spectral radius on a matrix.
We use `abs()`

which returns the magnitude of its (possibly
complex) argument.

// Scilab

function r = spectral_radius(m) // Return the spectral radius R of matrix M.

r = max(abs(spec(m))) endfunction

Finally, we write a Tela function which computes the Frobenius norm of a matrix.

// Tela

function x = frobenius(m) // Return the Frobenius norm X of matrix M. { x = sqrt(sum(abs(m)^2)) };

**Details:**

GNU/Octave's ``automagical'' function file loading works the following way:
if Octave runs into an undefined function name it searches the list of
directories specified by the built-in variable `LOADPATH`

for files
ending in .m that have the same base name as the undefined function; for
example, `x = my_square_root(2.0)`

looks for the file *
my_square_root.m* in the directories listed in `LOADPATH`

.

All code we have written thus far executes strictly top-to-bottom, we have not used any flow control statements such as conditionals or loops.

Before we manipulate the flow of control, we should look at logical expressions because the conditions used in conditionals and loops depend on them. Logical expressions are formed from (1.) numbers, (2.) comparisons, and (3.) logical expressions catenated with logical operators.

- Zero means logically false, any number not equal to zero means logically true, hence C-programmers should feel at home.
- The usual gang of comparison operators exist: less-than
``
`<`

'', less-or-equal ```<=`

'', greater-than ```>`

'', greater-or-equal ```>=`

'', and equal ```==`

''.**Differences:**The inequality operator varies quite a bit among the programs. (Octave cannot decide whether it feels like C, Smalltalk, or Pascal. Scilab wants to be Smalltalk and Pascal at the same time. :-)!= ~= <> # Octave ~= <> // Scilab != // Tela

- Complex logical expressions are formed with logical operators ``and'',
``or'' and ``not'' whose syntax is borrowed from C. However, each program uses
its own set of operators. Thus, we have to list some
**Differences:**and or not ---- ---- ---- & | ! ~ # Octave & | ~ // Scilab && || ! // Tela

We are all set now for the first conditional, the ```
if
```

-statement. Note that the parenthesis around the conditions are
mandatory (as they are in C). The `else`

-branches are optional in
any case.

# Octave // Scilab // Tela

if (cond) if cond then if (cond) { # then-body // then-body // then-body else else } else { # else-body // else-body // else-body endif end };

*cond* is a logical expression as described above.

`while`

-statements:

# Octave // Scilab // Tela

while (cond) while cond while (cond) { # body // body // body endwhile end };

Again, *cond* is a logical expression.

`for`

-statements in Octave and Scilab walk through the columns
of *expr* one by one. Most often *expr* will be a vector
generated with the range operator ```:`

'', like ```
for i =
1:10
```

. Tela's `for`

-statement is the same as C's.

# Octave // Scilab // Tela

for var = expr for var = expr for (init; cond; step) { # body // body // body endfor end };

Here come some examples which only show things we have discussed so far.

Octave

function n = catch22(x0) ## The famous catch-22 function: it is ## impossible to compute that it will ## stop for a specific input. Returns ## the number of loops.

n = 0 x = x0 while (x != 1) if (x - floor(x/2)*2 == 0) x = x / 2 else x = 3*x + 1 endif n = n + 1 endwhile endfunction

Scilab

function m = vandermonde(v) // Return the Vandermonde matrix M based on // vector V.

[rows, cols] = size(v) m = [] // empty matrix if rows < cols then for i = 0 : (cols-1) m = [m; v^i] end else for i = 0 : (rows-1) m = [m, v^i] end end endfunction

Tela

function vp = sieve(n) // Sieve of Erathostenes; returns vector of // all primes VP that are strictly less than // 2*N. 1 is not considered to be a prime // number in sieve(). { vp = #(); // empty vector if (n <= 2) { return };

vp = #(2); flags = ones(1, n + 1); for (i = 0; i <= n - 2; i = i + 1) { if (flags[i + 1]) { p = i + i + 3; vp = #(vp, p); for (j = p + i; j <= n; j = j + p) { flags[j + 1] = 0 } } } };

We have been using with the workbenches a lot. At some point we would like to call it a day, but we do not want to lose all of our work. Our functions are already stored in files. It is time to see how to make our data persist.

All three applications at least have one input/output (I/O) model that borrows heavily from the C programming language. This model allows close control of the items read or written. Often though, it is unnecessary to take direct control over the file format written. If variables must be saved just to be restored later, simplified I/O commands will do.

- Octave offers the most flexible solution with the
`save`

/`load`

command pair.save filename varname1 varname2 ... varnameN

saves the variables named

*varname1*,*varname2*, ...,*varnameN*in file*filename*. The complementaryload filename varname1 varname2 ... varnameN

command restores them from

*filename*. If`load`

is given no variable names, all variables form*filename*are loaded. Handing over names to`load`

selects only the named variables for loading.Note that the

`save`

and`load`

commands do not have parenthesis and their arguments are separated by spaces not commas. Filename and variable names are strings.save "model.oct-data" "prantl" "reynolds" "grashoff" load "model.oct-data" "reynolds"

By default

`load`

does not overwrite existing variables, but complain with an error if the user tries to do so. When it is save to discard of the values of existing variables, add option ```-force`

'' to`load`

, likeload -force "model.oct-data" "reynolds"

and variable

`reynolds`

will be loaded from file*model.oct-data*no matter whether it has existed before or not. - Scilab's simple I/O parallels that of GNU/Octave:
save(filename, var1, var2, ..., varN)

However, the variables

*var1*, ...,*varN*are not strings, but appear literally. This means that the name of a variable is not stored in the file. The association between name and contents is lost!The complementary function

load(filename, varname1, varname2, ..., varnameN)

restores the contents of

*filename*in the variables named*varname1*,*varname2*, ...*varnameN*. - Tela lets the users save her variables with the
save(filename, varname1, varname2, ..., varnameN)

function, preserving the association between variable name and variable contents. The complementary

load(filename)

function loads all variables stored in

*filename*. It is not possible to select specific variables.

As we use matrices so often, specialized functions exist to load and save whole matrices. Especially loading a matrix with a single command is convenient and efficient to read data from experiments or other programs.

Let us assume, we have the ASCII file *datafile.ascii* which
contains the lines

# run 271 # 2000-4-27 # # P/bar T/K R/Ohm # ====== ====== ====== 19.6 0.118352 0.893906e4 15.9846 0.1 0.253311e5 39.66 0.378377 0.678877e4 13.6 0.752707 0.00622945e4 12.4877 0.126462 0.61755e5

and sits in the current working directory. The file's five leading lines
are non-numeric. They are skipped by the workbenches, but possibly aid the
user in identifying her data. I have intentionally taken a data set which is
not neatly formatted, as are most data files. Matrix-loading functions split
the input at whitespace not at a specific column, thus they are happy with
*datafile.ascii*.

We load the data into GNU/Octave with

octave:1> data = load("datafile.ascii") data =

1.9600e+01 1.1835e-01 8.9391e+03 1.5985e+01 1.0000e-01 2.5331e+04 3.9660e+01 3.7838e-01 6.7888e+03 1.3600e+01 7.5271e-01 6.2294e+01 1.2488e+01 1.2646e-01 6.1755e+04

or into Scilab

-->data = fscanfMat("datafile.ascii") data =

! 19.6 0.118352 8939.06 ! ! 15.9846 0.1 25331.1 ! ! 39.66 0.378377 6788.77 ! ! 13.6 0.752707 62.2945 ! ! 12.4877 0.126462 61755. !

or into Tela

>data1 = import1("datafile.ascii") >data1 #( 19.6, 0.118352, 8939.06; 15.9846, 0.1, 25331.1; 39.66, 0.378377, 6788.77; 13.6, 0.752707, 62.2945; 12.4877, 0.126462, 61755)

In all three examples data will contain a 5-times-3 matrix with all the
values from *datafile.ascii*.

The complementary commands for saving a single matrix in ASCII format are

save("data.ascii", "data") # GNU/Octave fprintfMat("data.ascii", data, "%12.6g") // Scilab export_ASCII("data.ascii", data) // Tela

Note that Scilab's `fprintfMat()`

requires a third parameter
that defines the output format with a C-style template string.

Of course none of the above save commands writes the original header, the
lines starting with hash-symbols, of *datafile.ascii*. To write these,
we need the ``low-level'', C-like input/output functions, which featured in
each of the three workbenches.

For a precise control of the input and the output, C-like I/O models are offered. All three applications implement function

printf(format, ...)

Moreover, GNU/Octave and Tela follow the C naming scheme with their C-style file I/O:

handle = fopen(filename) fprintf(handle, format, ...) fclose(handle)

whereas Scilab prefixes these functions with an ```m`

'' instead
of an ```f`

''

handle = mopen(filename) mprintf(handle, format, ...) mclose(handle)

Whether the function is called `fprintf()`

or ```
mprintf()
```

, they work the same way.

*Next Month: Graphics, function plotting and data plotting.*

Copying license http://www.linuxgazette.net/copying.html

Published in Issue 70 of