xref: /phasta/phSolver/compressible/i3lu.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)

        subroutine i3LU (Diag, r, code)
c
c----------------------------------------------------------------------
c
c This routine performs a LU factorization/solve of a set of matrices,
c used for block diagonal preconditioning in the iterative driver.

c input:
c  Diag (nshg,nflow,nflow)  : block diagonal (symmetric storage)
c  r    (nshg,nflow)    : residual
c  code                 : operation code
c                           .eq. 'LU_Fact ', Cholesky Factor
c                           .eq. 'forward ', forward reduction
c                           .eq. 'backward', backward substitution
c                           .eq. 'product ', product Diag.r
c
c output:
c  Diag (nshg,nflow,nflow)  : LU decomp. of block diagonal
c  r    (nshg,nflow)    : reduced residual
c
c
c Note: the formulation used here is taken from Golub's "Matrix
c       Computations" Book (1984), pages 82-85 algorithm P5.1-3.
c
c       L and U overwrite Diag
c
c       The diagonal terms (i,i) of U are stored in inverted form.
c
c
c----------------------------------------------------------------------
c
        include "common.h"
c
        dimension Diag(nshg,nflow,nflow),        r(nshg,nflow)
c
        character*8 code

c
c.... perform LU decomposition with the Diagonal terms inverted
c
        if (code .eq. 'LU_Fact ') then
c
          Diag(:,1,1) = one   / Diag(:,1,1)
c
          Diag(:,2,1) = Diag(:,1,1) *  Diag(:,2,1)
          Diag(:,3,1) = Diag(:,1,1) *  Diag(:,3,1)
          Diag(:,4,1) = Diag(:,1,1) *  Diag(:,4,1)
          Diag(:,5,1) = Diag(:,1,1) *  Diag(:,5,1)
c
          Diag(:,2,2) = Diag(:,2,2) - Diag(:,2,1) *  Diag(:,1,2)
          Diag(:,2,3) = Diag(:,2,3) - Diag(:,2,1) *  Diag(:,1,3)
          Diag(:,2,4) = Diag(:,2,4) - Diag(:,2,1) *  Diag(:,1,4)
          Diag(:,2,5) = Diag(:,2,5) - Diag(:,2,1) *  Diag(:,1,5)
          Diag(:,2,2) = one   / Diag(:,2,2)
c
          Diag(:,3,2) = Diag(:,2,2)*(Diag(:,3,2)
     &                             - Diag(:,3,1) * Diag(:,1,2))
          Diag(:,4,2) = Diag(:,2,2)*(Diag(:,4,2)
     &                             - Diag(:,4,1) * Diag(:,1,2))
          Diag(:,5,2) = Diag(:,2,2)*(Diag(:,5,2)
     &                             - Diag(:,5,1) * Diag(:,1,2))
c

          Diag(:,3,3) = Diag(:,3,3) - Diag(:,3,1) * Diag(:,1,3)
     &                              - Diag(:,3,2) * Diag(:,2,3)
          Diag(:,3,4) = Diag(:,3,4) - Diag(:,3,1) * Diag(:,1,4)
     &                              - Diag(:,3,2) * Diag(:,2,4)
          Diag(:,3,5) = Diag(:,3,5) - Diag(:,3,1) * Diag(:,1,5)
     &                              - Diag(:,3,2) * Diag(:,2,5)
          Diag(:,3,3) = one   / Diag(:,3,3)
c
          Diag(:,4,3) = Diag(:,3,3) *(Diag(:,4,3)
     &                              - Diag(:,4,1) * Diag(:,1,3)
     &                              - Diag(:,4,2) * Diag(:,2,3))
          Diag(:,5,3) = Diag(:,3,3) *(Diag(:,5,3)
     &                              - Diag(:,5,1) * Diag(:,1,3)
     &                              - Diag(:,5,2) * Diag(:,2,3))
c
          Diag(:,4,4) = Diag(:,4,4) - Diag(:,4,1) * Diag(:,1,4)
     &                              - Diag(:,4,2) * Diag(:,2,4)
     &                              - Diag(:,4,3) * Diag(:,3,4)
          Diag(:,4,4) = one / Diag(:,4,4)
c
          Diag(:,5,4) = Diag(:,4,4) *(Diag(:,5,4)
     &                              - Diag(:,5,1) * Diag(:,1,4)
     &                              - Diag(:,5,2) * Diag(:,2,4)
     &                              - Diag(:,5,3) * Diag(:,3,4))
c
          Diag(:,5,5) = Diag(:,5,5) - Diag(:,5,1) * Diag(:,1,5)
     &                              - Diag(:,5,2) * Diag(:,2,5)
     &                              - Diag(:,5,3) * Diag(:,3,5)
     &                              - Diag(:,5,4) * Diag(:,4,5)
          Diag(:,5,5) = one / Diag(:,5,5)
c
c  INACCURATE NOW          flops = flops + 110*nshg
c
          return
        endif
c
c.... perform forward reduction
c
        if (code .eq. 'forward ') then
c
c         r(:,1) =  r(:,1) !no-op
          r(:,2) =  r(:,2) - Diag(:,2,1) *  r(:,1)
          r(:,3) =  r(:,3) - Diag(:,3,1) *  r(:,1)
     &                     - Diag(:,3,2) *  r(:,2)
          r(:,4) =  r(:,4) - Diag(:,4,1) *  r(:,1)
     &                     - Diag(:,4,2) *  r(:,2)
     &                     - Diag(:,4,3) *  r(:,3)
          r(:,5) =  r(:,5) - Diag(:,5,1) *  r(:,1)
     &                     - Diag(:,5,2) *  r(:,2)
     &                     - Diag(:,5,3) *  r(:,3)
     &                     - Diag(:,5,4) *  r(:,4)
c
c.... flop count
c
c  INACCURATE           flops = flops + 25*nshg
c
          return
        endif
c
c.... perform backward substitution
c
        if (code .eq. 'backward') then
c
          r(:,5) = Diag(:,5,5) *   r(:,5)
          r(:,4) = Diag(:,4,4) * ( r(:,4)
     &                         -   r(:,5) * Diag(:,4,5) )
          r(:,3) = Diag(:,3,3) * ( r(:,3)
     &                         -   r(:,5) * Diag(:,3,5)
     &                         -   r(:,4) * Diag(:,3,4) )
          r(:,2) = Diag(:,2,2) * ( r(:,2)
     &                         -   r(:,5) * Diag(:,2,5)
     &                         -   r(:,4) * Diag(:,2,4)
     &                         -   r(:,3) * Diag(:,2,3) )
          r(:,1) = Diag(:,1,1) * ( r(:,1)
     &                         -   r(:,5) * Diag(:,1,5)
     &                         -   r(:,4) * Diag(:,1,4)
     &                         -   r(:,3) * Diag(:,1,3)
     &                         -   r(:,2) * Diag(:,1,2) )
c
c.... flop count
c
          flops = flops + 25*nshg

          return
        endif
c
c.... perform product U.r
c
        if (code .eq. 'product ') then
c
          r(:,1) = r(:,1) / Diag(:,1,1) + r(:,2) * Diag(:,1,2) +
     &             r(:,3) * Diag(:,1,3) + r(:,4) * Diag(:,1,4) +
     &             r(:,5) * Diag(:,1,5)
          r(:,2) = r(:,2) / Diag(:,2,2) +
     &             r(:,3) * Diag(:,2,3) + r(:,4) * Diag(:,2,4) +
     &             r(:,5) * Diag(:,2,5)
          r(:,3) = r(:,3) / Diag(:,3,3) +
     &             r(:,4) * Diag(:,3,4) +
     &             r(:,5) * Diag(:,3,5)
          r(:,4) = r(:,4) / Diag(:,4,4) +
     &             r(:,5) * Diag(:,4,5)
          r(:,5) = r(:,5) / Diag(:,5,5)
c
c.... flop count
c
          flops = flops + 40*nshg

          return
        endif
c
        call error ('i3LU    ', code, 0)
c
c.... return
c
c
111     format(5(e14.7,2x))
        return
        end