precd4.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE precd4
!                   *****************
!
     &(x1,x2,a11,a12,a21,a22,
     & b1,b2,d1,d2,mesh,precon,prexsm,diadon)
!
!***********************************************************************
! BIEF   V7P1
!***********************************************************************
!
!brief    DIAGONAL PRECONDITIONING OF A SYSTEM A X = B
!+               (SEE EXPLANATIONS IN PRECDT).
!+
!+            A IS A 4-MATRIX BLOCK HERE.
!
!history  J-M HERVOUET (LNHE)
!+        06/07/2009
!+        V5P0
!+   First version (of the header...)
!
!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
!
!history  J-M HERVOUET (LNHE)
!+        08/12/2015
!+        V7P1
!+   Rebuilding the diagonals with MESH%IFAC in parallel in case of
!+   diagonal preconditioning and DIADON=FALSE.
!+   Correction of a bug in parallel with preconditioning 5.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| 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
!| D1             |<--| DIAGONAL MATRIX
!| D2             |<--| DIAGONAL MATRIX
!| DIADON         |-->| .TRUE. : DIAGONALS ARE GIVEN
!| MESH           |-->| MESH STRUCTURE
!| PRECON         |-->| CHOICE OF PRECONDITIONING
!| PREXSM         |-->| .TRUE. : PRECONDITIONING X1,X2 AND B1,B2
!| X1             |<->| FIRST INITIAL GUESS
!| X2             |-->| SECOND INITIAL GUESS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_precd4 => precd4
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: PRECON
      LOGICAL, INTENT(IN) :: PREXSM,DIADON
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!  VECTOR STRUCTURES
!
      TYPE(bief_obj), INTENT(INOUT) :: X1,X2,B1,B2,D1,D2
!
!-----------------------------------------------------------------------
!
!  MATRIX STRUCTURES
!
      TYPE(bief_obj), INTENT(INOUT) :: A11,A12,A21,A22
!
!-----------------------------------------------------------------------
!
!  MESH STRUCTURE
!
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I
!
!-----------------------------------------------------------------------
!
!  PREPARES THE DIAGONALS:
!
      IF(.NOT.diadon) THEN
!
!       COPY
!
        CALL os( 'X=Y     ' , x=d1 , y=a11%D )
        CALL os( 'X=Y     ' , x=d2 , y=a22%D )
!
!       PARALLEL MODE: COMPLETE DIAGONAL BEFORE GOING FURTHER
!
        IF(ncsize.GT.1) THEN
          CALL parcom(d1,2,mesh)
          CALL parcom(d2,2,mesh)
        ENDIF
!
!       POSSIBLY ABSOLUTE VALUES
!
        IF(precon.EQ.5) THEN
          CALL os( 'X=ABS(Y)' , x=d1 , y=d1 )
          CALL os( 'X=ABS(Y)' , x=d2 , y=d2 )
        ENDIF
!
!       SQUARE ROOT
!
        CALL os( 'X=SQR(Y)' , x=d1 , y=d1 )
        CALL os( 'X=SQR(Y)' , x=d2 , y=d2 )
!
!-----------------------------------------------------------------------
!                                                    -1
!  CHANGE OF VARIABLES (D1,D2 AND D3 ACTUALLY HOLD D1 ,...)
!
        IF(prexsm) THEN
          CALL os( 'X=XY    ' , x=x1 , y=d1 )
          CALL os( 'X=XY    ' , x=x2 , y=d2 )
        ENDIF
!
!-----------------------------------------------------------------------
!
!  COMPUTES THE INVERSE OF THE SQUARE ROOTS OF THE DIAGONALS
!  THIS GIVES BACK TRUE D1 AND D2 AND NOT D1 AND D2 INVERTED
!
        CALL os('X=1/Y   ', x=d1, y=d1, iopt=2,
     &                      infini=1.d0, zero=1.d-10)
        CALL os('X=1/Y   ', x=d2, y=d2, iopt=2,
     &                      infini=1.d0, zero=1.d-10)
!
      ELSE
!
!  CASE WHERE D IS GIVEN, CHANGE OF VARIABLES
!  CHANGE OF VARIABLE (D1,D2 REALLY HOLD D1,D2)
!
        IF(prexsm) THEN
          CALL os('X=Y/Z   ' , x=x1 , y=x1 , z=d1)
          CALL os('X=Y/Z   ' , x=x2 , y=x2 , z=d2)
        ENDIF
!
      ENDIF
!
!=======================================================================
! PRECONDITIONING OF A11 :
!=======================================================================
!
      CALL om('M=DMD   ', m=a11, d=d1, mesh=mesh)
!
!=======================================================================
! PRECONDITIONING OF A12 :
!=======================================================================
!
      CALL om('M=DM    ', m=a12, d=d1, mesh=mesh)
      CALL om('M=MD    ', m=a12, d=d2, mesh=mesh)
!
!=======================================================================
! PRECONDITIONING OF A21 :
!=======================================================================
!
      CALL om('M=DM    ', m=a21, d=d2, mesh=mesh)
      CALL om('M=MD    ', m=a21, d=d1, mesh=mesh)
!
!=======================================================================
! PRECONDITIONING OF A22 :
!=======================================================================
!
      CALL om('M=DMD   ', m=a22, d=d2, mesh=mesh)
!
!=======================================================================
!
!     CASES WHERE THE DIAGONALS ARE KNOWN
!     (VALID ONLY WITH ONE SINGLE DOMAIN)
!
      IF(ncsize.LE.1.OR.nptir.EQ.0) THEN
!
!       IF PRECON = 2 OR 3
        IF(2*(precon/2).EQ.precon.AND..NOT.diadon) THEN
          a11%TYPDIA='I'
          a22%TYPDIA='I'
        ELSEIF(3*(precon/3).EQ.precon.AND..NOT.diadon) THEN
          a11%TYPDIA='I'
          a22%TYPDIA='I'
          a12%TYPDIA='0'
          a21%TYPDIA='0'
        ENDIF
!
      ELSE
!
!       CASE OF DIAGONAL=IDENTITY, BUT ONLY AFTER ASSEMBLING
!
        IF((2*(precon/2).EQ.precon.OR.3*(precon/3).EQ.precon).AND.
     &                                                .NOT.diadon) THEN
!         HERE THE DIAGONAL IS REDONE WITH IFAC, SO THAT A MATRIX-VECTOR
!         PRODUCT WILL LEAD TO A SUM OF THE SAME NUMBERS (BUT POSSIBLY
!         NOT IN THE SAME ORDER). OTHERWISE IT IS NOT SURE THAT THE
!         ASSEMBLED DIAGONAL WOULD GIVE EXACT 1.D0.
          DO i=1,a11%D%DIM1
            a11%D%R(i)=mesh%IFAC%I(i)
          ENDDO
          DO i=1,a22%D%DIM1
            a22%D%R(i)=mesh%IFAC%I(i)
          ENDDO
        ENDIF
!
      ENDIF
!
!=======================================================================
!
! PRECONDITIONING OF THE SECOND MEMBER
!
      IF(prexsm) THEN
        CALL os( 'X=XY    ' , x=b1 , y=d1 )
        CALL os( 'X=XY    ' , x=b2 , y=d2 )
      ENDIF
!
!=======================================================================
!
      RETURN
      END