Basic types:
============

1. The floating point type: Numeric
2. The integer type:        INDEX

Always use these two types for numbers, unless you have a very good
reason to use something else.

The program should run with Numeric beeing either float or
double. With float it should be faster. If you change the INDEX type
you have to be sure to be consistent with MTL. At the moment it is
size_t.


General MTL information:
========================

For detailed information see the thesis mentioned on
http://www.lsc.nd.edu/research/mtl/publications.php3 (the document is
http://www.lsc.nd.edu/downloads/research/mtl/papers/thesis.ps.gz). Is
is easy to read and explains the MTL concept much better than is
possible here.

Assignments of vectors and matrices are shallow! That means if you
really want to make a duplicate of the contents, you have to use the
copy() algorithm. For copy, the target object (matrix or vector) must
have the right dimensions. Dimensions are not adjusted automatically,
you will just get a core dump.

MTL installation (for MTL-2.1.2-19):
====================================

1. Get MTL source

2. In the top directory, type ./configure

3. Go to the mtl subdirectory

4. In file mtl_config.h, comment out all lines which define PACKAGE
   and VERSION (4 lines altogether, because the #defines occur twice).

5. Go back to top directory, become root, type make install.


VEC: (description also valid for ARRAY!)
=======

There is no technical need to distinguish between ARRAY and VEC. At
the moment the two are separate, but one could simply use
ARRAY<Numeric> in place of VEC. Unfortunately, this nice symmetry
brakes down for matrices. Leaving the type open (using MAT<Numeric>
instad of MAT) does not work in this case, due to the design of MTL
(if I understand correcly).

1. Define a vector:   

VEC x(10);

1.a. Resize a vector:

There is a resize function:

x.resize(20);

This is at the moment heavily used in the program. However, I
discovered that sometimes resize is dangerous, in particular if you
resize a vector that is part of another vector. The better way to
resize is by assigning a vector of the desired size:

x = VEC(20) 


2. Element access:

x[i] = i;

3. Print vector to output stream:

print_vector(cout,x);

One could easily overload the << operator for VEC. In fact I tried
this, but discovered that then there is a loss of generality. With MTL
you have not just VEC, but also some derived types (for example for
subranges). Function print_vector (and print_all_matrix below) will
work for all these types because it uses genericity. It is not
possible to program a general solution for the << operator and having
separate ones for each subtype would be a mess.


4. Vector algorithms:

set(x,alpha);         // set all elements to alpha
scale(x,alpha);       // multiply all elements by alpha
s = sum(x);           // sum up all elements
s = one_norm(x);      // sum up absolute values
s = two_norm(x);      // sum up squares
s = infinity_norm(x); // return maximum absolute value
s = max(x);           // return maximum
s = min(x);           // return minimum
print_vector(os,x);   // print to stream os

5. Vector Vector algorithms:

The target vector must have the right size! If not you will get a core
dump. 

copy(x,y);            // copy contents of x to y
swap(x,y);            // exchange contents of x and y
ele_mult(x,y,z);      // elementwise multiplication of x and y
ele_div(x,y,z);       // elementwise division y / x
add(x,y);             // add x to y
s = dot(x,y);         // dot (scalar) product

6. Adaptor helper functions:

We use the add() algorithm to explain these:

add(scaled(x,alpha),y); // multiply x by alpha before adding to y
add(strided(x,alpha),y); // use steps of alpha instead of 1 inside x

These should be efficient, because they just present a different
`view' of the vector. For example, with scaled and add, the scaling is
done within the add algorithm, not before.

7. Vector associated types:

VEC::subrange_type     // the type of a subrange

8. Vector methods:

a(s,f)                 // access subrange s to f



ARRAY<type>:
===========

Is just like VEC, just that the type can be anything.


MAT:
=======
Conceptually a `container of containers'.
We use row major order. The outer container contains the rows of the
matrix. You use a `twoDiterator' to select the row and a
plain `iterator' to select a column. Put differently, each row of a
matrix behaves just like a VECTOR.

1. Define a matrix:   

MAT A(10,20);

2. Element access:

You could use:

A(i,j)  = 1.0;

THIS IS DANGEROUS, because in TNT the same notation meant a different
thing (1-based indices).

Recommended access:

A[i][j] = 1.0

This means the same thing in TNT, so it is safe. ARTS code should be
changed to this notation BEFORE switching to MTL.


3. Print matrix to output stream:

print_all_matrix(cout,A);

See also what I wrote for print_vector.


4. Select a row of a matrix:

print_vector(cout, A[3]);

Unfortunately, although A[3] behaves like a vector, you cannot assign
it to a VEC directly. The correct type for it is MAT::OneD. However,
you can use it in all vector algorithms (like print_vector above). So
you can copy it to a vector like this:

VEC x(A.ncols());
copy(A[3],x);

5. Matrix algorithms:

set(A,alpha);          	// set all elements of A to alpha
scale(A,alpha);        	// multiply all elements of A by alpha
set_diagonal(A,alpha); 	// set diagonal elements to alpha
transpose(A);          	// replace A by A transposed
print_all_matrix(os,A); // print to stream os

6. Matrix vector algorithms:

mult(A,x,y);            // y = A*x
mult(A,x,y,z);          // z = A*x + y
tri_solve(T,x);         // x = T^-1 * x

7. Matrix matrix algorithms:

copy(A,B);              // copy A to B
swap(A,B);              // exchange A and B
add(A,C);		// C = A+C
ele_mult(A,B,C);        // elementwise multiplication of A and B, stored in C
mult(A,B,C);            // C = A*B (matrix multiplication)
mult(A,B,C,E);          // E = A*B + C
tri_solve(?);           // Calculate inverse

8. Adaptor helper functions:

scaled(A,alpha)         // scale by factor alpha
trans(A)                // transpose
rows(A)                 // Access to rows for columns matrix
columns(A)              // Access to columns for row matrix (ARTS case)

9. Matrix associated types:

MAT::shape              // tag to describe shape (e.g., rectangular)
MAT::orientation        // in our case always row_tag
MAT::sparsity           // dense_tag or sparse_tag
MAT::OneD               // type of inner container (the individual
                           rows, columns, or diagonals)
MAT::submatrix_type     // for submatrices
MAT::value_type         // in our case Numeric

OneD and submatrix_type are very useful for row, column, or submatrix
access (see usage examples in demonstrate_linalg.cc). Sparsity should
be useful for the maybe_sparse types.

10. Matrix methods:

A.sub_matrix(r0,r1,c0,c1) // Selects submatrix from row r0 to r1 and
                             columns c0 to c1
A.nrows()                 // number of rows
A.ncols()                 // number of columns
A.nnz()                   // number of nonzero elements


SPARSE_MAT:
===========

You can use basically all algorithms and member functions of MAT.

1. Define:

SPARSE_MAT S(n,m);        // n,m are the dimensions, the number of
                             nonzero elements is 0 at first

2. Set elements by:

S(3,4) = 1;

3. Number of nonzero elements:

S.nnz()