Re: [petsc-dev] PDIPDM questions

[Abhyankar, Shrirang G via petsc-dev]

From: Pierre Jolivet <pierre.joli...@enseeiht.fr>
Date: Wednesday, September 16, 2020 at 12:59 AM
To: Barry Smith <bsm...@petsc.dev>
Cc: "Abhyankar, Shrirang G" <shrirang.abhyan...@pnnl.gov>, PETSc Development 
<petsc-dev@mcs.anl.gov>
Subject: Re: [petsc-dev] PDIPDM questions


On 15 Sep 2020, at 11:18 PM, Barry Smith 
<bsm...@petsc.dev<mailto:bsm...@petsc.dev>> wrote:


  I see now, the TaoSetup_PDIPM() code explicitly builds the large matrix by 
calling MatGetRow() on the smaller matrices, no matrix abstractions at all and 
assuming a global AIJ storage of the large matrix.

  I think this can be made more general by building a MatNest() from the 
smaller matrices which allows the smaller matrices to have whatever unique 
structure is appropriate for them. Of course, if needed, for example, for a 
direct solver the MatNest can convert itself to a global AIJ matrix and get the 
current result.

Before jumping into the PDIPM train, I thought a MatNest was used internally. 
Like you can solve [[A, b],[b^T, 0]] (b of dimension M x 1) with a MatNest + 
fieldsplit or just MatConvert() + LU as you mentioned.

The current implementation of PDIPM forms a monolithic KKT matrix in AIJ format 
and then passes that to the linear solver. We did try using fieldsplit by 
setting ISes for primal and dual variables. But, the performance of the solver 
using Fieldsplit was not that great. We’ll probably try it again after the 
current development of PDIPM.


   As Pierre rightly points out, MPIAIJ is just not a good format for many 
constrained optimization problems. This has been something seriously neglected 
in PETSc/TAO. I don't think any single format is ideal for such problems, 
likely hybrid formats are needed and this gets super tricky when one needs to 
use direct solvers such a Cholesky. If MUMPS has non-row oriented storage 
formats we could leverage that but I don't know if it does.

Just FYI, MUMPS is the most general (direct solver) when it comes to input 
storage. It handles non-sorted coordinate format, possibly with duplicated 
entries which are then summed among processes (so, in theory, it could be used 
to “invert" a MatIS/MatNest, without having to first convert the MatIS/MatNest 
to MatAIJ).

Thanks for the pointer on MUMPS’s non-sorted coordinate format.

Thanks,
Pierre

Thanks,
Shri


  Barry



  ierr = MatGetOwnershipRange(tao->hessian,&rjstart,NULL);CHKERRQ(ierr);
  for (i=0; i<pdipm->nx; i++){
    row = rstart + i;

    ierr = MatGetRow(tao->hessian,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
    proc = 0;
    for (j=0; j < nc; j++) {
      while (aj[j] >= cranges[proc+1]) proc++;
      col = aj[j] - cranges[proc] + Jranges[proc];
      ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
    }
    ierr = MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);

    if (pdipm->ng) {
      /* Insert grad g' */
      ierr = MatGetRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
      ierr = 
MatGetOwnershipRanges(tao->jacobian_equality,&ranges);CHKERRQ(ierr);
      proc = 0;
      for (j=0; j < nc; j++) {
        /* find row ownership of */
        while (aj[j] >= ranges[proc+1]) proc++;
        nx_all = rranges[proc+1] - rranges[proc];
        col = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
        ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
      }
      ierr = 
MatRestoreRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
    }

    /* Insert Jce_xfixed^T' */
    if (pdipm->nxfixed) {
      ierr = MatGetRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
      ierr = MatGetOwnershipRanges(pdipm->Jce_xfixed,&ranges);CHKERRQ(ierr);
      proc = 0;
      for (j=0; j < nc; j++) {
        /* find row ownership of */
        while (aj[j] >= ranges[proc+1]) proc++;
        nx_all = rranges[proc+1] - rranges[proc];
        col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + ng_all[proc];
        ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
      }
      ierr = 
MatRestoreRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
    }

    if (pdipm->nh) {
      /* Insert -grad h' */
      ierr = 
MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
      ierr = 
MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);CHKERRQ(ierr);
      proc = 0;
      for (j=0; j < nc; j++) {
        /* find row ownership of */
        while (aj[j] >= ranges[proc+1]) proc++;
        nx_all = rranges[proc+1] - rranges[proc];
        col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
        ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
      }
      ierr = 
MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
    }

    /* Insert Jci_xb^T' */
    ierr = MatGetRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
    ierr = MatGetOwnershipRanges(pdipm->Jci_xb,&ranges);CHKERRQ(ierr);
    proc = 0;
    for (j=0; j < nc; j++) {
      /* find row ownership of */
      while (aj[j] >= ranges[proc+1]) proc++;
      nx_all = rranges[proc+1] - rranges[proc];
      col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc] + 
nh_all[proc];
      ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
    }
    ierr = MatRestoreRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
  }

  /* 2nd Row block of KKT matrix: [grad Ce, 0, 0, 0] */
  if (pdipm->Ng) {
    ierr = 
MatGetOwnershipRange(tao->jacobian_equality,&rjstart,NULL);CHKERRQ(ierr);
    for (i=0; i < pdipm->ng; i++){
      row = rstart + pdipm->off_lambdae + i;

      ierr = 
MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
      proc = 0;
      for (j=0; j < nc; j++) {
        while (aj[j] >= cranges[proc+1]) proc++;
        col = aj[j] - cranges[proc] + Jranges[proc];
        ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr); /* grad g */
      }
      ierr = 
MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
    }
  }
  /* Jce_xfixed */
  if (pdipm->Nxfixed) {
    ierr = MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);CHKERRQ(ierr);
    for (i=0; i < (pdipm->nce - pdipm->ng); i++){
      row = rstart + pdipm->off_lambdae + pdipm->ng + i;

      ierr = 
MatGetRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);CHKERRQ(ierr);
      if (nc != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nc != 1");

      proc = 0;
      j    = 0;
      while (cols[j] >= cranges[proc+1]) proc++;
      col = cols[j] - cranges[proc] + Jranges[proc];
      ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
      ierr = 
MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);CHKERRQ(ierr);
    }
  }

  /* 3rd Row block of KKT matrix: [ gradCi, 0, 0, -I] */
  if (pdipm->Nh) {
    ierr = 
MatGetOwnershipRange(tao->jacobian_inequality,&rjstart,NULL);CHKERRQ(ierr);
    for (i=0; i < pdipm->nh; i++){
      row = rstart + pdipm->off_lambdai + i;

      ierr = 
MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
      proc = 0;
      for (j=0; j < nc; j++) {
        while (aj[j] >= cranges[proc+1]) proc++;
        col = aj[j] - cranges[proc] + Jranges[proc];
        ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr); /* grad h */
      }
      ierr = 
MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
    }
    /* -I */
    for (i=0; i < pdipm->nh; i++){
      row = rstart + pdipm->off_lambdai + i;
      col = rstart + pdipm->off_z + i;
      ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
    }
  }

  /* Jci_xb */
  ierr = MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,NULL);CHKERRQ(ierr);
  for (i=0; i < (pdipm->nci - pdipm->nh); i++){
    row = rstart + pdipm->off_lambdai + pdipm->nh + i;

    ierr = MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);CHKERRQ(ierr);
    if (nc != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nc != 1");
    proc = 0;
    for (j=0; j < nc; j++) {
      while (cols[j] >= cranges[proc+1]) proc++;
      col = cols[j] - cranges[proc] + Jranges[proc];
      ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
    }
    ierr = MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);CHKERRQ(ierr);
    /* -I */
    col = rstart + pdipm->off_z + pdipm->nh + i;
    ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
  }

  /* 4-th Row block of KKT matrix: Z and Ci */
  for (i=0; i < pdipm->nci; i++) {
    row     = rstart + pdipm->off_z + i;
    cols1[0] = rstart + pdipm->off_lambdai + i;
    cols1[1] = row;
    ierr = MatPreallocateSet(row,2,cols1,dnz,onz);CHKERRQ(ierr);
  }



On Sep 15, 2020, at 4:03 PM, Abhyankar, Shrirang G 
<shrirang.abhyan...@pnnl.gov<mailto:shrirang.abhyan...@pnnl.gov>> wrote:

There is some code in PDIPM that copies over matrix elements from user 
jacobians/hessians to the big KKT matrix. I am not sure how/if that will work 
with MatMPI_Column. We’ll have to try out and we need an example for that 😊

Thanks,
Shri


From: Barry Smith <bsm...@petsc.dev<mailto:bsm...@petsc.dev>>
Date: Tuesday, September 15, 2020 at 1:21 PM
To: Pierre Jolivet 
<pierre.joli...@enseeiht.fr<mailto:pierre.joli...@enseeiht.fr>>
Cc: "Abhyankar, Shrirang G" 
<shrirang.abhyan...@pnnl.gov<mailto:shrirang.abhyan...@pnnl.gov>>, PETSc 
Development <petsc-dev@mcs.anl.gov<mailto:petsc-dev@mcs.anl.gov>>
Subject: Re: [petsc-dev] PDIPDM questions


  Pierre,

    Based on my previous mail I am hoping that the PDIPM algorithm itself won't 
need a major refactorization to be scalable, only a custom matrix type is 
needed to store and compute with the  Hessian in a scalable way.

   Barry




On Sep 15, 2020, at 12:50 PM, Pierre Jolivet 
<pierre.joli...@enseeiht.fr<mailto:pierre.joli...@enseeiht.fr>> wrote:





On 15 Sep 2020, at 5:40 PM, Abhyankar, Shrirang G 
<shrirang.abhyan...@pnnl.gov<mailto:shrirang.abhyan...@pnnl.gov>> wrote:

Pierre,
   You are right. There are a few MatMultTransposeAdd that may need conforming 
layouts for the equality/inequality constraint vectors and equality/inequality 
constraint Jacobian matrices. I need to check if that’s the case. We only have 
ex1 example currently, we need to add more examples. We are currently working 
on making PDIPM robust and while doing it will work on adding another example.

Very naive question, but given that I have a single constraint, how do I split 
a 1 x N matrix column-wise? I thought it was not possible.

When setting the size of the constraint vector, you need to set the local size 
on one rank to 1 and all others to zero. For the Jacobian, the local row size 
on that rank will be 1 and all others to zero. The column layout for the 
Jacobian should follow the layout for vector x. So each rank will set the local 
column size of the Jacobian to local size of x.

That is assuming I don’t want x to follow the distribution of the Hessian, 
which is not my case.
Is there some plan to make PDIPM handle different layouts?
I hope I’m not the only one thinking that having a centralized Hessian when 
there is a single constraint is not scalable?

Thanks,
Pierre

Shri

On 15 Sep 2020, at 2:21 AM, Abhyankar, Shrirang G 
<shrirang.abhyan...@pnnl.gov<mailto:shrirang.abhyan...@pnnl.gov>> wrote:

Hello Pierre,
   PDIPM works in parallel so you can have distributed Hessian, Jacobians, 
constraints, variables, gradients in any layout you want.  If you are using a 
DM then you can have it generate the Hessian.

Could you please show an example where this is the case?
pdipm->x, which I’m assuming is a working vector, is both used as input for 
Hessian and Jacobian functions, e.g., 
https://gitlab.com/petsc/petsc/-/blob/master/src/tao/constrained/impls/ipm/pdipm.c#L369
 (Hessian) + 
https://gitlab.com/petsc/petsc/-/blob/master/src/tao/constrained/impls/ipm/pdipm.c#L473
 (Jacobian)
I thus doubt that it is possible to have different layouts?
In practice, I end up with the following error when I try this (2 processes, 
distributed Hessian with centralized Jacobian):
[1]PETSC ERROR: --------------------- Error Message 
--------------------------------------------------------------
[1]PETSC ERROR: Nonconforming object sizes
[1]PETSC ERROR: Vector wrong size 14172 for scatter 0 (scatter reverse and 
vector to != ctx from size)
[1]PETSC ERROR: #1 VecScatterBegin() line 96 in 
/Users/jolivet/Documents/repositories/petsc/src/vec/vscat/interface/vscatfce.c
[1]PETSC ERROR: #2 MatMultTransposeAdd_MPIAIJ() line 1223 in 
/Users/jolivet/Documents/repositories/petsc/src/mat/impls/aij/mpi/mpiaij.c
[1]PETSC ERROR: #3 MatMultTransposeAdd() line 2648 in 
/Users/jolivet/Documents/repositories/petsc/src/mat/interface/matrix.c
[0]PETSC ERROR: Nonconforming object sizes
[0]PETSC ERROR: Vector wrong size 13790 for scatter 27962 (scatter reverse and 
vector to != ctx from size)
[1]PETSC ERROR: #4 TaoSNESFunction_PDIPM() line 510 in 
/Users/jolivet/Documents/repositories/petsc/src/tao/constrained/impls/ipm/pdipm.c
[0]PETSC ERROR: #5 TaoSolve_PDIPM() line 712 in 
/Users/jolivet/Documents/repositories/petsc/src/tao/constrained/impls/ipm/pdipm.c
[1]PETSC ERROR: #6 TaoSolve() line 222 in 
/Users/jolivet/Documents/repositories/petsc/src/tao/interface/taosolver.c
[0]PETSC ERROR: #1 VecScatterBegin() line 96 in 
/Users/jolivet/Documents/repositories/petsc/src/vec/vscat/interface/vscatfce.c
[0]PETSC ERROR: #2 MatMultTransposeAdd_MPIAIJ() line 1223 in 
/Users/jolivet/Documents/repositories/petsc/src/mat/impls/aij/mpi/mpiaij.c
[0]PETSC ERROR: #3 MatMultTransposeAdd() line 2648 in 
/Users/jolivet/Documents/repositories/petsc/src/mat/interface/matrix.c
[0]PETSC ERROR: #4 TaoSNESFunction_PDIPM() line 510 in 
/Users/jolivet/Documents/repositories/petsc/src/tao/constrained/impls/ipm/pdipm.c
[0]PETSC ERROR: #5 TaoSolve_PDIPM() line 712 in 
/Users/jolivet/Documents/repositories/petsc/src/tao/constrained/impls/ipm/pdipm.c
[0]PETSC ERROR: #6 TaoSolve() line 222 in 
/Users/jolivet/Documents/repositories/petsc/src/tao/interface/taosolver.c

I think this can be reproduced by ex1.c by just distributing the Hessian 
instead of having it centralized on rank 0.




Ideally, you want to have the layout below to minimize movement of 
matrix/vector elements across ranks.

•         The layout of vectors x, bounds on x, and gradient is same.

•         The row layout of the equality/inequality Jacobian is same as the 
equality/inequality constraint vector layout.

•         The column layout of the equality/inequality Jacobian is same as that 
for x.

Very naive question, but given that I have a single constraint, how do I split 
a 1 x N matrix column-wise? I thought it was not possible.

Thanks,
Pierre





•         The row and column layout for the Hessian is same as x.

The tutorial example ex1 is extremely small (only 2 variables) so its 
implementation is very simplistic. I think, in parallel, it ships off 
constraints etc. to rank 0. It’s not an ideal example w.r.t demonstrating a 
parallel implementation. We aim to add more examples as we develop PDIPM. If 
you have an example to contribute then we would most welcome it and provide 
help on adding it.

Thanks,
Shri
From: petsc-dev 
<petsc-dev-boun...@mcs.anl.gov<mailto:petsc-dev-boun...@mcs.anl.gov>> on behalf 
of Pierre Jolivet 
<pierre.joli...@enseeiht.fr<mailto:pierre.joli...@enseeiht.fr>>
Date: Monday, September 14, 2020 at 1:52 PM
To: PETSc Development <petsc-dev@mcs.anl.gov<mailto:petsc-dev@mcs.anl.gov>>
Subject: [petsc-dev] PDIPDM questions

Hello,
In my quest to help users migrate from Ipopt to Tao, I’ve a new question.
When looking at src/tao/constrained/tutorials/ex1.c, it seems that almost 
everything is centralized on rank 0 (local sizes are 0 but on rank 0).
I’d like to have my Hessian distributed more naturally, as in (almost?) all 
other SNES/TS examples, but still keep the Jacobian of my equality constraint, 
which is of dimension 1 x N (N >> 1), centralized on rank 0.
Is this possible?
If not, is it possible to supply the transpose of the Jacobian, of dimension N 
x 1, which could then be distributed row-wise like the Hessian?
Or maybe use some trick to distribute a MatAIJ/MatDense of dimension 1 x N 
column-wise? Use a MatNest with as many blocks as processes?

So, just to sum up, how can I have a distributed Hessian with a Jacobian with a 
single row?

Thanks in advance for your help,
Pierre