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.

This has one very important consequence for types like

ARRAY<VECTOR> my_data;

After you create such an object, all outer
ARRAY elements will point to the same data! You have to initialize
manually, by doing the following:

for (INDEX i=0; i<my_data.size(); ++i)
   resize( my_data[i], n );

where n is the size of the inner vectors that you want. Sorry for the
inconvenience, but I see no way to fix this. In practise, this is no
problem, since you anyway have to adjust the size of the inner vectors
at some point (empty vectors are not useful). I made sure that the
resize function above does the right thing (see also about resizing
below). NEVER USE a.resize(i), ALWAYS USE resize(a,i).


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.

(PE: This didn't work for me. I did instead the following

2. Go to the mtl subdirectory

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

4. In the top directory, type ./configure

5. Become root, type make install.

i.e. I changed instead mtl_config.h.in. )




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

VECTORS are MTL based, ARRAYs are STL based. The most important
difference is that VECTOR assignments are shallow, whereas ARRAY
assignments are deep. Actually, STL and MTL based ARRAYs are rather
similar. Maybe we'll switch back to MTL based ARRAYs in the future,
but at the moment I see no advantage either way. 

The two types are meant to behave very similarly, except that ARRAY
does not offer special math and matrix/vector functions. If you find
inconsistencies, let me know.


1. Define a vector:   

VECTOR x(10);

you can initialize at the same time:

VECTOR x( 10, 0.0 );

1.a. Resize a vector or ARRAY:

resize(x,20);

This is the safe way to resize. Calling x.resize() is depreciated! I
did this function myself, internally it does a `x = VECTOR(20)'.

2. Element access:

x[i] = i;

3. Print vector to output stream:

print_vector(cout,x);

I have also overloaded the << operator for VECTOR (and
MATRIX). However, with MTL you have not just VECTOR, 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. So, use print_vector if << does not work.

4. Vector algorithms:

setto(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
add(x,y,z);           // z = x+y
add(x,y,z,w);         // w = x+y+z
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:

VECTOR::subrange_type     // the type of a subrange

8. Vector methods:

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

Note that a(2,3) will give you exactly one element, the one with index
2. Similarly a(2,4) will give you the element 2 and 3, and so on.



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

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


MATRIX:
=======
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:   

MATRIX A(10,20);

initialization at the same time is not possible. Use the setto() algorithm.

1.a. Resize a matrix:

resize(A,15,20);

I implemented this function myself, to make MATRIX more consistent
with ARRAY and VECTOR. Resizing will destroy the contents of the
matrix! 

This is the safe way to resize. Internally it does a `A = MATRIX(15,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). On the other hand, if you write template
functions that should work also with the more exotic matrix-like types
(for example trans(A) or A.sub_range()) then you must use this kind of
indexing since [][] will lead to errors. 

Recommended access:

A[i][j] = 1.0

This means the same thing in TNT, so it is safe. 

3. Print matrix to output stream:

print_all_matrix(cout,A);

In most cases you can also use the << operator. 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 VECTOR directly. The correct type for it is MATRIX::OneD. However,
you can use it in all vector algorithms (like print_vector above). So
you can copy it to a vector like this:

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

5. Matrix algorithms:

setto(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:

MATRIX::shape              // tag to describe shape (e.g., rectangular)
MATRIX::orientation        // in our case always row_tag
MATRIX::sparsity           // dense_tag or sparse_tag
MATRIX::OneD               // type of inner container (the individual
                           rows, columns, or diagonals)
MATRIX::submatrix_type     // for submatrices
MATRIX::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. Indices behave in the
			  // same way as for VECTOR subranges!
A.nrows()                 // number of rows
A.ncols()                 // number of columns
A.nnz()                   // number of nonzero elements


SPARSE_MATRIX:
===========

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

1. Define:

SPARSE_MATRIX 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()




// Special functions added for ARTS:

For VECTOR and MATRIX:

transf(x,sqrt,y) or y = transf(x,exp):

The generic transformation function. You can use this to apply any
mathematical function that takes only a single argument to all
elements of a VECTOR or MATRIX x. This replaces the old exp, sqrt,
etc. functions.

interp_lin_vector() and interp_lin_matrix():

These are meant to replace the various linear interpolation
functions. They are very flexible, because they are template
functions. For example, you can use interp_lin_matrix to interpolate
either the rows or the columns of the matrix, by using the trans()
adapter on the input and output matrix argument. (This is used inside
ARTS in one place, just do a grep for the function.)



// Common bugs:

1. Assigning a scalar to a vector as in:
   VECTOR a;
   a = 1;
   At runtime this will lead to a segmentation fault later on in the
   program. Unfortunately this bug is hard to locate. There could be
   still a few instances of it creeping through ARTS.

2. Calling a MTL algorithm without making sure that the arguments have
   compatible dimensions. There is no safeguard against this in MTL!

3. Matrix indices not properly adapted to 0-based.