%
% To start the document, use
%  \chapter{...}
% For lover level, sections use
%  \section{...}
%  \subsection{...}
%
\chapter{Vectors, matrices, tensors, and arrays}
%-------------------------------------------------------------------------
\label{sec:matpack}


%
% Document history, format:
%  \starthistory
%    date1 & text .... \\
%    date2 & text .... \\
%    ....
%  \stophistory
%
\starthistory
  030807 & Sparse added by Mattias Ekstr\"om. \\
  030109 & Documentation for using jokers without Range added by Stefan Buehler.\\
  020516 & Tensors added by Stefan Buehler.\\
  011018 & Created and written by Stefan Buehler.\\
\stophistory




%
% Introduction
%

This section describes how vectors and matrices are implemented in
ARTS and how they are used. Furthermore it describes how arrays of
arbitrary type can be constructed and used.


\section{Implementation files}
%-------------------------------------------------------------------------
\label{sec:matpack:files}

The \builtindoc{Matrix} and \builtindoc{Vector} classes described below reside in the files:
\begin{itemize}
\item \fileindex{matpackI.h}
\item \fileindex{make\_vector.h}
\item \fileindex{matpackI.cc}
\item \fileindex{make\_vector.cc}
\end{itemize}

Tensors of order 3 to 7 reside in the files:
\begin{itemize}
\item \fileindex{matpackIII.h}
\item \fileindex{matpackIV.h}
\item \fileindex{matpackV.h}
\item \fileindex{matpackVI.h}
\item \fileindex{matpackVII.h}
\end{itemize}

The template class \shortcode{Array} (also described below) is implemented 
in the files:
\begin{itemize}
\item \fileindex{array.h}
\item \fileindex{make\_array.h}
\end{itemize}

The \builtindoc{Sparse} class is described in the file:
\begin{itemize}
\item \fileindex{matpackII.h}
\end{itemize}

The file \fileindex{test\_matpack.cc} contains test cases and usage
examples. For Sparse there is a separate test file, 
\fileindex{test\_sparse.cc}.

\section{Vectors}
%-------------------------------------------------------------------------
\label{sec:matpack:vectors}

The class \typeindex{Vector} implements the mathematical concept of a
vector. (Surprise, surprise.) This means that:
\begin{itemize}
\item A Vector contains a list of floating point values of type \typeindex{Numeric}.
\item A Vector can be multiplied with another Vector (scalar product),
  or with a Matrix.
\item Sub-ranges of a Vector can easily be accessed, and used as if
  they were Vectors.
\item Resizing a Vector is expensive and should be avoided.
\end{itemize}

\subsection{Constructing a Vector}
%-------------------------------------------------------------------------
You can construct an object of class Vector in any of these ways:

\begin{code}
Vector a;         // Create empty Vector.
Vector b(3);      // Create Vector of length 3, if
                  // created like this it will contain
                  // arbitrary values.
Vector c(3,0.0);  // Create Vector of length 3, and
                  // fill it with 0.

Vector d=c;       // Make d a copy of c.

Vector e(1,5,1);  // 1, 2, 3, 4, 5
Vector f(1,5,.5); // 1, 1.5, 2, 2.5, 3
Vector g(5,5,-1); // 5, 4, 3, 2, 1
\end{code}

The last three examples all use the same constructor, which takes
the three arguments `start', `extent', and `stride'. It will create a
Vector containing `extent' elements, starting with `start', with a
step of `stride'.

There also exists a special sub-class of Vector that can be initialized
explicitly. This must be a special class in order to avoid ambiguities
with the standard constructors. Usage:

\begin{code}
MakeVector a(1.0,2.0,3.0);  // Creates a vector of length 3 
                            // containing the values 
                            // 1.0, 2.0, and 3.0.
\end{code}

You can use MakeVectors just like Vectors, except that the
constructors are different. Otherwise you can mix them freely with
Vectors. 

\subsection{VectorViews}
%-------------------------------------------------------------------------
\label{sec:vector_views}

An object of class \typeindex{VectorView} is, like the name says, just
another view on an existing Vector. It does not have its own
data. This has the important consequence that it cannot be resized,
since that would mess up the original Vector that the view is
referring to. You can create VectorViews from Vectors using the index
operator `[]', the class \typeindex{Range}, and the special \shortcode{joker}
object. Examples:

\begin{code}
MakeVector x(1,2,3,4,5,6,7); 
VectorView a = x;                 // Now a refers to the 
                                  // whole of x;
VectorView b = x[Range(joker)];   // Same effect.
VectorView c = x[Range(0,2)];     // Take 2 elements of x, 
                                  // starting at the 
                                  // beginning,
                                  // in this case: 1,2.
VectorView d = x[Range(0,3,2)];   // In this case: 1,3,5.
VectorView e = x[Range(3,joker)]; // In this case: 4,5,6,7.
\end{code}

As you can see, most useful ways to create VectorViews involve the
Range class. The general constructor to this class takes three
arguments, `start', `extent', and `stride'. This means that you will
select `extent' elements from the Vector, starting with index `start',
with a step-width of stride. Note that indices are 0-based, so 0
refers to the first element. The last argument, `stride', can be
omitted, in that case the default of 1 is assumed. As a special case,
`extent'==\shortcode{joker} means `to the end', and calling Range with only one
argument \shortcode{joker} means `all elements'.

Usually, you will not have to use VectorView explicitly, because you
can use expressions like:

\begin{code}
Vector a(1,5,1);                // a = 1,2,3,4,5
Vector b = a[Range(1,3)];       // b = 2,3,4
\end{code}

However, \shortcode{VectorView} and the related class
\typeindex{ConstVectorView} are extremely useful as the argument types of
functions operating on Vectors. You should define your functions like
this:

\begin{code}
void silly_function(VectorView a,      // Output argument
                    ConstVectorView b  // Input argument
                                       // (read only)
                   )
{
   // Do some silly stuff with a and b.
}
\end{code}

Note that there must not be any `\&' after VectorView or
ConstVectorView. In other words they have to be passed by value, not
by reference. This is ok, since they do not contain the actual
data, so that passing by value is efficient. Passing VectorViews by
reference is forbidden.

You should use these kind of arguments for all input Vectors, and also
for the output if you have a function that does not resize the output
Vector. This has the great advantage that you can call the function
with Vector sub-ranges, e.g., 
\begin{code}
Vector a(1,5,1);                  // a = 1,2,3,4,5
Vector b(3);                      // Set size of b.
silly_function(b,a[Range(0,3)]);  // Call fuction with
                                  // sub-range of a.
\end{code}

An exception to this rule are workspace methods, which use
conventional argument types \shortcode{const Vector\&} for input and
\shortcode{Vector\&} for output.

\subsection{What you can do with a Vector (or VectorView)}
%-------------------------------------------------------------------------

All examples below (except for the first) assume that \shortcode{a} is a
Vector, MakeVector, or VectorView.

\subsubsection{Resize (only for Vector, not for VectorView!):}
\begin{code}
a.resize(5);
\end{code}
This makes \shortcode{a} a 5 element vector. The new Vector is not
initialized (i.e., the contents will be unpredictable). Also, note
that the previous content will be completely lost. Appending to a
Vector is not possible.

\subsubsection{Get the number of elements:}
\begin{code}
cout << a.nelem();
\end{code}

\subsubsection{Sum up all elements:}
\begin{code}
cout << a.sum();
\end{code}

\subsubsection{Element access:}
\begin{code}
cout << a[3];   // Print 4th element.
a[0] = 3.5;     // Assign 3.5 to first element.
\end{code}

Note that we use 0-based indexing! Furthermore note that the operator
`[]' can be also used with \shortcode{Range}, as explained above.

\subsubsection{Copying Vectors:}
\begin{code}
Vector b;
b = a;
\end{code}
In this case the size of b will be adjusted to that of a
automatically. Maybe you have noticed that there is a way to formulate the 
example above in even shorter fashion:
\begin{code}
Vector b = a;
\end{code}
The result is exactly the same. Note, though, that in this case b is
\emph{constructed} from a, not copied (see section about constructing
Vectors above). 

\subsubsection{Copying in connection with views:}
This one is a bit tricky. Obviously, the size of views can not be
adjusted, because a view is just some selection of the underlying
object. The `=' operator in this case copies the
\emph{contents}, so the sizes of the left-hand and right-hand argument
must match. VectorView internally uses assertions to make sure of
this. So, if you get an assertion failure one reason could be that you
forgot to make the target the correct size. Here is an example:
\begin{code}
b[Range(5,5,-1)] = a[Range(3,5)];  // Copy 5 elements from 
                                   // a to b, reversing 
                                   // the order and starting
                                   // with index 3 in a.
\end{code}
Great, isn't it?

\subsubsection{Assigning a scalar:}
\begin{code}
a = 1.0;                        // Assign 1 to all elements.
\end{code}

\subsubsection{Mathematical operators:}
\begin{code}
Vector a(1,3,1), b(3,1); // a = 1,2,3; b = 1,1,1
a *= 2;                  // a = 2,4,6
                         // Similarly, /=, +=, -=
a += b;                  // a = 3,5,7
                         // Similarly, -=, *=, /=
a += a;                  // a = 6,10,14
                         // So a can appear on both sides.
\end{code}

All these operate element-wise.  Note, that there are no return
versions of these operators (i.e., expressions like \shortcode{b = a+1} are
not possible). This is again for efficiency reasons. It is currently
an active area of research in programming techniques how to make this
kind of expression efficient. None of the available solutions works,
so ARTS has to live without it.

\subsubsection{Maximum, minimum and mean:}
\begin{code}
cout << max(a);
cout << min(a);
cout << mean(a);
\end{code}

\subsubsection{Scalar product:}
\begin{code}
cout << a*a;
\end{code}

This is an exception to the rule not to have return versions of
operators. The reason is quite obvious: The return value is only a
scalar. 

\subsubsection{Arbitrary single-argument math functions:}
\begin{code}
Vector b(a.nelem());
transform(b,sin,a);  // b = sin(a)
transform(b,cos,b);  // b = sin(b)
                     // So b can appear on both sides.
\end{code}

The transform function operates on each element of \shortcode{a} with the
function you specify and puts the result in \shortcode{b}. Note that the
order of the arguments is swapped compared to the old function
\shortcode{trans} that we had in the pre-Matpack era.


\section{Matrices}
%-------------------------------------------------------------------------
\label{sec:matpack:matrices}

The class \typeindex{Matrix} implements the mathematical concept of a
matrix. (Who would have guessed this?) This means that:
\begin{itemize}
\item A Matrix contains floating point values of type \builtindoc{Numeric}.
\item The values are arranged in rows and columns and can be accessed
  by indices. The first index is the row, the second the column. In
  other words, we use \emph{row-major} order, similar to C, Matlab,
  and most math textbooks. Note, however, that some languages like
  FORTRAN and IDL use \emph{column-major} order.
\item A Matrix can be multiplied with a Vector, or with another
  Matrix.
\item A sub-range of a Matrix in both dimensions (submatrix) can
  easily be accessed, and used as if it was just a normal matrix.
\item Resizing a Matrix is expensive and should be avoided.
\end{itemize}

\subsection{Constructing a Matrix}
%-------------------------------------------------------------------------
You can construct an object of class Matrix in any of these ways:

\begin{code}
Matrix a;          // Create empty Matrix.
Matrix b(3,4);     // Create Matrix with 3 rows
                   // and 4 columns. When
                   // created like this it will contain
                   // arbitrary values.
Matrix c(3,4,0.0); // Similar, but
                   // fill it with 0.

Matrix d=c;        // Make d a copy of c.
\end{code}

That is all. More fancy constructors, like for Vector, do not exist
for Matrix. There is also no equivalent to the \shortcode{MakeVector}
class.

\subsection{MatrixViews}
%-------------------------------------------------------------------------

A \typeindex{MatrixView} is a view on an existing Matrix, in the same way
as a \shortcode{VectorView} is a view on an existing Vector. Like a
VectorView, a MatrixView cannot be resized and does not contain the
actual data. A view is generated by using Ranges:

\begin{code}
Matrix x(10,20);                  // Create 10x20 matrix.
MatrixView a = x;                 // Now a refers to the
                                  // whole of x;
MatrixView b = x(Range(joker),Range(joker));
                                  // Same effect.
MatrixView c = x[Range(0,2),Range(0,2)];
                                  // 2x2 sub-matrix.
\end{code}

You probably get the idea. Note that the second argument of Range gives
the number of elements to take, not the index of the last element. See
the section about Vectors for more examples how to use Range. You can
use \shortcode{joker}, and also the third argument of Range to select only
every nth row, or column, or reverse the order of the rows or columns.

In analogy to the Vector case, you should use the two classes
\shortcode{MatrixView} and \typeindex{ConstMatrixView} as function arguments.
Please refer to the discussion in the Vector section for details. As
in the case of VectorViews, all arguments of these types should be
passed by value, not by reference. Also, similar to the Vector case,
workspace methods are the exception, because they have to use the
conventional \shortcode{const Matrix\&} or \shortcode{Matrix\&} as input/output
arguments.

\subsection{What you can do with a Matrix (or MatrixView)}
%-------------------------------------------------------------------------

All examples below (except for the first) assume that \shortcode{a} is a
Matrix or MatrixView.

\subsubsection{Resize (only for Matrix, not for MatrixView!):}
\begin{code}
a.resize(5,10);
\end{code}
This makes \shortcode{a} a 5x10 Matrix (5 rows, 10 columns). The new Matrix
is not initialized (i.e., the contents will be unpredictable). Also,
note that the previous content will be completely lost.

\subsubsection{Get the number of rows or columns:}
\begin{code}
cout << a.nrows();
cout << a.ncols();
\end{code}

\subsubsection{Refer to a row or column:}
\begin{code}
Vector x = a(0,Range(joker));           // First row.
Vector y = a(Range(joker),a.ncols()-1); // Last column.
\end{code}
Of course, you can use more complicated Range expressions to refer to
only parts of a row or column. However, the case that you want all
elements of a given dimension is so
much more common than the more sophisticated uses of the
\shortcode{Range} class, that it is worth to introduce a simplified
notation for this case. Therefore, you are allowed to
omit the \shortcode{Range} and just write:
\begin{code}
Vector x = a(0,joker);           // First row.
Vector y = a(joker,a.ncols()-1); // Last column.
\end{code}

Technically, expressions of this kind return the type
\shortcode{VectorView}. This means, they can be used in all cases
where an object of that type is expected, for example with the
function defined in Section \ref{sec:vector_views}:

\begin{code}
silly_function(a(0,Range(joker)),
               a(1,Range(joker))); // Call silly_function
                                   // with first and
                                   // second row of a.
\end{code}



\subsubsection{Element access:}
\begin{code}
cout << a(3,4); // Print that element.
a(0,0) = 3.5;   // Assign 3.5 to the top-left element
\end{code}

Note that we use 0-based indexing! Furthermore note that the operator
`()' can be also used with one or two \shortcode{Range} arguments, as
explained above. To summarize:

\begin{itemize}
\item (\shortcode{Index},\shortcode{Index}) returns \shortcode{Numeric} (element access).
\item (\shortcode{Index},\shortcode{Range}) or (\shortcode{Range},\shortcode{Index}) returns
  \shortcode{VectorView} (row or column access).
\item (\shortcode{Range},\shortcode{Range}) returns \shortcode{MatrixView}
  (sub-matrix access).
\end{itemize}

You may find it unlogical, that Matrix uses `()' for indexing, whereas
Vector uses `[]'. However, using `[]' for Matrix is not possible,
since it can have only one argument. On the other hand, using `()' for
Vector element access seemed not a good idea, since that would break
with the established use of `[]' for element access in C and C++.

\subsubsection{Copying Matrices:}
\begin{code}
Matrix b;
b = a;
\end{code}

As in the case of Vectors, the `=' operator adjusts the size of the
target automatically.

\subsubsection{Copying in connection with views:}

As in the case of Vectors, the `=' operator copies only the
\emph{contents} for views, so the dimensions must match. An attempt to justify
this behavior has been made above in the Section about Vector. As for
Vector, you can use `=' with complicated expressions. Here is a more
elaborate example:

\begin{code}
b(Range(0,3),Range(0,4)) =
   a(Range(10,3),Range(3,4,-1)); // Copy a row 10-12,
                                 // column 0-3
                                 // to b row 0-2,
                                 // column 0-3, reversing
                                 // the order of columns.
\end{code}
Note that in this case the dimensions must match exactly, as explained
in the Section about Vector.

If you do not understand the use of Range here, refer to Section
\ref{sec:vector_views}. 

\subsubsection{Assigning a scalar:}
\begin{code}
a = 1.0;                        // Assign 1 to all elements.
\end{code}

\subsubsection{Mathematical operators:}

You can use the operators `+=', `-=', `*=', and `/=', which operate
element-vise, just as for Vector.

\subsubsection{Maximum, minimum and mean:}
\begin{code}
cout << max(a);
cout << min(a);
cout << mean(a);
\end{code}
These operations act on the complete matrix. Hence, the output is a scalar
value.

\subsubsection{Arbitrary single-argument math functions:}

The function \funcindex{transform} works just like for Vector.

\subsubsection{Transpose:}
\begin{code}
Matrix b = transpose(a); // Make b the transpose of a.
\end{code}

The function \funcindex{transpose} creates a MatrixView, for which rows and
columns are interchanged. Note, that only the way the data is accessed
is changed, not the data itself. So Matrix \shortcode{a} in the example
above is not changed. For this reason, transposing is very efficient.
You can use \shortcode{transpose(a)} instead of \shortcode{a} in any matrix
expression practically without additional cost. (This is not strictly
true, after all, the view has to be generated and passed. But that
cost should be negligible except for very small matrices.)

\subsubsection{Matrix multiplication:}
\begin{code}
// Matrix-Vector:
Vector b(a.nrows()), c(a.ncols());
mult(b,a,c);                           // b = a * c

// Matrix-Matrix:
Matrix d(a.nrows(),5), e(a.ncols(),5);
mult(d,a,e);                           // d = a * e
\end{code}

Note, that the result is put in the first argument, consistent with
the general ARTS policy, but different from the old MTL based
multiplication function. Furthermore note, that as you can see from
the first example, a Vector is always considered to be a 1-column
Matrix.

\textbf{Important: The matrices or vectors that you give for the three
arguments must not overlap, or you will get garbage.} In particular,
this means that
\begin{code}
mult(x,y,x);            // x = y * x FORBIDDEN!!!
\end{code}
does not work. No, even worse: It works, but it gives the wrong
result.  The reason for this behavior is that the result is
constructed in the first argument variable. If that is also an input
variable it will change while it is multiplied, which will lead to a
different result.  There is no efficient way to detect overlap, so the
only way to allow input and output arguments to be identical would be
to use another internal dummy variable to store the result. However,
this would be much less efficient.

Another thing: You can use transpose, of course. These two examples
should obviously give the same result:

\begin{code}
// Define b and c as in first example above.
mult(c,transpose(a),b);                      // c = a' * b

// Vector-Matrix:
mult(transpose(c),transpose(b),a);           // c' = b' * a
\end{code}


\section{Tensors}
%-------------------------------------------------------------------------
\label{sec:matpack:tensors}

ARTS has tensors with rank 3 to 7. They are called
\typeindex{Tensor3}, \typeindex{Tensor4}, \typeindex{Tensor5},
\typeindex{Tensor6}, \typeindex{Tensor7}, and work very much like
matrices, just with more dimensions. Some properties:
\begin{itemize}
\item A Tensor contains floating point values of type \builtindoc{Numeric}.
\item The \emph{rank} of a tensor means the number of dimensions, so a
  \builtindoc{Tensor4} has 4 dimensions. Tensors of different rank are
  different classes. That means, the rank is fixed at compile time and
  cannot be changed at runtime. We will use rank and dimension as
  synonyms. 
\item The different dimensions are named:
\begin{itemize}
\item Library
\item Vitrine
\item Shelf
\item Book
\item Page
\item Row
\item Column
\end{itemize}
For example, \shortcode{Tensor3 b(2,4,3)} defines a third order tensor with 2
pages, 4 rows, and 3 columns. Note that the column dimension is always
last. (Incidentally, a \builtindoc{Matrix} behaves exactly like a second
order tensor, except that it has some additional features.) 
\item A sub-range of a tensor in all dimensions (sub tensor) can
  easily be accessed, and used as if it was just a normal tensor.
\item More importantly, you can easily access lower dimensional
  `slices' of a tensor. 
\item Resizing a tensor is expensive and should be avoided.
\end{itemize}

\subsection{Constructing a tensor}
%-------------------------------------------------------------------------
You can construct an object of a tensor class like this:

\begin{code}
Tensor7 a;            // Create empty tensor of rank 7
Tensor3 b(2,4,3);     // 2 pages, 4 rows, 3 columns
Tensor3 c(2,4,3,0.0); // Similar, but
                      // fill it with 0.

Tensor3 d=c;          // Make d a copy of c.
\end{code}


\subsection{Tensor views}
%-------------------------------------------------------------------------

Tensor views work exactly like matrix and vector views. Example:
\begin{code}
Tensor4 a(10,20,5,4);
Tensor3View b = a(3,Range(1,3),joker,joker);
\end{code}
If you have read the previous sections carefully, it should be clear
what this expression does.

This is what is meant by slicing: You can easily create a view of a
tensor that picks out an object of lower dimension. Note that you can
use either an \builtindoc{Index} or a \shortcode{Range}
argument\footnote{Using just \shortcode{joker} is equivalent to using
  \shortcode{Range(joker)}, as explained in Section
  \ref{sec:matpack:matrices}.} for any of the dimensions. The
dimensionality of the result will adjust accordingly, as in the
example above.

Everything that was said about matrix and vector views holds also
here. In particular, please always use views as function arguments. 

\subsection{What you can do with a tensor (or tensor view)}
%-------------------------------------------------------------------------

All examples below (except for the first) assume that \shortcode{a} is a
Tensor7 or Tensor7View.

\subsubsection{Resize (only for tensors, not for views):}
\begin{code}
a.resize(5, 10, 4, 5, 3, 6, 8);
\end{code}
This makes \shortcode{a} the requested size. The new tensor is not
initialized (i.e., the contents will be unpredictable). Also, note
that the previous content will be completely lost.

\subsubsection{Get the extent of the various dimensions:}
\begin{code}
Index nl = a.nlibraries();
Index nv = a.nvitrines();
Index ns = a.nshelves();
Index nb = a.nbooks();
Index np = a.npages();
Index nr = a.nrows();
Index nc = a.ncols();
\end{code}

Which of these functions are available depends on the dimension of
your tensor. For example, \shortcode{nlibraries()} is only available for
\builtindoc{Tensor7}. Note, that I took care that the first letters of the
dimension names are unique, which is very convenient if you prefer short
names for your variables that refer so some dimension of a tensor.

\subsubsection{Slicing:}
\begin{code}
Vector x = a(0,2,1,8,3,4,joker);
// Select row 4
// on page 3 
// in book 8
// on shelf 1 
// in vitrine 2 
// in library 0
// and copy it to the Vector x.
\end{code}

Any \shortcode{Range} or \builtindoc{Index} expression is allowed in any of the
arguments, of course.

\subsubsection{Element access:}
\begin{code}
cout << a(3,4,0,0,0,0,0); // Print that element.
a(0,0,0,0,0,0,0) = 3.5;   // Assign 3.5 to this element.
\end{code}

\subsubsection{Copying tensors:}

Works exactly like copying matrices. Size of output argument is
adjusted for Tensors, but must already have the correct size for
TensorViews.

\subsubsection{Assigning a scalar:}
\begin{code}
a = 1.0;                        // Assign 1 to all elements.
\end{code}

\subsubsection{Mathematical operators:}

You can use the operators `+=', `-=', `*/', and `/=', which operate
element-vise, just as for Vector.

\subsubsection{Maximum and minimum:}
\begin{code}
cout << max(a);
cout << min(a);
\end{code}

\subsubsection{Arbitrary single-argument math functions:}

The function \funcindex{transform} works just like for Vector.

\subsection{Making things appear larger than they are}

Assume that you have written a function that performs some calculation
for a \builtindoc{Tensor5}:
\begin{code}
void my_function(Tensor5View x);
\end{code}

\noindent Can you call this function with a \builtindoc{Tensor4}? Yes, you can:
\begin{code}
Tensor4 a;                      
Tensor5View b = a;              // The extent of the first
                                // dimension of b will be 1.
my_function(b);                 // Call the function.
\end{code}

In general, you can always create a view that is one dimension bigger
than what you have. The leading dimension then has extent 1. There is
one important exception: If you interpret a \builtindoc{Vector} as a
\builtindoc{Matrix}, the trailing dimension will be 1, not the leading
dimension. This is necessary, because the vector has to act like a
column vector, so that matrix-vector products work in the normal way.
Of course you can use telescoping to blow up anything to
\builtindoc{Tensor7}: 

\begin{code}
Numeric b = 3.1415;             // Just any number here.
Tensor7View bt7 =
  Tensor6View(
     Tensor5View(
        Tensor4View(
           Tensor3View(
              MatrixView(
                 VectorView(b)
              )
           )
        )
     )
  );                 // All dimensions of bt7 will be 1!
\end{code}

This kind of conversion works also implicitly. So, instead of the
first example, you could have simply written:
\begin{code}
my_function(a);
\end{code}

The purpose of this feature is to avoid having to make special
versions of auxiliary functions for all different tensor
dimensions. Use it wisely! There is very little runtime overhead for
this, since the data itself is not copied.

\subsection{Summary}

This is all. In particular, we have no tensor products, since they
are not needed. So, tensors are mostly used to store things, notably
atmospheric field and so on.  There is a set of sophisticated
interpolation routines, though, which are described separately in
Chapter \ref{sec:interpolation}.


\section{Arrays}
%-------------------------------------------------------------------------
\label{sec:matpack:arrays}

The template class \typeindex{Array} can be used to make arrays out of
anything. I do not know a good definition for `array', but I guess
anybody who has written a computer program in any programming language
is familiar with the concept. Of course, it is rather similar to the
concept of a Vector, just missing all the mathematical functionality
like Matrix-Vector multiplication and sub-range access.

The implementation of our \shortcode{Array} class is based on the STL class
\shortcode{std::vector}, whereas the implementation of our \builtindoc{Vector}
class is done from scratch. So the two implementations are completely
independent. Nevertheless, I tried to make \shortcode{Array} behave
consistently with \builtindoc{Vector}, as much as possible. There are a number
of important differences, though, hopefully sufficiently explained in
this part. A short summary of important differences:

\begin{itemize}
\item An Array can contain elements of any type, whereas a Vector
  always contains elements of type Numeric.
\item No mathematical functionality for Array (no sub-ranges (nothing
  like VectorView); no +=, -=, *=, /=; no scalar product; no
  \shortcode{transform} function; no \shortcode{mult} function; no
  \shortcode{transpose} function).
\item On the other hand, resizing (for example adding to the end) of
  an Array is ok. (See the \shortcode{push\_back} method below.) It is still
  rather expensive, though, at least for large Arrays. 
\end{itemize}

\subsection{Constructing an Array}
%-------------------------------------------------------------------------
You can construct an object of an Array class like this:

\begin{code}
Array<Index>  a;        // Empty Array of class Index.

Array<String> b(5);     // String Array with 5
                        // elements. Without initialization, 
                        // elements contain random values.
Array<String> c(5,"x"); // The same, but fill with "x".

Array<Index>  d=a;      // Make d a copy of a;
\end{code}

There are already a lot of predefined Array classes. The naming
convention for them is: \typeindex{ArrayOfIndex}, \typeindex{ArrayOfString},
etc.. Normally you should use these predefined classes. But if you want
to define an Array of some uncommon type, you can do it with `$<>$',
as in the above examples.

As for Vector, there is a special sub-class of Array that can be
initialized explicitly. Usage:

\begin{code}
MakeArray<String> a("ARTS",
                    "is",
                    "great"); // Creates an array of String
                              // with these 3 elements.
\end{code}

\subsection{What you can do with an Array}
%-------------------------------------------------------------------------

All examples below assume that \shortcode{a} is an ArrayOfString.

\subsubsection{Resize:}
\begin{code}
a.resize(5);
\end{code}

This adjusts the size of \shortcode{a} to 5. Resizing is more efficiently
implemented than for Vector, but still expensive.

\subsubsection{Get the number of elements:}
\begin{code}
cout << a.nelem();  // Just as for Vector.
\end{code}

In particular, note that the return type of this method is
\builtindoc{Index}, just as for Vector. This is an extension compared to
std::vector, which just has a method \shortcode{size()} that returns the
positive integer type \shortcode{size\_t}.

\subsubsection{Element access:}
\begin{code}
cout << a[3];   // Print 4th element.
a[0] = "Hello"; // Assign string "Hello" to first element.
\end{code}

In other words, this works just like for Vector.

\subsubsection{Copying Arrays:}

This works also the same as for Vector. The size of the target must
match! In this respect, I have modified the behavior with respect to
the underlying std::vector, which has different copy semantics.

\subsubsection{Assigning a scalar of the base type:}
\begin{code}
a = "Hello";    // Assign string "Hello" to all elements.
\end{code}

\subsubsection{Append to the end:}
\begin{code}
a.push_back("Hello"); // Adds this new element at the
                      // end of a.
\end{code}

This can be an expensive operation, especially for large Arrays.
Therefore, use it with care. Actually, the \funcindex{push\_back} method
comes from the \shortcode{std::vector} class that Array is based on. You
can do a lot more with \shortcode{std::vector}, all of which also works
with \shortcode{Array}. However, to explain the Standard Template Library
is beyond the scope of this text. You can read about it in C++ or even
dedicated STL textbooks.

\section{Sparse matrices}
%-------------------------------------------------------------------------
\label{sec:matpack:sparse}

The class \typeindex{Sparse} implements the mathematical concept of a
matrix, same as Matrix does, but the data is stored in a
different manner. Sparse offers a memory saving storage when
most of the matrix is filled with zeros. This means that:
\begin{itemize}
\item A Sparse contains floating point values of type \builtindoc{Numeric}.
\item The values are arranged in rows and columns in the same ways as for
  ordinary matrices, in \emph{row-major} order.
\item A Sparse can be multiplied with a Vector, a Matrix or with another
  Sparse.
\item There exist no views for Sparse.
\item Resizing a Sparse is expensive and should be avoided.
\end{itemize}

To calculate the maximum number of non-zero elements for efficient storage,
take the product of number of columns and number of rows, subtract the
number of columns plus one and then divide by two,
\mbox{($nnz \leq 0.5\times (ncols\times nrows - (ncols+1)$).}

\subsection{Constructing a Sparse}
%-------------------------------------------------------------------------
You can construct an object of class Sparse in any of these ways:

\begin{code}
Sparse a;          // Create empty Sparse.
Sparse b(3,4);     // Create Sparse with 3 rows
                   // and 4 columns. When
                   // created like this it will
                   // contain only zeros, i.e.
                   // be an empty Sparse.

Sparse d=c;        // Make d a copy of c.
\end{code}

That is all. As for Matrix there exist no more fancy constructors, like an
equivalent to the \shortcode{MakeVector} class.

\subsection{What you can do with a Sparse}
%-------------------------------------------------------------------------
All examples below assume that \shortcode{a} is a Sparse.

\subsubsection{Identity matrix:}
\begin{code}
a.make_I(10,10);
\end{code}
This sets \shortcode{a} to be the identity matrix of size 10x10 (10 rows 
and 10 columns). Using this function is much faster than setting the 
diagonal elements to one by yourself. The number of rows and columns 
doesn't have to match each other. In the case that they don't, the rule 
that each row and column only will have only one position occupied by a one
is applied. That is, in the case there are more rows than columns, the last
rows will be empty.

\subsubsection{Resize:}
\begin{code}
a.resize(5,10);
\end{code}
This makes \shortcode{a} a 5x10 Sparse (5 rows, 10 columns). Note that the
previous content will be completely lost. The new Sparse will be empty.

\subsubsection{Get the number of rows, columns or non-zero elements:}
\begin{code}
cout << a.nrows();
cout << a.ncols();
cout << a.nnz();
\end{code}

\subsubsection{Element access:}
There are two different ways to access individual elements. One used for
read only and one for read and write. The distinction is necessary since
the read and write method creates elements if they don't already exist.
Note that we use 0-based indexing. For reading only use:
\begin{code}
cout << a.ro(3,4);  // Print that element. If it
                    // it doen't exist a zero will
                    // be printed.
cout << a(0,0);     // Short version of the above.
\end{code}

For reading and writing, such as assigning values to elements, use:
\begin{code}
a.rw(0,0) = 1.5;    // Assigns the value 1.5 to the
                    // first row and first column.
cout << a.rw(0,0);  // Also returns the value of the
                    // first row and first column,
                    // if the element doesn't exist
                    // it will be created and set
                    // to zero.
\end{code}

\subsubsection{Copying Matrices:}
\begin{code}
Sparse b;
b = a;
\end{code}

As in the case of Vectors, the `=' operator adjusts the size of the
target automatically.

\subsubsection{Transpose:} The function \shortcode{transpose} works a bit
differently for Sparse than for Vector and Matrix. This is due to the
fact that we don't have any views for Sparse. Thus,
\shortcode{transpose} for a Sparse creates a new Sparse variable that
contains the transpose of the original Sparse, whereas
\shortcode{transpose} for a Matrix just creates a transposed view of
the original Matrix.

The target variable for the transposed Sparse has to have the right
dimensions before the function is called.
\begin{code}
Sparse b(a.ncols(),a.nrows());
transpose(b,a);     // Make b the transpose of a.
                    // Note the argument order!
\end{code}

\subsubsection{Matrix multiplication:}
\begin{code}
// Sparse-Vector
Vector b(a.nrows()), c(a.ncols());
mult(b,a,c);        // b = a * c

// Sparse-Matrix
Matrix d(a.nrows(),5), e(a.ncols(),5);
mult(d,a,e);        // d = a * e

// Sparse-Sparse
Sparse f(a.nrows(),5), g(a.ncols(),5);
mult(f,a,g);        // f = a * g
\end{code}
The result is put in the first argument, consistent with the Matrix class.
Note that for the Sparse -- Matrix multiplication the output is a Matrix.
\textbf{Important: As for Matrix, the matrices or vectors
that you give for the three arguments must not overlap, or you will get
garbage.}




%%% Local Variables:
%%% mode: latex
%%% TeX-master: "uguide"
%%% End:

