text:
c=================================================================
c     Daniel R. Reynolds
c     SMU Mathematics
c     Math 4370/6370
c     7 February 2015
c=================================================================

      program DotProd
c-----------------------------------------------------------------
c     Example routine to compute the dot-product of two vectors.
c-----------------------------------------------------------------
c======= Inclusions ===========

c======= Declarations =========
      implicit none

      integer, parameter :: MAX_INTEGER_DIGITS = 10
      integer :: i, n
      character(MAX_INTEGER_DIGITS) :: n_string
      double precision, allocatable :: a(:), b(:)
      double precision :: sum, stime, ftime
      double precision :: alloctime, inittime, runtime

c======= Internals ============

c     get n from the command line
      call getarg(1, n_string)

c     ensure that an argument was passed in
      if ( n_string == '' ) then
         stop 'Error: function requires one argument (vector length)'
      endif

c     convert n_string to integer, and ensure it's positive
      read (n_string, *) n
      if ( n < 1 ) then
         stop 'Error: vector length must be greater than 0'
      endif

c     allocate the vectors
      call get_time(stime)
      allocate(a(n),b(n))
      call get_time(ftime)
      alloctime = ftime - stime

c     initialize the vectors
      call get_time(stime)
      do i=1,n
         a(i) = 1.d-3*(i)/n
         b(i) = 1.d-3*(n-i)/n
      enddo
      call get_time(ftime)
      inittime = ftime - stime

c     compute dot-product
      call get_time(stime)
      sum = 0.d0
      do i=1,n,1
         sum = sum + a(i)*b(i)
      enddo
      call get_time(ftime)
      runtime = ftime-stime

c     output computed value and runtime
      print *, ' vector length =',n
      print *, '   dot-product =',sum
      print *, '    alloc time =',alloctime
      print *, '     init time =',inittime
      print *, '      run time =',runtime

c     free vectors
      deallocate(a,b)

      end program DotProd
c=================================================================

ext: f


text:
      subroutine dkqqb_wz_g(p,msq)
      implicit none
C----Author R.K.Ellis June 2012
c----Matrix element for WZ production with radiation in the
C----hadronic decay of the Z
C----averaged over initial colours and spins
C----massless final state particles
c     q(-p1)+qbar(-p2)-->n(p3)+ebar(p4)+q'(p5)+bar{q'}(p6)+g(p7)
      include 'constants.f'
      include 'masses.f'
      include 'sprods_com.f'
      include 'zprods_com.f'
      include 'zerowidth.f'
      include 'ewcouple.f'
      include 'qcdcouple.f'
      include 'zcouple.f'
      include 'ewcharge.f'
      include 'ckm.f'
c      include 'anomcoup.f'
      include 'plabel.f'
      include 'srdiags.f'
      include 'nwz.f'
      include 'flqqb.f'
      double precision msq(-nf:nf,-nf:nf),p(mxpart,4),
     & fac,s567,cotw,qqbsq,qbqsq
      double complex prop12,prop34,prop567,AWZqbq(2,2),AWZqqb(2,2)
      double complex cprop,props,propw,propz
      double complex Fa12(2,2),Fb12(2,2),Fc12(2,2),Fd12(2,2)
      double complex Fa21(2,2),Fb21(2,2),Fc21(2,2),Fd21(2,2)
      double complex qqb(4,2,2),qbq(4,2,2)
      double complex cs_z(2,2),cs_g(2,2)
      integer j,k,hg,hq

      fac=aveqq*xn*gw**8*gsq*2d0*CF
      flqqb=2
      srdiags=.true.
      do j=-nf,nf
      do k=-nf,nf
c--set msq=0 to initalize
      msq(j,k)=0d0
      enddo
      enddo

c   We have --- f(p1) + f'(p2)--> nu(p3)+e^+(p4)+q(p5)+qbar(p6)+g(p7)
c   We have --- f(p1) + f'(p2)--> q(p3)+qbar(p4)+e^-(p5)+nu~(p6)+g(p7)

      call spinoru(7,p,za,zb)

      if     (zerowidth  .eqv. .true.) then
      write(6,*) 'dkqqb_ww_g.f:zerowidth should be false'
      stop
      elseif (zerowidth .neqv. .true.) then
c--   calculate propagators
      cotw=dsqrt((one-xw)/xw)
      s567=s(5,6)+s(5,7)+s(6,7)
      if     (zerowidth  .eqv. .true.) then
      prop12=s(1,2)/dcmplx(s(1,2)-wmass**2,wmass*wwidth)
      prop34=s(3,4)/dcmplx(s(3,4)-wmass**2,wmass*wwidth)
      prop567=s567/dcmplx(s567-zmass**2,zmass*zwidth)
      cprop=dcmplx(1d0)
      elseif (zerowidth .neqv. .true.) then
      prop12=dcmplx(s(1,2)/(s(1,2)-wmass**2))
      prop34=dcmplx(s(3,4)/(s(3,4)-wmass**2))
      prop567=dcmplx(s567/(s567-zmass**2))
      props=(s(1,2)-wmass**2)/dcmplx(s(1,2)-wmass**2,wmass*wwidth)
      propw=(s(3,4)-wmass**2)/dcmplx(s(3,4)-wmass**2,wmass*wwidth)
      propz=(s(5,6)-zmass**2)/dcmplx(s(5,6)-zmass**2,zmass*zwidth)
      cprop=props*propw*propz
      endif

c--- DEBUG to compare with Madgraph
      prop12=s(1,2)/dcmplx(s(1,2)-wmass**2,wmass*wwidth)
      prop34=s(3,4)/dcmplx(s(3,4)-wmass**2,wmass*wwidth)
      prop567=s567/dcmplx(s567-zmass**2,zmass*zwidth)
      cprop=dcmplx(1d0)
c--- DEBUG to compare with Madgraph

      endif

c-- couplings with or without photon pole
c---first index is quark helicity, second is weak isospin
      do j=1,2
      cs_z(1,j)=2d0*xw*l(j)*prop567
      cs_z(2,j)=2d0*xw*r(j)*prop567
      cs_g(1,j)=2d0*Q(j)*xw
      cs_g(2,j)=2d0*Q(j)*xw
      enddo

      if ((plabel(3) .eq. 'nl').and.(plabel(5) .eq. 'qj')) then
      call dkdrwz(1,2,3,4,5,6,7,Fa12,Fb12,Fc12,Fd12)
      call dkdrwz(2,1,3,4,5,6,7,Fa21,Fb21,Fc21,Fd21)
      elseif ((plabel(3) .eq. 'el').and.(plabel(5) .eq. 'qj')) then
      write(6,*) 'nwz',nwz
      call dkdrwz(1,2,3,4,5,6,7,Fb12,Fa12,Fd12,Fc12)
      call dkdrwz(2,1,3,4,5,6,7,Fb21,Fa21,Fd21,Fc21)
      endif

C----hq is the helicity of the 56 line
C----hg is the helicity of the gluon

      do hg=1,2
      do hq=1,2
      AWZqqb(hq,hg)=
     &  +prop34*(cs_z(hq,flqqb)*L(2)+cs_g(hq,flqqb)*Q(2))*Fa12(hq,hg) ! Coupling of strange line to up
     &  +prop34*(cs_z(hq,flqqb)*L(1)+cs_g(hq,flqqb)*Q(1))*Fb12(hq,hg) ! Coupling of strange line to down
     &  +prop12*prop34*(cotw*cs_z(hq,flqqb)+cs_g(hq,flqqb))*Fc12(hq,hg)

      AWZqbq(hq,hg)=
     &  +prop34*(cs_z(hq,flqqb)*L(2)+cs_g(hq,flqqb)*Q(2))*Fa21(hq,hg)
     &  +prop34*(cs_z(hq,flqqb)*L(1)+cs_g(hq,flqqb)*Q(1))*Fb21(hq,hg)
     &  +prop12*prop34*(cotw*cs_z(hq,flqqb)+cs_g(hq,flqqb))*Fc21(hq,hg)
      enddo
      enddo

      if (srdiags) then
c---we need supplementary diagrams for gauge invariance.

      if (nwz .eq. +1) then
      call dksrwz(1,2,3,4,5,6,7,qqb)   ! qqb
      call dksrwz(2,1,3,4,5,6,7,qbq)   ! qbq

      do hg=1,2
      do hq=1,2

      AWZqqb(hq,hg)=AWZqqb(hq,hg)
     & +prop12*qqb(1,hq,hg)*(cs_z(hq,flqqb)*ln+0d0*cs_g(hq,flqqb))
     & +prop12*qqb(2,hq,hg)*(cs_z(hq,flqqb)*le+qe*cs_g(hq,flqqb))
     & +prop12*prop34*qqb(2+flqqb,hq,hg)
      AWZqbq(hq,hg)=AWZqbq(hq,hg)
     & +prop12*qbq(1,hq,hg)*(cs_z(hq,flqqb)*ln+0d0*cs_g(hq,flqqb))
     & +prop12*qbq(2,hq,hg)*(cs_z(hq,flqqb)*le+qe*cs_g(hq,flqqb))
     & +prop12*prop34*qbq(2+flqqb,hq,hg)
      enddo
      enddo

      elseif (nwz .eq. -1) then
      call dksrwz(1,2,3,4,5,6,7,qqb)   ! qqb
      call dksrwz(2,1,3,4,5,6,7,qbq)   ! qbq
      do hg=1,2
      do hq=1,2
      AWZqqb(hq,hg)=AWZqqb(hq,hg)
     & +prop12*qqb(2,hq,hg)*(cs_z(hq,flqqb)*ln+0d0*cs_g(hq,flqqb))
     & +prop12*qqb(1,hq,hg)*(cs_z(hq,flqqb)*le+qe*cs_g(hq,flqqb))
     & +prop12*prop34*qqb(2+flqqb,hq,hg)
      AWZqbq(hq,hg)=AWZqbq(hq,hg)
     & +prop12*qbq(2,hq,hg)*(cs_z(hq,flqqb)*ln+0d0*cs_g(hq,flqqb))
     & +prop12*qbq(1,hq,hg)*(cs_z(hq,flqqb)*le+qe*cs_g(hq,flqqb))
     & +prop12*prop34*qbq(2+flqqb,hq,hg)
      enddo
      enddo
      endif
      endif

C-- Inclusion of width for W's a la Baur and Zeppenfeld with cprop.
      do hg=1,2
      do hq=1,2
      AWZqqb(hq,hg)=cprop*AWZqqb(hq,hg)
      AWZqbq(hq,hg)=cprop*AWZqbq(hq,hg)
      enddo
      enddo

      qqbsq=fac*
     & (+abs(AWZqqb(1,1))**2+abs(AWZqqb(1,2))**2
     &  +abs(AWZqqb(2,1))**2+abs(AWZqqb(2,2))**2)

      qbqsq=fac*
     & (+abs(AWZqbq(1,1))**2+abs(AWZqbq(1,2))**2
     &  +abs(AWZqbq(2,1))**2+abs(AWZqbq(2,2))**2)

      do j=-nf,nf
      do k=-nf,nf
c--set msq=0 to initalize
      msq(j,k)=0d0
          if ((j .gt. 0) .and. (k .lt. 0)) then
            msq(j,k)=Vsq(j,k)*qqbsq
          elseif ((j .lt. 0) .and. (k .gt. 0)) then
            msq(j,k)=Vsq(j,k)*qbqsq
          endif
      enddo
      enddo
      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE BD_USET ( CARD )

! Processes USET Bulk Data Cards. Reads and checks data and then write a record to file LINK1X for later processing.
! Each record in file LINK1X has:

!          USET_NAME, COMPJ, GRIDJ, DOFSET

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06, L1X
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, ECHO, FATAL_ERR, IERRFL, JCARD_LEN, JF, NUM_USET_RECORDS, WARN_ERR
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  BD_USET_BEGEND
      USE CONSTANTS_1, ONLY           :  ZERO
      USE PARAMS, ONLY                :  SUPWARN
      USE DOF_TABLES, ONLY            :  TSET_CHR_LEN

      USE BD_USET_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'BD_USET'
      CHARACTER(LEN=*),INTENT(IN)     :: CARD              ! A Bulk Data card
      CHARACTER(LEN=JCARD_LEN)        :: CHRFLD            ! A character field from the entry
      CHARACTER( 8*BYTE)              :: IP6TYP            ! An output from subr IP6CHK called herein
      CHARACTER(LEN=JCARD_LEN)        :: JCARD(10)         ! The 10 fields of characters making up CARD
      CHARACTER(LEN(JCARD))           :: JCARDO            ! An output from subr IP6CHK called herein
      CHARACTER( 2*BYTE)              :: USET_NAME         ! Name in field 2 of the USET entry

      INTEGER(LONG)                   :: COMPJ     = 0     ! Displ components constrained at GRIDJ
      INTEGER(LONG)                   :: GRIDJ     = 0     ! Grid ID on USET card
      INTEGER(LONG)                   :: I                 ! DO loop index
      INTEGER(LONG)                   :: IDUM              ! Dummy arg in subr IP6CHK not used herein
      INTEGER(LONG)                   :: JERR      = 0     ! A local error count
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = BD_USET_BEGEND

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
!  USET Bulk Data Card routine

!    FIELD   ITEM           ARRAY ELEMENT
!    -----   ------------   -------------
!     2      USET name (i.e. "U1", "U2")
!     3      Grid ID
!     4      Comp. numbers
!     5      Grid ID
!     6      Comp.numbers
!     7      Grid ID
!     8      Comp.numbers

! Initialize

      JERR = 0

!  Make JCARD from CARD

      CALL MKJCARD ( SUBR_NAME, CARD, JCARD )

! Read USET name

      CALL CHAR_FLD ( JCARD(2), JF(2), CHRFLD )            ! Read field 2: USET name
      IF (IERRFL(2) == 'N') THEN
         CALL LEFT_ADJ_BDFLD ( CHRFLD )
         USET_NAME = CHRFLD(1:2)
         IF ((USET_NAME /= 'U1      ') .AND. (USET_NAME /= 'U2      ')) THEN
            JERR = JERR + 1
            FATAL_ERR = FATAL_ERR + 1
            WRITE(ERR,1199) SUBR_NAME, USET_NAME
            WRITE(F06,1199) SUBR_NAME, USET_NAME
         ENDIF
      ENDIF

! Process data in fields 3-8

      DO I=3,7,2                                           ! There can be 3 sets of (Grid, Comp) on each card

         IF (JCARD(I)(1:) /= ' ') THEN

            CALL I4FLD ( JCARD(I), JF(I), GRIDJ )          ! Read Grid ID
                                                           ! Read displ components
            CALL IP6CHK ( JCARD(I+1), JCARDO, IP6TYP, IDUM )
            IF ((IP6TYP == 'COMP NOS') .OR. (IP6TYP == 'ZERO    ') .OR. (IP6TYP == 'BLANK   ')) THEN
               CALL I4FLD ( JCARDO, JF(I+1), COMPJ )
            ELSE
               JERR      = JERR + 1
               FATAL_ERR = FATAL_ERR + 1
               WRITE(ERR,1124) JF(I+1), JCARD(1), JCARD(2), JF(I+1), JCARD(I+1)
               WRITE(F06,1124) JF(I+1), JCARD(1), JCARD(2), JF(I+1), JCARD(I+1)
            ENDIF

            IF ((JERR == 0 ) .AND. (IERRFL(I) == 'N') .AND. (IERRFL(I+1) == 'N')) THEN
               NUM_USET_RECORDS = NUM_USET_RECORDS + 1     ! Incr count of number of entries written to file LINK1X

               WRITE(L1X) USET_NAME(1:2), COMPJ, GRIDJ, GRIDJ
            ENDIF

         ELSE                                              ! Field 3, 5 or 7 is blank

            IF (JCARD(I+1)(1:) /= ' ') THEN
               WARN_ERR = WARN_ERR + 1
               WRITE(ERR,101) CARD
               WRITE(ERR,1144) JCARD(1), JCARD(2), JF(I+1), JF(I)
               IF (SUPWARN == 'N') THEN
                  IF (ECHO == 'NONE  ') THEN
                     WRITE(F06,101) CARD
                  ENDIF
                  WRITE(F06,1144) JCARD(1), JCARD(2), JF(I+1), JF(I)
               ENDIF
            ENDIF

         ENDIF

      ENDDO

      CALL BD_IMBEDDED_BLANK   ( JCARD,0,3,0,5,0,7,0,0 )   ! Make sure that there are no imbedded blanks in fields 2,3,5,6,8,9
      CALL CARD_FLDS_NOT_BLANK ( JCARD,0,0,0,0,0,0,8,9 )
      CALL CRDERR ( CARD )                                 ! CRDERR prints errors found when reading fields

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
  101 FORMAT(A)

 1124 FORMAT(' *ERROR  1124: INVALID DOF NUMBER IN FIELD ',I3,' ON ',A,' ENTRY WITH ID = ',A                                       &
                    ,/,14X,' MUST BE A COMBINATION OF DIGITS 1-6. HOWEVER, FIELD ',I3, ' HAS: "',A,'"')

 1144 FORMAT(' *WARNING    : ON ',A,' SNAME = ',A,' FIELD ',I3,' IS IGNORED SINCE GRID FIELD ',I3,' IS BLANK.')

 1199 FORMAT(' *ERROR  1199: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' SNAME ON THE ABOVE USET ENTRY MUST BE "U1" OR "U2" BUT IS "',A,'"')

! **********************************************************************************************************************************

      END SUBROUTINE BD_USET


text:
!*************************************************************************
!* Dual Number Automatic Differentiation (DNAD) of Fortran Codes
!*----------------------------------------------------------------
!* COPYRIGHT (c) Wenbin Yu, All rights reserved, you are free to copy,
!* modify or translate this code to other languages such as c/c++. If
!* you find a bug please let me know through wenbinyu.heaven@gmail.com. If
!* you added new functions and want to share with others, please let me know too.
!* You are welcome to share your successful stories with us through
!* http://groups.google.com/group/hifi-comp.
!******************************************************************************
!* Simple Instruction:
!*---------------------
!* This is the procedure file: define the procedures (functions/subroutines)
!* needed for overloading intrinsic functions and operators. If a function
!* or operation is not defined (unlikely), you need to create a corresponding
!* interface in header file (DNADHeaders.f90) and a corresponding
!* function/subroutine in DNAD.f90.
!* If you are uncertain whether I did things in the correct way or not, you can
!* create a new module with the same definition of the type AD_D, but delete
!* all the interfaces and procedures first. Then compile it with your analysis
!* codes after you changed all the definitions of real numbers to be dual numbers.
!* The compiler will complain about some operations and functions are not defined.
!* Then you do it one by one by looking at my original module, copying the needed
!* interfaces/procedures into your module. This way, you can also get a leaner
!* model and use those only needed for your analysis code.
!*********************************************************************************
!* Acknowledgements
!*-------------------
!* The development of DNAD is supported, in part, by the Chief Scientist
!* Innovative Research Fund at AFRL/RB WPAFB, and by Department of Army
!* SBIR (Topic A08-022) through Advanced Dynamics Inc. The views and
!* conclusions contained herein are those of the authors and should not be
!* interpreted as necessarily representing the official policies or endorsement,
!* either expressed or implied, of the funding agency.
!**********************************************************************************
!* Citation
!*-------------------
!* Your citation of the following two papers are appreciated:
!* Yu, W. and Blair, M.: "DNAD, a Simple Tool for Automatic Differentiation of
!* Fortran Codes Using Dual Numbers," Computer Physics Communications, vol. 184,
!* 2013, pp. 1446-1452.
!*
!* Spall, R. and Yu, W.: "Imbedded Dual-Number Automatic Differentiation for CFD
!* Sensitivity Analysis," Journal of Fluids Engineering, vol. 135, 2013, 014501.
!**********************************************************************************
!* Basic idea
!*-------------------------------------
!* Carry out Automatic Differentiation using the arithmetic of dual numbers.
!* It is supposed to be more efficient and more accurate than complex-step
!* approximation. What one needs to do is to create a new data type according
!* to the rules of dual numbers, which is equivalent to create a quasi-complex
!* numbers. And overload all the intrinsic operators and functions for
!* such a type.
!* For simplicity and memoric convenience, we call a dual number as a dual (D),
!* a double precision real number as a real (R) and a single precision real
!* number as a single (S), and an integer number as an integer (I). We also
!* use M to denote a 2D array, V to denote a vector, M3 to denote a 3D arry.
!* By defining functions/subroutines to the ELEMENTAL, the inputs/outputs
!* can be automatically overloaded for vectors and matrices of the same shape
!**********************************************************************************
!* To AD a Fortran code use DNAD
!* Step 0: compile DNAD.f90 to be DNAD.o
!* Step 1: Replace all the definitions of real numbers with dual numbers
!*        For example
!*              replace REAL(8) :: x
!*               with  TYPE(AD_D) :: x
!*              replace REAL(8), PARAMETER:: ONE=1.0D0
!*               with TYPE(AD_D),PARAMETER::ONE=AD_D(1.0D0,0.D0)
!* Step 2: Insert USE Dual_Num_Auto_Diff right after Module/Function/Subroutine/Program
!*         statements used TYPE(AD_D)
!* Step 3: Change IO commands correspondingly if the code does not use free formatting
!*         read and write (can be automated by written some general-purpose utility subroutines)
!* Step 4: Recompile the source along with DNAD.o
!* The whole process can be automated, and even manually it only takes just a few minutes
!* for most real analysis codes, although step 3 is code dependent.
!*****************************************************************************************
!* Example
!* The analysis version of the code computing the area of a circle with input radius
!* PROGRAM CircleArea
!*    REAL(8),PARAMETER:: PI=3.141592653589793D0
!*    REAL(8):: radius, area
!*    READ(*,*) radius
!*    Area=PI*radius**2
!*    WRITE(*,*) "AREA=", Area
!* END PROGRAM CircleArea
!*Input: 5
!*Output: AREA=78.5398163397448
!*---------------------------------------------------------------------------------------
!* The AD version of the code computing the area of a circle and sensitivity with input radius
!*PROGRAM CircleArea
!*    USE DNAD
!*    TYPE (AD_D),PARAMETER:: PI=AD_D(3.141592653589793D0,0.D0)
!*    TYPE (AD_D):: radius,area
!*    READ(*,*) radius
!*    Area=PI*radius**2
!*    WRITE(*,*) "AREA=",Area
!*END PROGRAM CircleArea
!* Input 5, 1
!* Output: AREA=78.5398163397448, 31.4159265358979
!*--------------------------------------------------------
!*****************************************************************************************
!* Possible Mistakes to Avoid
!* 11/23/2011: DNAD should always compute the same functional value. However, in very rare situations,
!* sensitivities are not computed as expected. One such case has been identified with
!* dgemv in the blas package. The main reason is that dgemv skip some calculations due to
!* the fact of multiplying zero. However, sometimes, the functional value is zero, but
!* the sensitivity is not zero. In DNAD the comparision is done only between functional values
!* hence is sensitivity calculation is also avoid altogether. For example, the following piece of code
!*                IF( X( JX ).NE.ZERO )THEN
!*                  TEMP = ALPHA*X( JX )
!*                  DO 50, I = 1, M
!*                     Y( I ) = Y( I ) + TEMP*A( I, J )
!*   50             CONTINUE
!*              END IF
!* the solution is to either to rewrite dgemv or assign epsilon(0.0d0) to x(jx) if the functional value
!* while the corresponding derivatives are not zero.
!************************************************************************************************************

MODULE DNAD_D
    use mod_kinds,      only: rsingle,rdouble,ishort,ilong,rk
    use mod_constants,  only: ZERO, ONE

IMPLICIT NONE

    external DGEMM
    external DAXPY

! Dynamic allocation, so this is no longer used
!INTEGER(2), PUBLIC,PARAMETER:: NDV_AD=4 ! number of design variables

PRIVATE
! Imported kinds from common file for consistency
integer, parameter :: SNG_AD = rsingle
integer, parameter :: DBL_AD = rdouble
!INTEGER, PARAMETER:: SNG_AD=SELECTED_REAL_KIND(6)  ! single precision
!INTEGER, PARAMETER:: DBL_AD=SELECTED_REAL_KIND(15) ! double precision
REAL(DBL_AD)      ::negative_one=-1.0d0

    TYPE,PUBLIC:: AD_D
    !TYPE :: AD_D
                                    ! make this private will create difficulty to use the original write/read commands,
                                   ! hence x_ad_ and xp_ad_ are variables which can be accessed using D%x_ad_ and
                                    ! D%xp_ad_ in other units using this module in which D is defined as TYPE DUAL_NUM.
        !REAL(DBL_AD)              :: x_ad_      ! functional value
        !REAL(DBL_AD), allocatable :: xp_ad_(:)  ! derivative
        REAL(rk)              :: x_ad_      ! functional value
        REAL(rk), allocatable :: xp_ad_(:)  ! derivative

!    contains
!
!        final :: destructor

    END TYPE AD_D

    interface AD_D
        module procedure new_dual_num_s
        module procedure new_dual_num_l
    end interface

INCLUDE 'DNADHeaders_D.f90'

CONTAINS
!=======================
!
!   Contstructors
!
!=======================
    ! SHORT INTEGER
    type(AD_D) function new_dual_num_s(nd)
        integer(ishort), intent(in) :: nd
        integer(ishort)    :: ierr

        allocate(new_dual_num_s%xp_ad_(nd),stat=ierr)
        if (ierr /= 0) stop "Allocation error: DNAD_D - new_dual_num"

        new_dual_num_s%xp_ad_ = real(0.,DBL_AD)
    end function

    ! LONG INTEGER
    type(AD_D) function new_dual_num_l(nd)
        integer(ilong), intent(in) :: nd
        integer(ishort)    :: ierr

        allocate(new_dual_num_l%xp_ad_(nd),stat=ierr)
        if (ierr /= 0) stop "Allocation error: DNAD_D - new_dual_num"

        new_dual_num_l%xp_ad_ = real(0.,DBL_AD)
    end function

!*********Begin: functions/subroutines for overloading operators

!******* Begin: (=)
!---------------------
    !** dual=integer:     ! <u,up>=<a,0.D0>
    ELEMENTAL SUBROUTINE ASSIGN_DI_D(u,n)
         TYPE (AD_D), INTENT(INOUT) :: u
         INTEGER, INTENT(IN)::n

 !         u=AD_D(n,0.d0)    ! It is shown that this approach is much slower than the new approach
         u%x_ad_= REAL(n,DBL_AD) ! REAL(n,DBL_AD) is faster than let the code do the conversion.
         u%xp_ad_=0.0D0

    END SUBROUTINE ASSIGN_DI_D

    !** dual=real: <u,up>=<a,0.D0>
    ELEMENTAL SUBROUTINE ASSIGN_DR_D(u,n)
         TYPE (AD_D), INTENT(INOUT)::u
         REAL(DBL_AD), INTENT(IN)::n

          u%x_ad_=n
          u%xp_ad_=0.0D0

    END SUBROUTINE ASSIGN_DR_D

    !** dual=single
    ! <u,up>=<n,0.D0>
    ELEMENTAL SUBROUTINE ASSIGN_DS_D(u,n)
         TYPE (AD_D), INTENT(INOUT)::u
         REAL(SNG_AD), INTENT(IN)::n

         u%x_ad_=REAL(n,DBL_AD)
         u%xp_ad_=0.0D0

    END SUBROUTINE ASSIGN_DS_D

    !** integer=dual
    ! n=u%x_ad_
    ELEMENTAL SUBROUTINE ASSIGN_ID_D(n,u)
         TYPE (AD_D), INTENT(IN)::u
         INTEGER, INTENT(OUT)::n

         ! Wukie
         n=int(u%x_ad_)
         !n=u%x_ad_

    END SUBROUTINE ASSIGN_ID_D

!******* End: (=)
!---------------------

!******* Begin: (+)
!---------------------

    !dual=+<v,vp>
    !-----------------------------------------
    ELEMENTAL FUNCTION ADD_D_D(u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res=u ! It is faster than assigning component wise

    END FUNCTION ADD_D_D

    ! <u,up>+<v,vp>=<u+v,up+vp>
    !-----------------------------------------
    ELEMENTAL FUNCTION ADD_DD_D(u,v) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u,v
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(v%xp_ad_)))

         res%x_ad_  = u%x_ad_+v%x_ad_
         res%xp_ad_ = u%xp_ad_+v%xp_ad_

    END FUNCTION ADD_DD_D

    ! dual+integer
    ! <v,vp>+<n,0>=<n+v,vp>
    !-----------------------------------------
    ELEMENTAL FUNCTION ADD_DI_D(v,n) RESULT(res)
         TYPE (AD_D), INTENT(IN)::v
         INTEGER,INTENT(IN)::n
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(v%xp_ad_)))

         res%x_ad_  = REAL(n,DBL_AD)+v%x_ad_
         res%xp_ad_ = v%xp_ad_

    END FUNCTION ADD_DI_D

    ! dual +real
    ! <v,vp>+<n,0>=<n+v,vp>
    !-----------------------------------------
    ELEMENTAL FUNCTION ADD_DR_D(v,n) RESULT(res)
         TYPE (AD_D), INTENT(IN)::v
         REAL(DBL_AD),INTENT(IN)::n
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(v%xp_ad_)))

         res%x_ad_  = n+v%x_ad_
         res%xp_ad_ = v%xp_ad_

    END FUNCTION ADD_DR_D

    !-----------------------------------------
    ! dual+single
    ! <v,vp>+<n,0>=<n+v,vp>
    !-----------------------------------------
    ELEMENTAL FUNCTION ADD_DS_D(v,n) RESULT(res)
         TYPE (AD_D), INTENT(IN)::v
         REAL(SNG_AD),INTENT(IN)::n
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(v%xp_ad_)))

         res%x_ad_  = REAL(n,DBL_AD)+v%x_ad_
         res%xp_ad_ = v%xp_ad_

    END FUNCTION ADD_DS_D
    !-----------------------------------------

    !-----------------------------------------
    !an integer+ dual number
    ! <n,0>+<v,vp>=<n+v,vp>
    !-----------------------------------------
    ELEMENTAL FUNCTION ADD_ID_D(n,v) RESULT(res)
         INTEGER,INTENT(IN)::n
         TYPE (AD_D), INTENT(IN)::v
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(v%xp_ad_)))

         res%x_ad_  = REAL(n,DBL_AD)+v%x_ad_
         res%xp_ad_ = v%xp_ad_

    END FUNCTION ADD_ID_D

    !-----------------------------------------
    ! real + dual
    ! <n,0>+<v,vp>=<n+v,vp>
    !-----------------------------------------
    ELEMENTAL FUNCTION ADD_RD_D(n,v) RESULT(res)
         REAL(DBL_AD),INTENT(IN)::n
         TYPE (AD_D), INTENT(IN)::v
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(v%xp_ad_)))

         res%x_ad_  = n+v%x_ad_
         res%xp_ad_ = v%xp_ad_

    END FUNCTION ADD_RD_D

     ! single + dual
    ! <n,0>+<v,vp>=<n+v,vp>
    !-----------------------------------------
    ELEMENTAL FUNCTION ADD_SD_D(n,v) RESULT(res)
         REAL(SNG_AD),INTENT(IN)::n
         TYPE (AD_D), INTENT(IN)::v
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(v%xp_ad_)))

         res%x_ad_  = REAL(n,DBL_AD)+v%x_ad_
         res%xp_ad_ = v%xp_ad_

    END FUNCTION ADD_SD_D

!******* End: (+)
!---------------------

!******* Begin: (-)
!---------------------

    !-----------------------------------------
    ! negate a dual number
    !-------------------------------------------------
    ELEMENTAL FUNCTION MINUS_D_D(u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = -u%x_ad_
         res%xp_ad_= -u%xp_ad_

    END FUNCTION MINUS_D_D

    !-----------------------------------------
    ! <u,up>-<v,vp>=<u-v,up-vp>
    !-------------------------------------------------
    ELEMENTAL FUNCTION MINUS_DD_D(u,v) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u,v
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = u%x_ad_-v%x_ad_
         res%xp_ad_= u%xp_ad_-v%xp_ad_

    END FUNCTION MINUS_DD_D

    !-----------------------------------------
    ! dual number - integer
    ! <u,up>-<n,0>=<u-n,up>
    !-------------------------------------------------
    ELEMENTAL FUNCTION MINUS_DI_D(u,n) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         INTEGER,INTENT(IN)::n
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = u%x_ad_-REAL(n,DBL_AD)
         res%xp_ad_= u%xp_ad_

    END FUNCTION MINUS_DI_D

    !-----------------------------------------
    ! dual number - real
    ! <u,up>-<n,0>=<u-n,up>
    !-------------------------------------------------
    ELEMENTAL FUNCTION MINUS_DR_D(u,n) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         REAL(DBL_AD),INTENT(IN)::n
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = u%x_ad_-n
         res%xp_ad_= u%xp_ad_

    END FUNCTION MINUS_DR_D

    !-----------------------------------------
    ! dual number - single
    ! <u,up>-<n,0>=<u-n,up>
    !-------------------------------------------------
    ELEMENTAL FUNCTION MINUS_DS_D(u,n) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         REAL(SNG_AD),INTENT(IN)::n
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = u%x_ad_-REAL(n,DBL_AD)
         res%xp_ad_= u%xp_ad_

    END FUNCTION MINUS_DS_D

    !-----------------------------------------
    ! integer-dual number
    !-------------------------------------------------
    ELEMENTAL FUNCTION MINUS_ID_D(n,u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         INTEGER,INTENT(IN)::n
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = REAL(n,DBL_AD)-u%x_ad_
         res%xp_ad_= -u%xp_ad_

    END FUNCTION MINUS_ID_D

    !-----------------------------------------
    ! real-dual number
    !-------------------------------------------------
    ELEMENTAL FUNCTION MINUS_RD_D(n,u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         REAL(DBL_AD),INTENT(IN)::n
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = n-u%x_ad_
         res%xp_ad_= -u%xp_ad_

    END FUNCTION MINUS_RD_D

    !-----------------------------------------
    ! single-dual number
    ! <n,0>-<u,up>=<n-u,-up>
    !-------------------------------------------------
    ELEMENTAL FUNCTION MINUS_SD_D(n,u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         REAL(SNG_AD),INTENT(IN)::n
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ =REAL(n,DBL_AD) - u%x_ad_
         res%xp_ad_=- u%xp_ad_

    END FUNCTION MINUS_SD_D

!******* END: (-)
!---------------------

!******* BEGIN: (*)
!---------------------

    !-----------------------------------------
    ! <u,up>*<v,vp>=<u*v,up*v+u*vp>
    !----------------------------------------
    ELEMENTAL FUNCTION MULT_DD_D(u,v) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u,v
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = u%x_ad_*v%x_ad_
         res%xp_ad_= u%xp_ad_*v%x_ad_ + u%x_ad_*v%xp_ad_

    END FUNCTION MULT_DD_D

    !-----------------------------------------
    !  dual*integer
    ! <u,up>*<n,0>=<u*n,up*n>
    !----------------------------------------
    ELEMENTAL FUNCTION MULT_DI_D(u,n) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         INTEGER,INTENT(IN)::n
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = REAL(n,DBL_AD)*u%x_ad_
         res%xp_ad_= REAL(n,DBL_AD)*u%xp_ad_

    END FUNCTION MULT_DI_D

    !-----------------------------------------
    !  dual*real
    ! <u,up>*<n,0>=<u*n,up*n>
    !----------------------------------------
    ELEMENTAL FUNCTION MULT_DR_D(u,n) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         REAL(DBL_AD),INTENT(IN)::n
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = n*u%x_ad_
         res%xp_ad_= n*u%xp_ad_

    END FUNCTION MULT_DR_D

    !-----------------------------------------
    !  dual*single
    ! <u,up>*<n,0>=<u*n,up*n>
    !----------------------------------------
    ELEMENTAL FUNCTION MULT_DS_D(u,n) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         REAL(SNG_AD),INTENT(IN)::n
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = REAL(n,DBL_AD)*u%x_ad_
         res%xp_ad_= REAL(n,DBL_AD)*u%xp_ad_

    END FUNCTION MULT_DS_D

    !-----------------------------------------
    ! integer*dual
    ! <n,0>*<v,vp>=<n*v,n*vp>
    !----------------------------------------
    ELEMENTAL FUNCTION MULT_ID_D(n,v) RESULT(res)
         TYPE (AD_D), INTENT(IN)::v
         INTEGER,INTENT(IN)::n
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(v%xp_ad_)))

         res%x_ad_ = REAL(n,DBL_AD)*v%x_ad_
         res%xp_ad_= REAL(n,DBL_AD)*v%xp_ad_

    END FUNCTION MULT_ID_D

    !-----------------------------------------
    ! real* dual
    ! <n,0>*<u,up>=<u*n,up*n>
    !----------------------------------------
    ELEMENTAL FUNCTION MULT_RD_D(n,u) RESULT(res)
!    FUNCTION MULT_RD_D(n,u) RESULT(res)
         REAL(DBL_AD),  INTENT(IN)  ::n
         TYPE(AD_D),    INTENT(IN)  ::u
         TYPE(AD_D) ::res

         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = n*u%x_ad_
         res%xp_ad_= n*u%xp_ad_

    END FUNCTION MULT_RD_D

!
!    FUNCTION MULT_RD_D(n,u) RESULT(res)
!         REAL(DBL_AD),  INTENT(IN)  ::n
!         TYPE(AD_D),    INTENT(IN)  ::u(:)
!         TYPE(AD_D), allocatable                 ::res(:)
!
!        integer :: i
!!         allocate(res(size(u)))
!         res = u
!         res = 0._DBL_AD
!!         allocate(res%xp_ad_(size(u%xp_ad_)))
!
!         res%x_ad_ = n*u%x_ad_
!
!         do i = 1,size(res)
!            res(i)%xp_ad_= n*u(i)%xp_ad_
!        end do
!
!    END FUNCTION MULT_RD_D

    !-----------------------------------------
    ! MULTIPLY a dual number with REAL number
    ! <n,0>*<u,up>=<u*n,up*n>
    !----------------------------------------
    ELEMENTAL FUNCTION MULT_SD_D(n,u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         TYPE (AD_D)::res
         REAL(SNG_AD),INTENT(IN)::n
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = REAL(n,DBL_AD)*u%x_ad_
         res%xp_ad_= REAL(n,DBL_AD)*u%xp_ad_

    END FUNCTION MULT_SD_D

!******* END: (*)
!---------------------

!******* BEGIN: (/)
!---------------------

    !-----------------------------------------
    ! <u,up>/<v,vp>=<u/v,(up-u vp/v)/v>
    !----------------------------------------
    ELEMENTAL FUNCTION DIV_DD_D(u,v) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u,v
         REAL(DBL_AD)::tmp
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         tmp=1.D0/v%x_ad_
         res%x_ad_ = u%x_ad_*tmp
         res%xp_ad_ =(u%xp_ad_- res%x_ad_*v%xp_ad_)*tmp

    END FUNCTION DIV_DD_D

    !-----------------------------------------
    ! <u,up>/<n,0>=<u/n,up/n>
    !----------------------------------------
    ELEMENTAL FUNCTION DIV_DI_D(u,n) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         INTEGER,INTENT(IN)::n
         REAL(DBL_AD)::tmp
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         tmp=1.D0/REAL(n,DBL_AD)
         res%x_ad_ = u%x_ad_*tmp
         res%xp_ad_ =u%xp_ad_*tmp

    END FUNCTION DIV_DI_D

    !-----------------------------------------
    ! DIVIDE dual number with respect to real numbers
    ! <u,up>/<n,0>=<u/n,up/v>
    !----------------------------------------
    ELEMENTAL FUNCTION DIV_DR_D(u,n) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         REAL(DBL_AD),INTENT(IN)::n
         TYPE (AD_D):: res
         REAL(DBL_AD)::tmp
         allocate(res%xp_ad_(size(u%xp_ad_)))

         tmp=1.0D0/n
         res%x_ad_ = u%x_ad_*tmp
         res%xp_ad_ =u%xp_ad_*tmp

    END FUNCTION DIV_DR_D

    !-----------------------------------------
    ! DIVIDE dual number with respect to single
    ! <u,up>/<n,0>=<u/n,up/v>
    !----------------------------------------
    ELEMENTAL FUNCTION DIV_DS_D(u,n) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         REAL(SNG_AD),INTENT(IN)::n
         TYPE (AD_D):: res
         REAL(DBL_AD)::tmp
         allocate(res%xp_ad_(size(u%xp_ad_)))

         tmp=1.0D0/REAL(n,DBL_AD)
         res%x_ad_ = u%x_ad_*tmp
         res%xp_ad_ =u%xp_ad_*tmp

    END FUNCTION DIV_DS_D

    !-----------------------------------------
    ! integer/<v,vp>
    ! <n,0>/<v,vp>=<n/v,-n vp/v/v>=n*<1,-vp/v>/v
    !----------------------------------------
    ELEMENTAL FUNCTION DIV_ID_D(n,v) RESULT(res)
         INTEGER,INTENT(IN)::n
         TYPE (AD_D), INTENT(IN)::v
         REAL(DBL_AD)::tmp
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(v%xp_ad_)))

         tmp=1.D0/v%x_ad_
         res%x_ad_=REAL(n,DBL_AD)*tmp
         res%xp_ad_=-res%x_ad_*tmp*v%xp_ad_

    END FUNCTION DIV_ID_D

    !-----------------------------------------
    ! real/<v,vp>
    ! <n,0>/<v,vp>=<n/v,-n vp/v/v>=n*<1,-vp/v>/v
    !----------------------------------------
    ELEMENTAL FUNCTION DIV_RD_D(n,v) RESULT(res)
         REAL(DBL_AD),INTENT(IN)::n
         TYPE (AD_D), INTENT(IN)::v
         REAL(DBL_AD)::tmp
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(v%xp_ad_)))

         tmp=1.D0/v%x_ad_
         res%x_ad_=n*tmp
         res%xp_ad_=-res%x_ad_*tmp*v%xp_ad_

    END FUNCTION DIV_RD_D

    !-----------------------------------------
    ! single/<v,vp>
    ! <n,0>/<v,vp>=<n/v,-n vp/v/v>=n*<1,-vp/v>/v
    !----------------------------------------
    ELEMENTAL FUNCTION DIV_SD_D(n,v) RESULT(res)
         REAL(SNG_AD),INTENT(IN)::n
         TYPE (AD_D), INTENT(IN)::v
         REAL(DBL_AD)::tmp
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(v%xp_ad_)))

         tmp=1.D0/v%x_ad_
         res%x_ad_=REAL(n,DBL_AD)*tmp
         res%xp_ad_=-res%x_ad_*tmp*v%xp_ad_

    END FUNCTION DIV_SD_D

!******* END: (/)
!---------------------

!******* BEGIN: (**)
!---------------------

    !-----------------------------------------
    ! POWER dual numbers
    ! <u,up>^k=<u^k,k u^{k-1} up>
    !----------------------------------------
    ELEMENTAL FUNCTION POW_I_D(u,k) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         INTEGER,INTENT(IN)::k
         REAL(DBL_AD)::tmp
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         tmp=u%x_ad_**(real(k,rk)-1._rk)
         res%x_ad_ = u%x_ad_*tmp
         res%xp_ad_=REAL(k,DBL_AD)*tmp*u%xp_ad_

    END FUNCTION POW_I_D

    !-----------------------------------------
    ! POWER dual numbers
    ! <u,up>^k=<u^k,k u^{k-1} up>
    !----------------------------------------
    ELEMENTAL FUNCTION POW_R_D(u,k) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         REAL(DBL_AD),INTENT(IN)::k
         REAL(DBL_AD)::tmp
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         tmp=u%x_ad_**(k-1._rk)
         res%x_ad_ = u%x_ad_*tmp
         res%xp_ad_=k*tmp*u%xp_ad_

    END FUNCTION POW_R_D

    !-----------------------------------------
    ! POWER dual numbers
    ! <u,up>^k=<u^k,k u^{k-1} up>
    !----------------------------------------
    ELEMENTAL FUNCTION POW_S_D(u,k) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         REAL(SNG_AD),INTENT(IN)::k
         REAL(DBL_AD)::tmp
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         tmp=u%x_ad_**(real(k,rk)-1._rk)
         res%x_ad_ = u%x_ad_*tmp
         res%xp_ad_=k*tmp*u%xp_ad_

    END FUNCTION POW_S_D

    !-----------------------------------------
    ! POWER dual numbers to a dual power
    ! <u,up>^(v,vp)=<u^v,u^v (v/u*up+Log(u)*vp>
    !----------------------------------------
    ELEMENTAL FUNCTION POW_D_D(u,v) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         TYPE (AD_D), INTENT(IN)::v
         REAL(DBL_AD)::uf,vf
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))
         uf=u%x_ad_
         vf=v%x_ad_

         res%x_ad_ =uf**vf
         res%xp_ad_=res%x_ad_*(vf/uf*u%xp_ad_+LOG(uf)*v%xp_ad_)

    END FUNCTION POW_D_D
!******* END: (**)
!---------------------

!******* BEGIN: (==)
!---------------------
    !-----------------------------------------
    ! compare two dual numbers, simply compare
    ! the functional value.
    !----------------------------------------
    ELEMENTAL FUNCTION EQ_DD_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs,rhs
         LOGICAL::res

         res = (lhs%x_ad_ == rhs%x_ad_)

    END FUNCTION EQ_DD_D

   ! compare a dual with an integer
   !-----------------------------------------
    ELEMENTAL FUNCTION EQ_DI_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         INTEGER,INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ == REAL(rhs,DBL_AD))

    END FUNCTION EQ_DI_D

    !-----------------------------------------
    ! compare a dual number with a real number, simply compare
    ! the functional value.
    !----------------------------------------
    ELEMENTAL FUNCTION EQ_DR_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         REAL(DBL_AD),INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ == rhs)

    END FUNCTION EQ_DR_D

   ! compare a dual with a single
   !-----------------------------------------
    ELEMENTAL FUNCTION EQ_DS_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         REAL(SNG_AD),INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ == REAL(rhs,DBL_AD))

    END FUNCTION EQ_DS_D

    !-----------------------------------------
    ! compare an integer with a dual
    !----------------------------------------
    ELEMENTAL FUNCTION EQ_ID_D(lhs, rhs) RESULT(res)
         INTEGER,INTENT(IN)::lhs
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (REAL(lhs,DBL_AD)==rhs%x_ad_)

    END FUNCTION EQ_ID_D

    !-----------------------------------------
    ! compare a real with a dual
    !----------------------------------------
    ELEMENTAL FUNCTION EQ_RD_D(lhs, rhs) RESULT(res)
         REAL(DBL_AD),INTENT(IN)::lhs
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (lhs==rhs%x_ad_)

    END FUNCTION EQ_RD_D

    !-----------------------------------------
    ! compare a single with a dual
    !----------------------------------------
    ELEMENTAL FUNCTION EQ_SD_D(lhs, rhs) RESULT(res)
         REAL(SNG_AD),INTENT(IN)::lhs
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (REAL(lhs,DBL_AD)==rhs%x_ad_)

    END FUNCTION EQ_SD_D

!******* END: (==)
!---------------------

!******* BEGIN: (<=)
!---------------------
    !-----------------------------------------
    ! compare two dual numbers, simply compare
    ! the functional value.
    !----------------------------------------
    ELEMENTAL FUNCTION LE_DD_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs,rhs
         LOGICAL::res

         res = (lhs%x_ad_ <= rhs%x_ad_)

    END FUNCTION LE_DD_D

   ! compare a dual with an integer
   !-----------------------------------------
    ELEMENTAL FUNCTION LE_DI_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         INTEGER,INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ <= REAL(rhs,DBL_AD))

    END FUNCTION LE_DI_D

    !-----------------------------------------
    ! compare a dual number with a real number, simply compare
    ! the functional value.
    !----------------------------------------
    ELEMENTAL FUNCTION LE_DR_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         REAL(DBL_AD),INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ <= rhs)

    END FUNCTION LE_DR_D

  ! compare a dual with a single
   !----------------------------------------
    ELEMENTAL FUNCTION LE_DS_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         REAL(SNG_AD),INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ <= REAL(rhs,DBL_AD))

    END FUNCTION LE_DS_D

    !-----------------------------------------
    ! compare a dual number with an integer
    !----------------------------------------
    ELEMENTAL FUNCTION LE_ID_D(n, rhs) RESULT(res)
         INTEGER,INTENT(IN)::n
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (REAL(n,DBL_AD)<=rhs%x_ad_)

    END FUNCTION LE_ID_D

!-----------------------------------------
    ! compare a real with a dual
    !----------------------------------------
    ELEMENTAL FUNCTION LE_RD_D(lhs, rhs) RESULT(res)
         REAL(DBL_AD),INTENT(IN)::lhs
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (lhs<=rhs%x_ad_)

    END FUNCTION LE_RD_D

    !-----------------------------------------
    ! compare a single with a dual
    !----------------------------------------
    ELEMENTAL FUNCTION LE_SD_D(lhs, rhs) RESULT(res)
         REAL(SNG_AD),INTENT(IN)::lhs
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (REAL(lhs,DBL_AD)<=rhs%x_ad_)

    END FUNCTION LE_SD_D

!******* END: (<=)
!---------------------

!******* BEGIN: (<)
!---------------------
    !-----------------------------------------
    ! compare two dual numbers, simply compare
    ! the functional value.
    !----------------------------------------
    ELEMENTAL FUNCTION LT_DD_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs,rhs
         LOGICAL::res

         res = (lhs%x_ad_ < rhs%x_ad_)

    END FUNCTION LT_DD_D

   ! compare a dual with an integer
   !-----------------------------------------
    ELEMENTAL FUNCTION LT_DI_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         INTEGER,INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ < REAL(rhs,DBL_AD))

    END FUNCTION LT_DI_D

    !-----------------------------------------
    ! compare a dual number with a real number, simply compare
    ! the functional value.
    !----------------------------------------
    ELEMENTAL FUNCTION LT_DR_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         REAL(DBL_AD),INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ < rhs)

    END FUNCTION LT_DR_D

   ! compare a dual with a single
   !----------------------------------------
    ELEMENTAL FUNCTION LT_DS_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         REAL(SNG_AD),INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ < REAL(rhs,DBL_AD))

    END FUNCTION LT_DS_D

 !-----------------------------------------
    ! compare a dual number with an integer
    !----------------------------------------
    ELEMENTAL FUNCTION LT_ID_D(n, rhs) RESULT(res)
         INTEGER,INTENT(IN)::n
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (REAL(n,DBL_AD)<rhs%x_ad_)

    END FUNCTION LT_ID_D

    !-----------------------------------------
    ! compare a real with a dual
    !----------------------------------------
    ELEMENTAL FUNCTION LT_RD_D(lhs, rhs) RESULT(res)
         REAL(DBL_AD),INTENT(IN)::lhs
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (lhs < rhs%x_ad_)

    END FUNCTION LT_RD_D

    !-----------------------------------------
    ! compare a single with a dual
    !----------------------------------------
    ELEMENTAL FUNCTION LT_SD_D(lhs, rhs) RESULT(res)
         REAL(SNG_AD),INTENT(IN)::lhs
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (REAL(lhs,DBL_AD) < rhs%x_ad_)

    END FUNCTION LT_SD_D

!******* END: (<)
!---------------------

!******* BEGIN: (>=)
!---------------------
    !-----------------------------------------
    ! compare two dual numbers, simply compare
    ! the functional value.
    !----------------------------------------
    ELEMENTAL FUNCTION GE_DD_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs,rhs
         LOGICAL::res

         res = (lhs%x_ad_ >= rhs%x_ad_)

    END FUNCTION GE_DD_D

   ! compare a dual with an integer
   !-----------------------------------------
    ELEMENTAL FUNCTION GE_DI_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         INTEGER,INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ >= REAL(rhs,DBL_AD))

    END FUNCTION GE_DI_D

    !-----------------------------------------
    ! compare a dual number with a real number, simply compare
    ! the functional value.
    !----------------------------------------
    ELEMENTAL FUNCTION GE_DR_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         REAL(DBL_AD),INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ >= rhs)

    END FUNCTION GE_DR_D

   ! compare a dual with a single
   !----------------------------------------
    ELEMENTAL FUNCTION GE_DS_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         REAL(SNG_AD),INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ >=REAL(rhs,DBL_AD))

    END FUNCTION GE_DS_D
 !-----------------------------------------
    ! compare a dual number with an integer
    !----------------------------------------
    ELEMENTAL FUNCTION GE_ID_D(n, rhs) RESULT(res)
         INTEGER,INTENT(IN)::n
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (REAL(n,DBL_AD)>=rhs%x_ad_)

    END FUNCTION GE_ID_D

!-----------------------------------------
    ! compare a real with a dual
    !----------------------------------------
    ELEMENTAL FUNCTION GE_RD_D(lhs, rhs) RESULT(res)
         REAL(DBL_AD),INTENT(IN)::lhs
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (lhs>=rhs%x_ad_)

    END FUNCTION GE_RD_D

    !-----------------------------------------
    ! compare a single with a dual
    !----------------------------------------
    ELEMENTAL FUNCTION GE_SD_D(lhs, rhs) RESULT(res)
         REAL(SNG_AD),INTENT(IN)::lhs
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (REAL(lhs,DBL_AD)>=rhs%x_ad_)

    END FUNCTION GE_SD_D

!******* END: (>=)
!---------------------

!******* BEGIN: (>)
!---------------------
    !-----------------------------------------
    ! compare two dual numbers, simply compare
    ! the functional value.
    !----------------------------------------
    ELEMENTAL FUNCTION GT_DD_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs,rhs
         LOGICAL::res

         res = (lhs%x_ad_ > rhs%x_ad_)

    END FUNCTION GT_DD_D

   ! compare a dual with an integer
   !-----------------------------------------
    ELEMENTAL FUNCTION GT_DI_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         INTEGER,INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ > REAL(rhs,DBL_AD))

    END FUNCTION GT_DI_D

    !-----------------------------------------
    ! compare a dual number with a real number, simply compare
    ! the functional value.
    !----------------------------------------
    ELEMENTAL FUNCTION GT_DR_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         REAL(DBL_AD),INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ > rhs)

    END FUNCTION GT_DR_D

   ! compare a dual with a single
   !----------------------------------------
    ELEMENTAL FUNCTION GT_DS_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         REAL(SNG_AD),INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ >REAL(rhs,DBL_AD))

    END FUNCTION GT_DS_D

  !-----------------------------------------
    ! compare a dual number with an integer
    !----------------------------------------
    ELEMENTAL FUNCTION GT_ID_D(n, rhs) RESULT(res)
         INTEGER,INTENT(IN)::n
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (REAL(n,DBL_AD)>rhs%x_ad_)

    END FUNCTION GT_ID_D

!-----------------------------------------
    ! compare a real with a dual
    !----------------------------------------
    ELEMENTAL FUNCTION GT_RD_D(lhs, rhs) RESULT(res)
         REAL(DBL_AD),INTENT(IN)::lhs
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (lhs>rhs%x_ad_)

    END FUNCTION GT_RD_D

    !-----------------------------------------
    ! compare a single with a dual
    !----------------------------------------
    ELEMENTAL FUNCTION GT_SD_D(lhs, rhs) RESULT(res)
         REAL(SNG_AD),INTENT(IN)::lhs
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (REAL(lhs,DBL_AD)>rhs%x_ad_)

    END FUNCTION GT_SD_D

!******* END: (>)
!---------------------

!******* BEGIN: (/=)
!---------------------
    !-----------------------------------------
    ! compare two dual numbers, simply compare
    ! the functional value.
    !----------------------------------------
    ELEMENTAL FUNCTION NE_DD_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs,rhs
         LOGICAL::res

         res = (lhs%x_ad_ /= rhs%x_ad_)

    END FUNCTION NE_DD_D

   ! compare a dual with an integer
   !-----------------------------------------
    ELEMENTAL FUNCTION NE_DI_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         INTEGER,INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ /= REAL(rhs,DBL_AD))

    END FUNCTION NE_DI_D

    !-----------------------------------------
    ! compare a dual number with a real number, simply compare
    ! the functional value.
    !----------------------------------------
    ELEMENTAL FUNCTION NE_DR_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         REAL(DBL_AD),INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ /= rhs)

    END FUNCTION NE_DR_D

   ! compare a dual with a single
   !----------------------------------------
    ELEMENTAL FUNCTION NE_DS_D(lhs, rhs) RESULT(res)
         TYPE (AD_D), INTENT(IN):: lhs
         REAL(SNG_AD),INTENT(IN)::rhs
         LOGICAL::res

         res = (lhs%x_ad_ /= REAL(rhs,DBL_AD))

    END FUNCTION NE_DS_D

 !-----------------------------------------
    ! compare a dual number with an integer
    !----------------------------------------
    ELEMENTAL FUNCTION NE_ID_D(n, rhs) RESULT(res)
         INTEGER,INTENT(IN)::n
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (REAL(n,DBL_AD) /= rhs%x_ad_)

    END FUNCTION NE_ID_D

!-----------------------------------------
    ! compare a real with a dual
    !----------------------------------------
    ELEMENTAL FUNCTION NE_RD_D(lhs, rhs) RESULT(res)
         REAL(DBL_AD),INTENT(IN)::lhs
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (lhs/=rhs%x_ad_)

    END FUNCTION NE_RD_D

    !-----------------------------------------
    ! compare a single with a dual
    !----------------------------------------
    ELEMENTAL FUNCTION NE_SD_D(lhs, rhs) RESULT(res)
         REAL(SNG_AD),INTENT(IN)::lhs
         TYPE (AD_D), INTENT(IN):: rhs
         LOGICAL::res

         res = (REAL(lhs,DBL_AD)/=rhs%x_ad_)

    END FUNCTION NE_SD_D

!******* END: (/=)
!---------------------

    !---------------------------------------------------
    ! ABS of dual numbers
    ! ABS<u,up>=<ABS(u),up SIGN(u)>
    !----------------------------------------------------
    ELEMENTAL FUNCTION ABS_D_D(u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         IF(u%x_ad_>0) THEN
            res%x_ad_ = u%x_ad_
            res%xp_ad_= u%xp_ad_
         ELSE IF (u%x_ad_<0) THEN
            res%x_ad_ = -u%x_ad_
            res%xp_ad_= -u%xp_ad_
         ELSE
            res%x_ad_ = 0.0D0
!            res%xp_ad_= 0.0 !Set_NaN_D() ! Indicating an undefined derivative, however, for some codes it will cause problem.

            res%xp_ad_ = u%xp_ad_   !> NEED THIS TO ENFORCE ALLOCATION OF res%xp_ad_
            res%xp_ad_ = 0.0D0      !> THEN, SET TO ZERO

         ENDIF

    END FUNCTION ABS_D_D

    !-----------------------------------------
    ! ACOS of dual numbers
    ! ACOS<u,up>=<ACOS(u),-up/sqrt(1-(u%x_ad_)**2)>
    !----------------------------------------
    ELEMENTAL FUNCTION ACOS_D_D(u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         TYPE (AD_D)::res
         REAL(DBL_AD)::tmp
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = ACOS(u%x_ad_)
         IF(u%x_ad_==1.0D0.OR.u%x_ad_==-1.0D0) THEN
            res%xp_ad_=Set_NaN_D() ! Indicating an undefined derivative
         ELSE
            tmp= -1.0d0/SQRT(1.0D0-(u%x_ad_)**2)
            res%xp_ad_= u%xp_ad_*tmp
         ENDIF
    END FUNCTION ACOS_D_D

    !-----------------------------------------
    ! ASIN of dual numbers
    ! ASIN<u,up>=<ASIN(u),up 1/SQRT(1-u^2)>
    !----------------------------------------
    ELEMENTAL FUNCTION ASIN_D_D(u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         TYPE (AD_D)::res
         REAL (DBL_AD):: tmp
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = ASIN(u%x_ad_)
         IF(u%x_ad_==1.0D0.OR.u%x_ad_==-1.0D0) THEN
            res%xp_ad_=Set_NaN_D() ! Indicating an undefined derivative
         ELSE
            tmp= 1.0d0/SQRT(1.0D0-(u%x_ad_)**2)
            res%xp_ad_= u%xp_ad_*tmp
         ENDIF
    END FUNCTION ASIN_D_D

    !-----------------------------------------
    !
    ! @author   Nathan A. Wukie
    ! @date     2/23/2016
    !
    ! ATAN of dual numbers
    ! ATAN<u,up>=<ATAN(u), up/(1+(u%x_ad_)**2)>
    !----------------------------------------
    ELEMENTAL FUNCTION ATAN_D_D(u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         TYPE (AD_D)::res
         REAL(DBL_AD)::tmp
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = ATAN(u%x_ad_)

         tmp= 1.0d0/(1.0D0+(u%x_ad_)**2)
         res%xp_ad_= u%xp_ad_*tmp
    END FUNCTION ATAN_D_D

!    !-----------------------------------------
!    !
!    ! @author   Nathan A. Wukie
!    ! @date     2/23/2016
!    !
!    ! ATAN2 of dual numbers
!    ! ATAN2<u,up>=<ATAN2(u), up/(1+(u%x_ad_)**2)>
!    !----------------------------------------
!    ELEMENTAL FUNCTION ATAN2_D_D(u,v) RESULT(res)
!         TYPE (AD_D), INTENT(IN)::u
!         TYPE (AD_D), INTENT(IN)::v
!         TYPE (AD_D)::res
!         REAL(DBL_AD)::tmp
!         allocate(res%xp_ad_(size(u%xp_ad_)))
!
!         res%x_ad_ = ATAN2(u%x_ad_,v%x_ad_)
!
!         tmp= 1.0d0/(1.0D0+(u%x_ad_)**2)
!         res%xp_ad_= u%xp_ad_*tmp
!    END FUNCTION ATAN2_D_D
!

    !-----------------------------------------
    ! COS of dual numbers
    ! COS<u,up>=<COS(u),-up*SIN(u)>
    !----------------------------------------
    ELEMENTAL FUNCTION COS_D_D(u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         TYPE (AD_D)::res
         REAL(DBL_AD):: tmp
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = COS(u%x_ad_)
         tmp=-SIN(u%x_ad_)
         res%xp_ad_= u%xp_ad_*tmp

    END FUNCTION COS_D_D

    !-----------------------------------------
    ! DOT PRODUCT two dual number vectors
    ! <u,up>.<v,vp>=<u.v,u.vp+up.v>
    !----------------------------------------
!    FUNCTION DOT_PRODUCT_DD_D(u,v) RESULT(res)
!         TYPE (AD_D), INTENT(IN)::u(:),v(:)
!         TYPE (AD_D)::res
!!         INTEGER:: i
!!
!         stop "DOT_PRODUCT not implemented"
!!         res = AD_D(size(u(1)%xp_ad_))
!!         res = 0.
!!
!!         res%x_ad_ = DOT_PRODUCT(u%x_ad_,v%x_ad_)
!!
!!         DO i=1,size(u(1)%xp_ad_)
!!           res%xp_ad_(i) =DOT_PRODUCT(u%x_ad_,v%xp_ad_(i))+DOT_PRODUCT(u%xp_ad_(i),v%x_ad_)
!!         ENDDO
!    END FUNCTION DOT_PRODUCT_DD_D

    !-----------------------------------------
    ! EXPONENTIAL OF dual numbers
    ! EXP<u,up>=<EXP(u),up EXP(u)>
    !----------------------------------------
    ELEMENTAL FUNCTION EXP_D_D(u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         REAL(DBL_AD)::tmp
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         tmp=EXP(u%x_ad_)
         res%x_ad_ = tmp
         res%xp_ad_ =u%xp_ad_*tmp

    END FUNCTION EXP_D_D

    !-----------------------------------------
    ! EXPONENTIAL OF dual numbers
    ! EXP<u,up>=<EXP(u),up EXP(u)>
    !----------------------------------------
    ELEMENTAL FUNCTION INT_D_D(u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         REAL(DBL_AD)::tmp
         INTEGER::res

         tmp=u%x_ad_
         res = INT(tmp)

    END FUNCTION INT_D_D

    !-----------------------------------------
    ! LOG OF dual numbers,defined for u%x>0 only
    ! the error control should be done in the original code
    ! in other words, if u%x<=0, it is not possible to obtain LOG.
    ! LOG<u,up>=<LOG(u),up/u>
    !----------------------------------------
    ELEMENTAL FUNCTION LOG_D_D(u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         TYPE (AD_D)::res
         REAL(DBL_AD)::tmp
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = LOG(u%x_ad_)
         tmp=1.0d0/u%x_ad_
         res%xp_ad_ =u%xp_ad_*tmp

    END FUNCTION LOG_D_D

    !-----------------------------------------
    ! LOG OF dual numbers,defined for u%x>0 only
    ! the error control should be done in the original code
    ! in other words, if u%x<=0, it is not possible to obtain LOG.
    ! LOG<u,up>=<LOG(u),up/u>
    !----------------------------------------
    ELEMENTAL FUNCTION LOG10_D_D(u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         TYPE (AD_D)::res
         REAL(DBL_AD)::tmp
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = LOG10(u%x_ad_)
         tmp=1.0d0/u%x_ad_/LOG(10.0D0)
         res%xp_ad_ =u%xp_ad_*tmp

    END FUNCTION LOG10_D_D

    !-----------------------------------------
    ! MULTIPLY two dual number matrices
    !
    ! <u,up>.<v,vp>=<u.v,up.v+u.vp>
    !----------------------------------------
    function MATMUL_DD_D(u,v) result(res)
         type(AD_D), intent(in) :: u(:,:)
         type(AD_D), intent(in) :: v(:,:)

         type(AD_D)     :: res(size(u,1),size(v,2))
         real(DBL_AD)   :: xp_ad_u(size(u,1),size(u,2)), xp_ad_v(size(v,1),size(v,2)), xp_ad(size(res,1),size(res,2))
         integer        :: i,j,k

         do k = 1,size(v,2)
            do j = 1,size(u,1)
                allocate(res(j,k)%xp_ad_(size(u(1,1)%xp_ad_)))
            end do
        end do

         res%x_ad_ = matmul(u%x_ad_,v%x_ad_)
         DO i=1,size(u(1,1)%xp_ad_)
            ! Assemble derivative components for U
            do k = 1,size(u,2)
                do j = 1,size(u,1)
                    xp_ad_u(j,k) = u(j,k)%xp_ad_(i)
                end do
            end do

            ! Assemble derivative components for V
            do k = 1,size(v,2)
                do j = 1,size(v,1)
                    xp_ad_v(j,k) = v(j,k)%xp_ad_(i)
                end do
            end do

            xp_ad = matmul(xp_ad_u,v%x_ad_) + matmul(u%x_ad_,xp_ad_v)
!            res%xp_ad_(i)= MATMUL(xp_ad_u,v%x_ad_) + MATMUL(u%x_ad_,xp_ad_v)

            do k = 1,size(res,2)
                do j = 1,size(res,1)
                    res(j,k)%xp_ad_(i) = xp_ad(j,k)
                end do
            end do
         ENDDO

!         DO i=1,NDV_AD
!            res%xp_ad_(i)= MATMUL(u%xp_ad_(i),v%x_ad_) + MATMUL(u%x_ad_,v%xp_ad_(i))
!         ENDDO
    end function MATMUL_DD_D
    !**************************************************************************

    !-----------------------------------------
    ! MULTIPLY a dual number matrix with a dual number
    ! vector
    !
    ! <u,up>.<v,vp>=<u.v,up.v+u.vp>
    !----------------------------------------
    function MATMUL_DV_D(u,v) result(res)
         type(AD_D), intent(in)::u(:,:),v(:)

         type(AD_D)     :: res(size(u,1))
         real(DBL_AD)   :: xp_ad_u(size(u,1),size(u,2)), xp_ad_v(size(v)), xp_ad(size(res))
         integer        :: i,j,k

         do j = 1,size(u,1)
                allocate(res(j)%xp_ad_(size(u(1,1)%xp_ad_)))
         end do

         res%x_ad_ = matmul(u%x_ad_,v%x_ad_)
         DO i=1,size(u(1,1)%xp_ad_)
            ! Assemble derivative components for U
            do k = 1,size(u,2)
                do j = 1,size(u,1)
                    xp_ad_u(j,k) = u(j,k)%xp_ad_(i)
                end do
            end do

            ! Assemble derivative components for V
            do j = 1,size(v)
                    xp_ad_v(j) = v(j)%xp_ad_(i)
            end do

           xp_ad = matmul(xp_ad_u,v%x_ad_) + matmul(u%x_ad_,xp_ad_v)

           do j = 1,size(res)
               res(j)%xp_ad_(i) = xp_ad(j)
           end do
         ENDDO

!         DO i=1,NDV_AD
!            res%xp_ad_(i)= MATMUL(u%xp_ad_(i),v%x_ad_) + MATMUL(u%x_ad_,v%xp_ad_(i))
!         ENDDO
    end function MATMUL_DV_D

!    !-----------------------------------------
!    ! MULTIPLY Matrix(Real)-Vector(Dual)
!    !
!    ! <u,up>.<v,vp>=<u.v,up.v+u.vp>
!    !----------------------------------------
!    FUNCTION MATMUL_MV_RD_D(u,v) RESULT(res)
!        !REAL(DBL_AD), intent(in)     :: u(:,:)
!        REAL(rk), intent(in)     :: u(:,:)
!        TYPE(AD_D), intent(in)       :: v(:)
!!        TYPE(AD_D)                      :: res(size(v))
!        TYPE(AD_D)                      :: res(size(u,1))  !> Causes memory leak in ifort 15.0.3. Doesn't seem to be deallocating everything correctly
!!        TYPE(AD_D), allocatable, dimension(:) :: res        !> Declaring as allocatable to fix memory leak
!!        REAL(DBL_AD)                    :: xp_ad_v(size(v)), &
!!                                           xp_ad(size(u,1))
!!        REAL(rk), dimension(size(v),size(v(1)%xp_ad_))      :: xp_ad_vm
!        REAL(rk), dimension(size(v(1)%xp_ad_), size(v))     :: xp_ad_vm_tr
!        REAL(rk), dimension(size(u,1),size(v(1)%xp_ad_))    :: res_xp_m
!!        REAL(rk), dimension(size(v(1)%xp_ad_), size(u,1))   :: res_xp_m_tr
!        real(rk), dimension(size(v))   :: tmp
!        real(rk), dimension(size(u,1)) :: tmp2
!        INTEGER:: i,j
!        integer(4) :: nrows_a, ncols_a, ncols_b
!
!
!
!!        allocate(res(size(u,1)))
!        do j = 1,size(res)
!               allocate(res(j)%xp_ad_(size(v(1)%xp_ad_)))
!        end do
!
!        !
!        ! Standard matrix multiplication of function values
!        !
!
!        ! Broken with GCC 6.1.0
!!        res%x_ad_ = MATMUL(u,v%x_ad_)
!
!        ! Fix for GCC 6.1.0
!        tmp = v%x_ad_
!        tmp2 = matmul(u,tmp)
!        res%x_ad_ = tmp2
!
!!        !
!!        ! Assemble derivative components as a matrix
!!        !
!!        do i = 1,size(v)
!!            xp_ad_vm(i,:) = v(i)%xp_ad_
!!        end do
!
!        !
!        ! Assemble transposed matrix instead to improve memory access
!        !
!        do i = 1,size(v)
!            xp_ad_vm_tr(:,i) = v(i)%xp_ad_
!        end do
!
!
!
!
!        !
!        ! Multiply derivatives
!        !
!        !res_xp_m = matmul(u,xp_ad_vm)
!!        res_xp_m = matmul(u,transpose(xp_ad_vm_tr))
!        nrows_a = size(u,1)
!        ncols_a = size(u,2)
!        ncols_b = size(xp_ad_vm_tr,1) ! since it is getting transposed
!        call DGEMM('N', 'T', nrows_a, ncols_b, ncols_a, ONE, u, nrows_a, xp_ad_vm_tr, ncols_b, ZERO, res_xp_m, nrows_a)
!
!
!
!        !
!        ! Distribute derivatives
!        !
!        do i = 1,size(res)
!            res(i)%xp_ad_ = res_xp_m(i,:)
!        end do
!
!
!    END FUNCTION MATMUL_MV_RD_D

    !-----------------------------------------
    ! MULTIPLY Matrix(Real)-Vector(Dual)
    !
    ! <u,up>.<v,vp>=<u.v,up.v+u.vp>
    !----------------------------------------
    function MATMUL_MV_RD_D(u,v) result(res)
        real(rk),   intent(in)          :: u(:,:)
        type(AD_D), intent(in)          :: v(:)

        type(AD_D)                      :: res(size(u,1))  !> Causes memory leak in ifort 15.0.3. Doesn't seem to be deallocating everything correctly
        real(rk), dimension(size(v))    :: tmp
        real(rk), dimension(size(u,1))  :: tmp2
        integer                         :: i,j,vecsize

        do j = 1,size(res)
               allocate(res(j)%xp_ad_(size(v(1)%xp_ad_)))
               res(j)%xp_ad_ = ZERO
        end do

        !
        ! Standard matrix multiplication of function values
        !

        ! Broken with GCC 6.1.0
        ! res%x_ad_ = MATMUL(u,v%x_ad_)

        ! Fix for GCC 6.1.0
        tmp = v%x_ad_
        tmp2 = matmul(u,tmp)
        res%x_ad_ = tmp2

        !
        ! Multiply derivatives and accumulate
        !
        vecsize = size(v(1)%xp_ad_)
        do j = 1,size(u,2)
            do i = 1,size(u,1)
!                res(i)%xp_ad_ = res(i)%xp_ad_ + u(i,j)*v(j)%xp_ad_
                call DAXPY(vecsize,u(i,j),v(j)%xp_ad_,1,res(i)%xp_ad_,1)    ! Noticable improvement
            end do
        end do

    end function MATMUL_MV_RD_D
    !*****************************************************************************************

    !-----------------------------------------
    ! MULTIPLY a dual vector with a  dual matrix
    !
    ! <u,up>.<v,vp>=<u.v,up.v+u.vp>
    !----------------------------------------
    FUNCTION MATMUL_VD_D(u,v) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u(:),v(:,:)
         TYPE (AD_D)::res(SIZE(v,2))
         REAL(DBL_AD) :: xp_ad_u(size(u)), xp_ad_v(size(v,1),size(v,2)), xp_ad(size(res))
         INTEGER::i,j,k

        ! Disable vector-matrix multiplication
        stop "Vector-Matrix multiplication disabled"

         do j = 1,size(v,2)
                allocate(res(j)%xp_ad_(size(u(1)%xp_ad_)))
         end do

         res%x_ad_= MATMUL(u%x_ad_,v%x_ad_)
         DO i=1,size(u(1)%xp_ad_)
            ! Assemble derivative components for U
            do j = 1,size(u)
                xp_ad_u(j) = u(j)%xp_ad_(i)
            end do

            ! Assemble derivative components for V
            do k = 1,size(v,2)
                do j = 1,size(v,1)
                    xp_ad_v(j,k) = v(j,k)%xp_ad_(i)
                end do
            end do

            xp_ad = MATMUL(xp_ad_u,v%x_ad_) + MATMUL(u%x_ad_,xp_ad_v)

            do j = 1,size(res)
                res(j)%xp_ad_(i) = xp_ad(j)
            end do
!            res%xp_ad_(i)= MATMUL(xp_ad_u,v%x_ad_) + MATMUL(u%x_ad_,xp_ad_v)
         ENDDO

!         DO i=1,size(u(1)%xp_ad_)
!            res%xp_ad_(i)= MATMUL(u%xp_ad_(i),v%x_ad_) + MATMUL(u%x_ad_,v%xp_ad_(i))
!         ENDDO
    END FUNCTION MATMUL_VD_D

    !-----------------------------------------
    ! Obtain the max of 2 to 5 dual numbers
    !----------------------------------------
    ELEMENTAL FUNCTION MAX_DD_D(val1, val2, val3, val4,val5) RESULT(res)
        TYPE (AD_D), INTENT(IN)::val1, val2
        TYPE (AD_D), INTENT(IN),OPTIONAL:: val3, val4,val5
        TYPE (AD_D)::res
        allocate(res%xp_ad_(size(val1%xp_ad_)))

        IF (val1%x_ad_ > val2%x_ad_) THEN
            res = val1
        ELSE
            res = val2
        ENDIF
        IF(PRESENT(val3))THEN
           IF(res%x_ad_<val3%x_ad_) res=val3
        ENDIF
        IF(PRESENT(val4))THEN
           IF(res%x_ad_<val4%x_ad_) res=val4
        ENDIF
        IF(PRESENT(val5))THEN
           IF(res%x_ad_<val5%x_ad_) res=val5
        ENDIF

    END FUNCTION MAX_DD_D

    !-----------------------------------------
    ! Obtain the max of a dual number and an integer
    !----------------------------------------
    ELEMENTAL FUNCTION MAX_DI_D(u, n) RESULT(res)
        TYPE (AD_D), INTENT(IN)::u
        INTEGER,INTENT(IN)::n
        TYPE (AD_D)::res
        allocate(res%xp_ad_(size(u%xp_ad_)))

        IF (u%x_ad_ > n) THEN
            res = u
        ELSE
            res = n
        ENDIF

    END FUNCTION MAX_DI_D

    !-----------------------------------------
    ! Obtain the max of a dual number and a real number
    !----------------------------------------
    ELEMENTAL FUNCTION MAX_DR_D(u, n) RESULT(res)
        TYPE (AD_D), INTENT(IN)::u
        REAL(DBL_AD),INTENT(IN)::n
        TYPE (AD_D)::res
        allocate(res%xp_ad_(size(u%xp_ad_)))

        IF (u%x_ad_ > n) THEN
            res = u
        ELSE
            res = n
        ENDIF

    END FUNCTION MAX_DR_D

 !-----------------------------------------
    ! Obtain the max of a dual number and a single
    !----------------------------------------
    ELEMENTAL FUNCTION MAX_DS_D(u, n) RESULT(res)
        TYPE (AD_D), INTENT(IN)::u
        REAL(SNG_AD),INTENT(IN)::n
        TYPE (AD_D)::res
        allocate(res%xp_ad_(size(u%xp_ad_)))

        IF (u%x_ad_ > n) THEN
            res = u
        ELSE
            res = n
        ENDIF

    END FUNCTION MAX_DS_D

    !---------------------------------------------------
    ! Obtain the max of a real and a dual
    ! note the real argument is renamed as r to avoid
    ! the ambiguity due to positional association using
    ! keywords when functions/subroutines are called
    !---------------------------------------------------
     ELEMENTAL FUNCTION MAX_RD_D(r,u) RESULT(res)
        REAL(DBL_AD),INTENT(IN)::r
        TYPE (AD_D), INTENT(IN)::u
        TYPE (AD_D)::res
        allocate(res%xp_ad_(size(u%xp_ad_)))

        IF (u%x_ad_ > r) THEN
            res = u
        ELSE
            res = r
        ENDIF

    END FUNCTION MAX_RD_D

    !-----------------------------------------
    ! Obtain the max value of vector u
    !----------------------------------------
    FUNCTION MAXVAL_D_D(u) RESULT(res)
        TYPE (AD_D), INTENT(IN)::u(:)
        INTEGER::iloc(1)
        TYPE (AD_D)::res
        allocate(res%xp_ad_(size(u(1)%xp_ad_)))

        stop "MAXVAL - might not work properly"

        iloc=MAXLOC(u%x_ad_)
        res=u(iloc(1))

    END FUNCTION MAXVAL_D_D

    !-----------------------------------------
    ! Obtain the min of 2 to 4 dual numbers
    !----------------------------------------
    ELEMENTAL FUNCTION MIN_DD_D(val1, val2, val3, val4) RESULT(res)
        TYPE (AD_D), INTENT(IN)::val1, val2
        TYPE (AD_D), INTENT(IN),OPTIONAL:: val3, val4
        TYPE (AD_D)::res

        IF (val1%x_ad_ < val2%x_ad_) THEN
            res = val1
        ELSE
            res = val2
        ENDIF
        IF(PRESENT(val3))THEN
           IF(res%x_ad_>val3%x_ad_) res=val3
        ENDIF
        IF(PRESENT(val4))THEN
           IF(res%x_ad_>val4%x_ad_) res=val4
        ENDIF

    END FUNCTION MIN_DD_D

    !-----------------------------------------
    ! Obtain the min of a dual number and a single
    !----------------------------------------
    ELEMENTAL FUNCTION MIN_DR_D(u, n) RESULT(res)
        TYPE (AD_D), INTENT(IN)::u
        REAL(DBL_AD),INTENT(IN)::n
        TYPE (AD_D)::res

        IF (u%x_ad_ < n) THEN
            res = u
        ELSE
            res = n
        ENDIF

    END FUNCTION MIN_DR_D

    !-----------------------------------------
    ! Obtain the min of a dual number and a single
    !----------------------------------------
    ELEMENTAL FUNCTION MIN_DS_D(u, n) RESULT(res)
        TYPE (AD_D), INTENT(IN)::u
        REAL(SNG_AD),INTENT(IN)::n
        TYPE (AD_D)::res

        IF (u%x_ad_ < n) THEN
            res = u
        ELSE
            res = n
        ENDIF

    END FUNCTION MIN_DS_D

    !-----------------------------------------
    ! Obtain the min value of vector u
    !----------------------------------------
    FUNCTION MINVAL_D_D(u) RESULT(res)
        TYPE (AD_D), INTENT(IN)::u(:)
        INTEGER::iloc(1)
        TYPE (AD_D)::res

        stop "MINVAL - might not work properly"

        iloc=MINLOC(u%x_ad_)
        res=u(iloc(1))

    END FUNCTION MINVAL_D_D

    !------------------------------------------------------
    !Returns the nearest integer to u%x, ELEMENTAL
    !------------------------------------------------------
    ELEMENTAL FUNCTION NINT_D_D(u) RESULT(res)
        TYPE (AD_D), INTENT(IN):: u
        INTEGER::res

        res=NINT(u%x_ad_)

    END FUNCTION NINT_D_D

    !----------------------------------------------------------------
    ! SIGN(a,b) with two dual numbers as inputs,
    ! the result will be |a| if b%x>=0, -|a| if b%x<0,ELEMENTAL
    !----------------------------------------------------------------
    ELEMENTAL FUNCTION SIGN_DD_D(val1, val2) RESULT(res)
        TYPE (AD_D), INTENT(IN) :: val1, val2
        TYPE (AD_D)::res

        IF (val2%x_ad_ < 0.D0) THEN
            res = -ABS(val1)
        ELSE
            res =  ABS(val1)
        ENDIF

     END FUNCTION SIGN_DD_D

    !----------------------------------------------------------------
    ! SIGN(a,b) with one real and one dual number as inputs,
    ! the result will be |a| if b%x>=0, -|a| if b%x<0,ELEMENTAL
    !----------------------------------------------------------------
    ELEMENTAL FUNCTION SIGN_RD_D(val1, val2) RESULT(res)
        REAL(DBL_AD),INTENT(IN)::val1
        TYPE (AD_D), INTENT(IN) :: val2
        TYPE (AD_D)::res
        allocate(res%xp_ad_(size(val2%xp_ad_)))

        IF (val2%x_ad_ < 0.D0) THEN
            res = -ABS(val1)
        ELSE
            res = ABS(val1)
        ENDIF

     END FUNCTION SIGN_RD_D

    !-----------------------------------------
    ! SIN of dual numbers
    ! SIN<u,up>=<SIN(u),up COS(u)>
    !----------------------------------------
    ELEMENTAL FUNCTION SIN_D_D(u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         TYPE (AD_D)::res
         REAL (DBL_AD):: tmp
         allocate(res%xp_ad_(size(u%xp_ad_)))

         res%x_ad_ = SIN(u%x_ad_)
         tmp=COS(u%x_ad_)
         res%xp_ad_= u%xp_ad_*tmp

    END FUNCTION SIN_D_D

    !-----------------------------------------
    ! SQRT of dual numbers
    ! SQRT<u,up>=<SQRT(u),up/2/sqrt(u) >
    !----------------------------------------
    !IMPURE ELEMENTAL FUNCTION SQRT_D_D(u) RESULT(res)
    ELEMENTAL FUNCTION SQRT_D_D(u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u
         REAL(DBL_AD):: tmp
         TYPE (AD_D)::res
         allocate(res%xp_ad_(size(u%xp_ad_)))

         tmp=SQRT(u%x_ad_)

         res%x_ad_ = tmp
         IF(tmp==0.0d0) THEN
            res%xp_ad_ = Set_NaN_D()
            !res%xp_ad_ = 0.0d0
         ELSE
            tmp= 0.5D0/tmp
            res%xp_ad_ = u%xp_ad_*tmp
         ENDIF

    END FUNCTION SQRT_D_D

    !-----------------------------------------
    ! SUM OF A DUAL ARRAY
    !----------------------------------------
    FUNCTION SUM_D_D(u) RESULT(res)
         TYPE (AD_D), INTENT(IN)::u(:)
         TYPE (AD_D)::res
         INTEGER:: i,j
         REAL(DBL_AD):: tmp
         allocate(res%xp_ad_(size(u(1)%xp_ad_)))

         res%x_ad_ = SUM(u%x_ad_)
         DO i=1,size(u(1)%xp_ad_)
         tmp = REAL(0.0,DBL_AD)

            DO j=1,size(u)
                tmp = tmp + u(j)%xp_ad_(i)
            ENDDO
            res%xp_ad_(i) = tmp
         ENDDO

    END FUNCTION SUM_D_D

    ELEMENTAL FUNCTION Set_NaN_D() RESULT(res)
         REAL(DBL_AD)::res

         res=SQRT(negative_one)

    END FUNCTION Set_NaN_D

!    subroutine destructor(self)
!        type(AD_D), intent(inout) :: self
!
!        if (allocated(self%xp_ad_)) deallocate(self%xp_ad_)
!
!
!    end subroutine

END MODULE  DNAD_D

max_stars_repo_name: wanglican/ChiDG

max_issues_repo_name: haohb/ChiDG

max_forks_repo_name: haohb/ChiDG


text: http://www.dateline.ucdavis.edu/dl_detail.lasso?id7745 Lon Hatamiya among his other accomplishments appeared in an ad regarding November 2006 Election/Measure H Measure H in the November 2006 Election. He was an Executive in Residence for the Graduate School of Management.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      FUNCTION SQUARES(X)
      X = X^2
      PRINT X
      GOTO 5
10    CONTINUE
      RETURN X
      X = 20
      PRINT X
      PRINT SQUARES(5)
      PRINT X
      STOP
5     X = 10
      GOTO 10
      END

max_stars_repo_name: ElliotAlexander/EDSAC-FORTRAN-Compiler

max_issues_repo_name: ElliotAlexander/EDSAC-FORTRAN-Compiler

max_forks_repo_name: ElliotAlexander/EDSAC-FORTRAN-Compiler