ovd.f

Go to the documentation of this file.

!                   **************
                    SUBROUTINE ovd
!                   **************
!
     & ( op , x , y , z , c , npoin , iopt , d , eps )
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    OPERATIONS ON VECTORS INCLUDING DIVISIONS
!+                DIVISION BY 0 CAN BE TESTED.
!+
!+            IN THE EVENT OF A DIVIDE CHECK, CAN EITHER STOP THE
!+                PROGRAM OR SET THE RESULT OF THE OPERATION TO
!+                A VALUE: D.
!code
!+   OP IS A STRING OF 8 CHARACTERS, WHICH INDICATES THE OPERATION TO BE
!+   PERFORMED ON VECTORS X,Y AND Z AND CONSTANT C.
!+
!+   THE RESULT IS VECTOR X.
!+
!+   OP = 'X=1/Y   '     :  COPIES INVERSE OF Y IN X
!+   OP = 'X=C/Y   '     :  DIVIDES C BY Y
!+   OP = 'X=Y/Z   '     :  DIVIDES Y BY Z
!+   OP = 'X=CY/Z  '     :  DIVIDES C.Y BY Z
!+   OP = 'X=CXY/Z '     :  DIVIDES C.X.Y BY Z
!+   OP = 'X=X+CY/Z'     :  ADDS C.Y/Z TO X
!
!warning  DIVIDE OPERATIONS INTERNALLY TAKE CARE OF DIVISIONS BY 0.
!+            SUCCESSFUL EXIT OF OVD IS THEREFORE NOT A PROOF THAT Y
!+            OR Z NEVER ARE 0
!
!history  J-M HERVOUET (LNH)    ; F  LEPEINTRE (LNH)
!+        26/11/93
!+        V5P2
!+
!
!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
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| C              |-->| A GIVEN CONSTANT
!| D              |-->| A DIAGONAL MATRIX
!| EPS            |-->| THRESHOLD TO AVOID DIVISIONS BY ZERO
!| IOPT           |-->| OPTION FOR DIVISIONS BY ZERO
!|                |   | 1: NO TEST DONE (WILL CRASH IF DIVISION BY 0.).
!|                |   | 2: INFINITE TERMS REPLACED BY CONSTANT INFINI.
!|                |   | 3: STOP IF DIVISION BY ZERO.
!|                |   | 4: DIVISIONS BY 0. REPLACED BY DIVISIONS/ZERO
!|                |   |    ZERO BEING AN OPTIONAL ARGUMENT
!| NPOIN          |-->| SIZE OF VECTORS
!| OP             |-->| STRING INDICATING THE OPERATION TO BE DONE
!| X              |<--| RESULTING VECTOR
!| Y              |-->| TO BE USED IN THE OPERATION
!| Z              |-->| TO BE USED IN THE OPERATION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER,          INTENT(IN)    :: NPOIN,IOPT
      DOUBLE PRECISION, INTENT(INOUT) :: X(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: Y(npoin),Z(npoin),C,D,EPS
      CHARACTER(LEN=8), INTENT(IN)    :: OP
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I
!
!-----------------------------------------------------------------------
!
      IF(op(1:8).EQ.'X=1/Y   ') THEN
!
        IF(iopt.EQ.1) THEN
!
        DO i=1,npoin
            x(i) = 1.d0/y(i)
        ENDDO ! I
!
        ELSEIF(iopt.EQ.2) THEN
!
        DO i=1,npoin
!
          IF (abs(y(i)).GT.eps) THEN
            x(i) = 1.d0/y(i)
          ELSE
            x(i) = d
          ENDIF
!
        ENDDO ! I
!
        ELSEIF(iopt.EQ.3) THEN
!
        DO i=1,npoin
!
          IF (abs(y(i)).GT.eps) THEN
            x(i) = 1.d0/y(i)
          ELSE
            WRITE(lu,2000) i,op,eps
            CALL plante(1)
            stop
          ENDIF
!
        ENDDO ! I
!
        ELSEIF(iopt.EQ.4) THEN
!
        DO i=1,npoin
!
          IF (abs(y(i)).GT.eps) THEN
            x(i) = 1.d0/y(i)
          ELSEIF (y(i).GE.0.d0) THEN
            x(i) =  1.d0/eps
          ELSE
            x(i) = -1.d0/eps
          ENDIF
!
        ENDDO ! I
!
        ENDIF
!
!-----------------------------------------------------------------------
!
      ELSEIF(op(1:8).EQ.'X=C/Y   ') THEN
!
        IF(iopt.EQ.1) THEN
!
        DO i=1,npoin
            x(i) = c/y(i)
        ENDDO ! I
!
        ELSEIF(iopt.EQ.2) THEN
!
        DO i=1,npoin
!
          IF (abs(y(i)).GT.eps) THEN
            x(i) = c/y(i)
          ELSE
            x(i) = d
          ENDIF
!
        ENDDO ! I
!
        ELSEIF(iopt.EQ.3) THEN
!
        DO i=1,npoin
!
          IF (abs(y(i)).GT.eps) THEN
            x(i) = c/y(i)
          ELSE
            WRITE(lu,2000) i,op,eps
            CALL plante(1)
            stop
          ENDIF
!
        ENDDO ! I
!
        ELSEIF(iopt.EQ.4) THEN
!
        DO i=1,npoin
!
          IF (abs(y(i)).GT.eps) THEN
            x(i) = c/y(i)
          ELSEIF (y(i).GE.0.d0) THEN
            x(i) =  1.d0/eps
          ELSE
            x(i) = -1.d0/eps
          ENDIF
!
        ENDDO ! I
!
        ENDIF
!
!-----------------------------------------------------------------------
!
      ELSEIF(op(1:8).EQ.'X=Y/Z   ') THEN
!
        IF(iopt.EQ.1) THEN
!
        DO i=1,npoin
            x(i) = y(i) / z(i)
        ENDDO ! I
!
        ELSEIF(iopt.EQ.2) THEN
!
        DO i=1,npoin
!
          IF (abs(z(i)).GT.eps) THEN
            x(i) = y(i) / z(i)
          ELSE
            x(i) = d
          ENDIF
!
        ENDDO ! I
!
        ELSEIF(iopt.EQ.3) THEN
!
        DO i=1,npoin
!
          IF (abs(z(i)).GT.eps) THEN
            x(i) = y(i) / z(i)
          ELSE
            WRITE(lu,2000) i,op,eps
            CALL plante(1)
            stop
          ENDIF
!
        ENDDO ! I
        ELSEIF(iopt.EQ.4) THEN
!
        DO i=1,npoin
!
          IF (abs(z(i)).GT.eps) THEN
            x(i) = y(i) / z(i)
          ELSEIF (abs(y(i)).LT.eps) THEN
            x(i) = d
          ELSEIF (z(i).GE.0.d0) THEN
            x(i) =  y(i) / eps
          ELSE
            x(i) = -y(i) / eps
          ENDIF
!
        ENDDO ! I
!
        ENDIF
!
!-----------------------------------------------------------------------
!
      ELSEIF(op(1:8).EQ.'X=CY/Z  ') THEN
!
        IF(iopt.EQ.1) THEN
!
        DO i=1,npoin
            x(i) = c*y(i) / z(i)
        ENDDO ! I
!
        ELSEIF(iopt.EQ.2) THEN
!
        DO i=1,npoin
!
          IF (abs(z(i)).GT.eps) THEN
            x(i) = c*y(i) / z(i)
          ELSE
            x(i) = d
          ENDIF
!
        ENDDO ! I
!
        ELSEIF(iopt.EQ.3) THEN
!
        DO i=1,npoin
!
          IF (abs(z(i)).GT.eps) THEN
            x(i) = c*y(i) / z(i)
          ELSE
            WRITE(lu,2000) i,op,eps
            CALL plante(1)
            stop
          ENDIF
!
        ENDDO ! I
!
        ELSEIF(iopt.EQ.4) THEN
!
        DO i=1,npoin
!
          IF (abs(z(i)).GT.eps) THEN
            x(i) = c*y(i) / z(i)
          ELSEIF (abs(c*y(i)).LT.eps) THEN
            x(i) = d
          ELSEIF (z(i).GE.0.d0) THEN
            x(i) =  c*y(i) / eps
          ELSE
            x(i) = -c*y(i) / eps
          ENDIF
!
        ENDDO ! I
!
        ENDIF
!
!-----------------------------------------------------------------------
!
      ELSEIF(op(1:8).EQ.'X=CXY/Z ') THEN
!
        IF(iopt.EQ.1) THEN
!
        DO i=1,npoin
            x(i) = c*x(i)*y(i) / z(i)
        ENDDO ! I
!
        ELSEIF(iopt.EQ.2) THEN
!
        DO i=1,npoin
!
          IF (abs(z(i)).GT.eps) THEN
            x(i) = c*x(i)*y(i) / z(i)
          ELSE
            x(i) = d
          ENDIF
!
        ENDDO ! I
!
        ELSEIF(iopt.EQ.3) THEN
!
        DO i=1,npoin
!
          IF (abs(z(i)).GT.eps) THEN
            x(i) = c*x(i)*y(i) / z(i)
          ELSE
            WRITE(lu,2000) i,op,eps
            CALL plante(1)
            stop
          ENDIF
!
        ENDDO ! I
!
        ELSEIF(iopt.EQ.4) THEN
!
        DO i=1,npoin
!
          IF (abs(z(i)).GT.eps) THEN
            x(i) = c*x(i)*y(i) / z(i)
          ELSEIF (abs(c*x(i)*y(i)).LT.eps) THEN
            x(i) = d
          ELSEIF (z(i).GE.0.d0) THEN
            x(i) =  c*y(i) / eps
          ELSE
            x(i) = -c*y(i) / eps
          ENDIF
!
        ENDDO ! I
!
!
        ENDIF
!
!-----------------------------------------------------------------------
!
      ELSEIF(op(1:8).EQ.'X=X+CY/Z') THEN
!
        IF(iopt.EQ.1) THEN
!
        DO i=1,npoin
            x(i) = x(i) + c * y(i) / z(i)
        ENDDO ! I
!
        ELSEIF(iopt.EQ.2) THEN
!
        DO i=1,npoin
!
          IF (abs(z(i)).GT.eps) THEN
            x(i) = x(i) + c * y(i) / z(i)
          ELSE
            x(i) = d
          ENDIF
!
        ENDDO ! I
!
        ELSEIF(iopt.EQ.3) THEN
!
        DO i=1,npoin
!
          IF (abs(z(i)).GT.eps) THEN
            x(i) = x(i) + c * y(i) / z(i)
          ELSE
            WRITE(lu,2000) i,op,eps
            CALL plante(1)
            stop
          ENDIF
!
        ENDDO ! I
!
        ELSEIF(iopt.EQ.4) THEN
!
        DO i=1,npoin
!
          IF (abs(z(i)).GT.eps) THEN
            x(i) = x(i) + c*y(i) / z(i)
          ELSEIF (abs(c*y(i)).LT.eps) THEN
            x(i) = d
          ELSEIF (z(i).GE.0.d0) THEN
            x(i) = x(i) + c*y(i) / eps
          ELSE
            x(i) = x(i) - c*y(i) / eps
          ENDIF
!
        ENDDO ! I
!
        ENDIF
!
!-----------------------------------------------------------------------
!
      ELSE
!
        WRITE(lu,4000) op
        CALL plante(1)
        stop
!
      ENDIF
!
!-----------------------------------------------------------------------
!
2000  FORMAT(1x,'OVD (BIEF) : DIVIDE BY ZERO AT POINT ',1i6,
     &          ' FOR OPERATION ',a8,/,1x,
     &          'THE CRITERION IS ',g16.7)
4000  FORMAT(1x,'OVD (BIEF) : UNKNOWN OPERATION: ',a8)
!
!-----------------------------------------------------------------------
!
      RETURN
      END