A sub-matrix can be defined with any instance of a dense or sparse matrix. Arbitrary rows and columns may be selected.
Given a base matrix A, a sub-matrix is declared with a list of rows row_list and a list of columns column_list:
Matrix<double> A(4, 6); Vector<int> row_list(2); row_list(0) = 1; row_list(1) = 2; Vector<int> column_list(3); column_list(0) = 0; column_list(1) = 1; column_list(2) = 5; SubMatrix<Matrix<double> > SubA(A, row_list, column_list);
The template argument of SubMatrix must be the exact type of the base matrix A.
The row and column lists can be any index list, providing the indexes refer to an existing row or column in the base matrix. The same index may appear several times in a list, and the lists do not need to be sorted. As a consequence, the sub-matrices might be used for a more general purpose than simply extracting a sub-set of the rows and columns. For example, one may swap two columns and duplicate a line:
Matrix<double> A(3, 3); A.Fill(); cout << "Complete matrix:" << endl; A.Print(); Vector<int> row_list(4); row_list.Fill(); row_list(3) = 2; Vector<int> column_list(3); column_list(0) = 1; column_list(1) = 0; column_list(2) = 2; SubMatrix<Matrix<double> > SubA(A, row_list, column_list); cout << "Sub-matrix:" << endl; SubA.Print();
will result in
Complete matrix: 0 1 2 3 4 5 6 7 8 Sub-matrix: 1 0 2 4 3 5 7 6 8 7 6 8
A sub-matrix is a light structure. It essentially stores the list of rows and columns, along with a pointer to the base matrix. When an element of the sub-matrix is accessed, the sub-matrix performs an indirection to the corresponding element in the base matrix. Hence, if the sub-matrix is modified, so will be the base matrix. This also means that the sub-matrix should not be manipulated once the base matrix is destroyed.
A sub-matrix is viewed as a matrix by Seldon. Therefore, all generic functions that are implemented in C++ can be called, which includes:
Warning: In the case of sub-matrices, these functions will not call Blas, Lapack or another external library. As a consequence, they may be slow.
#define SELDON_DEBUG_LEVEL_4
#include "SeldonLib.hxx"
using namespace Seldon;
int main()
{
TRY;
/*** Full matrix ***/
Matrix<double> A(4, 6);
A.Fill();
cout << "Complete matrix:" << endl;
A.Print();
Vector<int> row_list(2);
row_list(0) = 1;
row_list(1) = 2;
Vector<int> column_list(3);
column_list(0) = 0;
column_list(1) = 1;
column_list(2) = 5;
SubMatrix<Matrix<double> > SubA(A, row_list, column_list);
cout << "Sub-matrix:" << endl;
SubA.Print();
// Basic operations are supported, but they are slow (Blas/Lapack will not
// be called).
Vector<double> X(3), Y(2);
X.Fill();
cout << "Multiplied by X = [" << X << "]:" << endl;
Mlt(SubA, X, Y);
Y.Print();
/*** Symmetric matrix ***/
Matrix<double, General, ColSymPacked> B(4);
B.Fill();
cout << "\nComplete matrix:" << endl;
B.Print();
row_list(0) = 1;
row_list(1) = 3;
column_list(0) = 0;
column_list(1) = 2;
column_list(2) = 3;
SubMatrix<Matrix<double, General, ColSymPacked> >
SubB(B, row_list, column_list);
cout << "Sub-matrix (no more symmetric):" << endl;
SubB.Print();
// Assignments in the sub-matrix.
for (int i = 0; i < 2; i++)
for (int j = 0; j < 3; j++)
SubB(i, j) = -1.;
// 'B' will remain symmetric.
cout << "Complete matrix after the sub-matrix is filled with -1:" << endl;
B.Print();
END;
return 0;
}
Output:
Complete matrix: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 Sub-matrix: 6 7 11 12 13 17 Multiplied by X = [0 1 2]: 29 47 Complete matrix: 0 1 3 6 1 2 4 7 3 4 5 8 6 7 8 9 Sub-matrix (no more symmetric): 1 4 7 6 8 9 Complete matrix after the sub-matrix is filled with -1: 0 -1 3 -1 -1 2 -1 -1 3 -1 5 -1 -1 -1 -1 -1