prebd4.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE prebd4
!                   *****************
!
     &(x1,x2,a11,a12,a21,a22,b1,b2,d11,d12,d21,d22,
     & mesh,prexsm,diadon)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    BLOCK-DIAGONAL PRECONDITIONING OF A SYSTEM A X = B.
!
!history  J.M. HERVOUET (LNH)
!+        23/12/94
!+        V5P1
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A11            |<->| TERM (1,1) OF MATRIX
!| A12            |<->| TERM (1,2) OF MATRIX
!| A21            |<->| TERM (2,1) OF MATRIX
!| A22            |<->| TERM (2,2) OF MATRIX
!| B1             |<->| FIRST RIGHT-HAND SIDE
!| B2             |<->| SECOND RIGHT-HAND SIDE
!| D11            |<--| DIAGONAL MATRIX
!| D12            |<--| DIAGONAL MATRIX
!| D21            |<--| DIAGONAL MATRIX
!| D22            |<--| DIAGONAL MATRIX
!| DIADON         |-->| .TRUE. : DIAGONALS ARE GIVEN
!| MESH           |-->| MESH STRUCTURE
!| PREXSM         |-->| .TRUE. : PRECONDITIONING X1,X2 AND B1,B2
!| X1             |<->| FIRST INITIAL GUESS
!| X2             |-->| SECOND INITIAL GUESS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_prebd4 => prebd4
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      LOGICAL, INTENT(IN) :: PREXSM,DIADON
!
!-----------------------------------------------------------------------
!
!  VECTOR STRUCTURES
!
      TYPE(bief_obj), INTENT(INOUT) :: X1,B2,D11,D12,D21,D22
      TYPE(bief_obj), INTENT(IN)    :: X2,B1
!
!-----------------------------------------------------------------------
!
!  MATRIX STRUCTURES
!
      TYPE(bief_obj), INTENT(INOUT) :: A11,A12,A21,A22
!
!-----------------------------------------------------------------------
!
!  MESH STRUCTURE
!
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,NPOIN1,NPOIN2
!
!-----------------------------------------------------------------------
!
      npoin1 = x1%DIM1
      npoin2 = x2%DIM1
!
      IF(npoin2.NE.npoin1) THEN
        WRITE(lu,200)
200     FORMAT(1x,'PREBD4 (BIEF) : RECTANGULAR MATRICES',/,1x,
     &  'BLOCK-DIAGONAL PRECONDITIONING IMPOSSIBLE IN THIS CASE')
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
!  PREPARES THE DIAGONALS:
!
      IF(.NOT.diadon) THEN
!
        CALL os('X=Y     ', x=d11, y=a11%D)
        CALL os('X=Y     ', x=d12, y=a12%D)
        CALL os('X=Y     ', x=d21, y=a21%D)
        CALL os('X=Y     ', x=d22, y=a22%D)
!
      ENDIF
!
!-----------------------------------------------------------------------
!
!  L D U FACTORISATION OF THE DIAGONAL BLOCK:
!
!     ONLY D11 INVERTED IS NOW USED
      CALL os('X=1/Y   ', x=d11, y=d11)
!
      DO i = 1,npoin1
!
        d21%R(i) = d21%R(i) * d11%R(i)
        d22%R(i) = d22%R(i) - d21%R(i) * d12%R(i)
        d12%R(i) = d12%R(i) * d11%R(i)
!
      ENDDO
!
!-----------------------------------------------------------------------
!
! CHANGE OF VARIABLES:
!
      IF(prexsm) THEN
!
        CALL os('X=X+YZ  ', x=x1, y=x2, z=d12)
!
      ENDIF
!
!  COMPUTES THE SQUARE ROOT
!  INVERTS D11,D22,D33
!  (THEY ARE ONLY USED IN THIS FORM FROM NOW ON)
!
!     INVERSION OF D11 ALREADY PERFORMED
      CALL os('X=1/Y   ', x=d22, y=d22)
      CALL os('X=SQR(Y)', x=d11, y=d11)
      CALL os('X=SQR(Y)', x=d22, y=d22)
!
!=======================================================================
! MULTIPLIES A ON THE LEFT BY L INVERTED
!=======================================================================
!
! A21 :
      CALL om('M=M-DN  ', m=a21, n=a11, d=d21, mesh=mesh)
! A22 :
      CALL om('M=M-DN  ', m=a22, n=a12, d=d21, mesh=mesh)
!
!=======================================================================
! MULTIPLIES A ON THE RIGHT BY U INVERTED
!=======================================================================
!
! A12 :
      CALL om('M=M-ND  ', m=a12, n=a11, d=d12, mesh=mesh)
! A22 :
      CALL om('M=M-ND  ', m=a22, n=a21, d=d12, mesh=mesh)
!
!-----------------------------------------------------------------------
!
! NEW SECOND MEMBER
!
      IF(prexsm) THEN
!
      DO i = 1,npoin1
        b2%R(i) = b2%R(i) - d21%R(i) * b1%R(i)
      ENDDO
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END