[petsc-users] [EXTERNAL] Re: Example code of linear solve of a block matrix system in Fortran

[Satish Balay]

Its best to keep the discussion on the list [or petsc-maint] - as others would have answers to some of these qestions.

Wrt -info - yeah this is the correct option [my mistake in suggesting -log_info]

Wrt understanding performance -log_view should help

Satish


On Thu, 7 May 2020, Thomas S. Chyczewski wrote:

> Satish,
> 
> Thanks for your help. Sorry for not being able to figure it out on my own. I had a little trouble following the discussion in the manual. I got the block version working now and the linear solver time is cut in half compared to the monolithic version for a given set of parameters. I have experimented with a number of solver and preconditioner options as well as solver convergence criteria. But I'm wondering if there are any other parameters I should be playing with. I ask because I also had some trouble following the PCFIELDSPLIT discussion in the manual and I'm wondering if the default is the best option for me.
> 
> The -log_info option isn't available in the Fortran version, so I couldn't check the inode information as you suggested. However, below is the output when I run with the -info option. I know that having no mallocs during MatSetValues is good, but that's about it.
> 
> Thanks,
> Tom
> 
> [0] PetscGetHostName(): Rejecting domainname, likely is NIS arl19814.(none)
> [0] petscinitialize_internal(): (Fortran):PETSc successfully started: procs 1
> [0] PetscGetHostName(): Rejecting domainname, likely is NIS arl19814.(none)
> [0] petscinitialize_internal(): Running on machine: arl19814
> [0] PetscCommDuplicate(): Duplicating a communicator 2 2 max tags = 100000000
> [0] MatAssemblyEnd_SeqBAIJ(): Matrix size: 374400 X 374400, block size 3; storage space: 11199240 unneeded, 5606640 used
> [0] MatAssemblyEnd_SeqBAIJ(): Number of mallocs during MatSetValues is 0
> [0] MatAssemblyEnd_SeqBAIJ(): Most nonzeros blocks in any row is 5
> [0] MatCheckCompressedRow(): Found the ratio (num_zerorows 0)/(num_localrows 124800) < 0.6. Do not use CompressedRow routines.
> [0] PCSetUp(): Setting up PC for first time
> [0] PetscCommDuplicate(): Duplicating a communicator 1 3 max tags = 100000000
> [0] PetscCommDuplicate(): Using internal PETSc communicator 1 3
> [0] PCSetUp(): Leaving PC with identical preconditioner since operator is unchanged
> [0] PCSetUp(): Leaving PC with identical preconditioner since operator is unchanged
> [0] PCSetUp(): Leaving PC with identical preconditioner since operator is unchanged
> [0] PCSetUp(): Leaving PC with identical preconditioner since operator is unchanged
> [0] PCSetUp(): Leaving PC with identical preconditioner since operator is unchanged
> [0] PCSetUp(): Leaving PC with identical preconditioner since operator is unchanged
> [0] PCSetUp(): Leaving PC with identical preconditioner since operator is unchanged
> [0] PCSetUp(): Leaving PC with identical preconditioner since operator is unchanged
> [0] PCSetUp(): Leaving PC with identical preconditioner since operator is unchanged
> 
> 
> 
> -----Original Message-----
> From: Satish Balay <balay at mcs.anl.gov> 
> Sent: Wednesday, May 6, 2020 4:05 PM
> To: Thomas S. Chyczewski <tsc109 at arl.psu.edu>
> Cc: petsc-users <petsc-users at mcs.anl.gov>
> Subject: Re: [petsc-users] [EXTERNAL] Re: Example code of linear solve of a block matrix system in Fortran
> 
> You can use https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSetValuesBlocked.html
> 
> MatSetValues() will also work - but MatSetValuesBlocked() is more efficient wrt BAIJ
> 
> Satish
> 
> On Wed, 6 May 2020, Thomas S. Chyczewski wrote:
> 
> > Thanks Satish. I have seen that page and can create a block matrix. It's not clear to me how to fill it and use it in a Fortran code.
> > 
> > -----Original Message-----
> > From: Satish Balay <balay at mcs.anl.gov> 
> > Sent: Wednesday, May 6, 2020 3:43 PM
> > To: Thomas S. Chyczewski <tsc109 at arl.psu.edu>
> > Cc: petsc-users at mcs.anl.gov
> > Subject: [EXTERNAL] Re: [petsc-users] Example code of linear solve of a block matrix system in Fortran
> > 
> > What you are looking for is:
> > 
> > https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateBAIJ.html
> > 
> > indoe is optimization for AIJ type - where you might not have a simple block size like 3 that you have.
> > 
> > BAIJ will perform better than AIJ/with-inode
> > 
> > To verify indoe usage - you can run the code with '-log_info' option - and do a 'grep inode'
> > 
> > Satish
> > 
> > On Wed, 6 May 2020, Thomas S. Chyczewski wrote:
> > 
> > > 
> > > All,
> > > 
> > > I'm relatively new to PETSc and have relied pretty heavily on the example codes included in the distribution to figure out the finer points of using the PETSc library that I couldn't deduce from the manual. One thing I can't figure out is how to solve block matrix systems and I couldn't find an example in Fortran. I'm writing a 2D incompressible CFD solver so I have a 3x3 block Imax*Jmax system. The closest I've come to finding an example is ex19.c in the snes directory, but that is in c and for the nonlinear solver.
> > > 
> > > I have been able to run PETSc but unwrapping the block matrix into a monolithic system. But the manual says "Block matrices represent an important class of problems in numerical linear algebra and offer the possibility of far more efficient iterative solvers than just treating the entire matrix as black box." However, in the FAQs I saw a comment that PETSc scans the AIJ matrices for rows that have the same column layout and can deduce if it's a block system and use the more efficient solvers. I also saw in the archives for this email list a thread where it seems workaround for building fields in a Fortran code is discussed ("Back to struct in Fortran to represent field with dof > 1"), so I'm beginning to suspect building a block system in Fortran might not be straight forward.
> > > 
> > > All that being said, my questions:
> > > 
> > > Is there a significant advantage to building the block system as opposed to the analogous monolithic system if PETSc can figure out that it's a block system? Can you confirm that PETSc does figure this out?
> > > If there is an advantage to loading the matrix as a block matrix, is there an example Fortran code that builds and solves a linear block system?
> > > 
> > > Thanks,
> > > Tom C
> > > 
> > > 
> > 
>