Thank you, Matt. It took me a while but I understand how to work with
submatrices.
(Below is simple stand-alone code which assembles two submatrices into a matrix
of 5 x 5, which might be of use to someone else)
Q1) Looks like I don’t need to call MatrixAssembleBegin/End of on submatrices?
Q2) Though I have preallocated A for 5 x 5 – why do I still need to call
MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
else I get an error while assembling using MatSetValuesLocal()?
Q3) Am I calling MatRestoreLocalSubmatrix() at the right place and what does it
do?
Q4) When assembling a large submatrix, can I call MatGetOwershipRange on the
submatrix?
Best,
Karthik.
#include <petsc.h>
Mat A; /* linear system matrix */
PetscInt i,j,m = 3,n = 3, P = 5, Q = 5;
PetscErrorCode ierr;
PetscScalar v;
// create matrix A
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetType(A, MATMPIAIJ); CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,P,Q);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);CHKERRQ(ierr);
//ierr = MatSeqAIJSetPreallocation(A,5,NULL);CHKERRQ(ierr);
// local indices is always 0 1 2
PetscInt c = 3, global_col_indices[] = {0, 1, 2};
PetscInt r = 2, global_row_indices[] = {3, 4};
ISLocalToGlobalMapping col_mapping, row_mapping;
ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, c, global_col_indices,
PETSC_COPY_VALUES, &col_mapping);
ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, r, global_row_indices,
PETSC_COPY_VALUES, &row_mapping);
MatSetLocalToGlobalMapping(A,row_mapping, col_mapping);
for (i=0; i<m; i++) {
for (j=0; j<n; j++) {
v = i*n + j + 1;
std::cout << "inserting value " << v << " into position " << i << ",
" << j << std::endl;
MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES); CHKERRQ(ierr);
}
}
// assemble matrix
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
// Test MatGetLocalSubMatrix
Mat SubA, SubB;
IS row_isx, col_isx;
PetscInt row_isx_vals[2];
row_isx_vals[0] = 0;
row_isx_vals[1] = 1;
PetscInt col_isx_vals[3];
col_isx_vals[0] = 0;
col_isx_vals[1] = 1;
col_isx_vals[2] = 2;
CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD, 3, col_isx_vals,
PETSC_COPY_VALUES, &col_isx));
CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD, 2, row_isx_vals,
PETSC_COPY_VALUES, &row_isx));
CHKERRQ(MatGetLocalSubMatrix(A, row_isx, col_isx, &SubA));
// MatSetLocalToGlobalMapping(A,row_mapping, row_mapping);
MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
int irow = 0, icol = 0;
PetscScalar val = 10;
MatSetValuesLocal(SubA,1,&irow,1,&icol,&val,INSERT_VALUES);
// ierr = MatAssemblyBegin(SubA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
// ierr = MatAssemblyEnd(SubA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
MatSetLocalToGlobalMapping(A,col_mapping, row_mapping);
CHKERRQ(MatGetLocalSubMatrix(A, col_isx, row_isx, &SubB));
irow = 0, icol = 0;
MatSetValuesLocal(SubB,1,&irow,1,&icol,&val,INSERT_VALUES);
irow = 2, icol = 0;
MatSetValuesLocal(SubB,1,&irow,1,&icol,&val,INSERT_VALUES);
irow = 2, icol = 1;
MatSetValuesLocal(SubB,1,&irow,1,&icol,&val,INSERT_VALUES);
// ierr = MatAssemblyBegin(SubB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
// ierr = MatAssemblyEnd(SubB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
MatRestoreLocalSubMatrix(A, row_isx, col_isx, &SubA);
MatRestoreLocalSubMatrix(A, col_isx, row_isx, &SubB);
// cleanup
ierr = MatDestroy(&A); CHKERRQ(ierr);
ierr = PetscFinalize();
From: Matthew Knepley <[email protected]>
Date: Thursday, 20 April 2023 at 11:37
To: Chockalingam, Karthikeyan (STFC,DL,HC) <[email protected]>
Cc: [email protected] <[email protected]>
Subject: Re: [petsc-users] question about MatSetLocalToGlobalMapping
On Thu, Apr 20, 2023 at 6:13 AM Karthikeyan Chockalingam - STFC UKRI via
petsc-users <[email protected]<mailto:[email protected]>> wrote:
Hello,
I created a new thread, thought would it be more appropriate (and is a
continuation of my previous post). I want to construct the below K matrix
(which is composed of submatrices)
K = [A P^T
P 0]
Where K is of type MatMPIAIJ. I first constructed the top left [A] using
MatSetValues().
Now, I would like to construct the bottom left [p] and top right [p^T] using
MatSetValuesLocal().
To use MatSetValuesLocal(), I first have to create a local-to-global mapping
using ISLocalToGlobalMappingCreate. I have created two mapping row_mapping and
column_mapping.
I do not understand why they are not the same map. Maybe I was unclear before.
It looks like you have two fields, say phi and lambda, where lambda is a
Lagrange multiplier imposing some constraint. Then you get a saddle point like
this. You can imagine matrices
(phi, phi) --> A
(phi, lambda) --> P^T
(lambda, phi) --> P
So you make a L2G map for the phi field and the lambda field. Oh, you are
calling them row and col map, but they are my phi and lambda
maps. I do not like the row and col names since in P they reverse.
Q1) At what point should I declare MatSetLocalToGlobalMapping – is it just
before I use MatSetValuesLocal()?
Okay, it is good you are asking this because my thinking was somewhat confused.
I think the precise steps are:
1) Create the large saddle point matrix K
1a) We must call
https://petsc.org/main/manualpages/Mat/MatSetLocalToGlobalMapping/ on it. In
the simplest case, this just maps
the local rows numbers [0, Nrows) to the global rows numbers
[rowStart, rowStart + Nrows).
2) To form each piece:
2a) Extract that block using
https://petsc.org/main/manualpages/Mat/MatGetLocalSubMatrix/
This gives back a Mat object that you subsequently restore using
https://petsc.org/main/manualpages/Mat/MatRestoreLocalSubMatrix/
2b) Insert values using
https://petsc.org/main/manualpages/Mat/MatSetValuesLocal/
The local indices used for insertion here are indices relative to
the block itself, and the L2G map for this matrix
has been rewritten to insert into that block in the larger matrix.
Thus this looks like just calling MatSetValuesLocal()
on the smaller matrix block, but inserts correctly into the larger
matrix.
Therefore, the code you write code in 2) could work equally well making the
large matrix from 1), or independent smaller matrix blocks.
Does this make sense?
Thanks,
Matt
I will use MatSetLocalToGlobalMapping(K, row_mapping, column_mapping) to build
the bottom left [P].
Q2) Can now I reset the mapping as MatSetLocalToGlobalMapping(K,
column_mapping, row_mapping) to build the top right [P^T]?
Many thanks!
Kind regards,
Karthik.
--
What most experimenters take for granted before they begin their experiments is
infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>
