PRACE NPT

      Subroutine sym7pre( m, n1, n2, n3, a, b, fac )

! ---------------------------------------------------------------------

! --- sym7pre computes the inverse of D as used in the preconditioners.

!     Also the elements of A and RHS b are adjusted.

!     Note that Inv(D) is approximate as no elements of A of

!     other processors are used in order to keep communication

!     within reasonable bounds.

! ---------------------------------------------------------------------

      Use         numerics

      Use         mpi_module

      Use         floptime

      Implicit    None

      Integer  :: m, n1 ,n2 ,n3

      Real(l_) :: a(m,0:3), b(n1*n2*n3), fac(n1*n2*n3)

      Integer  :: i, inm, in12, ntot, n12

! ---------------------------------------------------------------------

      n12  = n1*n2

      ntot = n12*n3

      fac(lb) = 1.0_l_/a(1,0)

      Do i = 2, m

         inm    = i - 1

         fac(glb+i) = a(i,0) - fac(glb+inm)*a(inm,1)*a(inm,1)

      End Do

      Do i = n1 + 1, m

         inm    = i - n1

         fac(glb+i) = fac(glb+i) - fac(glb+inm)*a(inm,2)*a(inm,2)

      End Do

      Do i = n12 + 1, m

         inm    = i - n12

         fac(glb+i) = fac(glb+i) - fac(glb+inm)*a(inm,3)*a(inm,3)

      End Do

      fac(lb:gub) = 1.0_l_/fac(lb:gub)

! ---------------------------------------------------------------------

! --- D-inverse is now complete, scale A and b.

      a(1:m,0) = a(1:m,0)*fac(lb:gub)

      fac(lb:gub) = Sqrt( fac(lb:gub) )

      b(lb:gub)   = b(lb:gub)*fac(lb:gub)

      Do i = 2, m

         inm      = i - 1

         a(inm,1) = fac(glb+i)*fac(glb+inm)*a(inm,1)

      End Do

      Do i = n1 + 1, m

         inm      = i - n1

         a(inm,2) = fac(glb+i)*fac(glb+inm)*a(inm,2)

      End Do

      Do i = n12 + 1, m

         inm     = i - n12

         a(inm,3) = fac(glb+i)*fac(glb+inm)*a(inm,3)

      End Do

      Call MPI_Allgatherv( fac(lb), m, rtyp, fac, sizes, offset, rtyp,

     &                     comm, ierr )

      Call MPI_Allgatherv( b(lb), m, rtyp, b, sizes, offset, rtyp,

     &                     comm, ierr )

      flops = flops + 25*m - 6*( n1 + n12 ) - 5 ! <-- No. of flops.

! ---------------------------------------------------------------------

! --- Keep fac for backscaling the solution x.

! ---------------------------------------------------------------------

      End Subroutine sym7pre