[petsc-users] Problem with MatZeroRowsColumnsIS()

[Barry Smith]

    Given that Jed wrote MatZeroRowsColumns_MPIAIJ() it is unlikely to be wrong.

  You wrote

Calculating the norm2 of the residuals defined above in each case gives:
MatZeroRowsIS() 1cpu:  norm(res,2) =      0
MatZeroRowsIS() 4cpu:  norm(res,2) =      0
MatZeroRowsColumnsIS() 1cpu:  norm(res,2) =    1.6880e-10
MatZeroRowsColumnsIS() 4cpu:  norm(res,2) =    7.3786e+06

  why do you conclude this is wrong? MatZeroRowsColumnsIS() IS suppose to change the right hand side in a way different than MatZeroRowsIS(). 

   Explanation. For simplicity reorder the matrix rows/columns so that zeroed ones come last and the matrix is symmetric. Then you have 

        (  A   B )   (x_A)  =   (b_A)
        (  B   D )   (x_B)       (b_B)

   with MatZeroRows the new system is     

       (  A   B )   (x_A)  =   (b_A)
        (  0   I  )   (x_B)       (x_B)

   it has the same solution as the original problem with the give x_B

  with MatZeroRowsColumns the new system is 

      (  A   0 )   (x_A)  =   (b_A)  - B*x_B
      (  0   I  )   (x_B)       (x_B)

note the right hand side needs to be changed so that the new problem has the same solution.

   Barry


On May 9, 2014, at 9:50 AM, Püsök, Adina-Erika <puesoek at uni-mainz.de> wrote:

> Yes, I tested the implementation with both MatZeroRowsIS() and MatZeroRowsColumnsIS(). But first, I will be more explicit about the problem I was set to solve: 
> 
> We have a Dirichlet block of size (L,W,H) and centered (xc,yc,zc), which is much smaller than the model domain, and we set Vx = Vpush, Vy=0 within the block (Vz is let free for easier convergence). 
> As I said before, since the code does not have a monolithic matrix, but 4 submatrices (VV VP; PV PP), and the rhs has 2 sub vectors rhs=(f; g), my approach is to modify only (VV, VP, f) for the Dirichlet BC.
> 
> The way I tested the implementation:
> 1) Output (VV, VP, f, Dirichlet dofs) - unmodified (no Dirichlet BC)
> 2) Output (VV, VP, f, Dirichlet dofs) - a) modified with MatZeroRowsIS(), 
>        - b) modified with MatZeroRowsColumnsIS() -> S_PETSc
> Again, the only difference between a) and b) is:
> //
> ierr = MatZeroRowsColumnsIS(VV_MAT,isx,v_vv,x_push,rhs);  CHKERRQ(ierr);
>    // ierr = MatZeroRowsColumnsIS(VV_MAT,isy,v_vv,x_push,rhs);  CHKERRQ(ierr);
> 
> ierr = MatZeroRowsIS(VV_MAT,isx,v_vv,x_push,rhs);  CHKERRQ(ierr);
> ierr = MatZeroRowsIS(VV_MAT,isy,v_vv,x_push,rhs);  CHKERRQ(ierr);
> 3) Read them in Matlab and perform the exact same operations on the unmodified matrices and f vector. -> S_Matlab
> 4) Compare S_PETSc with S_Matlab. If the implementation is correct, they should be equal (VV, VP, f).
> 5) Check for 1 cpu and 4 cpus.
> 
> Now to answer your questions:
> 
> a,b,d) Yes, matrix modification is done correctly (check the spy diagrams below) in all cases:  MatZeroRowsIS() and MatZeroRowsColumnsIS() on 1 and 4 cpus.
> 
> I should have said that in the piece of code above:
>  v_vv = 1.0;
>  v_vp = 0.0;
> The vector x_push is a duplicate of rhs, with zero elements except the values for the Dirichlet dofs. 
> 
> c) The rhs is a different matter. With MatZeroRows() there is no problem. The rhs is equivalent with the one in Matlab, sequential and parallel.
>     However, with MatZeroRowsColumns(), the residual contains nonzero elements, and in parallel the nonzero pattern is even bigger (1 cpu - 63, 4 cpu - 554). But if you look carefully, the values of the nonzero residuals are very small < +/- 1e-10.
> So, I did a tolerance filter:
> 
> tol = 1e-10;
> res = f_petsc - f_mod_matlab;
> for i=1:length(res)
>     if abs(res(i))>0 & abs(res(i))<tol
>         res(i)=0;
>     end
> end
> 
> and then the f_petsc and f_mod_matlab are equivalent on 1 and 4 cpus (figure 5). So it seems that MatZeroRowsColumnsIS() might give some nonzero residuals. 
> 
> Calculating the norm2 of the residuals defined above in each case gives:
> MatZeroRowsIS() 1cpu:  norm(res,2) =      0
> MatZeroRowsIS() 4cpu:  norm(res,2) =      0
> MatZeroRowsColumnsIS() 1cpu:  norm(res,2) =    1.6880e-10
> MatZeroRowsColumnsIS() 4cpu:  norm(res,2) =    7.3786e+06
> 
> Since this is purely a problem of matrix and vector assembly/manipulation, I think the nonzero residuals of the rhs with MatZeroRowsColumnsIS() give the parallel artefacts that I showed last time.
> If you need the raw data and the matlab scripts that I used for testing for your consideration, please let me know. 
> 
> Thanks,
> Adina
> 
> When performing the manual operations on the unmodified matrices and rhs vector in Matlab, I took into account: 
> - matlab indexing = petsc indexing +1;
> - the vectors written to file for matlab (PETSC_VIEWER_BINARY_MATLAB) have the natural ordering, rather than the petsc ordering. On 1 cpu, they are equivalent, but on 4 cpus, the Dirichlet BC indices had to be converted to natural indices in order to perform the correct operations on the rhs.
> 
> <spy_zerorows_1cpu.png>
> <spy_zerorows_4cpu.png>
> <spy_zerorowscol_1cpu.png>
> <spy_zerorowscol_4cpu.png>
> <residual_tol.png>
> 
> On May 6, 2014, at 4:22 PM, Matthew Knepley wrote:
> 
>> On Tue, May 6, 2014 at 7:23 AM, Püsök, Adina-Erika <puesoek at uni-mainz.de> wrote:
>> Hello!
>> 
>> I was trying to implement some internal Dirichlet boundary conditions into an aij matrix of the form:  A=(  VV  VP; PV PP ). The idea was to create an internal block (let's say Dirichlet block) that moves with constant velocity within the domain (i.e. check all the dofs within the block and set the values accordingly to the desired motion).
>> 
>> Ideally, this means to zero the rows and columns in VV, VP, PV corresponding to the dirichlet dofs and modify the corresponding rhs values. However, since we have submatrices and not a monolithic matrix A,  we can choose to modify only VV and PV matrices. 
>> The global indices of the velocity points within the Dirichlet block are contained in the arrays rowid_array. 
>> 
>> What I want to point out is that the function MatZeroRowsColumnsIS() seems to create parallel artefacts, compared to MatZeroRowsIS() when run on more than 1 processor. Moreover, the results on 1 cpu are identical. 
>> See below the results of the test (the Dirichlet block is outlined in white) and the piece of the code involved where the 1) - 2) parts are the only difference. 
>> 
>> I am assuming that you are showing the result of solving the equations. It would be more useful, and presumably just as easy
>> to say:
>> 
>>   a) Are the correct rows zeroed out?
>> 
>>   b) Is the diagonal element correct?
>> 
>>   c) Is the rhs value correct?
>> 
>>   d) Are the columns zeroed correctly?
>> 
>> If we know where the problem is, its easier to fix. For example, if the rhs values are
>> correct and the rows are zeroed, then something is wrong with the solution procedure.
>> Since ZeroRows() works and ZeroRowsColumns() does not, this is a distinct possibility.
>> 
>>   Thanks,
>> 
>>      Matt
>>  
>> Thanks,
>> Adina Pusok
>> 
>> // Create an IS required by MatZeroRows()
>> ierr = ISCreateGeneral(PETSC_COMM_WORLD,numRowsx,rowidx_array,PETSC_COPY_VALUES,&isx);  CHKERRQ(ierr);
>> ierr = ISCreateGeneral(PETSC_COMM_WORLD,numRowsy,rowidy_array,PETSC_COPY_VALUES,&isy);  CHKERRQ(ierr);
>> ierr = ISCreateGeneral(PETSC_COMM_WORLD,numRowsz,rowidz_array,PETSC_COPY_VALUES,&isz);  CHKERRQ(ierr);
>> 
>> 1) /*
>> ierr = MatZeroRowsColumnsIS(VV_MAT,isx,v_vv,x_push,rhs);  CHKERRQ(ierr);
>> ierr = MatZeroRowsColumnsIS(VV_MAT,isy,v_vv,x_push,rhs);  CHKERRQ(ierr);
>> ierr = MatZeroRowsColumnsIS(VV_MAT,isz,v_vv,x_push,rhs);  CHKERRQ(ierr);*/
>> 
>> 2) ierr = MatZeroRowsIS(VV_MAT,isx,v_vv,x_push,rhs);  CHKERRQ(ierr);
>> ierr = MatZeroRowsIS(VV_MAT,isy,v_vv,x_push,rhs);  CHKERRQ(ierr);
>> ierr = MatZeroRowsIS(VV_MAT,isz,v_vv,x_push,rhs);  CHKERRQ(ierr);
>> 
>> ierr = MatZeroRowsIS(VP_MAT,isx,v_vp,PETSC_NULL,PETSC_NULL);  CHKERRQ(ierr);
>> ierr = MatZeroRowsIS(VP_MAT,isy,v_vp,PETSC_NULL,PETSC_NULL);  CHKERRQ(ierr);
>> ierr = MatZeroRowsIS(VP_MAT,isz,v_vp,PETSC_NULL,PETSC_NULL);  CHKERRQ(ierr);
>> 
>> ierr = ISDestroy(&isx); CHKERRQ(ierr);
>> ierr = ISDestroy(&isy); CHKERRQ(ierr);
>> ierr = ISDestroy(&isz); CHKERRQ(ierr);
>> 
>> 
>> Results (velocity) with MatZeroRowsColumnsIS().
>> 1cpu<r01_1cpu_rows_columns.png> 4cpu<r01_rows_columns.png>
>> 
>> Results (velocity) with MatZeroRowsIS():
>> 1cpu<r01_1cpu_rows.png> 4cpu<r01_rows.png>
>> 
>> 
>> 
>> -- 
>> 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
>