Linear algebra

The main categories of objects for doing linear algebra are

• VectDse: dense vectors,

• MatrDse: dense matrices,

• MatrStr: structured matrices,

• MatrSps: sparse matrices. 0

Vectors

If one want to construct vectors with double coefficients, one can use for instance the definition:

 typedef VectDse<double> vect_t;

This is equivalent to

 typedef VectDse<double,linalg::rep1d<double> > vect_t;

the default value for the second argument of the parameterised type VectDse being linalg::rep1d<C>. It specifies the container type linalg::rep1d<C> which is

• a one-dimensional array of coefficients of type C,

• with a field corresponding to its size.

But one can also use the containor std::vector<C> of the standard library or any other convenient container, which provides the same functions. See the file synaps/linalg/VectDse.h for more details on the implementation of this view.

A vector can be constructed and manipulated as follows:

 double u[]={1,-2,1}; vect_t V(3,u);
 vect_t W= V*2, U; U -= W+3*V;

#include <iostream>
#include <synaps/linalg/lapack.h>
#include <synaps/linalg/MatrDse.h>
#include <synaps/linalg/Eigen.h>
#include <synaps/linalg/Svd.h>

typedef VectDse<double> vect_t;
typedef MatrDse<double,lapack::rep2d<double> > matr_t;

int main(int argc, char** argv) 
{
  using std::cout;
  using std::endl;
  double u[]={1,-2,1}; vect_t V(3,u);

  vect_t W= V*2, U;
  U -= W+3*V;

  cout<<U<<endl;
  cout<<V<<endl;
  cout<<W<<endl;
}

Matrices

For the definition of dense matrices with double coefficients, you can use the definition:

 typedef MatrDse<double> matr_t;

Here the default container is linalg::rep2d<C> which extends linalg::rep1d<C>, by specifying the number of rows and columns If instead, we use the definition

 typedef MatrDse<double,lapack::rep2d<double> > matr_t;

that is the container lapack::rep2d<C>, we are able to use the routines of the LAPACK library. The internal representation is also a one-dimensional array with a number of rows, a number of columns, a size bigger than the product of these two last numbers and a way to access to the element of index (i, j). In order to be able to use this container, the following include instruction should be given:

#include <synaps/linalg/lapack.h>

More details on the implementation of the matrix interface can be found in the header file synaps/linalg/MatrDse.h.

Here are some illustrations of its use:

 matr_t A(3,3,"1 2 3 4 5 6 7 8 9 10",ByRow());

The matrix A is defined from a string of coefficients. It is a 3x3 matrix, filed by row. It is printed as:

The usual operations on matrices and vectors are available from this type:

 using std::cout; using std::endl;
 vect_t W = A*V; cout << W <<endl;

 matr_t B = A*A; cout << B << endl; 

 matr_t C(-A); C += B + A*2; cout <<C<< endl; 

Numerical linear algebra

Standard routines from numerical linear algebra are available on matrices with lapack containers, one the corresponding header files:

#include <synaps/linalg/Svd.h>

#include <synaps/linalg/Eigen.h>

are included. Some implementation (such as Svd) are also available for generic matrices.

 cout<< Svd(A)<<endl; cout << Eigenval(A) << endl; 

Here the type of the eigenvalues is complex<double>. They are printed with real and imaginary part between parentheses.

#include <synaps/init.h>
#ifdef SYNAPS_WITH_LAPACK
#include <iostream>
#include <synaps/linalg/lapack.h>
#include <synaps/linalg/MatrDse.h>
#include <synaps/linalg/Eigen.h>
#include <synaps/linalg/Svd.h>

typedef VectDse<double> vect_t;
typedef MatrDse<double,lapack::rep2d<double> > matr_t;

int main(int argc, char** argv) 
{
  using std::cout;
  using std::endl;
  matr_t A(3,3,"1 2 3 4 5 6 7 8 9",ByRow());
  cout<<A<<endl;
  cout<< Svd(A)<<endl;
  cout<< Eigenval(A)<<endl;
}
#else
int main(int argc, char** argv) {}
#endif //SYNAPS_WITH_LAPACK

Exercise: Estimate the mean time for computing the svd of a nxn matrix by iterating the computation on N random matrices with coefficients between -1 and 1 (using the routine 1-(2.0*rand())/RAND MAX of stdlib.h). Compare it with the representation using lapack routine. Compare with different compilation options. To count the time, we will use the instructions {{Clock c; c.start(); ...; c.stop(); c.time();}} provided by the file synaps/base/Clock.h.

Views on vectors and matrices

The linear algebra classes also offers local views on subobjects. A row (or a column) of a matrix, for instance, is a view on the corresponding elements of the matrix and share the properties of a vector. Indeed, it is a VectDse with a special container and can be manipulated as such. The types of the columns and the rows are accesible as

 MatrDse<double>::row_type
 MatrDse<double>::col_type

The characteristic of these views is that performing operations on the vectors will affect the object the view is looking at. For instance,

 Col(A,1) = Col(A,0);

transforms the matrix

into:

Similarly

 Row(A,1) -= Row(A,0)*4;

yields:

These functions Col and Row produce vector views on the matrix. They can be interprated as references on the internal data. Remember that the memory space of the views Row(A,i), Col(A,j) are not duplicated but shared with the matrix A. Such a vector view can be copied into a standard vector:

 vect_t V1 = Col(A,1)*2, V2= Row(A,0)+ Row(A,1);
 vect_t W = V1+V2;

yields

#include <iostream>
#include <synaps/linalg/MatrDse.h>
#include <synaps/linalg/lapack.h>
#include <synaps/linalg/Eigen.h>
#include <synaps/linalg/Svd.h>

typedef VectDse<double> vect_t;
typedef MatrDse<double,lapack::rep2d<double> > matr_t;

int main(int argc, char** argv) 
{
  using std::cout;
  using std::endl;
  matr_t A(3,3,"1 2 3 4 5 6 7 8 9",ByRow());
  std::cout<<A<<std::endl;

  col(A,1) = col(A,0);
  std::cout<<A<<std::endl;

  row(A,1) -= row(A,0)*4;
  std::cout<<A<<std::endl;

//   vect_t V1 = col(A,1)*2; 
//   std::cout<<V1<<std::endl;
  vect_t V2= row(A,0)+ row(A,1);
  std::cout<<V2<<std::endl;

  vect_t W = V2+V2;
  std::cout<<W<<std::endl;
}