precd1.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE precd1
!                   *****************
!
     &(x,a,b,d,mesh,precon,prexsm,diadon)
!
!***********************************************************************
! BIEF   V7P1
!***********************************************************************
!
!brief    DIAGONAL PRECONDITIONING OF A SYSTEM A X = B
!+               (SEE EXPLANATIONS IN PRECDT).
!+
!+            A IS A SIMPLE MATRIX 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 diagonal with MESH%IFAC in parallel in case of
!+   diagonal preconditioning and DIADON=FALSE.
!+   Correction of a bug in parallel with preconditioning 5.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A              |-->| BLOCK OF MATRICES
!| B              |-->| BLOCK OF RIGHT-HAND SIZES
!| D              |<--| BLOCK OF DIAGONALS
!| DIADON         |-->| .TRUE. : DIAGONALS ARE GIVEN
!| MESH           |-->| MESH STRUCTURE
!| PRECON         |-->| CHOICE OF PRECONDITIONING
!| PREXSM         |-->| .TRUE. : PRECONDITIONING X AND B
!| X              |<->| BLOCK OF UNKNOWN VECTORS IN THE SYSTEM
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_precd1 => precd1
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: PRECON
      LOGICAL, INTENT(IN) :: PREXSM,DIADON
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!  VECTOR STRUCTURES
!
      TYPE(bief_obj), INTENT(INOUT) :: X,B,D
!
!-----------------------------------------------------------------------
!
!  MATRIX STRUCTURES
!
      TYPE(bief_obj), INTENT(INOUT) :: A
!
!-----------------------------------------------------------------------
!
!  MESH STRUCTURE
!
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I
!
!-----------------------------------------------------------------------
!
!  PREPARES THE DIAGONALS:
!
      IF(.NOT.diadon) THEN
!
        CALL os( 'X=Y     ' , x=d , y=a%D )
!
!       PARALLEL MODE: COMPLETE DIAGONAL BEFORE GOING FURTHER
!
        IF(ncsize.GT.1) THEN
          CALL parcom(d,2,mesh)
        ENDIF
!
!       POSSIBLY TAKING THE ABSOLUTE VALUES
!
        IF(precon.EQ.5) THEN
          CALL os( 'X=ABS(Y)' , x=d , y=d )
        ENDIF
!
!       TAKING THE SQUARE ROOT
!
        CALL os( 'X=SQR(Y)' , x=d , y=d )
!
!-----------------------------------------------------------------------
!                                         -1
!  CHANGE OF VARIABLES (D ACTUALLY HOLDS D  )
!
        IF(prexsm) CALL os( 'X=XY    ' , x=x , y=d )
!
!-----------------------------------------------------------------------
!
!  COMPUTES THE INVERSE OF THE SQUARE ROOTS OF THE DIAGONALS
!  THIS GIVES BACK TRUE D AND NOT D INVERTED
!
        CALL os('X=1/Y   ', x=d, y=d, iopt=2, infini=1.d0, zero=1.d-10)
!
      ELSE
!
!  CASE WHERE D IS GIVEN, CHANGE OF VARIABLES
!  CHANGE OF VARIABLE (D REALLY HOLDS D)
!
        IF(prexsm) THEN
          CALL os('X=Y/Z   ' , x=x , y=x , z=d)
        ENDIF
!
      ENDIF
!
!=======================================================================
! PRECONDITIONING OF A:
!=======================================================================
!
      CALL om('M=DMD   ', m=a, d=d, mesh=mesh)
!     IF PRECON = 2 OR 3
      IF((2*(precon/2).EQ.precon.OR.3*(precon/3).EQ.precon).AND.
     &                                                 .NOT.diadon) THEN
!       VALID ONLY WITH ONE SINGLE DOMAIN
        IF(ncsize.LE.1.OR.nptir.EQ.0) THEN
          a%TYPDIA='I'
        ELSE
!         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,a%D%DIM1
            a%D%R(i)=mesh%IFAC%I(i)
          ENDDO
        ENDIF
      ENDIF
!
!=======================================================================
!
! PRECONDITIONING OF THE SECOND MEMBER
!
      IF(prexsm) CALL os( 'X=XY    ' , x=b , y=d )
!
!=======================================================================
!
      RETURN
      END