> On 15 Sep 2020, at 11:18 PM, Barry Smith <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.
> 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,
Pierre
> 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
>>>>
>>>> <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
>>>>
>>>> <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
