Distributed vectors and matrices
Distributed vectors and matrices are implemented in Seldon. Unitary tests are present in files src/Program/Unit/Algebra/distributed_vector.cc, src/Program/Unit/Algebra/distributed_matrix_test.cc and src/Program/Unit/Algebra/distributed_solver.cc for respectively distributed vectors, distributed matrices and distributed solvers. In Montjoie, distributed block-diagonal matrices are implemented. Below, we show a basic example with these matrices
// for the constructor, you use the local number of rows
DistributedMatrix<double, General, BlockDiagRow> A(n, n);
// then you need to initialise the arrays used for parallel functions (e.g. matrix-vector product)
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
n, 1, &SharingProc, &SharingRows, MPI::COMM_WORLD);
// example with two processors
// processor 0 owns rows 0, 2, 5, 6, 8
// processor 1 owns rows 1, 2, 3, 4, 6, 7
// block 0 : (0, 1, 4)
// block 1 : (2, 3, 6)
// block 2 : (5, 7, 8)
Vector<IVect> pattern(3), proc(3);
for (int k = 0; k < 3; k++)
{
pattern(k).Reallocate(3);
proc(k).Reallocate(3);
}
if (MPI::COMM_WORLD.Get_rank() == 0)
{
// block 0 : (0, 1, 4)
// rows 1 and 4 are owned by processor 1, we put global numbers
pattern(0)(0) = 0; pattern(0)(1) = 1; pattern(0)(2) = 4;
proc(0)(0) = 0; proc(0)(1) = 1; proc(0)(2) = 1;
// block 1 : (2, 3, 6)
// local numbers of rows 2, 6 : 1, 3
pattern(1)(0) = 1; pattern(1)(1) = 3; pattern(1)(2) = 3;
proc(1)(0) = 0; proc(1)(1) = 1; proc(1)(2) = 0;
// block 2 : (5, 7, 8)
// local numbers of rows 5, 8 : 2, 4
pattern(2)(0) = 2; pattern(2)(1) = 7; pattern(2)(2) = 4;
proc(2)(0) = 0; proc(2)(1) = 1; proc(2)(2) = 0;
}
else
{
// block 0 : (0, 1, 4)
// local numbers of rows 1, 4 : 0 3
pattern(0)(0) = 0; pattern(0)(1) = 0; pattern(0)(2) = 3;
proc(0)(0) = 0; proc(0)(1) = 1; proc(0)(2) = 1;
// block 1 : (2, 3, 6)
// local numbers of rows 2, 3, 6 : 1, 2, 4
pattern(1)(0) = 1; pattern(1)(1) = 2; pattern(1)(2) = 4;
proc(1)(0) = 1; proc(1)(1) = 1; proc(1)(2) = 1;
// block 2 : (5, 7, 8)
// local number of row 7 : 5
pattern(2)(0) = 5; pattern(2)(1) = 5; pattern(2)(2) = 8;
proc(2)(0) = 0; proc(2)(1) = 1; proc(2)(2) = 0;
}
A.SetPattern(pattern, proc);
// once the pattern has been defined, you can add non-zero entries to the matrix
A.AddInteraction(1, 3, -2.5);
// if the column is distant => AddDistantInteraction
A.AddDistantInteraction(0, 4, 1, 0.8);
// another solution is to convert a distributed sparse matrix into a block-diagonal one
// in that case SetPattern is not needed, you just call ConvertToBlockDiagonal
// once the matrix is filled, you can compute its inverse
GetInverse(A);
Methods for block-diagonal distributed matrices :
| Clear | clears the matrix |
| Init | sets pointers to the arrays containing global row numbers and overlapped rows |
| GetCommunicator | returns the MPI communicator associated with the matrix |
| GetGlobalM | returns the global number of rows |
| GetNodlScalar | returns the number of rows for a scalar unknown |
| GetNbScalarUnknowns | returns the number of scalar unknowns |
| GetProcNumber | returns a pointer to the array containing the processor number associated with each row. |
| GetGlobalRowNumber | returns local to global numbering |
| GetOverlapRowNumber | returns the array containing the numbers of rows already handled by an another processor |
| GetOverlapProcNumber | returns the array containing the processor numbers of rows already handled by an another processor |
| GetProcessorSharingRows | returns the list of processors that share rows with the current processor |
| GetSharingRowNumbers | returns the list of rows shared with each processor |
| SetPattern | provides the structure of the block-diagonal matrix |
| AddDistantInteraction | adds a value to A(i, j) where i is local and j global |
| AddRowDistantInteraction | adds a value to A(i, j) where i is global and j local |
| IndexGlobal | returns the global row number of row j of block i |
| WriteText | writes the distributed matrix on several files (one per processor) in ascii format |
| IsReadyForMltAdd | returns true if the structure is ready to perform a matrix vector without preparation |
| Assemble | assembles the block-diagonal matrix |
Functions for distributed block-diagonal matrices :
| Mlt | performs a matrix-vector product |
| MltAdd | performs a matrix-vector product |
| Add | adds a distributed matrix to another one |
| ConvertToSparse | converts a block-diagonal matrix to a sparse matrix |
| ConvertToBlockDiagonal | converts a sparse matrix to a block-diagonal matrix |
| Copy | Copies/Converts a distributed matrix into another one |
| GetInverse | replaces a block-diagonal matrix by its inverse |
| GetCholesky | replaces a block-diagonal matrix by its Cholesky factorisation |
| SolveCholesky | solves L x = b or L^T x = b for a block-diagonal matrix |
| MltCholesky | performs product with Cholesky factor L (or its transpose) |
| EraseRow | erases rows in the distributed matrix |
| EraseCol | erases columns in the distributed matrix |
GetCommunicator
Syntax
MPI::Comm& GetCommunicator();
This method returns the MPI communicator used to distribute the matrix.
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrixInline.cxx
GetGlobalM
Syntax
int GetGlobalM() const;
This function returns the global number of rows.
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// after filling nglob, global_row, overlap_row and overlap_proc
// calling Init in order to provide these datas
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
Nvol, nb_u, &ProcShared, &SharedRows, MPI::COMM_WORLD);
// you can access to nglob by calling GetGlobalM
nglob = A.GetGlobalM();
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
GetNodlScalar
Syntax
int GetNodlScalar() const;
This method returns the local number of rows for each scalar components. Several components can be specified when Init is called.
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// after filling nglob, global_row, overlap_row and overlap_proc
// calling Init in order to provide these datas
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
Nvol, nb_u, &ProcShared, &SharedRows, MPI::COMM_WORLD);
// you can access to Nvol by calling GetNodlScalar
Nvol = A.GetNodlScalar();
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
GetNbScalarUnknowns
Syntax
int GetNbScalarUnknowns() const;
This method returns the local number of scalar components. Several components can be specified when Init is called.
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// after filling nglob, global_row, overlap_row and overlap_proc
// calling Init in order to provide these datas
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
Nvol, nb_u, &ProcShared, &SharedRows, MPI::COMM_WORLD);
// you can access to nb_u by calling GetNbScalarUnknowns
nb_u = A.GetNbScalarUnknowns();
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
AddDistantInteraction
Syntax
void AddDistantInteraction(int i, int jglob, int proc, const T& val);
This member function adds val for the local row i, and the global row jglob, proc being the processor that treats the global row jglob.
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// after filling nglob, global_row, overlap_row and overlap_proc
// calling Init in order to provide these datas
A.Init(nglob, &global_row, &overlap_row, &overlap_proc,
Nvol, nb_u, &ProcShared, &SharedRows, MPI::COMM_WORLD);
// for local interaction, you have to use AddInteraction
A.AddInteraction(i, j, val);
// when the column is distant (ie located on another processor), you have to use AddDistantInteraction
// jglob : global column number
// proc : distant processor
// val : value to add
A.AddDistantInteraction(i, jglob, proc, val);
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
AddRowDistantInteraction
Syntax
void AddRowDistantInteraction(int iglob, int j, int proc, const T& val);
This member function adds val for the global row iglob, and the local column j, proc being the processor that treats the global row iglob.
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// after filling nglob, global_row, overlap_row and overlap_proc
// calling Init in order to provide these datas
A.Init(nglob, &global_row, &overlap_row, &overlap_proc,
Nvol, nb_u, &ProcShared, &SharedRows, MPI::COMM_WORLD);
// for a local interaction, you have to use AddInteraction
A.AddInteraction(i, j, val);
// when the column is distant (ie located on another processor), you have to use AddDistantInteraction
// i : local row number
// jglob : global column number
// proc : distant processor
// val : value to add
// for the global matrix, it means that A_{global_row(i), jglob} is incremented with val
A.AddDistantInteraction(i, jglob, proc, val);
// when the row is distant, you have to use AddRowDistantInteraction
// iglob : global row number
// j : local column number
// proc : distant processor
// val : value to add
// for the global matrix, it means that A_{iglob, global_row(j)} is incremented with val
A.AddRowDistantInteraction(iglob, j, proc, val);
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
IsReadyForMltAdd
Syntax
bool IsReadyForMltAdd() const
This member function returns true if the structure is ready to perform a matrix-vector product. This function is not really useful since if the user performs a matrix-vector product with a structure that is not ready, the program will create all the arrays needed to perform the matrix-vector product. The function IsReadyForMltAdd can be used for informative purposes, for example to avoid taking into account the cost of constructing all these arrays (they induce MPI communications) if the structure was not ready.
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// after filling nglob, global_row, overlap_row and overlap_proc
// calling Init in order to provide these datas
A.Init(nglob, &global_row, &overlap_row, &overlap_proc,
Nvol, nb_u, &ProcShared, &SharedRows, MPI::COMM_WORLD);
// for a local interaction, you have to use AddInteraction
A.AddInteraction(i, j, val);
// when the column is distant (ie located on another processor), you have to use AddDistantInteraction
// i : local row number
// jglob : global column number
// proc : distant processor
// val : value to add
// for the global matrix, it means that A_{global_row(i), jglob} is incremented with val
A.AddDistantInteraction(i, jglob, proc, val);
// when the row is distant, you have to use AddRowDistantInteraction
// iglob : global row number
// j : local column number
// proc : distant processor
// val : value to add
// for the global matrix, it means that A_{iglob, global_row(j)} is incremented with val
A.AddRowDistantInteraction(iglob, j, proc, val);
// the structure will be said "ready" for a matrix-vector operation
// if a first matrix-vector product has been performed
bool struct_ready = A.IsReadyForMltAdd();
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
Clear
Syntax
void Clear();
This function clears the matrix.
Example :
// default constructor DistributedMatrix<double, General, BlockDiagRow> A; // the matrix is erased A.Clear();
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
GetGlobalRowNumber
Syntax
IVect& GetGlobalRowNumber() const;
This function returns the global row numbers (local to global numbering).
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// after filling nglob, global_row, overlap_row and overlap_proc
// calling Init in order to provide these datas
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
n, 1, &SharingProc, &SharingRows, MPI::COMM_WORLD);
// you can access to global_row by calling GetGlobalRowNumber
IVect& global_row = A.GetGlobalRowNumber();
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
GetOverlapRowNumber
Syntax
IVect& GetOverlapRowNumber() const;
This function returns the overlapped row numbers, i.e. rows that are already treated by another processor.
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// after filling nglob, global_row, overlap_row and overlap_proc
// calling Init in order to provide these datas
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
n, 1, &SharingProc, &SharingRows, MPI::COMM_WORLD);
// you can access to overlap_row by calling GetOverlapRowNumber
IVect& overlap_row = A.GetOverlapRowNumber();
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
GetOverlapProcNumber
Syntax
IVect& GetOverlapProcNumber() const;
This function returns the processor numbers associated with overlapped row numbers, i.e. rows that are already treated by another processor.
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// after filling nglob, global_row, overlap_row and overlap_proc
// calling Init in order to provide these datas
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
n, 1, &SharingProc, &SharingRows, MPI::COMM_WORLD);
// you can access to original_proc by calling GetOverlapProcNumber
IVect& original_proc = A.GetOverlaProcNumber();
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
GetProcessorSharingRows
Syntax
IVect& GetProcessorSharingRows() const;
This function returns the processor numbers that share rows with the current processor.
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// after filling nglob, global_row, overlap_row and overlap_proc
// calling Init in order to provide these datas
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
n, 1, &SharingProc, &SharingRows, MPI::COMM_WORLD);
// you can access to SharingProc by calling GetProcessorSharingRows
IVect& SharingProc = A.GetProcessorSharingRows();
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
GetSharingRowNumbers
Syntax
Vector<IVect>& GetSharingRowNumbers() const;
This function returns the row numbers that are shared with other processors.
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// after filling nglob, global_row, overlap_row and overlap_proc
// calling Init in order to provide these datas
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
n, 1, &SharingProc, &SharingRows, MPI::COMM_WORLD);
// you can access to SharingRows by calling GetSharingRowNumbers
IVect& SharingProc = A.GetSharingRowNumbers();
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
WriteText
Syntax
void WriteText(string) const; void WriteText(ostream&) const;
This member function writes the matrix on files or output streams in text format. Each processor writes its part of the matrix on a different file (ending with _P0.dat, _P1.dat, etc).
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// you can write the matrix in files
// the matrix here will be written in files mat_P0.dat, mat_P1.dat, etc
A.WriteText("mat.dat");
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
MltAdd
Syntax
void MltAdd(const T& alpha, const DistributedMatrix<T>& A, const Vector<T>& x, const T& beta, Vector<T>& y, bool assemble); void MltAdd(const T& alpha, const SeldonTranspose&, const DistributedMatrix<T>& A, const Vector<T>& x, const T& beta, Vector<T>& y, bool assemble);
This function can perform matrix-vector product (with the original matrix or with its transpose). The matrix-matrix product is currently not implemented. If the communicator contains only one processor, the sequential function MltAdd (with the class Matrix) will be called, otherwise an error will be displayed during the execution.
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// global row numbers are provided with Init (glob_number)
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
Nvol, nb_u, &ProcShared, &SharedRows, MPI::COMM_WORLD);
// once the matrix is constructed, matrix-vector products can be performed
// each vector is stored as a usual Seldon vector (class Vector)
// each processor is assumed to store the values of the vector for the
// global rows given during the construction of the matrix (here, the array glob_number)
Vector<double> x(n), y(n), z;
// MltAdd will assume that the values of x for rows that are shared between processors are the same
// which is here not the case when calling FillRand
x.FillRand();
y.FillRand();
z = y;
// a solution is to call AssembleVector to ensure this property
AssembleVector(x, MPI::SUM, ProcShared, SharedRows, MPI::COMM_WORLD, Nvol, nb_u, 23);
double alpha = 2.5, beta = 0.4;
// matrix-vector product y = beta y + alpha A x
// the default value of the last argument (assemble) is true
// the result y is assembled and will contain the same result as a sequential matrix-vector product,
// the values being distributed in the different processors
MltAdd(alpha, A, x, beta, y);
// assemble = true, induces a final call to the function AssembleVector (to add values of shared rows)
// if the user wants to perform a matrix-vector product without performing this final assembling step, he puts false in the last optional argument
MltAdd(alpha, A, x, beta, z, false);
// and perform the assembling phase later
AssembleVector(z, MPI::SUM, ProcShared, SharedRows, MPI::COMM_WORLD, Nvol, nb_u, 23);
// MltAdd can be used to perform a matrix-vector with the transpose matrix or the transpose conjugate
MltAdd(alpha, SeldonTrans, A, x, beta, y);
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
Add
Syntax
void Add(const T& alpha, const DistributedMatrix<T>& A, DistributedMatrix<T>& B);
This function adds a distributed matrix multiplied by a scalar to another distributed matrix. The two matrices are assumed to have the same global row numbers.
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// global row numbers are provided with Init (glob_number)
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
Nvol, nb_u, &ProcShared, &SharedRows, MPI::COMM_WORLD);
// a way to have the same global row numbers
// is to use Init with the matrix A
DistributedMatrix<double, General, BlockDiagRow> B;
B.Init(A);
double alpha = 2.1;
// we compute B = B + alpha A
Add(alpha, A, B);
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
Copy
Syntax
void Copy(const DistributedMatrix<T>& A, DistributedMatrix<T>&);
This function converts a distributed sparse matrix into another one. The format can be different (RowSparse, ColSparse, ArrayRowSparse, etc), the local part of the matrix is converted by calling the appropriate Copy function (located in the file Matrix_Conversions.cxx), the distributed part is always stored with the same format.
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
EraseRow
Syntax
void EraseRow(IVect& num, DistributedMatrix& A);
This function erases some rows (local row numbers are provided in num) of the matrix A.
Example :
// default constructor DistributedMatrix<double, General, BlockDiagRow> A; // after filling nglob, global_row, overlap_row and overlap_proc // calling Init in order to provide these datas A.Init(nglob, &global_row, &overlap_row, &overlap_proc, MPI::COMM_WORLD); // then you can erase some rows IVect num(2); num(0) = 2; num(1) = 5; EraseRow(num, A);
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
EraseCol
Syntax
void EraseCol(IVect& num, DistributedMatrix& A);
This function erases some columns (local column numbers are provided in num) of the matrix A.
Example :
// default constructor DistributedMatrix<double, General, BlockDiagRow> A; // after filling nglob, global_row, overlap_row and overlap_proc // calling Init in order to provide these datas A.Init(nglob, &global_row, &overlap_row, &overlap_proc, MPI::COMM_WORLD); // then you can erase some columns IVect num(2); num(0) = 4; num(1) = 9; EraseCol(num, A);
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
Init
Syntax
void Init(int nglob, IVect* global_row, IVect* overlap_row, IVect* overlap_proc,
int Nvol, int nb_u, IVect* SharingProc, Vector<IVect>* SharingRows, MPI::Comm& comm);
void Init(const DistributedMatrix<T>&);
The function Init must be called after the constructor in order to provide the global row numbers (array global_row), the rows already treated by another processor (array overlap_row), the processors that treat originally these rows (array overlap_proc) and the MPI communicator. Nvol is the number of rows for a scalar component, nb_u the number of scalar components, SharingProc stores the processors that share rows with current processor, and SharingRows the local row numbers that are shared with processor SharingProc(i). These numbers are assumed to be sorted such that shared row numbers with processor j in processor i correspond to row numbers shared with processor i in processor j.
Example :
// on each processor, you specify rows that are already treated by another processor
// for example if the processor 0 handles rows [0, 3, 5, 6] and
// processor 1 the rows [1 2, 4, 5, 7], the row 5 is already treated by
// processor 0. Of course, if each row is treated by a processor and only one
// this array should be left empty
IVect OverlapRow;
int nglob = 8;
int n = 4;
if (MPI::COMM_WORLD.Get_rank() == 1)
{
n = 5;
OverlapRow.Reallocate(1);
// be careful because OverlapRow stores local numbers
// here the global row 5 has a local number equal to 3 on processor 1
OverlapRow(0) = 3;
}
// for distributed matrices
// we also need to know the global row numbers
// and for each overlapped row, the original processor
IVect glob_number, original_proc;
// for shared rows, all the row numbers shared with processors
IVect SharingProc(1); Vector<IVect> SharingRows(1);
if (MPI::COMM_WORLD.Get_rank() == 0)
{
glob_number.Reallocate(4);
glob_number(0) = 0;
glob_number(1) = 3;
glob_number(2) = 5;
glob_number(3) = 6;
// global row 5 is shared with processor 1 (local number = 2)
SharingProc(0) = 1;
SharingRows(0).PushBack(2);
}
else
{
glob_number.Reallocate(5);
glob_number(0) = 1;
glob_number(1) = 2;
glob_number(2) = 4;
glob_number(3) = 5;
glob_number(4) = 7;
// row 5 is already handled by processor 0
original_proc.Reallocate(1);
original_proc(0) = 0;
// here SharingRows is similar to OverlapRow because all shared dofs of this processor are already treated by processor 0
// global row 5 is shared with processor 0 (local number = 3)
SharingProc(0) = 0;
SharingRows(0).PushBack(3);
}
// for the constructor, you use the local number of rows
DistributedMatrix<double, General, BlockDiagRow> A(n, n);
// then you need to initialise the arrays used for parallel functions (e.g. matrix-vector product)
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
n, 1, &SharingProc, &SharingRows, MPI::COMM_WORLD);
Another example does not construct these arrays but calls the method Init of a sparse distributed matrix (not block-diagonal) to construct them automatically.
// local number of rows
int n = 4;
if (MPI::COMM_WORLD.Get_rank() == 1)
n = 5;
// each processor constructs its global numbers
// rows [0, 3, 5, 6] for processor 0
// rows [1 2, 4, 5, 7] for processor 1
// you notice that global row 5 is here shared by the two processors
IVect glob_number(n);
if (MPI::COMM_WORLD.Get_rank() == 0)
{
glob_number(0) = 0;
glob_number(1) = 3;
glob_number(2) = 5;
glob_number(3) = 6;
}
else
{
glob_number(0) = 1;
glob_number(1) = 2;
glob_number(2) = 4;
glob_number(3) = 5;
glob_number(4) = 7;
}
// for the constructor, you use the local number of rows
DistributedMatrix<double, General, ArrayRowSparse> Asp(n, n);
DistributedMatrix<double, General, BlockDiagRow> A(n, n);
// arrays overlap_row, original_proc, SharingProc, SharingRows
// are output arrays of Init here, they are constructed with the given global numbers
IVect overlap_row, original_proc;
IVect SharingProc;
Vector<IVect> SharingRows;
// these arrays are constructed by Init of a distributed sparse matrix (Asp here)
Asp.Init(glob_number, overlap_row, original_proc, SharingProc, SharingRows, MPI::COMM_WORLD);
// once these arrays have been constructed, you can use them for other distributed matrices that have the same distribution.
DistributedMatrix<double, General, BlockDiagRow> B(n, n);
B.Init(nglob, &glob_number, &OverlapRow, &original_proc,
n, 1, &SharingProc, &SharingRows, MPI::COMM_WORLD);
// another solution is to call Init with the original distributed matrix
A.Init(Asp);
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
GetProcNumber
Syntax
int* GetProcNumber() const
This method returns a pointer storing the processor number for each "row". Here, we take into account all the rows that are coupled with the rows owned by the current processor.
Example :
// for the constructor, you use the local number of rows
DistributedMatrix<double, General, BlockDiagRow> A(n, n);
// then you need to initialise the arrays used for parallel functions (e.g. matrix-vector product)
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
n, 1, &SharingProc, &SharingRows, MPI::COMM_WORLD);
A.SetPattern(pattern, proc);
int* proc_row = A.GetProcNumber();
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
SetPattern
Syntax
void SetPattern(const Vector<IVect>& NumDof_Blocks); void SetPattern(const Vector<IVect>& NumDof_Blocks, const Vector<IVect>& ProcDofs);
This method specifies the profile of the block-diagonal matrix. In the first syntax, the blocks are not shared between processors, and the user provides the local row numbers for each block. In the second syntax, the blocks may be shared between processors, the second argument gives the processor number for each row of the block. If the processor is distant (different from the current processor), the row number is the global row number, otherwise it is the local row number.
Example :
// for the constructor, you use the local number of rows
DistributedMatrix<double, General, BlockDiagRow> A(n, n);
// then you need to initialise the arrays used for parallel functions (e.g. matrix-vector product)
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
n, 1, &SharingProc, &SharingRows, MPI::COMM_WORLD);
// example with two processors
// processor 0 owns rows 0, 2, 5, 6, 8
// processor 1 owns rows 1, 2, 3, 4, 6, 7
// block 0 : (0, 1, 4)
// block 1 : (2, 3, 6)
// block 2 : (5, 7, 8)
Vector<IVect> pattern(3), proc(3);
for (int k = 0; k < 3; k++)
{
pattern(k).Reallocate(3);
proc(k).Reallocate(3);
}
if (MPI::COMM_WORLD.Get_rank() == 0)
{
// block 0 : (0, 1, 4)
// rows 1 and 4 are owned by processor 1, we put global numbers
pattern(0)(0) = 0; pattern(0)(1) = 1; pattern(0)(2) = 4;
proc(0)(0) = 0; proc(0)(1) = 1; proc(0)(2) = 1;
// block 1 : (2, 3, 6)
// local numbers of rows 2, 6 : 1, 3
pattern(1)(0) = 1; pattern(1)(1) = 3; pattern(1)(2) = 3;
proc(1)(0) = 0; proc(1)(1) = 1; proc(1)(2) = 0;
// block 2 : (5, 7, 8)
// local numbers of rows 5, 8 : 2, 4
pattern(2)(0) = 2; pattern(2)(1) = 7; pattern(2)(2) = 4;
proc(2)(0) = 0; proc(2)(1) = 1; proc(2)(2) = 0;
}
else
{
// block 0 : (0, 1, 4)
// local numbers of rows 1, 4 : 0 3
pattern(0)(0) = 0; pattern(0)(1) = 0; pattern(0)(2) = 3;
proc(0)(0) = 0; proc(0)(1) = 1; proc(0)(2) = 1;
// block 1 : (2, 3, 6)
// local numbers of rows 2, 3, 6 : 1, 2, 4
pattern(1)(0) = 1; pattern(1)(1) = 2; pattern(1)(2) = 4;
proc(1)(0) = 1; proc(1)(1) = 1; proc(1)(2) = 1;
// block 2 : (5, 7, 8)
// local number of row 7 : 5
pattern(2)(0) = 5; pattern(2)(1) = 5; pattern(2)(2) = 8;
proc(2)(0) = 0; proc(2)(1) = 1; proc(2)(2) = 0;
}
A.SetPattern(pattern, proc);
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
IndexGlobal
Syntax
int IndexGlobal(int i, int j)
This method returns the global row number of row j of block i.
Example :
// for the constructor, you use the local number of rows
DistributedMatrix<double, General, BlockDiagRow> A(n, n);
// then you need to initialise the arrays used for parallel functions (e.g. matrix-vector product)
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
n, 1, &SharingProc, &SharingRows, MPI::COMM_WORLD);
// example with two processors
// processor 0 owns rows 0, 2, 5, 6, 8
// processor 1 owns rows 1, 2, 3, 4, 6, 7
// block 0 : (0, 1, 4)
// block 1 : (2, 3, 6)
// block 2 : (5, 7, 8)
Vector<IVect> pattern(3), proc(3);
for (int k = 0; k < 3; k++)
{
pattern(k).Reallocate(3);
proc(k).Reallocate(3);
}
if (MPI::COMM_WORLD.Get_rank() == 0)
{
// block 0 : (0, 1, 4)
// rows 1 and 4 are owned by processor 1, we put global numbers
pattern(0)(0) = 0; pattern(0)(1) = 1; pattern(0)(2) = 4;
proc(0)(0) = 0; proc(0)(1) = 1; proc(0)(2) = 1;
// block 1 : (2, 3, 6)
// local numbers of rows 2, 6 : 1, 3
pattern(1)(0) = 1; pattern(1)(1) = 3; pattern(1)(2) = 3;
proc(1)(0) = 0; proc(1)(1) = 1; proc(1)(2) = 0;
// block 2 : (5, 7, 8)
// local numbers of rows 5, 8 : 2, 4
pattern(2)(0) = 2; pattern(2)(1) = 7; pattern(2)(2) = 4;
proc(2)(0) = 0; proc(2)(1) = 1; proc(2)(2) = 0;
}
else
{
// block 0 : (0, 1, 4)
// local numbers of rows 1, 4 : 0 3
pattern(0)(0) = 0; pattern(0)(1) = 0; pattern(0)(2) = 3;
proc(0)(0) = 0; proc(0)(1) = 1; proc(0)(2) = 1;
// block 1 : (2, 3, 6)
// local numbers of rows 2, 3, 6 : 1, 2, 4
pattern(1)(0) = 1; pattern(1)(1) = 2; pattern(1)(2) = 4;
proc(1)(0) = 1; proc(1)(1) = 1; proc(1)(2) = 1;
// block 2 : (5, 7, 8)
// local number of row 7 : 5
pattern(2)(0) = 5; pattern(2)(1) = 5; pattern(2)(2) = 8;
proc(2)(0) = 0; proc(2)(1) = 1; proc(2)(2) = 0;
}
A.SetPattern(pattern, proc);
if (MPI::COMM_WORLD.Get_rank() == 0)
{
// it should return 6 (global row number of row 2 of block 1)
cout << A.IndexGlobal(1, 2) << endl;
}
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
Assemble
Syntax
void Assemble()
This method assembles the block-diagonal matrix (values shared between processors are added).
Example :
// for the constructor, you use the local number of rows
DistributedMatrix<double, General, BlockDiagRow> A(n, n);
// then you need to initialise the arrays used for parallel functions (e.g. matrix-vector product)
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
n, 1, &SharingProc, &SharingRows, MPI::COMM_WORLD);
// example with two processors
// processor 0 owns rows 0, 2, 5, 6, 8
// processor 1 owns rows 1, 2, 3, 4, 6, 7
// block 0 : (0, 1, 4)
// block 1 : (2, 3, 6)
// block 2 : (5, 7, 8)
Vector<IVect> pattern(3), proc(3);
for (int k = 0; k < 3; k++)
{
pattern(k).Reallocate(3);
proc(k).Reallocate(3);
}
if (MPI::COMM_WORLD.Get_rank() == 0)
{
// block 0 : (0, 1, 4)
// rows 1 and 4 are owned by processor 1, we put global numbers
pattern(0)(0) = 0; pattern(0)(1) = 1; pattern(0)(2) = 4;
proc(0)(0) = 0; proc(0)(1) = 1; proc(0)(2) = 1;
// block 1 : (2, 3, 6)
// local numbers of rows 2, 6 : 1, 3
pattern(1)(0) = 1; pattern(1)(1) = 3; pattern(1)(2) = 3;
proc(1)(0) = 0; proc(1)(1) = 1; proc(1)(2) = 0;
// block 2 : (5, 7, 8)
// local numbers of rows 5, 8 : 2, 4
pattern(2)(0) = 2; pattern(2)(1) = 7; pattern(2)(2) = 4;
proc(2)(0) = 0; proc(2)(1) = 1; proc(2)(2) = 0;
}
else
{
// block 0 : (0, 1, 4)
// local numbers of rows 1, 4 : 0 3
pattern(0)(0) = 0; pattern(0)(1) = 0; pattern(0)(2) = 3;
proc(0)(0) = 0; proc(0)(1) = 1; proc(0)(2) = 1;
// block 1 : (2, 3, 6)
// local numbers of rows 2, 3, 6 : 1, 2, 4
pattern(1)(0) = 1; pattern(1)(1) = 2; pattern(1)(2) = 4;
proc(1)(0) = 1; proc(1)(1) = 1; proc(1)(2) = 1;
// block 2 : (5, 7, 8)
// local number of row 7 : 5
pattern(2)(0) = 5; pattern(2)(1) = 5; pattern(2)(2) = 8;
proc(2)(0) = 0; proc(2)(1) = 1; proc(2)(2) = 0;
}
A.SetPattern(pattern, proc);
// values can be added to A with AddInteraction, AddDistantInteraction or AddRowDistantInteraction
A.AddInteraction(0, 1, 2.3);
// then the matrix can be assembled
A.Assemble();
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
Mlt
Syntax
void Mlt(const DistributedMatrix<T>& A, const Vector<T>& x, Vector<T>& y, bool assemble); void Mlt(const SeldonTranspose& trans, const DistributedMatrix<T>& A, const Vector<T>& x, Vector<T>& y, bool assemble); void Mlt(const T& alpha, DistributedMatrix<T>& A); void Mlt(const DistributedMatrix<T>& A, Vector<T>& x, bool assemble);
This function can perform matrix-vector product (with the original matrix or with its transpose) and can be used to multiply the matrix by a scalar. The matrix-matrix product is currently not implemented. If the communicator contains only one processor, the sequential function Mlt (with the class Matrix) will be called, otherwise an error will be displayed during the execution.
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// global row numbers are provided with Init (glob_number)
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
Nvol, nb_u, &ProcShared, &SharedRows, MPI::COMM_WORLD);
// SetPattern is called to provide the profile of A
A.SetPattern(NumDofs, ProcDofs);
// for local interaction, you can use AddInteraction
A.AddInteraction(i, j, val);
// when the column is distant, you can use AddDistantInteraction
A.AddDistantInteraction(i, jglob, proc, val);
// when the row is distant, you can use AddRowDistantInteraction
A.AddRowDistantInteraction(iglob, j, proc, val);
// once the matrix is constructed, matrix-vector products can be performed
// each vector is stored as a usual Seldon vector (class Vector)
// each processor is assumed to store the values of the vector for the
// global rows given during the construction of the matrix (here, the array glob_number)
Vector<double> x(n), y(n), z;
// Mlt will assume that the values of x for rows that are shared between processors are the same
// which is here not the case when calling FillRand
x.FillRand();
y.FillRand();
z = y;
// a solution is to call AssembleVector to ensure this property
AssembleVector(x, MPI::SUM, ProcShared, SharedRows, MPI::COMM_WORLD, Nvol, nb_u, 23);
// classical matrix-vector product y = A x
// the default value of the last argument (assemble) is true
// the result y is assembled and will contain the same result as a sequential matrix-vector product,
// the values being distributed in the different processors
Mlt(A, x, y);
// assemble = true, induces a final call to the function AssembleVector (to add values of shared rows)
// if the user wants to perform a matrix-vector product without performing this final assembling step, he puts false in the last optional argument
Mlt(A, x, z, false);
// and perform the assembling phase later
AssembleVector(z, MPI::SUM, ProcShared, SharedRows, MPI::COMM_WORLD, Nvol, nb_u, 23);
// Mlt can be used to perform a matrix-vector with the transpose matrix or the transpose conjugate
Mlt(SeldonTrans, A, x, y);
// Finally Mlt can be used to multiply a distributed matrix by a scalar coefficient
Mlt(-2.1, A);
// the last syntax performs y = A x, on input x is provided and overwritten with y
y = x;
Mlt(A, y);
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
ConvertToSparse
Syntax
void ConvertToSparse(const DistributedMatrix<T>& A, DistributedMatrix<T>& B);
This function converts a block-diagonal matrix to a sparse distributed matrix. It is equivalent to call the function Copy (which handles other conversions as well).
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// global row numbers are provided with Init (glob_number)
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
Nvol, nb_u, &ProcShared, &SharedRows, MPI::COMM_WORLD);
// SetPattern is called to provide the profile of A
A.SetPattern(NumDofs, ProcDofs);
// for local interaction, you can use AddInteraction
A.AddInteraction(i, j, val);
// when the column is distant, you can use AddDistantInteraction
A.AddDistantInteraction(i, jglob, proc, val);
// when the row is distant, you can use AddRowDistantInteraction
A.AddRowDistantInteraction(iglob, j, proc, val);
// once the matrix is constructed, the matrix can be converted
// with Copy or ConvertToSparse
DistributedMatrix<double, General, ArrayRowSparse> B;
ConvertToSparse(A, B);
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
ConvertToBlockDiagonal
Syntax
void ConvertToBlockDiagonal(const DistributedMatrix<T>& A, DistributedMatrix<T>& B);
This function converts a sparse distributed matrix to a block-diagonal matrix. It is equivalent to call the function Copy (which handles other conversions as well).
Example :
// default constructor
DistributedMatrix<double, General, ArrayRowSparse> A;
// global row numbers are provided with Init (glob_number)
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
Nvol, nb_u, &ProcShared, &SharedRows, MPI::COMM_WORLD);
// A is allocated with the loca number of rows
A.Reallocate(n, n);
// for local interaction, you can use AddInteraction
A.AddInteraction(i, j, val);
// when the column is distant, you can use AddDistantInteraction
A.AddDistantInteraction(i, jglob, proc, val);
// when the row is distant, you can use AddRowDistantInteraction
A.AddRowDistantInteraction(iglob, j, proc, val);
// once the matrix is constructed, the matrix can be converted
// with Copy or ConvertToBlockDiagonal
DistributedMatrix<double, General, BlockDiagRow> B;
ConvertToBlockDiagonal(A, B);
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
GetInverse
Syntax
void GetInverse(DistributedMatrix<T>& B);
This function replaces a block-diagonal matrix by its inverse.
Example :
// default constructor
DistributedMatrix<double, General, BlockDiagRow> A;
// global row numbers are provided with Init (glob_number)
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
Nvol, nb_u, &ProcShared, &SharedRows, MPI::COMM_WORLD);
// SetPattern is called to provide the profile of A
A.SetPattern(NumDofs, ProcDofs);
// for local interaction, you can use AddInteraction
A.AddInteraction(i, j, val);
// when the column is distant, you can use AddDistantInteraction
A.AddDistantInteraction(i, jglob, proc, val);
// when the row is distant, you can use AddRowDistantInteraction
A.AddRowDistantInteraction(iglob, j, proc, val);
// once the matrix is constructed, the inverse can be computed
GetInverse(A);
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
GetCholesky
Syntax
void GetCholesky(DistributedMatrix<T>& B);
This function replaces a block-diagonal matrix by its Cholesky factorization. It is only available for BlockDiagRowSym storage, since the matrix is assumed to be symmetric definite positive.
Example :
// default constructor
DistributedMatrix<double, Symmetric, BlockDiagRowSym> A;
// global row numbers are provided with Init (glob_number)
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
Nvol, nb_u, &ProcShared, &SharedRows, MPI::COMM_WORLD);
// SetPattern is called to provide the profile of A
A.SetPattern(NumDofs, ProcDofs);
// for local interaction, you can use AddInteraction
A.AddInteraction(i, j, val);
// when the column is distant, you can use AddDistantInteraction
A.AddDistantInteraction(i, jglob, proc, val);
// when the row is distant, you can use AddRowDistantInteraction
A.AddRowDistantInteraction(iglob, j, proc, val);
// once the matrix is constructed, the Cholesky factorization can be computed
GetCholesky(A);
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
SolveCholesky
Syntax
void SolveCholesky(const SeldonTranspose&, const DistributedMatrix<T>& A, Vector<T>& x);
This function solves system L x = b or L^T x = b, where A = L L^T (Cholesky factorization). GetCholesky is assumed to have been called previously such that A stores the Cholesky factor L.
Example :
// default constructor
DistributedMatrix<double, Symmetric, BlockDiagRowSym> A;
// global row numbers are provided with Init (glob_number)
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
Nvol, nb_u, &ProcShared, &SharedRows, MPI::COMM_WORLD);
// SetPattern is called to provide the profile of A
A.SetPattern(NumDofs, ProcDofs);
// for local interaction, you can use AddInteraction
A.AddInteraction(i, j, val);
// when the column is distant, you can use AddDistantInteraction
A.AddDistantInteraction(i, jglob, proc, val);
// when the row is distant, you can use AddRowDistantInteraction
A.AddRowDistantInteraction(iglob, j, proc, val);
// once the matrix is constructed, the Cholesky factorization can be computed
GetCholesky(A);
// and used to solve L L^T x = b
Vector<T> x, b(n);
b.FillRand();
x = b;
SolveCholesky(SeldonNoTrans, A, x);
SolveCholesky(SeldonTrans, A, x);
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx
MltCholesky
Syntax
void MltCholesky(const SeldonTranspose&, const DistributedMatrix<T>& A, Vector<T>& x);
This function computes matrix-vector product y = L x or y = L^T x, where A = L L^T (Cholesky factorization). GetCholesky is assumed to have been called previously such that A stores the Cholesky factor L.
Example :
// default constructor
DistributedMatrix<double, Symmetric, BlockDiagRowSym> A;
// global row numbers are provided with Init (glob_number)
A.Init(nglob, &glob_number, &OverlapRow, &original_proc,
Nvol, nb_u, &ProcShared, &SharedRows, MPI::COMM_WORLD);
// SetPattern is called to provide the profile of A
A.SetPattern(NumDofs, ProcDofs);
// for local interaction, you can use AddInteraction
A.AddInteraction(i, j, val);
// when the column is distant, you can use AddDistantInteraction
A.AddDistantInteraction(i, jglob, proc, val);
// when the row is distant, you can use AddRowDistantInteraction
A.AddRowDistantInteraction(iglob, j, proc, val);
// once the matrix is constructed, the Cholesky factorization can be computed
GetCholesky(A);
// and used to compute y = L L^T x
Vector<T> y, x(n);
x.FillRand();
y = x;
MltCholesky(SeldonTrans, A, y);
MltCholesky(SeldonNoTrans, A, y);
Location :
DistributedBlockDiagonalMatrix.hxx
DistributedBlockDiagonalMatrix.cxx