[petsc-users] Setting up a matrix for Lagrange multiplier
Karthikeyan Chockalingam - STFC UKRI
karthikeyan.chockalingam at stfc.ac.uk
Tue Apr 25 04:38:51 CDT 2023
Thank you for your response. I will come back to the error later (for now I have it working with -pc_fieldsplit_detect_saddle_point.)
I would like to understand how I need to index while running in parallel.
(Q1) Index set (IS) for the ISLocalToGlobalMapping.
for (int i = 0; i < N; i++)
vec_col_index[i] = i;
for (int i = 0; i < C; i++)
vec_shifted_row_index[i] = i + N;
ISLocalToGlobalMapping phi_mapping, lamda_mapping;
ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, C, vec_shifted_row_index, PETSC_COPY_VALUES, &lamda_mapping);
ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, N, vec_col_index, PETSC_COPY_VALUES, &phi_mapping);
(i) I believe the above index set is created by all parallel processes – why then I would have to offset it with the first index of this process?
Since the mapping will always hold true and we would only use a subset of this index set while populating the submatrices. I can see that I might have an over-allocated index set than necessary but I can’t see why it can be wrong.
(ii) What would be the first index of this process anyway, since at this point I haven’t declared any parallel sub-matrix or sub-vector?
(Q2) Index set (IS) while creating the submatrix
for (int i = 0; i < C; i++)
vec_row_index[i] = i;
Mat SubP;
IS row_isx, col_isx;
ISCreateGeneral(PETSC_COMM_WORLD, C, vec_row_index, PETSC_COPY_VALUES, &row_isx);
ISCreateGeneral(PETSC_COMM_WORLD, N, vec_col_index, PETSC_COPY_VALUES, &col_isx);
MatSetOption(A[level], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
MatSetLocalToGlobalMapping(A, lamda_mapping, phi_mapping);
MatGetLocalSubMatrix(A, row_isx, col_isx, &SubP);
I believe I will have to request a submatrix SubP of global size C x N on all processes – right? Or do you think I can still break row_isx and col_isx into parallel parts?
(Q3) Populating the submatrix
const PetscInt *nindices;
PetscInt count;
ISGetIndices(isout[level], &nindices);
//isout is the global index of the nodes to be constrained.
PetscScalar val = 1;
for (int irow=0; irow < n; irow++)
{
int icol = nindices[irow];
MatSetValuesLocal(SubP,1,&irow,1,&icol,&val,ADD_VALUES);
}
I do think if I have the MatGetOwnerShipRange of SubP; I will be able to parallelize the loop. Is MatGetOwnerShipRange also available for submatrices as well?
If I have to parallelize the index sets maybe I have to reorder the above steps, can you please make a suggestion?
Kind regards,
Karthik.
From: Matthew Knepley <knepley at gmail.com>
Date: Monday, 24 April 2023 at 19:24
To: Chockalingam, Karthikeyan (STFC,DL,HC) <karthikeyan.chockalingam at stfc.ac.uk>
Cc: petsc-users at mcs.anl.gov <petsc-users at mcs.anl.gov>
Subject: Re: [petsc-users] Setting up a matrix for Lagrange multiplier
On Mon, Apr 24, 2023 at 1:58 PM Karthikeyan Chockalingam - STFC UKRI <karthikeyan.chockalingam at stfc.ac.uk<mailto:karthikeyan.chockalingam at stfc.ac.uk>> wrote:
Changing the names to NULL produced the same error.
This is just to make the code simpler.
Don’t I have to give the fields a name?
They get the default name, which is the numbering you used.
The problem solution didn’t change from one time step to the next.
[1;31m[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0;39m[0;49m[0]PETSC ERROR: No support for this operation for this object type
[0]PETSC ERROR: Unsupported viewer gmres
You are not completely catching the error, because it should print out the entire options database on termination.
Here it says you gave "gmres" as a Viewer, which is incorrect, but I cannot see all the options you used.
Thanks,
Matt
[0]PETSC ERROR: WARNING! There are option(s) set that were not used! Could be the program crashed before they were used or a spelling mistake, etc!
[0]PETSC ERROR: Option left: name:-fieldsplit_1_ksp_type value: preonly source: code
[0]PETSC ERROR: Option left: name:-fieldsplit_1_ksp_view value: gmres source: code
[0]PETSC ERROR: Option left: name:-fieldsplit_1_pc_type value: lu source: code
[0]PETSC ERROR: Option left: name:-options_left (no value) source: code
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.18.4-529-g995ec06f92 GIT Date: 2023-02-03 18:41:48 +0000
[0]PETSC ERROR: /Users/karthikeyan.chockalingam/AMReX/amrFEM/build/Debug/amrFEM on a named HC20210312 by karthikeyan.chockalingam Mon Apr 24 18:51:25 2023
[0]PETSC ERROR: Configure options --with-debugging=0 --prefix=/Users/karthikeyan.chockalingam/AMReX/petsc --download-fblaslapack=yes --download-scalapack=yes --download-mumps=yes --with-hypre-dir=/Users/karthikeyan.chockalingam/AMReX/hypre/src/hypre
[0]PETSC ERROR: #1 PetscOptionsGetViewer() at /Users/karthikeyan.chockalingam/AMReX/SRC_PKG/petsc/src/sys/classes/viewer/interface/viewreg.c:309
[0]PETSC ERROR: #2 KSPSetFromOptions() at /Users/karthikeyan.chockalingam/AMReX/SRC_PKG/petsc/src/ksp/ksp/interface/itcl.c:522
[0]PETSC ERROR: #3 PCSetUp_FieldSplit() at /Users/karthikeyan.chockalingam/AMReX/SRC_PKG/petsc/src/ksp/pc/impls/fieldsplit/fieldsplit.c:1054
[0]PETSC ERROR: #4 PCSetUp() at /Users/karthikeyan.chockalingam/AMReX/SRC_PKG/petsc/src/ksp/pc/interface/precon.c:994
[0]PETSC ERROR: #5 KSPSetUp() at /Users/karthikeyan.chockalingam/AMReX/SRC_PKG/petsc/src/ksp/ksp/interface/itfunc.c:405
[0]PETSC ERROR: #6 KSPSolve_Private() at /Users/karthikeyan.chockalingam/AMReX/SRC_PKG/petsc/src/ksp/ksp/interface/itfunc.c:824
[0]PETSC ERROR: #7 KSPSolve() at /Users/karthikeyan.chockalingam/AMReX/SRC_PKG/petsc/src/ksp/ksp/interface/itfunc.c:1070
End of program
solve time 0.03416619 seconds
It will be a bit difficult for me to produce a stand-alone code.
Thank you for your help.
Best,
Karthik.
From: Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>>
Date: Monday, 24 April 2023 at 18:42
To: Chockalingam, Karthikeyan (STFC,DL,HC) <karthikeyan.chockalingam at stfc.ac.uk<mailto:karthikeyan.chockalingam at stfc.ac.uk>>
Cc: petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov> <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>>
Subject: Re: [petsc-users] Setting up a matrix for Lagrange multiplier
On Mon, Apr 24, 2023 at 1:37 PM Karthikeyan Chockalingam - STFC UKRI <karthikeyan.chockalingam at stfc.ac.uk<mailto:karthikeyan.chockalingam at stfc.ac.uk>> wrote:
Great to know there is an example now.
Yes, -pc_fieldsplit_defect_saddle_point worked.
But I would like to understand what I did wrong (and why it didn’t work).
Can you show the error? Or give something self-contained that I can run.
After reading many posted posts, I needed to create PCFieldSplitSetIS to split the fields.
Yes. I would give NULL for the name here.
I set the first N indices to the field phi and the next C indices (starting from N) to the field lambda.
This looks fine. However, these are global indices, so this would not work in parallel. You would have to offset by
the first index on this process.
Thanks,
Matt
Please tell me what I did wrong and how I can fix the lines which are now commented?
PetscErrorCode ierr;
KSP ksp;
PC pc;
KSPCreate(PETSC_COMM_WORLD, &ksp);
KSPSetType(ksp, KSPGMRES);
KSPSetOperators(ksp, A[level], A[level]);
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCFIELDSPLIT);CHKERRQ(ierr);
……..
for (int i = 0; i < N; i++)
vec_zero_field_index[i] = i;
PetscInt * vec_first_field_index = NULL;
PetscMalloc(n * sizeof(PetscInt), &vec_first_field_index);
for (int i = 0; i < C; i++)
vec_first_field_index[i] = N + i;
IS zero_field_isx, first_field_isx;
CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD, N, vec_zero_field_index, PETSC_COPY_VALUES, &zero_field_isx));
CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD, C, vec_first_field_index, PETSC_COPY_VALUES, &first_field_isx));
ierr = PCFieldSplitSetIS(pc,"0",zero_field_isx); CHKERRQ(ierr);
ierr = PCFieldSplitSetIS(pc,"1",first_field_isx); CHKERRQ(ierr);
ierr = PetscOptionsSetValue(NULL,"-ksp_type", "fgmres"); CHKERRQ(ierr);
ierr = PetscOptionsSetValue(NULL,"-pc_type", "fieldsplit"); CHKERRQ(ierr);
/*
ierr = PetscOptionsSetValue(NULL,"-fieldsplit_0_ksp_type", "gmres"); CHKERRQ(ierr);
ierr = PetscOptionsSetValue(NULL,"-fieldsplit_0_pc_type", "jacobi"); CHKERRQ(ierr);
ierr = PetscOptionsSetValue(NULL,"-fieldsplit_1_ksp_type", "preonly"); CHKERRQ(ierr);
ierr = PetscOptionsSetValue(NULL,"-fieldsplit_1_pc_type", "lu"); CHKERRQ(ierr);*/
ierr = PetscOptionsSetValue(NULL, "-pc_fieldsplit_detect_saddle_point", NULL);CHKERRQ(ierr);
ierr = PetscOptionsSetValue(NULL, "-options_left", NULL);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
KSPSolve(ksp, b, x);
Best,
Karthik.
From: Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>>
Date: Monday, 24 April 2023 at 17:57
To: Chockalingam, Karthikeyan (STFC,DL,HC) <karthikeyan.chockalingam at stfc.ac.uk<mailto:karthikeyan.chockalingam at stfc.ac.uk>>
Cc: petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov> <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>>
Subject: Re: [petsc-users] Setting up a matrix for Lagrange multiplier
On Mon, Apr 24, 2023 at 10:22 AM Karthikeyan Chockalingam - STFC UKRI <karthikeyan.chockalingam at stfc.ac.uk<mailto:karthikeyan.chockalingam at stfc.ac.uk>> wrote:
Hello,
I was able to construct the below K matrix (using submatrices P and P^T), which is of type MATAIJ
K = [A P^T
P 0]
and solved them using a direct solver.
I modified your example to create either AIJ or Nest matrices, but use the same assembly code:
https://gitlab.com/petsc/petsc/-/merge_requests/6368
However, I was reading online that this is a saddle point problem and I should be employing PCFIELDSPLIT.
Since I have one monolithic matrix K, I was not sure how to split the fields.
With this particular matrix, you can use
-pc_fieldsplit_detect_saddle_point
and it will split it automatically.
Thanks,
Matt
Best regards,
Karthik.
From: Chockalingam, Karthikeyan (STFC,DL,HC) <karthikeyan.chockalingam at stfc.ac.uk<mailto:karthikeyan.chockalingam at stfc.ac.uk>>
Date: Wednesday, 19 April 2023 at 17:52
To: Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>>
Cc: petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov> <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>>
Subject: Re: [petsc-users] Setting up a matrix for Lagrange multiplier
I have declared the mapping
ISLocalToGlobalMapping mapping;
ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, n, nindices, PETSC_COPY_VALUES, &mapping);
But when I use MatSetValuesLocal(), how do I know the above mapping is employed because it is not one of the parameters passed to the function?
Thank you.
Kind regards,
Karthik.
From: Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>>
Date: Tuesday, 18 April 2023 at 16:21
To: Chockalingam, Karthikeyan (STFC,DL,HC) <karthikeyan.chockalingam at stfc.ac.uk<mailto:karthikeyan.chockalingam at stfc.ac.uk>>
Cc: petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov> <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>>
Subject: Re: [petsc-users] Setting up a matrix for Lagrange multiplier
On Tue, Apr 18, 2023 at 11:16 AM Karthikeyan Chockalingam - STFC UKRI <karthikeyan.chockalingam at stfc.ac.uk<mailto:karthikeyan.chockalingam at stfc.ac.uk>> wrote:
Thank you for your response. I spend some time understanding how
MatSetValuesLocal and ISLocalToGlobalMappingCreate work.
You can look at SNES ex28 where we do this with DMCOMPOSITE.
Q1) Will the matrix K be of type MATMPIAIJ or MATIS?
K = [A P^T
P 0]
I assume MPIAIJ since IS is only used for Neumann-Neumann decompositions.
Q2) Can I use both MatSetValues() to MatSetValuesLocal() to populate K? Since I have already used MatSetValues() to construct A.
You can, and there would be no changes in serial if K is exactly the upper left block, but in parallel global indices would change.
Q3) What are the advantages of using MatSetValuesLocal()? Is it that I can construct P directly using local indies and map the entrees to the global index in K?
You have a monolithic K, so that you can use sparse direct solvers to check things. THis is impossible with separate storage.
Q4) I probably don’t have to construct an independent P matrix
You wouldn't in this case.
Thanks,
Matt
Best regards,
Karthik.
From: Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>>
Date: Tuesday, 18 April 2023 at 11:08
To: Chockalingam, Karthikeyan (STFC,DL,HC) <karthikeyan.chockalingam at stfc.ac.uk<mailto:karthikeyan.chockalingam at stfc.ac.uk>>
Cc: petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov> <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>>
Subject: Re: [petsc-users] Setting up a matrix for Lagrange multiplier
On Tue, Apr 18, 2023 at 5:24 AM Karthikeyan Chockalingam - STFC UKRI via petsc-users <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>> wrote:
Hello,
I'm solving a problem using the Lagrange multiplier, the matrix has the form
K = [A P^T
P 0]
I am familiar with constructing K using MATMPIAIJ. However, I would like to know if had [A], can I augment it with [P], [P^T] and [0] of type MATMPIAIJ? Likewise for vectors as well.
Can you please point me to the right resource, if it is a common operation in PETSc?
You can do this at least 2 ways:
1) Assemble you submatrices directly into the larger matrix by constructing local-to-global maps for the emplacement. so that you do
not change your assembly code, except to change MatSetValues() to MatSetValuesLocal(). This is usually preferable.
2) Use MATNEST and VecNEST to put pointers to submatrices and subvectors directly in.
Thanks,
Matt
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/>
--
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/>
--
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/>
--
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/>
--
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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230425/3ded7cd7/attachment-0001.html>
More information about the petsc-users
mailing list