xref: /phasta/phSolver/compressible/elmgmr.f (revision 712d3df0b59ebebaaeaea358162c8d2c043c6e08)

        subroutine ElmGMRe (y,         ac,        x,
     &                     shp,       shgl,      iBC,
     &                     BC,        shpb,      shglb,
     &                     res,       rmes,      BDiag,
     &                     iper,      ilwork,    EGmass, rerr)
c
c----------------------------------------------------------------------
c
c This routine computes the LHS mass matrix, the RHS residual
c vector, and the preconditioning matrix, for use with the GMRES
c solver.
c
c Zdenek Johan, Winter 1991.      (Fortran 90)
c Chris Whiting, Winter 1998.     (Matrix EBE-GMRES)
c----------------------------------------------------------------------
c
        use pointer_data
        use timedataC
c
        include "common.h"
        include "mpif.h"
c
        dimension y(nshg,ndof),         ac(nshg,ndof),
     &            x(numnp,nsd),
     &            iBC(nshg),
     &            BC(nshg,ndofBC),
     &            res(nshg,nflow),
     &            rmes(nshg,nflow),      BDiag(nshg,nflow,nflow),
     &            iper(nshg),           EGmass(numel,nedof,nedof)
c
        dimension shp(MAXTOP,maxsh,MAXQPT),
     &            shgl(MAXTOP,nsd,maxsh,MAXQPT),
     &            shpb(MAXTOP,maxsh,MAXQPT),
     &            shglb(MAXTOP,nsd,maxsh,MAXQPT)
c
        dimension qres(nshg, idflx),     rmass(nshg)
c
        dimension ilwork(nlwork)

        real*8 Bdiagvec(nshg,nflow), rerr(nshg,10)

        real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:)
        real*8, allocatable :: tmpshpb(:,:), tmpshglb(:,:,:)

	ttim(80) = ttim(80) - secs(0.0)
c
c.... set up the timer
c

        call timer ('Elm_Form')
c
c.... -------------------->   interior elements   <--------------------
c
c.... set up parameters
c
        ires   = 1
c
        if (idiff==1 .or. idiff==3 .or. isurf==1) then ! global reconstruction
                                                       ! of qdiff
c
c loop over element blocks for the global reconstruction
c of the diffusive flux vector, q, and lumped mass matrix, rmass
c
        qres = zero
        rmass = zero

        do iblk = 1, nelblk
c
c.... set up the parameters
c
          nenl   = lcblk(5,iblk)   ! no. of vertices per element
          iel    = lcblk(1,iblk)
          lelCat = lcblk(2,iblk)
          lcsyst = lcblk(3,iblk)
          iorder = lcblk(4,iblk)
          nenl   = lcblk(5,iblk)   ! no. of vertices per element
          nshl   = lcblk(10,iblk)
          mattyp = lcblk(7,iblk)
          ndofl  = lcblk(8,iblk)
          nsymdl = lcblk(9,iblk)
          npro   = lcblk(1,iblk+1) - iel
          ngauss = nint(lcsyst)
c
c.... compute and assemble diffusive flux vector residual, qres,
c     and lumped mass matrix, rmass

          allocate (tmpshp(nshl,MAXQPT))
          allocate (tmpshgl(nsd,nshl,MAXQPT))

          tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:)
          tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:)

          call AsIq (y,                x,
     &               tmpshp,
     &               tmpshgl,
     &               mien(iblk)%p,     mxmudmi(iblk)%p,
     &               qres,
     &               rmass)

          deallocate ( tmpshp )
          deallocate ( tmpshgl )
       enddo

c
c.... take care of periodic boundary conditions
c

       call qpbc( rmass, qres, iBC,  iper, ilwork )
c
      endif                     ! computation of global diffusive flux
c
c.... loop over element blocks to compute element residuals
c
c
c.... initialize the arrays
c
        res    = zero
        rmes   = zero ! to avoid trap_uninitialized
        if (lhs. eq. 1)   EGmass = zero
        if (iprec .ne. 0) BDiag = zero
        flxID = zero

c
c.... loop over the element-blocks
c
        do iblk = 1, nelblk
c
c.... set up the parameters
c
          iblkts = iblk          ! used in timeseries
          nenl   = lcblk(5,iblk) ! no. of vertices per element
          iel    = lcblk(1,iblk)
          lelCat = lcblk(2,iblk)
          lcsyst = lcblk(3,iblk)
          iorder = lcblk(4,iblk)
          nenl   = lcblk(5,iblk)   ! no. of vertices per element
          nshl   = lcblk(10,iblk)
          mattyp = lcblk(7,iblk)
          ndofl  = lcblk(8,iblk)
          nsymdl = lcblk(9,iblk)
          npro   = lcblk(1,iblk+1) - iel
          inum   = iel + npro - 1
          ngauss = nint(lcsyst)
c
c.... compute and assemble the residual and tangent matrix
c

          allocate (tmpshp(nshl,MAXQPT))
          allocate (tmpshgl(nsd,nshl,MAXQPT))

          tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:)
          tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:)

          call AsIGMR (y,                   ac,
     &                 x,                   mxmudmi(iblk)%p,
     &                 tmpshp,
     &                 tmpshgl,             mien(iblk)%p,
     &                 mmat(iblk)%p,        res,
     &                 rmes,                BDiag,
     &                 qres,                EGmass(iel:inum,:,:),
     &                 rerr)
c
c.... satisfy the BC's on the implicit LHS
c
          call bc3LHS (iBC,                  BC,  mien(iblk)%p,
     &                 EGmass(iel:inum,:,:)  )

          deallocate ( tmpshp )
          deallocate ( tmpshgl )
c
c.... end of interior element loop
c
       enddo
c
c.... -------------------->   boundary elements   <--------------------
c
c.... loop over the boundary elements
c
        do iblk = 1, nelblb
c
c.... set up the parameters
c
          iel    = lcblkb(1,iblk)
          lelCat = lcblkb(2,iblk)
          lcsyst = lcblkb(3,iblk)
          iorder = lcblkb(4,iblk)
          nenl   = lcblkb(5,iblk)  ! no. of vertices per element
          nenbl  = lcblkb(6,iblk)  ! no. of vertices per bdry. face
          mattyp = lcblkb(7,iblk)
          ndofl  = lcblkb(8,iblk)
          nshl   = lcblkb(9,iblk)
          nshlb  = lcblkb(10,iblk)
          npro   = lcblkb(1,iblk+1) - iel
          if(lcsyst.eq.3) lcsyst=nenbl
          ngaussb = nintb(lcsyst)

c
c.... compute and assemble the residuals corresponding to the
c     boundary integral
c

          allocate (tmpshpb(nshl,MAXQPT))
          allocate (tmpshglb(nsd,nshl,MAXQPT))

          tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:)
          tmpshglb(:,1:nshl,:) = shglb(lcsyst,:,1:nshl,:)

          call AsBMFG (y,                       x,
     &                 tmpshpb,                 tmpshglb,
     &                 mienb(iblk)%p,           mmatb(iblk)%p,
     &                 miBCB(iblk)%p,           mBCB(iblk)%p,
     &                 res,                     rmes)

          deallocate (tmpshpb)
          deallocate (tmpshglb)
c
c.... end of boundary element loop
c
        enddo
c
      ttim(80) = ttim(80) + secs(0.0)
c
c before the commu we need to rotate the residual vector for axisymmetric
c boundary conditions (so that off processor periodicity is a dof add instead
c of a dof combination).  Take care of all nodes now so periodicity, like
c commu is a simple dof add.
c
c      if(iabc==1)               !are there any axisym bc's
c     &     call rotabc(res(1,2), iBC, BC, nflow,  'in ')
      if(iabc==1) then               !are there any axisym bc's
          call rotabc(res(1,2), iBC,  'in ')
c          Bdiagvec(:,1)=BDiag(:,1,1)
c          Bdiagvec(:,2)=BDiag(:,2,2)
c          Bdiagvec(:,3)=BDiag(:,3,3)
c          Bdiagvec(:,4)=BDiag(:,4,4)
c          Bdiagvec(:,5)=BDiag(:,5,5)
c          call rotabc(Bdiagvec(1,2), iBC, BC, 2,  'in ')
c          BDiag(:,:,:)=zero
c          BDiag(:,1,1)=Bdiagvec(:,1)
c          BDiag(:,2,2)=Bdiagvec(:,2)
c          BDiag(:,3,3)=Bdiagvec(:,3)
c          BDiag(:,4,4)=Bdiagvec(:,4)
c          BDiag(:,5,5)=Bdiagvec(:,5)
       endif

c.... -------------------->   communications <-------------------------
c

      if (numpe > 1) then
        call commu (res  , ilwork, nflow  , 'in ')

        call MPI_BARRIER (MPI_COMM_WORLD,ierr)

        if(iprec .ne. 0) call commu (BDiag, ilwork, nflow*nflow, 'in ')
      endif

c
c.... ---------------------->   post processing  <----------------------
c
c.... satisfy the BCs on the residual
c
      call bc3Res (y,  iBC,  BC,  res,  iper, ilwork)
c
c.... satisfy the BCs on the block-diagonal preconditioner
c
      if (iprec .ne. 0) then
         call bc3BDg (y,  iBC,  BC,  BDiag, iper, ilwork)
      endif
c
c.... return
c
      call timer ('Back    ')
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccc       SPARSE
c_______________________________________________________________

        subroutine ElmGMRs (y,         ac,        x,
     &                     shp,       shgl,      iBC,
     &                     BC,        shpb,      shglb,
     &                     res,       rmes,      BDiag,
     &                     iper,      ilwork,    lhsK,
     &                     col,       row,       rerr)
c
c----------------------------------------------------------------------
c
c This routine computes the LHS mass matrix, the RHS residual
c vector, and the preconditioning matrix, for use with the GMRES
c solver.
c
c Zdenek Johan, Winter 1991.      (Fortran 90)
c Chris Whiting, Winter 1998.     (Matrix EBE-GMRES)
c----------------------------------------------------------------------
c
        use pointer_data
        use timedataC
c
        include "common.h"
        include "mpif.h"
c
        integer col(nshg+1), row(nnz*nshg)
        real*8 lhsK(nflow*nflow,nnz_tot)

        dimension y(nshg,ndof),         ac(nshg,ndof),
     &            x(numnp,nsd),
     &            iBC(nshg),
     &            BC(nshg,ndofBC),
     &            res(nshg,nflow),
     &            rmes(nshg,nflow),      BDiag(nshg,nflow,nflow),
     &            iper(nshg)
c
        dimension shp(MAXTOP,maxsh,MAXQPT),
     &            shgl(MAXTOP,nsd,maxsh,MAXQPT),
     &            shpb(MAXTOP,maxsh,MAXQPT),
     &            shglb(MAXTOP,nsd,maxsh,MAXQPT)
c
        dimension qres(nshg, idflx),     rmass(nshg)
c
        dimension ilwork(nlwork)

        real*8 Bdiagvec(nshg,nflow), rerr(nshg,10)

        real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:)
        real*8, allocatable :: tmpshpb(:,:), tmpshglb(:,:,:)
        real*8, allocatable :: EGmass(:,:,:)
        ttim(80) = ttim(80) - secs(0.0)
c
c.... set up the timer
c

        call timer ('Elm_Form')
c
c.... -------------------->   interior elements   <--------------------
c
c.... set up parameters
c
        ires   = 1
c
        if (idiff==1 .or. idiff==3 .or. isurf==1) then ! global reconstruction
                                                       ! of qdiff
c
c loop over element blocks for the global reconstruction
c of the diffusive flux vector, q, and lumped mass matrix, rmass
c
        qres = zero
        rmass = zero

        do iblk = 1, nelblk
c
c.... set up the parameters
c
          nenl   = lcblk(5,iblk)   ! no. of vertices per element
          iel    = lcblk(1,iblk)
          lelCat = lcblk(2,iblk)
          lcsyst = lcblk(3,iblk)
          iorder = lcblk(4,iblk)
          nshl   = lcblk(10,iblk)
          mattyp = lcblk(7,iblk)
          ndofl  = lcblk(8,iblk)
          nsymdl = lcblk(9,iblk)
          npro   = lcblk(1,iblk+1) - iel
          ngauss = nint(lcsyst)
c
c.... compute and assemble diffusive flux vector residual, qres,
c     and lumped mass matrix, rmass

          allocate (tmpshp(nshl,MAXQPT))
          allocate (tmpshgl(nsd,nshl,MAXQPT))

          tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:)
          tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:)

          call AsIq (y,                x,
     &               tmpshp,
     &               tmpshgl,
     &               mien(iblk)%p,     mxmudmi(iblk)%p,
     &               qres,
     &               rmass)

          deallocate ( tmpshp )
          deallocate ( tmpshgl )
       enddo

c
c.... take care of periodic boundary conditions
c

       call qpbc( rmass, qres, iBC, iper, ilwork )
c
      endif                     ! computation of global diffusive flux
c
c.... loop over element blocks to compute element residuals
c
c
c.... initialize the arrays
c
        res    = zero
        rmes   = zero ! to avoid trap_uninitialized
        if (lhs. eq. 1)   lhsK = zero
        if (iprec .ne. 0) BDiag = zero
        flxID = zero
c
c.... loop over the element-blocks
c
        do iblk = 1, nelblk
c
c.... set up the parameters
c
          iblkts = iblk          ! used in timeseries
          nenl   = lcblk(5,iblk) ! no. of vertices per element
          iel    = lcblk(1,iblk)
          lelCat = lcblk(2,iblk)
          lcsyst = lcblk(3,iblk)
          iorder = lcblk(4,iblk)
          nenl   = lcblk(5,iblk)   ! no. of vertices per element
          nshl   = lcblk(10,iblk)
          mattyp = lcblk(7,iblk)
          ndofl  = lcblk(8,iblk)
          nsymdl = lcblk(9,iblk)
          npro   = lcblk(1,iblk+1) - iel
          inum   = iel + npro - 1
          ngauss = nint(lcsyst)
c
c.... compute and assemble the residual and tangent matrix
c

          if(lhs.eq.1) then
             allocate (EGmass(npro,nedof,nedof))
             EGmass = zero
          else
             allocate (EGmass(1,1,1))
          endif

          allocate (tmpshp(nshl,MAXQPT))
          allocate (tmpshgl(nsd,nshl,MAXQPT))
          tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:)
          tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:)

          call AsIGMR (y,                   ac,
     &                 x,                   mxmudmi(iblk)%p,
     &                 tmpshp,
     &                 tmpshgl,             mien(iblk)%p,
     &                 mmat(iblk)%p,        res,
     &                 rmes,                BDiag,
     &                 qres,                EGmass,
     &                 rerr )
          if(lhs.eq.1) then
c
c.... satisfy the BC's on the implicit LHS
c
             call bc3LHS (iBC,                  BC,  mien(iblk)%p,
     &                    EGmass  )

c
c.... Fill-up the global sparse LHS mass matrix
c
             call fillsparseC( mien(iblk)%p, EGmass,
     1                        lhsK, row, col)
          endif
c
          deallocate ( EGmass )
          deallocate ( tmpshp )
          deallocate ( tmpshgl )
c
c.... end of interior element loop
c
       enddo
!ifdef DEBUG !Nicholas Mati
!        call write_debug(myrank, 'res-afterAsIGMR'//char(0),
!     &                           'res'//char(0), res,
!     &                           'd'//char(0), nshg, nflow, lstep)
!        call write_debug(myrank, 'y-afterAsIGMR'//char(0),
!     &                           'y'//char(0), y,
!     &                           'd'//char(0), nshg, ndof, lstep)
!endif //DEBUG
c
c.... -------------------->   boundary elements   <--------------------
c
c.... loop over the boundary elements
c
        do iblk = 1, nelblb
c
c.... set up the parameters
c
          iel    = lcblkb(1,iblk)
          lelCat = lcblkb(2,iblk)
          lcsyst = lcblkb(3,iblk)
          iorder = lcblkb(4,iblk)
          nenl   = lcblkb(5,iblk)  ! no. of vertices per element
          nenbl  = lcblkb(6,iblk)  ! no. of vertices per bdry. face
          mattyp = lcblkb(7,iblk)
          ndofl  = lcblkb(8,iblk)
          nshl   = lcblkb(9,iblk)
          nshlb  = lcblkb(10,iblk)
          npro   = lcblkb(1,iblk+1) - iel
          if(lcsyst.eq.3) lcsyst=nenbl
          ngaussb = nintb(lcsyst)

c
c.... compute and assemble the residuals corresponding to the
c     boundary integral
c

          allocate (tmpshpb(nshl,MAXQPT))
          allocate (tmpshglb(nsd,nshl,MAXQPT))
          if(lhs.eq.1 .and. iLHScond >= 1) then
             allocate (EGmass(npro,nshl,nshl))
             EGmass = zero
          else
             allocate (EGmass(1,1,1))
          endif

          tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:)
          tmpshglb(:,1:nshl,:) = shglb(lcsyst,:,1:nshl,:)

          call AsBMFG (y,                       x,
     &                 tmpshpb,                 tmpshglb,
     &                 mienb(iblk)%p,           mmatb(iblk)%p,
     &                 miBCB(iblk)%p,           mBCB(iblk)%p,
     &                 res,                     rmes,
     &                 EGmass)
          if(lhs == 1 .and. iLHScond > 0) then
            call fillSparseC_BC(mienb(iblk)%p, EGmass,
     &                   lhsk, row, col)
          endif

          deallocate (EGmass)
          deallocate (tmpshpb)
          deallocate (tmpshglb)
        enddo   !end of boundary element loop

!ifdef DEBUG !Nicholas Mati
!        call write_debug(myrank, 'res-afterAsBMFG'//char(0),
!     &                           'res'//char(0), res,
!     &                           'd'//char(0), nshg, nflow, lstep)
!        call MPI_ABORT(MPI_COMM_WORLD)
!endif //DEBUG


c
      ttim(80) = ttim(80) + secs(0.0)
c
c before the commu we need to rotate the residual vector for axisymmetric
c boundary conditions (so that off processor periodicity is a dof add instead
c of a dof combination).  Take care of all nodes now so periodicity, like
c commu is a simple dof add.
c
      if(iabc==1) then               !are there any axisym bc's
          call rotabc(res(1,2), iBC,  'in ')
c          Bdiagvec(:,1)=BDiag(:,1,1)
c          Bdiagvec(:,2)=BDiag(:,2,2)
c          Bdiagvec(:,3)=BDiag(:,3,3)
c          Bdiagvec(:,4)=BDiag(:,4,4)
c          Bdiagvec(:,5)=BDiag(:,5,5)
c          call rotabc(Bdiagvec(1,2), iBC,  'in ')
c          BDiag(:,:,:)=zero
c          BDiag(:,1,1)=Bdiagvec(:,1)
c          BDiag(:,2,2)=Bdiagvec(:,2)
c          BDiag(:,3,3)=Bdiagvec(:,3)
c          BDiag(:,4,4)=Bdiagvec(:,4)
c          BDiag(:,5,5)=Bdiagvec(:,5)
       endif

c.... -------------------->   communications <-------------------------
c

      if (numpe > 1) then
        call commu (res  , ilwork, nflow  , 'in ')

        call MPI_BARRIER (MPI_COMM_WORLD,ierr)

        if(iprec .ne. 0) call commu (BDiag, ilwork, nflow*nflow, 'in ')
      endif

c
c.... ---------------------->   post processing  <----------------------
c
c.... satisfy the BCs on the residual
c
      call bc3Res (y,  iBC,  BC,  res,  iper, ilwork)
c
c.... satisfy the BCs on the block-diagonal preconditioner
c
c  This code fragment would extract the "on processor diagonal
c      blocks". lhsK alread has the BC's applied to it (using BC3lhs),
c      though it was on an ebe basis. For now, we forgo this and still
c      form BDiag before BC3lhs, leaving the need to still apply BC's
c      below.  Keep in mind that if we used the code fragment below we
c      would still need to make BDiag periodic since BC3lhs did not do
c      that part.
c
      if (iprec .ne. 0) then
c$$$         do iaa=1,nshg
c$$$            k = sparseloc( row(col(iaa)), col(iaa+1)-colm(iaa), iaa )
c$$$     &       + colm(iaa)-1
c$$$            do idof=1,nflow
c$$$               do jdof=1,nflow
c$$$                  idx=idof+(jdof-1)*nflow
c$$$                  BDiag(iaa,idof,jdof)=lhsK(idx,k)
c$$$               enddo
c$$$            enddo
c$$$         enddo
         call bc3BDg (y,  iBC,  BC,  BDiag, iper, ilwork)
      endif
c
c.... return
c
      call timer ('Back    ')
      return
      end
c
c

c
        subroutine ElmGMRSclr(y,      ac,
     &                        x,      elDw,
     &                        shp,    shgl,   iBC,
     &                        BC,     shpb,   shglb,
     &                        rest,   rmest,  Diag,
     &                        iper,   ilwork, EGmasst)
c
c----------------------------------------------------------------------
c
c This routine computes the LHS mass matrix, the RHS residual
c vector, and the preconditioning matrix, for use with the GMRES
c solver.
c
c Zdenek Johan, Winter 1991.      (Fortran 90)
c Chris Whiting, Winter 1998.     (Matrix EBE-GMRES)
c----------------------------------------------------------------------
c
        use pointer_data
c
        include "common.h"
        include "mpif.h"
c
        dimension y(nshg,ndof),         ac(nshg,ndof),
     &            x(numnp,nsd),
     &            iBC(nshg),
     &            BC(nshg,ndofBC),
     &            rest(nshg),           Diag(nshg),
     &            rmest(nshg),          BDiag(nshg,nflow,nflow),
     &            iper(nshg),           EGmasst(numel,nshape,nshape)
c
        dimension shp(MAXTOP,maxsh,MAXQPT),
     &            shgl(MAXTOP,nsd,maxsh,MAXQPT),
     &            shpb(MAXTOP,maxsh,MAXQPT),
     &            shglb(MAXTOP,nsd,maxsh,MAXQPT)
c
        dimension qrest(nshg),          rmasst(nshg)
c
        dimension ilwork(nlwork)
        dimension elDw(numel)
c
        real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:)
        real*8, allocatable :: tmpshpb(:,:), tmpshglb(:,:,:)
        real*8, allocatable :: elDwl(:)
c
	ttim(80) = ttim(80) - tmr()
c
c.... set up the timer
c

        call timer ('Elm_Form')
c
c.... -------------------->   interior elements   <--------------------
c
c.... set up parameters
c
        intrul = intg  (1,itseq)
        intind = intpt (intrul)
c
        ires   = 1

c       if (idiff>=1) then  ! global reconstruction of qdiff
c
c loop over element blocks for the global reconstruction
c of the diffusive flux vector, q, and lumped mass matrix, rmass
c
c        qrest = zero
c        rmasst = zero
c
c        do iblk = 1, nelblk
c
c.... set up the parameters
c
c          iel    = lcblk(1,iblk)
c          lelCat = lcblk(2,iblk)
c          lcsyst = lcblk(3,iblk)
c          iorder = lcblk(4,iblk)
c          nenl   = lcblk(5,iblk)   ! no. of vertices per element
c          mattyp = lcblk(7,iblk)
c          ndofl  = lcblk(8,iblk)
c          nsymdl = lcblk(9,iblk)
c          npro   = lcblk(1,iblk+1) - iel
c
c          nintg  = numQpt (nsd,intrul,lcsyst)
c
c
c.... compute and assemble diffusive flux vector residual, qres,
c     and lumped mass matrix, rmass
c
c          call AsIq (y,                x,
c     &               shp(1,intind,lelCat),
c     &               shgl(1,intind,lelCat),
c     &               mien(iblk)%p,    mxmudmi(iblk)%p,
c     &               qres,           rmass )
c
c       enddo
c
c
c.... compute qi for node A, i.e., qres <-- qres/rmass
c
c       if (numpe > 1) then
c          call commu (qres  , ilwork, (ndof-1)*nsd  , 'in ')
c          call commu (rmass , ilwork,  1            , 'in ')
c       endif
c
c.... take care of periodic boundary conditions
c
c       call qpbc( rmass, qres, iBC, iper )
c
c       rmass = one/rmass
c
c       do i=1, (nflow-1)*nsd
c          qres(:,i) = rmass*qres(:,i)
c       enddo
c
c       if(numpe > 1) then
c          call commu (qres, ilwork, (nflow-1)*nsd, 'out')
c       endif
c
c      endif                     ! computation of global diffusive flux
c
c.... loop over element blocks to compute element residuals
c
c
c.... initialize the arrays
c
        rest    = zero
        rmest   = zero ! to avoid trap_uninitialized
        if (lhs .eq. 1)   EGmasst = zero
        if (iprec. ne. 0)   Diag  = zero
c
c.... loop over the element-blocks
c
        do iblk = 1, nelblk
c
c
          nenl   = lcblk(5,iblk) ! no. of vertices per element
          iel    = lcblk(1,iblk)
          lelCat = lcblk(2,iblk)
          lcsyst = lcblk(3,iblk)
          iorder = lcblk(4,iblk)
          nenl   = lcblk(5,iblk)   ! no. of vertices per element
          nshl   = lcblk(10,iblk)
          mattyp = lcblk(7,iblk)
          ndofl  = lcblk(8,iblk)
          nsymdl = lcblk(9,iblk)
          npro   = lcblk(1,iblk+1) - iel
          inum   = iel + npro - 1
          ngauss = nint(lcsyst)
c
c.... compute and assemble the residual and tangent matrix
c
          allocate (tmpshp(nshl,MAXQPT))
          allocate (tmpshgl(nsd,nshl,MAXQPT))

          tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:)
          tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:)

          allocate (elDwl(npro))
c
          call AsIGMRSclr(y,
     &                    ac,
     &                    x,               elDwl,
     &                    tmpshp,          tmpshgl,
     &                    mien(iblk)%p,
     &                    mmat(iblk)%p,    rest,
     &                    rmest,
     &                    qrest,           EGmasst(iel:inum,:,:),
     &                    Diag )
c

c.... satisfy the BC's on the implicit LHS
c
          call bc3LHSSclr (iBC, mien(iblk)%p, EGmasst(iel:inum,:,:) )
c
          elDw(iel:inum)=elDwl(1:npro)
          deallocate ( elDwl )
          deallocate ( tmpshp )
          deallocate ( tmpshgl )
c.... end of interior element loop
c
       enddo
c
c.... -------------------->   boundary elements   <--------------------
c
c.... set up parameters
c
        intrul = intg   (2,itseq)
        intind = intptb (intrul)
c
c.... loop over the boundary elements
c
        do iblk = 1, nelblb
c
c.... set up the parameters
c
          iel    = lcblkb(1,iblk)
          lelCat = lcblkb(2,iblk)
          lcsyst = lcblkb(3,iblk)
          iorder = lcblkb(4,iblk)
          nenl   = lcblkb(5,iblk)  ! no. of vertices per element
          nenbl  = lcblkb(6,iblk)  ! no. of vertices per bdry. face
          mattyp = lcblkb(7,iblk)
          ndofl  = lcblkb(8,iblk)
          nshl   = lcblkb(9,iblk)
          nshlb  = lcblkb(10,iblk)
          npro   = lcblkb(1,iblk+1) - iel
          if(lcsyst.eq.3) lcsyst=nenbl
          ngaussb = nintb(lcsyst)
c
c.... compute and assemble the residuals corresponding to the
c     boundary integral
c

          allocate (tmpshpb(nshl,MAXQPT))
          allocate (tmpshglb(nsd,nshl,MAXQPT))

          tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:)
          tmpshglb(:,1:nshl,:) = shglb(lcsyst,:,1:nshl,:)
c
          call AsBMFGSclr (y,                  x,
     &                     tmpshpb,
     &                     tmpshglb,
     &                     mienb(iblk)%p,      mmatb(iblk)%p,
     &                     miBCB(iblk)%p,      mBCB(iblk)%p,
     &                     rest,               rmest)
c
          deallocate ( tmpshpb )
          deallocate ( tmpshglb )

c.... end of boundary element loop
c
        enddo


      ttim(80) = ttim(80) + tmr()
c
c.... -------------------->   communications <-------------------------
c

      if (numpe > 1) then
        call commu (rest  , ilwork, 1  , 'in ')

        call MPI_BARRIER (MPI_COMM_WORLD,ierr)

        if(iprec .ne. 0) call commu (Diag, ilwork, 1, 'in ')
      endif

c
c.... ---------------------->   post processing  <----------------------
c
c.... satisfy the BCs on the residual
c
      call bc3ResSclr (y,  iBC,  BC,  rest,  iper, ilwork)
c
c.... satisfy the BCs on the preconditioner
c
      call bc3BDgSclr (iBC, Diag, iper, ilwork)
c
c.... return
c
      call timer ('Back    ')
      return
      end