prebd9.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE prebd9
!                   *****************
!
     &(x1,x2,x3,a11,a12,a13,a21,a22,a23,a31,a32,a33,
     & b1,b2,b3,d11,d12,d13,d21,d22,d23,d31,d32,d33,
     & 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
!| ...            |<->| ...
!| A33            |<->| TERM (3,3) OF MATRIX
!| B1             |<->| FIRST RIGHT-HAND SIDE
!| B2             |<->| SECOND RIGHT-HAND SIDE
!| B3             |<->| THIRD RIGHT-HAND SIDE
!| D11            |<--| DIAGONAL MATRIX
!| ...            |<--| ...
!| D33            |<--| 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
!| X3             |-->| THIRD INITIAL GUESS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_prebd9 => prebd9
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      LOGICAL, INTENT(IN) :: PREXSM,DIADON
!
!-----------------------------------------------------------------------
!
!  VECTOR STRUCTURES
!
      TYPE(bief_obj), INTENT(IN)    :: X3,B1
      TYPE(bief_obj), INTENT(INOUT) :: X1,X2,B2,B3
      TYPE(bief_obj), INTENT(INOUT) :: D11,D12,D13,D21,D22,D23
      TYPE(bief_obj), INTENT(INOUT) :: D31,D32,D33
!
!-----------------------------------------------------------------------
!
!  MATRIX STRUCTURES
!
      TYPE(bief_obj), INTENT(INOUT) :: A11,A12,A13,A21,A22
      TYPE(bief_obj), INTENT(INOUT) :: A23,A31,A32,A33
!
!-----------------------------------------------------------------------
!
!  MESH STRUCTURE
!
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,NPOIN1,NPOIN2,NPOIN3
!
!-----------------------------------------------------------------------
!
      npoin1 = x1%DIM1
      npoin2 = x2%DIM1
      npoin3 = x3%DIM1
!
      IF(npoin2.NE.npoin1.AND.npoin3.NE.npoin1) THEN
        WRITE(lu,200)
200     FORMAT(1x,'PREBD9 (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=d13, y=a13%D)
        CALL os('X=Y     ', x=d21, y=a21%D)
        CALL os('X=Y     ', x=d22, y=a22%D)
        CALL os('X=Y     ', x=d23, y=a23%D)
        CALL os('X=Y     ', x=d31, y=a31%D)
        CALL os('X=Y     ', x=d32, y=a32%D)
        CALL os('X=Y     ', x=d33, y=a33%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)
        d31%R(i) =  d31%R(i) * d11%R(i)
        d22%R(i) =  d22%R(i) - d21%R(i) * d12%R(i)
!
      ENDDO
!
!     ONLY D22 INVERTED IS NOW USED
      CALL os('X=1/Y   ', x=d22, y=d22)
!
      DO i = 1,npoin1
!
        d32%R(i) = (d32%R(i) - d31%R(i) * d12%R(i)) * d22%R(i)
        d23%R(i) =  d23%R(i) - d21%R(i) * d13%R(i)
        d33%R(i) =  d33%R(i)
     &             -d31%R(i)*d13%R(i)-d32%R(i)*d23%R(i)
        d12%R(i) =  d12%R(i) * d11%R(i)
        d13%R(i) =  d13%R(i) * d11%R(i)
        d23%R(i) =  d23%R(i) * d22%R(i)
!
      ENDDO
!
!-----------------------------------------------------------------------
!
! CHANGE OF VARIABLES:
!
      IF(prexsm) THEN
!
        CALL os('X=X+YZ  ', x=x1, y=x2, z=d12)
        CALL os('X=X+YZ  ', x=x1, y=x3, z=d13)
        CALL os('X=X+YZ  ', x=x2, y=x3, z=d23)
!
      ENDIF
!
!  COMPUTES THE SQUARE ROOT
!  INVERTS D11,D22,D33
!  (THEY ARE ONLY USED IN THIS FORM FROM NOW ON)
!
!     INVERSION OF D11 ALREADY PERFORMED
!     INVERSION OF D22 ALREADY PERFORMED
      CALL os('X=1/Y   ', x=d33, y=d33)
      CALL os('X=SQR(Y)', x=d11, y=d11)
      CALL os('X=SQR(Y)', x=d22, y=d22)
      CALL os('X=SQR(Y)', x=d33, y=d33)
!
!=======================================================================
! 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)
!  A23 :
      CALL om('M=M-DN  ', m=a23, n=a13, d=d21, mesh=mesh)
!  A31 :
      CALL om('M=M-DN  ', m=a31, n=a11, d=d31, mesh=mesh)
      CALL om('M=M-DN  ', m=a31, n=a21, d=d32, mesh=mesh)
!  A32 :
      CALL om('M=M-DN  ', m=a32, n=a12, d=d31, mesh=mesh)
      CALL om('M=M-DN  ', m=a32, n=a22, d=d32, mesh=mesh)
!  A33 :
      CALL om('M=M-DN  ', m=a33, n=a13, d=d31, mesh=mesh)
      CALL om('M=M-DN  ', m=a33, n=a23, d=d32, 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)
!  A32 :
      CALL om('M=M-ND  ', m=a32, n=a31, d=d12, mesh=mesh)
!  A13 :
      CALL om('M=M-ND  ', m=a13, n=a11, d=d13, mesh=mesh)
      CALL om('M=M-ND  ', m=a13, n=a12, d=d23, mesh=mesh)
!  A23 :
      CALL om('M=M-ND  ', m=a23, n=a21, d=d13, mesh=mesh)
      CALL om('M=M-ND  ', m=a23, n=a22, d=d23, mesh=mesh)
!  A33 :
      CALL om('M=M-ND  ', m=a33, n=a31, d=d13, mesh=mesh)
      CALL om('M=M-ND  ', m=a33, n=a32, d=d23, mesh=mesh)
!
!-----------------------------------------------------------------------
!
! NEW SECOND MEMBER
!
      IF(prexsm) THEN
!
      DO i = 1,npoin1
        b2%R(i) = b2%R(i)-d21%R(i)*b1%R(i)
        b3%R(i) = b3%R(i)-d31%R(i)*b1%R(i)-d32%R(i)*b2%R(i)
      ENDDO
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END