text:
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc                                                                    cc
cc                 mnspro : sproot test program                       cc
cc      application : to find the intersection of a planar            cc
cc      cubic spline curve   x = sx(u)   y = sy(u)   with             cc
cc      a straight line   alfa*x + beta*y = gamma                     cc
cc                                                                    cc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real t(13),c(26),zero(20),sp(40),cc(13)
      integer i,idim,ier,is,i1,i2,j,k,k1,l1,l2,m,mest,n,nc,nk1
      real alfa,beta,gamma,per
c  we have a planar curve
      idim = 2
c  we have a cubic spline curve.
      k = 3
      k1 = k+1
c  set up the dimension information
      nc = 26
      mest = 20
c  n denotes the total number of knots.
      n = 13
c  set up the knots of the spline curve
      t(4) = 0.
      t(5) = 0.2e0
      t(6) = 0.3e0
      t(7) = 0.5e0
      t(8) = 0.6e0
      t(9) = 0.7e0
      t(10) = 0.1e+01
c  fetch the b-spline coefficients for sx(u)
      c(1) = 0.1e+01
      c(2) = 0.3e+01
      c(3) = 0.4e+01
      c(4) = 0.5e+01
      c(5) = 0.3e+01
      c(6) = -0.1e+01
c  fetch the b-spline coefficients for sy(u)
      c(14) = 0.1e+01
      c(15) = 0.2e+01
      c(16) = -0.3e+01
      c(17) = 0.2e+01
      c(18) = 0.1e+01
      c(19) = 0.4e+01
c  we have a closed curve.
c  incorporate the boundary conditions for periodic splines
      per = t(10)-t(4)
      do 10 i=1,3
c  the boundary knots
        t(i) = t(i+6)-per
        t(i+10) = t(i+4)+per
c  the boundary coefficients
        c(i+6) = c(i)
        c(i+19) = c(i+13)
  10  continue
c  print the data of the spline curve.
      write(6,900) k
      write(6,905)
      write(6,910) (t(i),i=1,n)
      write(6,915)
      nk1 = n-k1
      write(6,920) (c(i),i=1,nk1)
      write(6,925)
      i1 = n+1
      i2 = n+nk1
      write(6,920) (c(i),i=i1,i2)
c  loop for the different lines.
      do 200 is=1,5
        go to (110,120,130,140,150),is
c  fetch the parameters of the straight line.
 110    alfa = 0.
        beta = 0.1e+01
        gamma = 0.
        go to 160
 120    alfa = 0.1e+01
        beta = 0.
        go to 160
 130    beta = -0.1e+01
        go to 160
 140    alfa = 0.4e0
        beta = 0.3e0
        gamma = 0.12e+01
        go to 160
 150    beta = 0.4e0
        gamma = 0.
c  print the parameters of the straight line.
 160    write(6,930) alfa,beta,gamma
c  calculate the coefficients of s(u) = sx(u)*alfa + sy(u)*beta - gamma
        do 170 i=1,nk1
          j = i+n
          cc(i) = alfa*c(i)+beta*c(j)-gamma
 170    continue
c  find the zeros of s(u)
        call sproot(t,n,cc,zero,mest,m,ier)
        write(6,935) m
        if(m.eq.0) go to 200
c  find the intersection points
        call curev(idim,t,n,c,nc,k,zero,m,sp,nc,ier)
c  print the intersection points
        write(6,940)
        l2 = 0
        do 180 i=1,m
          l1 = l2+1
          l2 = l1+1
          write(6,945) i,zero(i),sp(l1),sp(l2)
 180    continue
 200  continue
      stop
c  format statements.
 900  format(31h0degree of the spline curve k =,i2)
 905  format(1x,21hposition of the knots)
 910  format(5x,7f6.1)
 915  format(1x,30hb-spline coefficients of sx(u))
 920  format(5x,14f5.0)
 925  format(1x,30hb-spline coefficients of sy(u))
 930  format(18h0intersection with,f6.1,5h *x +,f5.1,5h *y =,f5.1)
 935  format(1x,33hnumber of intersection points m =,i3)
 940  format(6x,1hi,7x,4hu(i),5x,8hsx(u(i)),4x,8hsy(u(i)))
 945  format(1x,i6,3f12.5)
      end

max_issues_repo_name: tbeu/ApproxSpline


text:
module Utils
    use Const
    use omp_lib
    implicit none

contains

    pure function norm2(x) result(a)
        real(mpc), dimension(:), intent(in) :: x
        real(mpc) :: a
        a = sqrt(sum(x**2))
    end function norm2

    pure function zeron(n) result(z)
        integer, intent(in) :: n
        real(mpc), dimension(n) :: z
        z = 0.0_mpc
    end function

    function identity(x) result(x1)
        real(mpc), intent(in), dimension(:) :: x
        real(mpc), dimension(size(x)) :: x1
        x1 = x
    end function

    function dyad_product(a,b) result(c)
        real(mpc), intent(in), dimension(:) :: a,b
        real(mpc), dimension(size(a),size(b)) :: c
        c = spread(a,dim=2,ncopies=size(b)) * spread(b,dim=1,ncopies=size(a))
    end function dyad_product

    pure function symplprod(a,b) result(c)
        real(mpc), intent(in), dimension(:) :: a
        real(mpc), intent(in), dimension(size(a)) :: b
        real(mpc), dimension(size(a)) :: c
        integer :: i,j,k,n
        n = size(a)
        c = 0
        do i = 1, n
            do j = 1, n
                do k = 1, n
                    c(i) = c(i) + asimeps(i,j,k)*a(j)*b(k)
                end do
            end do
        end do
    contains
        ! вроде должно работать
        pure function asimeps(i,j,k) result(e)
            integer, intent(in) :: i,j,k
            integer :: e
            e = (k-j)*(j-i)*(k-i)
            e = e/abs(e)
        end function asimeps
    end function

    subroutine reverse_cols(a)
        real(mpc), dimension(:, :), intent(inout) :: a
!         real(mpc), dimension(size(a,1), size(a, 2)) :: a1
        integer :: i, m
        m = size(a, 1)
        forall (i = 1:m/2) a(i, :) = a(m+1-i, :)
    end subroutine reverse_cols

    subroutine swap(a, i, j, a_cols)
        real(mpc), dimension(:, :) :: a
        integer :: i, j,k
        logical :: a_cols, cols
        intent(in) ::  i, j, a_cols
        optional :: a_cols
        real(mpc) :: tmp
        cols = .FALSE.
        if(present(a_cols)) cols = a_cols

        if (cols) then
            ! вроде присваивает разным кускам массива
            do k = 1, size(a,2)
                tmp = a(i,k)
                a(i,k) = a(j,k)
                a(j,k) = tmp
            end do
        else
            do k = 1, size(a,1)
                tmp = a(k,i)
                a(k,i) = a(k,j)
                a(k,j) = tmp
            end do
        end if

    end subroutine swap

    ! joins 2 arrays size by side
    function joinarr_sbs(a, b) result(c)
        real(mpc), intent(in), dimension(:)         :: a, b
        real(mpc), dimension(max(size(a), size(b)), 2) :: c
        integer :: i

        do i = 1, size(c,1)
            c(i,1) = a(i)
            c(i,2) = b(i)
        end do
    end function joinarr_sbs

end module Utils

max_stars_repo_name: taxus-d/integrators

max_issues_repo_name: taxus-d/integrators

max_forks_repo_name: taxus-d/integrators


text:
c { dg-do compile }
c
c  g77 used to warn for this case
c 19990218-1.f: In program `test':
c 19990218-1.f:13:
c             double precision function fun(a,b)
c                                       1
c  19990218-1.f:23: (continued):
c             c=fun(a,b)
c               2
c  Global name `fun' at (2) has different type at (1) [info -f g77 M GLOBALS]
c
        double precision function fun(a,b)
        double precision a,b
        print*,'in sub: a,b=',a,b
        fun=a*b
        print*,'in sub: fun=',fun
        return
        end
        program test
        double precision a,b,c
        data a,b/1.0d-46,1.0d0/
        c=fun(a,b) ! { dg-error "Return type mismatch of function" }
        print*,'in main: fun=',c
        end

ext: f

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
*DECK SNSQE
      SUBROUTINE SNSQE (FCN, JAC, IOPT, N, X, FVEC, TOL, NPRINT, INFO,
     +   WA, LWA)
C***BEGIN PROLOGUE  SNSQE
C***PURPOSE  An easy-to-use code to find a zero of a system of N
C            nonlinear functions in N variables by a modification of
C            the Powell hybrid method.
C***LIBRARY   SLATEC
C***CATEGORY  F2A
C***TYPE      SINGLE PRECISION (SNSQE-S, DNSQE-D)
C***KEYWORDS  EASY-TO-USE, NONLINEAR SQUARE SYSTEM,
C             POWELL HYBRID METHOD, ZEROS
C***AUTHOR  Hiebert, K. L., (SNLA)
C***DESCRIPTION
C
C 1. Purpose.
C
C
C       The purpose of SNSQE is to find a zero of a system of N non-
C       linear functions in N variables by a modification of the Powell
C       hybrid method.  This is done by using the more general nonlinear
C       equation solver SNSQ.  The user must provide a subroutine which
C       calculates the functions.  The user has the option of either to
C       provide a subroutine which calculates the Jacobian or to let the
C       code calculate it by a forward-difference approximation.  This
C       code is the combination of the MINPACK codes (Argonne) HYBRD1
C       and HYBRJ1.
C
C
C 2. Subroutine and Type Statements.
C
C       SUBROUTINE SNSQE(FCN,JAC,IOPT,N,X,FVEC,TOL,NPRINT,INFO,
C      *                  WA,LWA)
C       INTEGER IOPT,N,NPRINT,INFO,LWA
C       REAL TOL
C       REAL X(N),FVEC(N),WA(LWA)
C       EXTERNAL FCN,JAC
C
C
C 3. Parameters.
C
C       Parameters designated as input parameters must be specified on
C       entry to SNSQE and are not changed on exit, while parameters
C       designated as output parameters need not be specified on entry
C       and are set to appropriate values on exit from SNSQE.
C
C       FCN is the name of the user-supplied subroutine which calculates
C         the functions.  FCN must be declared in an EXTERNAL statement
C         in the user calling program, and should be written as follows.
C
C         SUBROUTINE FCN(N,X,FVEC,IFLAG)
C         INTEGER N,IFLAG
C         REAL X(N),FVEC(N)
C         ----------
C         Calculate the functions at X and
C         return this vector in FVEC.
C         ----------
C         RETURN
C         END
C
C         The value of IFLAG should not be changed by FCN unless the
C         user wants to terminate execution of SNSQE.  In this case, set
C         IFLAG to a negative integer.
C
C       JAC is the name of the user-supplied subroutine which calculates
C         the Jacobian.  If IOPT=1, then JAC must be declared in an
C         EXTERNAL statement in the user calling program, and should be
C         written as follows.
C
C         SUBROUTINE JAC(N,X,FVEC,FJAC,LDFJAC,IFLAG)
C         INTEGER N,LDFJAC,IFLAG
C         REAL X(N),FVEC(N),FJAC(LDFJAC,N)
C         ----------
C         Calculate the Jacobian at X and return this
C         matrix in FJAC.  FVEC contains the function
C         values at X and should not be altered.
C         ----------
C         RETURN
C         END
C
C         The value of IFLAG should not be changed by JAC unless the
C         user wants to terminate execution of SNSQE.  In this case, set
C         IFLAG to a negative integer.
C
C         If IOPT=2, JAC can be ignored (treat it as a dummy argument).
C
C       IOPT is an input variable which specifies how the Jacobian will
C         be calculated.  If IOPT=1, then the user must supply the
C         Jacobian through the subroutine JAC.  If IOPT=2, then the
C         code will approximate the Jacobian by forward-differencing.
C
C       N is a positive integer input variable set to the number of
C         functions and variables.
C
C       X is an array of length N.  On input, X must contain an initial
C         estimate of the solution vector.  On output, X contains the
C         final estimate of the solution vector.
C
C       FVEC is an output array of length N which contains the functions
C         evaluated at the output X.
C
C       TOL is a non-negative input variable.  Termination occurs when
C         the algorithm estimates that the relative error between X and
C         the solution is at most TOL.  Section 4 contains more details
C         about TOL.
C
C       NPRINT is an integer input variable that enables controlled
C         printing of iterates if it is positive.  In this case, FCN is
C         called with IFLAG = 0 at the beginning of the first iteration
C         and every NPRINT iteration thereafter and immediately prior
C         to return, with X and FVEC available for printing. Appropriate
C         print statements must be added to FCN (see example). If NPRINT
C         is not positive, no special calls of FCN with IFLAG = 0 are
C         made.
C
C       INFO is an integer output variable.  If the user has terminated
C         execution, INFO is set to the (negative) value of IFLAG.  See
C         description of FCN and JAC. Otherwise, INFO is set as follows.
C
C         INFO = 0  improper input parameters.
C
C         INFO = 1  algorithm estimates that the relative error between
C                   X and the solution is at most TOL.
C
C         INFO = 2  number of calls to FCN has reached or exceeded
C                   100*(N+1) for IOPT=1 or 200*(N+1) for IOPT=2.
C
C         INFO = 3  TOL is too small.  No further improvement in the
C                   approximate solution X is possible.
C
C         INFO = 4  iteration is not making good progress.
C
C         Sections 4 and 5 contain more details about INFO.
C
C       WA is a work array of length LWA.
C
C       LWA is a positive integer input variable not less than
C         (3*N**2+13*N))/2.
C
C
C 4. Successful Completion.
C
C       The accuracy of SNSQE is controlled by the convergence parame-
C       ter TOL.  This parameter is used in a test which makes a compar-
C       ison between the approximation X and a solution XSOL.  SNSQE
C       terminates when the test is satisfied.  If TOL is less than the
C       machine precision (as defined by the function R1MACH(4)), then
C       SNSQE attempts only to satisfy the test defined by the machine
C       precision.  Further progress is not usually possible.  Unless
C       high precision solutions are required, the recommended value
C       for TOL is the square root of the machine precision.
C
C       The test assumes that the functions are reasonably well behaved,
C       and, if the Jacobian is supplied by the user, that the functions
C       and the Jacobian  coded consistently.  If these conditions
C       are not satisfied, SNSQE may incorrectly indicate convergence.
C       The coding of the Jacobian can be checked by the subroutine
C       CHKDER.  If the Jacobian is coded correctly or IOPT=2, then
C       the validity of the answer can be checked, for example, by
C       rerunning SNSQE with a tighter tolerance.
C
C       Convergence Test.  If ENORM(Z) denotes the Euclidean norm of a
C         vector Z, then this test attempts to guarantee that
C
C               ENORM(X-XSOL) .LE.  TOL*ENORM(XSOL).
C
C         If this condition is satisfied with TOL = 10**(-K), then the
C         larger components of X have K significant decimal digits and
C         INFO is set to 1.  There is a danger that the smaller compo-
C         nents of X may have large relative errors, but the fast rate
C         of convergence of SNSQE usually avoids this possibility.
C
C
C 5. Unsuccessful Completion.
C
C       Unsuccessful termination of SNSQE can be due to improper input
C       parameters, arithmetic interrupts, an excessive number of func-
C       tion evaluations, errors in the functions, or lack of good prog-
C       ress.
C
C       Improper Input Parameters.  INFO is set to 0 if IOPT .LT. 1, or
C         IOPT .GT. 2, or N .LE. 0, or TOL .LT. 0.E0, or
C         LWA .LT. (3*N**2+13*N)/2.
C
C       Arithmetic Interrupts.  If these interrupts occur in the FCN
C         subroutine during an early stage of the computation, they may
C         be caused by an unacceptable choice of X by SNSQE.  In this
C         case, it may be possible to remedy the situation by not evalu-
C         ating the functions here, but instead setting the components
C         of FVEC to numbers that exceed those in the initial FVEC.
C
C       Excessive Number of Function Evaluations.  If the number of
C         calls to FCN reaches 100*(N+1) for IOPT=1 or 200*(N+1) for
C         IOPT=2, then this indicates that the routine is converging
C         very slowly as measured by the progress of FVEC, and INFO is
C         set to 2.  This situation should be unusual because, as
C         indicated below, lack of good progress is usually diagnosed
C         earlier by SNSQE, causing termination with INFO = 4.
C
C       Errors in the Functions.  When IOPT=2, the choice of step length
C         in the forward-difference approximation to the Jacobian
C         assumes that the relative errors in the functions are of the
C         order of the machine precision.  If this is not the case,
C         SNSQE may fail (usually with INFO = 4).  The user should
C         then either use SNSQ and set the step length or use IOPT=1
C         and supply the Jacobian.
C
C       Lack of Good Progress.  SNSQE searches for a zero of the system
C         by minimizing the sum of the squares of the functions.  In so
C         doing, it can become trapped in a region where the minimum
C         does not correspond to a zero of the system and, in this situ-
C         ation, the iteration eventually fails to make good progress.
C         In particular, this will happen if the system does not have a
C         zero.  If the system has a zero, rerunning SNSQE from a dif-
C         ferent starting point may be helpful.
C
C
C 6. Characteristics of the Algorithm.
C
C       SNSQE is a modification of the Powell hybrid method.  Two of
C       its main characteristics involve the choice of the correction as
C       a convex combination of the Newton and scaled gradient direc-
C       tions, and the updating of the Jacobian by the rank-1 method of
C       Broyden.  The choice of the correction guarantees (under reason-
C       able conditions) global convergence for starting points far from
C       the solution and a fast rate of convergence.  The Jacobian is
C       calculated at the starting point by either the user-supplied
C       subroutine or a forward-difference approximation, but it is not
C       recalculated until the rank-1 method fails to produce satis-
C       factory progress.
C
C       Timing.  The time required by SNSQE to solve a given problem
C         depends on N, the behavior of the functions, the accuracy
C         requested, and the starting point.  The number of arithmetic
C         operations needed by SNSQE is about 11.5*(N**2) to process
C         each evaluation of the functions (call to FCN) and 1.3*(N**3)
C         to process each evaluation of the Jacobian (call to JAC,
C         if IOPT = 1).  Unless FCN and JAC can be evaluated quickly,
C         the timing of SNSQE will be strongly influenced by the time
C         spent in FCN and JAC.
C
C       Storage.  SNSQE requires (3*N**2 + 17*N)/2 single precision
C         storage locations, in addition to the storage required by the
C         program.  There are no internally declared storage arrays.
C
C
C 7. Example.
C
C       The problem is to determine the values of X(1), X(2), ..., X(9),
C       which solve the system of tridiagonal equations
C
C       (3-2*X(1))*X(1)           -2*X(2)                   = -1
C               -X(I-1) + (3-2*X(I))*X(I)         -2*X(I+1) = -1, I=2-8
C                                   -X(8) + (3-2*X(9))*X(9) = -1
C
C       **********
C
C       PROGRAM TEST
C C
C C     Driver for SNSQE example.
C C
C       INTEGER J,N,IOPT,NPRINT,INFO,LWA,NWRITE
C       REAL TOL,FNORM
C       REAL X(9),FVEC(9),WA(180)
C       REAL ENORM,R1MACH
C       EXTERNAL FCN
C       DATA NWRITE /6/
C C
C       IOPT = 2
C       N = 9
C C
C C     The following starting values provide a rough solution.
C C
C       DO 10 J = 1, 9
C          X(J) = -1.E0
C    10    CONTINUE
C
C       LWA = 180
C       NPRINT = 0
C C
C C     Set TOL to the square root of the machine precision.
C C     Unless high precision solutions are required,
C C     this is the recommended setting.
C C
C       TOL = SQRT(R1MACH(4))
C C
C       CALL SNSQE(FCN,JAC,IOPT,N,X,FVEC,TOL,NPRINT,INFO,WA,LWA)
C       FNORM = ENORM(N,FVEC)
C       WRITE (NWRITE,1000) FNORM,INFO,(X(J),J=1,N)
C       STOP
C  1000 FORMAT (5X,' FINAL L2 NORM OF THE RESIDUALS',E15.7 //
C      *        5X,' EXIT PARAMETER',16X,I10 //
C      *        5X,' FINAL APPROXIMATE SOLUTION' // (5X,3E15.7))
C       END
C       SUBROUTINE FCN(N,X,FVEC,IFLAG)
C       INTEGER N,IFLAG
C       REAL X(N),FVEC(N)
C       INTEGER K
C       REAL ONE,TEMP,TEMP1,TEMP2,THREE,TWO,ZERO
C       DATA ZERO,ONE,TWO,THREE /0.E0,1.E0,2.E0,3.E0/
C C
C       DO 10 K = 1, N
C          TEMP = (THREE - TWO*X(K))*X(K)
C          TEMP1 = ZERO
C          IF (K .NE. 1) TEMP1 = X(K-1)
C          TEMP2 = ZERO
C          IF (K .NE. N) TEMP2 = X(K+1)
C          FVEC(K) = TEMP - TEMP1 - TWO*TEMP2 + ONE
C    10    CONTINUE
C       RETURN
C       END
C
C       Results obtained with different compilers or machines
C       may be slightly different.
C
C       FINAL L2 NORM OF THE RESIDUALS  0.1192636E-07
C
C       EXIT PARAMETER                         1
C
C       FINAL APPROXIMATE SOLUTION
C
C       -0.5706545E+00 -0.6816283E+00 -0.7017325E+00
C       -0.7042129E+00 -0.7013690E+00 -0.6918656E+00
C       -0.6657920E+00 -0.5960342E+00 -0.4164121E+00
C
C***REFERENCES  M. J. D. Powell, A hybrid method for nonlinear equa-
C                 tions. In Numerical Methods for Nonlinear Algebraic
C                 Equations, P. Rabinowitz, Editor.  Gordon and Breach,
C                 1988.
C***ROUTINES CALLED  SNSQ, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   800301  DATE WRITTEN
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  SNSQE
      INTEGER IOPT,N,NPRINT,INFO,LWA
      REAL TOL
      REAL X(*),FVEC(*),WA(LWA)
      EXTERNAL FCN, JAC
      INTEGER INDEX,J,LR,MAXFEV,ML,MODE,MU,NFEV,NJEV
      REAL EPSFCN,FACTOR,ONE,XTOL,ZERO
      SAVE FACTOR, ONE, ZERO
      DATA FACTOR,ONE,ZERO /1.0E2,1.0E0,0.0E0/
C***FIRST EXECUTABLE STATEMENT  SNSQE
      INFO = 0
C
C     CHECK THE INPUT PARAMETERS FOR ERRORS.
C
      IF (IOPT .LT. 1 .OR. IOPT .GT. 2 .OR. N .LE. 0
     1    .OR. TOL .LT. ZERO .OR. LWA .LT. (3*N**2 +13*N)/2)
     2   GO TO 20
C
C     CALL SNSQ.
C
      MAXFEV = 100*(N + 1)
      IF (IOPT .EQ. 2) MAXFEV = 2*MAXFEV
      XTOL = TOL
      ML = N - 1
      MU = N - 1
      EPSFCN = ZERO
      MODE = 2
      DO 10 J = 1, N
         WA(J) = ONE
   10    CONTINUE
      LR = (N*(N + 1))/2
      INDEX=6*N+LR
      CALL SNSQ(FCN,JAC,IOPT,N,X,FVEC,WA(INDEX+1),N,XTOL,MAXFEV,ML,MU,
     1           EPSFCN,WA(1),MODE,FACTOR,NPRINT,INFO,NFEV,NJEV,
     2           WA(6*N+1),LR,WA(N+1),WA(2*N+1),WA(3*N+1),WA(4*N+1),
     3           WA(5*N+1))
      IF (INFO .EQ. 5) INFO = 4
   20 CONTINUE
      IF (INFO .EQ. 0) CALL XERMSG ('SLATEC', 'SNSQE',
     +   'INVALID INPUT PARAMETER.', 2, 1)
      RETURN
C
C     LAST CARD OF SUBROUTINE SNSQE.
C
      END


text:
*----------------------------------------------------------------------*
!>     allocate a new slot for an operator structure in op_info
!>     the list is extended and the array is updated
!>     @param[in] label name of the new operator
!>     @param[inout] op_info information struct for all operators
*----------------------------------------------------------------------*
      subroutine add_operator(label,op_info)
*----------------------------------------------------------------------*

      implicit none

      include 'opdim.h'
      include 'mdef_operator_info.h'

      type(operator_info), intent(inout), target ::
     &     op_info
      character ::
     &     label*(*)

      type(operator_list), pointer ::
     &     list_pnt

      integer ::
     &     iops

      list_pnt => op_info%op_list
      ! advance to end of operator list:
      do while (associated(list_pnt%next))
        list_pnt => list_pnt%next
      end do
      ! is last entry already in use?
      if (associated(list_pnt%op)) then
        allocate(list_pnt%next)
        list_pnt%next%prev => list_pnt
        list_pnt => list_pnt%next
        list_pnt%next => null()
      end if
      allocate (list_pnt%op)

      ! add user-supplied label
      if (len_trim(label).gt.len_opname)
     &   call quit(1,'add_operator','name too long: "'
     &     //trim(label)//'"')
      list_pnt%op%name(1:len_opname) = ' '
      list_pnt%op%name = trim(label)
      list_pnt%op%assoc_list(1:2*len_opname+2) = ' '

      ! init all pointers
      list_pnt%op%ihpvca_occ => null()
      list_pnt%op%ica_occ => null()
      list_pnt%op%igasca_restr => null()

      ! increment counter
      op_info%nops = op_info%nops+1

      ! update operator array
      call update_op_arr(op_info)

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
!
! PP_MD
!
! Module containing forward, tangent-linear and adjoint subroutines
! for partial pressure to/from mass density units conversion.
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 24-Sep-2009
!                       paul.vandelst@noaa.gov
!

MODULE PP_MD

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Modules used
  USE Type_Kinds, ONLY: fp
  USE Profile_Utility_Parameters, ONLY: ZERO, ONE, &
                                        MAX_N_MOLECULAR_SPECIES, &
                                        R0, MOLECULAR_WEIGHT, &
                                        PA_TO_HPA, HPA_TO_PA
  ! Disable all implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  PUBLIC :: PP_to_MD, PP_to_MD_TL, PP_to_MD_AD
  PUBLIC :: MD_to_PP, MD_to_PP_TL, MD_to_PP_AD

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &

CONTAINS

!------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       PP_to_MD
!
! PURPOSE:
!       Elemental subroutine to convert gas concentrations from partial
!       pressure to mass density.
!
! CALLING SEQUENCE:
!       CALL PP_to_MD( Temperature              , &  ! Input
!                      Partial_Pressure         , &  ! Input
!                      Mass_Density             , &  ! Output
!                      Molecule_ID = Molecule_ID  )  ! Optional Input
!
! INPUTS:
!       Temperature:       Atmospheric temperature
!                          UNITS:      Kelvin, K
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Scalar or any rank
!                          ATTRIBUTES: INTENT(IN)
!
!       Partial_Pressure:  Partial pressure of gas.
!                          UNITS:      hectoPascals, hPa
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Same as Temperature
!                          ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Mass_Density:      Mass density of gas.
!                          Set to zero for invalid input.
!                          UNITS:      g/m^3
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Same as input Temperature
!                          ATTRIBUTES: INTENT(OUT)
!
! OPTIONAL INPUTS:
!       Molecule_ID:       HITRAN molecular designation identifying the
!                          molecule for which the concentration units
!                          conversion is required. If not specified, the
!                          default value is that for water vapor.
!                          Valid values are:
!                            1: H2O       9: SO2      17: HI       25: H2O2
!                            2: CO2      10: NO2      18: ClO      26: C2H2
!                            3: O3       11: NH3      19: OCS      27: C2H6
!                            4: N2O      12: HNO3     20: H2CO     28: PH3
!                            5: CO       13: OH       21: HOCl     29: COF2
!                            6: CH4      14: HF       22: N2       30: SF6
!                            7: O2       15: HCl      23: HCN      31: H2S
!                            8: NO       16: HBr      24: CH3Cl    32: HCOOH
!                          UNITS:      N/A
!                          TYPE:       INTEGER
!                          DIMENSION:  Scalar
!                          ATTRIBUTES: OPTIONAL, INTENT(IN)
!
! PROCEDURE:
!       The ideal gas law is
!
!         p.V = n.R.T
!
!       where,
!         p = Pressure
!         V = volume
!         n = number of moles of the gas,
!         R = molar gas constant,
!         T = Temperature.
!
!       This can be written in terms of the number of molecules, N,
!
!                 N
!         p.V = ---- . R.T                  ..............(1)
!                NA
!
!       where NA = Avogadro's constant (number of molecules in 1 mole.)
!
!       Eqn.(1) can be recast to provide the number of molecules in the
!       volume, the number density, nd,
!
!               N     p.NA
!         nd = --- = ------  molecules/m^3  ..............(2)
!               V     R.T
!
!       The mass density is related to the number density by the simple
!       relation,
!
!                    MW
!         md = nd . ----                    ..............(3)
!                    NA
!
!       Substituting equation (2) into (3) gives,
!
!               p.MW
!         md = ------  g/m^3                ..............(4)
!               R.T
!
!       The value of p0 used in this routine is expressed in pascals (Pa) whereas
!       the input partial Pressure is expected to be hPa (more common unit). Thus
!       there is a factor of 100 to include,
!
!               100.p.MW
!         md = ----------  g/m^3
!                 R.T
!
! UNITS ANALYSIS:
!       p  :  kg.m^-1.s^-2  (Pascals)
!       MW :  g.mol^-1
!       R  :  J.mol^-1.K^-1 == kg.m^2.s^-2.mol^-1.K^-1
!       T  :  K
!
!                  kg       g      s^2.mol.K     1
!         md == ------- . ----- . ----------- . ---
!                m.s^2     mol      kg.m^2       K
!
!                 g
!            == -----
!                m^3
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 08-Nov-2004
!                       paul.vandelst@noaa.gov
!:sdoc-:
!------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE PP_to_MD( &
    Temperature     , &  ! Input
    Partial_Pressure, &  ! Input
    Mass_Density    , &  ! Output
    Molecule_ID       )  ! Optional Input
    ! Arguments
    REAL(fp),          INTENT(IN)  :: Temperature
    REAL(fp),          INTENT(IN)  :: Partial_Pressure
    REAL(fp),          INTENT(OUT) :: Mass_Density
    INTEGER, OPTIONAL, INTENT(IN)  :: Molecule_ID
    ! Local variables
    INTEGER :: Id

    ! Setup
    Mass_Density = ZERO
    IF ( Temperature < ZERO .OR. Partial_Pressure < ZERO ) RETURN
    IF ( PRESENT( Molecule_ID ) ) THEN
      IF ( Molecule_ID < 1 .OR. Molecule_ID > MAX_N_MOLECULAR_SPECIES ) RETURN
      Id = Molecule_ID
    ELSE
      Id = 1  ! Default value is for water vapor
    END IF

    ! Convert partial pressure to mass density
    Mass_Density = HPA_TO_PA * Partial_Pressure * MOLECULAR_WEIGHT(Id) / (Temperature * R0)

  END SUBROUTINE PP_to_MD

!------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       PP_to_MD_TL
!
! PURPOSE:
!       Tangent-linear form of elemental subroutine to convert gas
!       concentrations from partial pressure to mass density.
!
! CALLING SEQUENCE:
!       CALL PP_to_MD_TL( Temperature              , &  ! FWD Input
!                         Partial_Pressure         , &  ! FWD Input
!                         Temperature_TL           , &  ! TL  Input
!                         Partial_Pressure_TL      , &  ! TL  Input
!                         Mass_Density_TL          , &  ! TL  Output
!                         Molecule_ID = Molecule_ID  )  ! Optional Input
!
! INPUTS:
!       Temperature:          Atmospheric temperature
!                             UNITS:      Kelvin, K
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Scalar or any rank
!                             ATTRIBUTES: INTENT(IN)
!
!       Partial_Pressure:     Partial pressure of gas.
!                             UNITS:      hectoPascals, hPa
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as Temperature
!                             ATTRIBUTES: INTENT(IN)
!
!       Temperature_TL:       Tangent-linear atmospheric temperature
!                             UNITS:      Kelvin, K
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as Temperature
!                             ATTRIBUTES: INTENT(IN)
!
!       Partial_Pressure_TL:  Tangent-linear partial pressure of gas.
!                             UNITS:      hectoPascals, hPa
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as Temperature
!                             ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Mass_Density_TL:      Tangent-linear mass density of gas.
!                             Set to zero for invalid input.
!                             UNITS:      g/m^3
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as input Temperature
!                             ATTRIBUTES: INTENT(OUT)
!
! OPTIONAL INPUTS:
!       Molecule_ID:          HITRAN molecular designation identifying the
!                             molecule for which the concentration units
!                             conversion is required. If not specified, the
!                             default value is that for water vapor.
!                             Valid values are:
!                               1: H2O       9: SO2      17: HI       25: H2O2
!                               2: CO2      10: NO2      18: ClO      26: C2H2
!                               3: O3       11: NH3      19: OCS      27: C2H6
!                               4: N2O      12: HNO3     20: H2CO     28: PH3
!                               5: CO       13: OH       21: HOCl     29: COF2
!                               6: CH4      14: HF       22: N2       30: SF6
!                               7: O2       15: HCl      23: HCN      31: H2S
!                               8: NO       16: HBr      24: CH3Cl    32: HCOOH
!                             UNITS:      N/A
!                             TYPE:       INTEGER
!                             DIMENSION:  Scalar
!                             ATTRIBUTES: OPTIONAL, INTENT(IN)
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 24-Sep-2009
!                       paul.vandelst@noaa.gov
!:sdoc-:
!------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE PP_to_MD_TL( &
    Temperature        , &  ! FWD Input
    Partial_Pressure   , &  ! FWD Input
    Temperature_TL     , &  ! TL  Input
    Partial_Pressure_TL, &  ! TL  Input
    Mass_Density_TL    , &  ! TL  Output
    Molecule_ID          )  ! Optional Input
    ! Arguments
    REAL(fp),          INTENT(IN)  :: Temperature
    REAL(fp),          INTENT(IN)  :: Partial_Pressure
    REAL(fp),          INTENT(IN)  :: Temperature_TL
    REAL(fp),          INTENT(IN)  :: Partial_Pressure_TL
    REAL(fp),          INTENT(OUT) :: Mass_Density_TL
    INTEGER, OPTIONAL, INTENT(IN)  :: Molecule_ID
    ! Local variables
    INTEGER :: Id
    REAL(fp) :: c

    ! Setup
    Mass_Density_TL = ZERO
    IF ( Temperature < ZERO .OR. Partial_Pressure < ZERO ) RETURN
    IF ( PRESENT( Molecule_ID ) ) THEN
      IF ( Molecule_ID < 1 .OR. Molecule_ID > MAX_N_MOLECULAR_SPECIES ) RETURN
      Id = Molecule_ID
    ELSE
      Id = 1  ! Default value is for water vapor
    END IF

    ! Tangent-linear form of partial pressure to mass density conversion
    c = HPA_TO_PA * MOLECULAR_WEIGHT(Id) / R0
    Mass_Density_TL = c * ((Temperature * Partial_Pressure_TL) - &
                           (Partial_Pressure * Temperature_TL)) / Temperature**2

  END SUBROUTINE PP_to_MD_TL

!------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       PP_to_MD_AD
!
! PURPOSE:
!       Adjoint form of elemental subroutine to convert gas
!       concentrations from partial pressure to mass density.
!
! CALLING SEQUENCE:
!       CALL PP_to_MD_AD( Temperature              , &  ! FWD Input
!                         Partial_Pressure         , &  ! FWD Input
!                         Mass_Density_AD          , &  ! AD  Input
!                         Temperature_AD           , &  ! AD  Output
!                         Partial_Pressure_AD      , &  ! AD  Output
!                         Molecule_ID = Molecule_ID  )  ! Optional Input
!
! INPUTS:
!       Temperature:          Atmospheric temperature
!                             UNITS:      Kelvin, K
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Scalar or any rank
!                             ATTRIBUTES: INTENT(IN)
!
!       Partial_Pressure:     Partial pressure of gas.
!                             UNITS:      hectoPascals, hPa
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as Temperature
!                             ATTRIBUTES: INTENT(IN)
!
!       Mass_Density_AD:      Adjoint mass density of gas.
!                             *** SET TO ZERO ON EXIT ***
!                             UNITS:      g/m^3
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as Temperature
!                             ATTRIBUTES: INTENT(IN OUT)
!
! OUTPUTS:
!       Temperature_AD:       Adjoint atmospheric temperature
!                             *** MUST HAVE VALUE ON ENTRY ***
!                             UNITS:      Kelvin, K
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as input Temperature
!                             ATTRIBUTES: INTENT(IN OUT)
!
!       Partial_Pressure_AD:  Adjoint partial pressure of gas.
!                             *** MUST HAVE VALUE ON ENTRY ***
!                             UNITS:      hectoPascals, hPa
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as input Temperature
!                             ATTRIBUTES: INTENT(IN OUT)
!
! OPTIONAL INPUTS:
!       Molecule_ID:          HITRAN molecular designation identifying the
!                             molecule for which the concentration units
!                             conversion is required. If not specified, the
!                             default value is that for water vapor.
!                             Valid values are:
!                               1: H2O       9: SO2      17: HI       25: H2O2
!                               2: CO2      10: NO2      18: ClO      26: C2H2
!                               3: O3       11: NH3      19: OCS      27: C2H6
!                               4: N2O      12: HNO3     20: H2CO     28: PH3
!                               5: CO       13: OH       21: HOCl     29: COF2
!                               6: CH4      14: HF       22: N2       30: SF6
!                               7: O2       15: HCl      23: HCN      31: H2S
!                               8: NO       16: HBr      24: CH3Cl    32: HCOOH
!                             UNITS:      N/A
!                             TYPE:       INTEGER
!                             DIMENSION:  Scalar
!                             ATTRIBUTES: OPTIONAL, INTENT(IN)
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 24-Sep-2009
!                       paul.vandelst@noaa.gov
!:sdoc-:
!------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE PP_to_MD_AD( &
    Temperature        , &  ! FWD Input
    Partial_Pressure   , &  ! FWD Input
    Mass_Density_AD    , &  ! AD  Input
    Temperature_AD     , &  ! AD  Output
    Partial_Pressure_AD, &  ! AD  Output
    Molecule_ID          )  ! Optional Input
    ! Arguments
    REAL(fp),          INTENT(IN)     :: Temperature
    REAL(fp),          INTENT(IN)     :: Partial_Pressure
    REAL(fp),          INTENT(IN OUT) :: Mass_Density_AD
    REAL(fp),          INTENT(IN OUT) :: Temperature_AD
    REAL(fp),          INTENT(IN OUT) :: Partial_Pressure_AD
    INTEGER, OPTIONAL, INTENT(IN)     :: Molecule_ID
    ! Local variables
    INTEGER :: Id
    REAL(fp) :: c

    ! Setup
    IF ( Temperature < ZERO .OR. Partial_Pressure < ZERO ) THEN
      Mass_Density_AD     = ZERO
      Temperature_AD      = ZERO
      Partial_Pressure_AD = ZERO
      RETURN
    END IF
    IF ( PRESENT( Molecule_ID ) ) THEN
      IF ( Molecule_ID < 1 .OR. Molecule_ID > MAX_N_MOLECULAR_SPECIES ) THEN
        Mass_Density_AD     = ZERO
        Temperature_AD      = ZERO
        Partial_Pressure_AD = ZERO
        RETURN
      END IF
      Id = Molecule_ID
    ELSE
      Id = 1  ! Default value is for water vapor
    END IF

    ! Adjoint form of partial pressure to mass density conversion
    c = HPA_TO_PA * MOLECULAR_WEIGHT(Id) / R0
    Temperature_AD = Temperature_AD - (c * Partial_Pressure * Mass_Density_AD / Temperature**2)
    Partial_Pressure_AD = Partial_Pressure_AD + (c * Mass_Density_AD / Temperature)
    Mass_Density_AD = ZERO

  END SUBROUTINE PP_to_MD_AD

!------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       MD_to_PP
!
! PURPOSE:
!       Elemental subroutine to convert gas concentration mass density
!       to partial pressure.
!
! CALLING SEQUENCE:
!       CALL MD_to_PP( Temperature              , &  ! Input
!                      Mass_Density             , &  ! Input
!                      Partial_Pressure         , &  ! Output
!                      Molecule_ID = Molecule_ID  )  ! Optional Input
!
! INPUTS:
!       Temperature:       Atmospheric temperature
!                          UNITS:      Kelvin, K
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Scalar or any rank
!                          ATTRIBUTES: INTENT(IN)
!
!       Mass_Density:      Mass density of gas.
!                          UNITS:      g/m^3
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Same as Temperature
!                          ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Partial_Pressure:  Partial pressure of gas.
!                          Set to zero for invalid input.
!                          UNITS:      hectoPascals, hPa
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Same as input Temperature
!                          ATTRIBUTES: INTENT(OUT)
!
! OPTIONAL INPUTS:
!       Molecule_ID:       HITRAN molecular designation identifying the
!                          molecule for which the concentration units
!                          conversion is required. If not specified, the
!                          default value is that for water vapor.
!                          Valid values are:
!                            1: H2O       9: SO2      17: HI       25: H2O2
!                            2: CO2      10: NO2      18: ClO      26: C2H2
!                            3: O3       11: NH3      19: OCS      27: C2H6
!                            4: N2O      12: HNO3     20: H2CO     28: PH3
!                            5: CO       13: OH       21: HOCl     29: COF2
!                            6: CH4      14: HF       22: N2       30: SF6
!                            7: O2       15: HCl      23: HCN      31: H2S
!                            8: NO       16: HBr      24: CH3Cl    32: HCOOH
!                          UNITS:      N/A
!                          TYPE:       INTEGER
!                          DIMENSION:  Scalar
!                          ATTRIBUTES: OPTIONAL, INTENT(IN)
!
! PROCEDURE:
!       The ideal gas law is
!
!         p.V = n.R.T
!
!       where,
!         p = Pressure
!         V = volume
!         n = number of moles of the gas,
!         R = molar gas constant,
!         T = Temperature.
!
!       This can be written in terms of the number of molecules, N,
!
!                 N
!         p.V = ---- . R.T                  ..............(1)
!                NA
!
!       where NA = Avogadro's constant (number of molecules in 1 mole.)
!
!       Eqn.(1) can be recast to provide the number of molecules in the
!       volume, the number density, nd,
!
!               N     p.NA
!         nd = --- = ------  molecules/m^3  ..............(2)
!               V     R.T
!
!       The mass density is related to the number density by the simple
!       relation,
!
!                    MW
!         md = nd . ----                    ..............(3)
!                    NA
!
!       Substituting equation (2) into (3) gives,
!
!               p.MW
!         md = ------  g/m^3                ..............(4)
!               R.T
!
!       and inverting equation (4) gives,
!
!              md.R.T
!         p = -------- Pa
!                MW
!
!       The pressure result above is determined in units of Pascals, so there
!       is a factor of 100 to include to return the pressure in units of
!       hectoPascals,
!
!              md.R.T
!         p = -------- hPa
!              100.MW
!
! UNITS ANALYSIS:
!       md :  g.m^-3
!       R  :  J.mol^-1.K^-1 == kg.m^2.s^-2.mol^-1.K^-1
!       T  :  K
!       MW :  g.mol^-1
!
!                g        kg.m^2          mol
!         p == ----- . ----------- . K . -----
!               m^3     s^2.mol.K          g
!
!                 kg
!           == -------
!               m.s^2
!
!               kg.m.s^-2      F
!           == ----------- == --- == pressure.
!                  m^2         A
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 08-Nov-2004
!                       paul.vandelst@noaa.gov
!:sdoc-:
!------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE MD_to_PP( &
    Temperature     , &  ! Input
    Mass_Density    , &  ! Input
    Partial_Pressure, &  ! Output
    Molecule_ID       )  ! Optional Input
    ! Arguments
    REAL(fp),          INTENT(IN)  :: Temperature
    REAL(fp),          INTENT(IN)  :: Mass_Density
    REAL(fp),          INTENT(OUT) :: Partial_Pressure
    INTEGER, OPTIONAL, INTENT(IN)  :: Molecule_ID
    ! Local variables
    INTEGER :: Id

    ! Setup
    Partial_Pressure = ZERO
    IF ( Temperature < ZERO .OR. Mass_Density < ZERO ) RETURN
    IF ( PRESENT( Molecule_ID ) ) THEN
      IF ( Molecule_ID < 1 .OR. Molecule_ID > MAX_N_MOLECULAR_SPECIES ) RETURN
      Id = Molecule_ID
    ELSE
      Id = 1  ! Default value is for water vapor
    END IF

    ! Convert mass density to partial pressure
    Partial_Pressure = PA_TO_HPA * Mass_Density * R0 * Temperature / MOLECULAR_WEIGHT(Id)

  END SUBROUTINE MD_to_PP

!------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       MD_to_PP_TL
!
! PURPOSE:
!       Tangent-linear form of elemental subroutine to convert gas
!       concentration mass density to partial pressure.
!
! CALLING SEQUENCE:
!       CALL MD_to_PP_TL( Temperature              , &  ! FWD Input
!                         Mass_Density             , &  ! FWD Input
!                         Temperature_TL           , &  ! TL  Input
!                         Mass_Density_TL          , &  ! TL  Input
!                         Partial_Pressure_TL      , &  ! TL  Output
!                         Molecule_ID = Molecule_ID  )  ! Optional Input
!
! INPUTS:
!       Temperature:          Atmospheric temperature
!                             UNITS:      Kelvin, K
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Scalar or any rank
!                             ATTRIBUTES: INTENT(IN)
!
!       Mass_Density:         Mass density of gas.
!                             UNITS:      g/m^3
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as Temperature
!                             ATTRIBUTES: INTENT(IN)
!
!       Temperature_TL:       Tangent-linear atmospheric temperature
!                             UNITS:      Kelvin, K
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as Temperature
!                             ATTRIBUTES: INTENT(IN)
!
!       Mass_Density_TL:      Tangent-linear mass density of gas.
!                             UNITS:      g/m^3
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as Temperature
!                             ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Partial_Pressure_TL:  Tangent-linear partial pressure of gas.
!                             Set to zero for invalid input.
!                             UNITS:      hectoPascals, hPa
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as input Temperature
!                             ATTRIBUTES: INTENT(OUT)
!
! OPTIONAL INPUTS:
!       Molecule_ID:          HITRAN molecular designation identifying the
!                             molecule for which the concentration units
!                             conversion is required. If not specified, the
!                             default value is that for water vapor.
!                             Valid values are:
!                               1: H2O       9: SO2      17: HI       25: H2O2
!                               2: CO2      10: NO2      18: ClO      26: C2H2
!                               3: O3       11: NH3      19: OCS      27: C2H6
!                               4: N2O      12: HNO3     20: H2CO     28: PH3
!                               5: CO       13: OH       21: HOCl     29: COF2
!                               6: CH4      14: HF       22: N2       30: SF6
!                               7: O2       15: HCl      23: HCN      31: H2S
!                               8: NO       16: HBr      24: CH3Cl    32: HCOOH
!                             UNITS:      N/A
!                             TYPE:       INTEGER
!                             DIMENSION:  Scalar
!                             ATTRIBUTES: OPTIONAL, INTENT(IN)
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 24-Sep-2009
!                       paul.vandelst@noaa.gov
!:sdoc-:
!------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE MD_to_PP_TL( &
    Temperature        , &  ! FWD Input
    Mass_Density       , &  ! FWD Input
    Temperature_TL     , &  ! TL  Input
    Mass_Density_TL    , &  ! TL  Input
    Partial_Pressure_TL, &  ! TL  Output
    Molecule_ID          )  ! Optional Input
    ! Arguments
    REAL(fp),          INTENT(IN)  :: Temperature
    REAL(fp),          INTENT(IN)  :: Mass_Density
    REAL(fp),          INTENT(IN)  :: Temperature_TL
    REAL(fp),          INTENT(IN)  :: Mass_Density_TL
    REAL(fp),          INTENT(OUT) :: Partial_Pressure_TL
    INTEGER, OPTIONAL, INTENT(IN)  :: Molecule_ID
    ! Local variables
    INTEGER :: Id
    REAL(fp) :: c

    ! Setup
    Partial_Pressure_TL = ZERO
    IF ( Temperature < ZERO .OR. Mass_Density < ZERO ) RETURN
    IF ( PRESENT( Molecule_ID ) ) THEN
      IF ( Molecule_ID < 1 .OR. Molecule_ID > MAX_N_MOLECULAR_SPECIES ) RETURN
      Id = Molecule_ID
    ELSE
      Id = 1  ! Default value is for water vapor
    END IF

    ! Tangent-linear form of mass density to partial pressure conversion
    c = PA_TO_HPA * R0 / MOLECULAR_WEIGHT(Id)
    Partial_Pressure_TL = c * ((Mass_Density * Temperature_TL) + &
                               (Temperature * Mass_Density_TL) )

  END SUBROUTINE MD_to_PP_TL

!------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       MD_to_PP_AD
!
! PURPOSE:
!       Adjoint form of elemental subroutine to convert gas
!       concentration mass density to partial pressure.
!
! CALLING SEQUENCE:
!       CALL MD_to_PP_AD( Temperature              , &  ! FWD Input
!                         Mass_Density             , &  ! FWD Input
!                         Partial_Pressure_AD      , &  ! AD  Input
!                         Temperature_AD           , &  ! AD  Output
!                         Mass_Density_AD          , &  ! AD  Output
!                         Molecule_ID = Molecule_ID  )  ! Optional Input
!
! INPUTS:
!       Temperature:          Atmospheric temperature
!                             UNITS:      Kelvin, K
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Scalar or any rank
!                             ATTRIBUTES: INTENT(IN)
!
!       Mass_Density:         Mass density of gas.
!                             UNITS:      g/m^3
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as Temperature
!                             ATTRIBUTES: INTENT(IN)
!
!       Partial_Pressure_AD:  Adjoint partial pressure of gas.
!                             *** SET TO ZERO ON EXIT ***
!                             UNITS:      hectoPascals, hPa
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as Temperature
!                             ATTRIBUTES: INTENT(IN OUT)
!
! OUTPUTS:
!       Temperature_AD:       Adjoint atmospheric temperature
!                             *** MUST HAVE VALUE ON ENTRY ***
!                             UNITS:      Kelvin, K
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as input Temperature
!                             ATTRIBUTES: INTENT(IN OUT)
!
!       Mass_Density_AD:      Adjoint mass density of gas.
!                             *** MUST HAVE VALUE ON ENTRY ***
!                             UNITS:      g/m^3
!                             TYPE:       REAL(fp)
!                             DIMENSION:  Same as input Temperature
!                             ATTRIBUTES: INTENT(IN OUT)
!
! OPTIONAL INPUTS:
!       Molecule_ID:          HITRAN molecular designation identifying the
!                             molecule for which the concentration units
!                             conversion is required. If not specified, the
!                             default value is that for water vapor.
!                             Valid values are:
!                               1: H2O       9: SO2      17: HI       25: H2O2
!                               2: CO2      10: NO2      18: ClO      26: C2H2
!                               3: O3       11: NH3      19: OCS      27: C2H6
!                               4: N2O      12: HNO3     20: H2CO     28: PH3
!                               5: CO       13: OH       21: HOCl     29: COF2
!                               6: CH4      14: HF       22: N2       30: SF6
!                               7: O2       15: HCl      23: HCN      31: H2S
!                               8: NO       16: HBr      24: CH3Cl    32: HCOOH
!                             UNITS:      N/A
!                             TYPE:       INTEGER
!                             DIMENSION:  Scalar
!                             ATTRIBUTES: OPTIONAL, INTENT(IN)
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 24-Sep-2009
!                       paul.vandelst@noaa.gov
!:sdoc-:
!------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE MD_to_PP_AD( &
    Temperature        , &  ! FWD Input
    Mass_Density       , &  ! FWD Input
    Partial_Pressure_AD, &  ! AD  Input
    Temperature_AD     , &  ! AD  Output
    Mass_Density_AD    , &  ! AD  Output
    Molecule_ID          )  ! Optional Input
    ! Arguments
    REAL(fp),          INTENT(IN)     :: Temperature
    REAL(fp),          INTENT(IN)     :: Mass_Density
    REAL(fp),          INTENT(IN OUT) :: Partial_Pressure_AD
    REAL(fp),          INTENT(IN OUT) :: Temperature_AD
    REAL(fp),          INTENT(IN OUT) :: Mass_Density_AD
    INTEGER, OPTIONAL, INTENT(IN)     :: Molecule_ID
    ! Local variables
    INTEGER :: Id
    REAL(fp) :: c

    ! Setup
    IF ( Temperature < ZERO .OR. Mass_Density < ZERO ) THEN
      Partial_Pressure_AD = ZERO
      Temperature_AD      = ZERO
      Mass_Density_AD     = ZERO
      RETURN
    END IF
    IF ( PRESENT( Molecule_ID ) ) THEN
      IF ( Molecule_ID < 1 .OR. Molecule_ID > MAX_N_MOLECULAR_SPECIES ) THEN
        Partial_Pressure_AD = ZERO
        Temperature_AD      = ZERO
        Mass_Density_AD     = ZERO
        RETURN
      END IF
      Id = Molecule_ID
    ELSE
      Id = 1  ! Default value is for water vapor
    END IF

    ! Adjoint form of mass density to partial pressure conversion
    c = PA_TO_HPA * R0 / MOLECULAR_WEIGHT(Id)
    Mass_Density_AD = Mass_Density_AD + (c * Temperature  * Partial_Pressure_AD)
    Temperature_AD  = Temperature_AD  + (c * Mass_Density * Partial_Pressure_AD)
    Partial_Pressure_AD = ZERO

  END SUBROUTINE MD_to_PP_AD

END MODULE PP_MD

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
module persistent_object_m
    use attribute_value_pair_m
    use sqlite

    implicit none
    private

    type, public :: persistent_object_t
        private
            type(SQLITE_DATABASE), public :: db_m
            character(len=80) :: db_name_m
            integer, public :: id_m

        contains
            procedure, public, pass(this) :: set_db_name, &
                                             open_database, &
                                             close_database, &
                                             create_new, &
                                             update_existing, &
                                             delete_existing, &
                                             delete_by_id, &
                                             map_from_data

    end type persistent_object_t

contains
    subroutine set_db_name(this, db_name)
        class(persistent_object_t), intent(inout) :: this
        character(len=*), intent(in) :: db_name

        this%db_name_m = db_name
    end subroutine set_db_name

    subroutine open_database(this)
        class(persistent_object_t), intent(inout) :: this

        call sqlite3_open(this%db_name_m, this%db_m)
    end subroutine open_database

    subroutine close_database(this)
        class(persistent_object_t), intent(inout) :: this

        call sqlite3_close(this%db_m)
    end subroutine close_database

    integer function create_new(this)
        class(persistent_object_t), intent(inout) :: this
    end function create_new

    subroutine update_existing(this)
        class(persistent_object_t), intent(inout) :: this
    end subroutine update_existing

    subroutine delete_existing(this)
        class(persistent_object_t), intent(inout) :: this
    end subroutine delete_existing

    subroutine delete_by_id(this, id)
        class(persistent_object_t), intent(inout) :: this
        integer, intent(in) :: id

        character(len=80) :: query

        write (query, '(a,i5,a)') 'delete from student where id=', id, ';'

        call this%set_db_name('../cgi-data/students.db')
        call this%open_database()
        call sqlite3_do(this%db_m, query)
        call this%close_database()
    end subroutine delete_by_id

    subroutine map_from_data(this, the_data)
        class(persistent_object_t), intent(inout) :: this
        class(attribute_value_pair_t), dimension(:), intent(in) :: the_data

    end subroutine map_from_data

end module persistent_object_m

max_stars_repo_name: bceverly/FABS

max_issues_repo_name: bceverly/FABS

max_forks_repo_name: bceverly/FABS


subroutine remeshing(nrefl)
common/ray_param/ds_ini,ds_segm_min,bend_min0,bend_max0
common/ray/ nodes,xray(1000),yray(1000)
real xtmp(1000),ytmp(1000)

nmic0=200
dsmic0=ds_ini/nmic0

xtmp(1)=xray(1)
ytmp(1)=yray(1)

scur=0
ntmp=1
nrefl_new=0
do inode=1,nodes-1
	x1=xray(inode)
	x2=xray(inode+1)
	y1=yray(inode)
	y2=yray(inode+1)
	ds=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
	nmic=int_best(ds/dsmic0)
	dsmic=ds/nmic
	dss=0
	do imic=1,nmic
		scur=scur+dsmic
		if(inode+1.eq.nrefl) then
			dss=dss+dsmic
			if(ds-dss.lt.ds_ini/2.) then
				ntmp=ntmp+1
				xtmp(ntmp)=x2
				ytmp(ntmp)=y2
				nrefl_new=ntmp
				yrefl=yyy_surf(x2)
				!write(*,*)' reflected:',ntmp,x2,y2,yrefl

				exit
			end if
		end if
		!write(*,*)imic,scur
		if(scur.ge.ds_ini) then
			ntmp=ntmp+1
			sss=dsmic*imic
			xtmp(ntmp)=x1+((x2-x1)/(ds))*sss
			ytmp(ntmp)=y1+((y2-y1)/(ds))*sss
			!write(*,*)' ntmp=',ntmp,' x=',xtmp(ntmp),' y=',ytmp(ntmp)
			scur=0
		end if
	end do

end do

ntmp=ntmp+1
xtmp(ntmp)=xray(nodes)
ytmp(ntmp)=yray(nodes)
!write(*,*)' ntmp=',ntmp

!open(23,file='ray_tmp.dat')
!do i=1,ntmp
	!write(23,*)xtmp(i),ytmp(i)
!end do
!close(23)

!write(*,*)' nodes=',nodes,' ntmp=',ntmp

nodes=ntmp
xray=xtmp
yray=ytmp
!write(*,*)' nodes=',nodes,' ntmp=',ntmp,' nrefl_new=',nrefl_new
nrefl=nrefl_new

return
end


text:
{ -- NO$GBA logging functions
  --
  -- Copyright (C) 2015 Ulrich Hecht
  --
  -- This software is provided 'as-is', without any express or implied
  -- warranty.  In no event will the authors be held liable for any damages
  -- arising from the use of this software.
  --
  -- Permission is granted to anyone to use this software for any purpose,
  -- including commercial applications, and to alter it and redistribute it
  -- freely, subject to the following restrictions:
  --
  -- 1. The origin of this software must not be misrepresented; you must not
  --    claim that you wrote the original software. If you use this software
  --    in a product, an acknowledgment in the product documentation would be
  --    appreciated but is not required.
  -- 2. Altered source versions must be plainly marked as such, and must not be
  --    misrepresented as being the original software.
  -- 3. This notice may not be removed or altered from any source distribution.
  --
  -- Ulrich Hecht
  -- ulrich.hecht@gmail.com }

\ send a string to the NO$GBA console
code log ( a -- )
	\ NO$GBA only allows inline debug strings, so
	\ we construct a suitable code snippet at runtime
	\ in EWRAM.
	\ XXX: This address should not be fixed!
	$203b000 a LITERAL

	$e1a0c00c v0 LITERAL	\ mov r12, r12
	a 0@ v0 str,
	$ea00003e v0 LITERAL	\ b   pc+256
	a 4 #( v0 str,
	$6464 v0 LITERAL	\ magic number
	a 8 #( v0 str,

	tos 1 (# v0 ldrb,	\ string length
	12 ## a v2 add,		\ start of string in EWRAM

	\ copy string to the inside of our code fragment
	l: loop
		tos 1 (# v1 ldrb,
		v2 1 (# v1 strb,
		1 ## v0 v0 s! sub,
		loop ne? b,

	\ zero-terminate
	0 ## v0 mov,
	v2 1 (# v0 strb,

	260 ## a v2 add,
	$e12fff1e v0 LITERAL	\ ret
	v2 0@ v0 str,

	\ call our code fragment
	lr push
	pc lr mov,
	a bx,
	lr pop

	tos pop
	ret
end-code

max_stars_repo_name: uli/dragonbasic

max_issues_repo_name: funkygallo/dragonbasic

max_forks_repo_name: funkygallo/dragonbasic


text:
      subroutine ainvg (res, adda, neq, t, y, ydot, miter,
     1                   ml, mu, pw, ipvt, ier )
clll. optimize
      external res, adda
      integer neq, miter, ml, mu, ipvt, ier
      integer i, lenpw, mlp1, nrowpw
      double precision t, y, ydot, pw
      dimension y(1), ydot(1), pw(1), ipvt(1)
c-----------------------------------------------------------------------
c this subroutine computes the initial value
c of the vector ydot satisfying
c     a * ydot = g(t,y)
c when a is nonsingular.  it is called by lsodi for
c initialization only, when istate = 0 .
c ainvg returns an error flag ier..
c   ier  =  0  means ainvg was successful.
c   ier .ge. 2 means res returned an error flag ires = ier.
c   ier .lt. 0 means the a-matrix was found to be singular.
c-----------------------------------------------------------------------
c
      if (miter .ge. 4)  go to 100
c
c full matrix case -----------------------------------------------------
c
      lenpw = neq*neq
      do 10  i = 1, lenpw
   10    pw(i) = 0.0d0
c
      ier = 1
      call res ( neq, t, y, pw, ydot, ier )
      if (ier .gt. 1) return
c
      call adda ( neq, t, y, 0, 0, pw, neq )
      call dgefa ( pw, neq, neq, ipvt, ier )
      if (ier .eq. 0) go to 20
         ier = -ier
         return
   20 call dgesl ( pw, neq, neq, ipvt, ydot, 0 )
      return
c
c band matrix case -----------------------------------------------------
c
  100 continue
      nrowpw = 2*ml + mu + 1
      lenpw = neq * nrowpw
      do 110  i = 1, lenpw
  110    pw(i) = 0.0d0
c
      ier = 1
      call res ( neq, t, y, pw, ydot, ier )
      if (ier .gt. 1) return
c
      mlp1 = ml + 1
      call adda ( neq, t, y, ml, mu, pw(mlp1), nrowpw )
      call dgbfa ( pw, nrowpw, neq, ml, mu, ipvt, ier )
      if (ier .eq. 0) go to 120
         ier = -ier
         return
  120 call dgbsl ( pw, nrowpw, neq, ml, mu, ipvt, ydot, 0 )
      return
c-------------------- end of subroutine ainvg --------------------------
      end

max_stars_repo_name: lesserwhirls/scipy-cwt

max_issues_repo_name: lesserwhirls/scipy-cwt

max_forks_repo_name: lesserwhirls/scipy-cwt


text:
!
! The Laboratory of Algorithms
!
! The MIT License
!
! Copyright 2011-2015 Andrey Pudov.
!
! 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.
!
! 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.
!

module MArrays

    implicit none
    private

    type, public :: TArrays
    contains
        procedure, nopass :: fillWithSequence1d
        procedure, nopass :: fillWithSequence2d

        procedure, nopass :: fillWithDirtySequence
        procedure, nopass :: fillWithRandom
        procedure, nopass :: fillWithInversedSequence

        procedure, nopass :: increaseWhenRequiredAllocatable
        procedure, nopass :: increaseWhenRequiredAllocatable2d
        procedure, nopass :: increaseWhenRequiredPointer

        procedure, nopass :: resizeAllocatable
        procedure, nopass :: resizeAllocatable2d

        procedure, nopass :: isSorted

        procedure, nopass :: print1d
        procedure, nopass :: print2d

        procedure, nopass :: swap1d

        generic :: fillWithSequence     => fillWithSequence1d, fillWithSequence2d
        generic :: increaseWhenRequired => increaseWhenRequiredAllocatable,   &
                                           increaseWhenRequiredAllocatable2d, &
                                           increaseWhenRequiredPointer

        generic :: resize => resizeAllocatable, resizeAllocatable2d

        generic :: print => print1d, print2d
        generic :: swap  => swap1d
    end type
contains
    !
    ! Fills the array with sequence from one to the last index.
    !
    subroutine fillWithSequence1d(array)
        integer, dimension(:), intent(in out) :: array
        integer index

        array = (/ (index, index = 1, size(array)) /)
    end subroutine

    subroutine fillWithSequence2d(array)
        integer, dimension(:,:), intent(in out) :: array
        integer index

        array = reshape((/ (index, index = 1, size(array)) /), shape(array))
    end subroutine

    !
    ! Fills the array with sequence from one to the last index
    ! and ten percent ramdom values.
    !
    subroutine fillWithDirtySequence(array)
        integer, dimension(:), intent(in out) :: array

        integer index
        integer step

        array = (/ (index, index = 1, size(array)) /)
        step  = size(array) / (size(array) * 0.1)

        do index = step, size(array), step
            array(index) =  array(index) * 2
        end do
    end subroutine

    !
    ! Fills the array with random numbers between zero and maximal
    ! value of type of the array.
    !
    subroutine fillWithRandom(array)
        integer, dimension(:), intent(in out) :: array
        real, dimension(size(array)) :: temporary

        call random_seed()
        call random_number(temporary)

        array = floor(temporary * huge(array))
    end subroutine

    !
    ! Fills the array with sequence from last index of array to the one.
    !
    subroutine fillWithInversedSequence(array)
        integer, dimension(:), intent(in out) :: array
        integer index

        array = (/ (index, index = size(array), 1, -1) /)
    end subroutine

    !
    ! Increase the size of given array if limit is larger than size of array.
    !
    subroutine increaseWhenRequiredAllocatable(array, limit)
        integer, dimension(:), allocatable, intent(in out) :: array
        integer, intent(in) :: limit

        integer, dimension(:), allocatable :: temporary_array
        integer length

        ! increase array size when required
        if (limit > size(array)) then
            length = size(array) * 3 / 2
            allocate(temporary_array(length))
            temporary_array = 0
            temporary_array(1:size(array)) = array
            deallocate(array)
            call move_alloc(temporary_array, array)
        end if
    end subroutine

    subroutine increaseWhenRequiredAllocatable2d(array, limit1, limit2)
        integer, dimension(:,:), allocatable, intent(in out) :: array
        integer, intent(in) :: limit1
        integer, intent(in) :: limit2

        integer, dimension(:,:), allocatable :: temporary_array
        integer length1
        integer length2

        ! increase array size when required
        if ((limit1 > size(array, 1)) .or. (limit2 > size(array, 2))) then
            !print *, limit1, limit2
            length1 = size(array, 1) * 3 / 2
            length2 = size(array, 2) * 3 / 2

            allocate(temporary_array(length1, length2))
            temporary_array = 0
            temporary_array(1:size(array, 1), 1:size(array, 2)) = array
            deallocate(array)
            call move_alloc(temporary_array, array)
        end if
    end subroutine

    subroutine increaseWhenRequiredPointer(array, limit)
        integer, dimension(:), pointer, intent(in out) :: array
        integer, intent(in) :: limit

        integer, dimension(:), pointer :: temporary_array
        integer length

        ! increase array size when required
        if (limit > size(array)) then
            length = size(array) * 3 / 2
            allocate(temporary_array(length))
            temporary_array = 0
            temporary_array(1:size(array)) = array
            deallocate(array)
            array => temporary_array
        end if
    end subroutine

    !
    ! Reallocate array to a given length.
    !
    subroutine resizeAllocatable(array, length)
        integer, dimension(:), allocatable, intent(in out) :: array
        integer, intent(in) :: length

        integer, dimension(:), allocatable :: temporary_array

        if (length /= size(array)) then
            allocate(temporary_array(length))
            temporary_array = 0
            temporary_array(1:length) = array(1:length)
            deallocate(array)
            call move_alloc(temporary_array, array)
        end if
    end subroutine

    subroutine resizeAllocatable2d(array, length1, length2)
        integer, dimension(:,:), allocatable, intent(in out) :: array
        integer, intent(in) :: length1
        integer, intent(in) :: length2

        integer, dimension(:,:), allocatable :: temporary_array

        if ((length1 /= size(array, 1)) .or. (length2 /= size(array, 2))) then
            allocate(temporary_array(length1, length2))
            temporary_array = 0
            temporary_array(1:length1, 1:length2) = array(1:length1, 1:length2)
            deallocate(array)
            call move_alloc(temporary_array, array)
        end if
    end subroutine

    !
    ! Validates given array to be sorted.
    !
    function isSorted(array) result(sorted)
        integer, dimension(:), intent(in) :: array
        logical                           :: sorted

        integer index

        do index = 2, size(array)
            if (array(index - 1) > array(index)) then
                sorted =  .false.
                return
            end if
         end do

         sorted = .true.
     end function

     !
     ! Prints given array to console.
     !
     subroutine print1d(array)
         integer, dimension(:), intent(in) :: array

         character(len=80) :: count
         character(len=80) :: length
         integer           :: maximum
         integer           :: value

         maximum = (maxval(array))
         if (maximum == 0) then
             value = 2
         else if ((mod(maximum, 10) == 0) .or. (maximum == 1)) then
             value = ceiling(log10(real(maximum))) + 2
         else
             value = ceiling(log10(real(maximum))) + 1
         end if

         write(count,  '(i0.0)') size(array)
         write(length, '(i0.0)') value

         print '(' // trim(count) // 'i' // trim(length) // ')', array
     end subroutine

     subroutine print2d(array)
         integer, dimension(:,:), intent(in) :: array

         character(len=80) :: count
         character(len=80) :: length
         integer           :: maximum
         integer           :: value

         maximum = (maxval(array))
         if (maximum == 0) then
             value = 2
         else if ((mod(maximum, 10) == 0) .or. (maximum == 1)) then
             value = ceiling(log10(real(maximum))) + 2
         else
             value = ceiling(log10(real(maximum))) + 1
         end if

         write(count,  '(i0.0)') size(array, 2)
         write(length, '(i0.0)') value

         print '(' // trim(count) // 'i' // trim(length) // ')', transpose(array)
     end subroutine

     subroutine swap1d(array, index, jndex)
         integer, dimension(:), intent(in out) :: array
         integer, intent(in) :: index
         integer, intent(in) :: jndex

         integer buffer

         buffer = array(index)
         array(index) = array(jndex)
         array(jndex) = buffer
     end subroutine
end module

max_stars_repo_name: andreypudov/Algorithms

max_issues_repo_name: jlokimlin/Algorithms

max_forks_repo_name: jlokimlin/Algorithms


text:
MODULE FUN_xyEXPAND
!!#### PURPOSE
!! Calculate the expansion of R^2 polygon into R^3 polygon or
!! dalculate the expansion of R^2 point to R^3 point.

!!#### EXTERNAL KINDS
USE KND_IntrinsicTypes,ONLY: KIND_Rsp,KIND_Rdp !!((01-A-KND_IntrinsicTypes.f90))

!!#### DEFAULT IMPLICIT
IMPLICIT NONE

!!#### DEFAULT ACCESS
PRIVATE

!!#### PROCEDURE OVERLOADING
INTERFACE xyEXPAND_Pg
 MODULE PROCEDURE xyEXPAND_Pg_Rsp
 MODULE PROCEDURE xyEXPAND_Pg_Rdp
END INTERFACE

INTERFACE xyEXPAND_P
 MODULE PROCEDURE xyEXPAND_P_Rsp
 MODULE PROCEDURE xyEXPAND_P_Rdp
END INTERFACE

!!#### PUBLIC ACCESS LIST
PUBLIC :: xyEXPAND_Pg
PUBLIC :: xyEXPAND_P

CONTAINS

PURE FUNCTION xyEXPAND_Pg_Rsp( N , Pg , M32_transform ) RESULT(xyzPg)
USE KND_IntrinsicTypes,ONLY: KIND_R=>KIND_Rsp  !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "05-B-FUN_xyEXPAND_Pg.f90.hdr"
!!--begin--
INCLUDE "05-B-FUN_xyEXPAND_Pg.f90.bdy"
!!--end--
END FUNCTION

PURE FUNCTION xyEXPAND_Pg_Rdp( N , Pg , M32_transform ) RESULT(xyzPg)
USE KND_IntrinsicTypes,ONLY: KIND_R=>KIND_Rdp  !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "05-B-FUN_xyEXPAND_Pg.f90.hdr"
!!--begin--
INCLUDE "05-B-FUN_xyEXPAND_Pg.f90.bdy"
!!--end--
END FUNCTION

PURE FUNCTION xyEXPAND_P_Rsp( P , M32_transform ) RESULT(xyzP)
USE KND_IntrinsicTypes,ONLY: KIND_R=>KIND_Rsp  !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "05-B-FUN_xyEXPAND_P.f90.hdr"
!!--begin--
INCLUDE "05-B-FUN_xyEXPAND_P.f90.bdy"
!!--end--
END FUNCTION

PURE FUNCTION xyEXPAND_P_Rdp( P , M32_transform ) RESULT(xyzP)
USE KND_IntrinsicTypes,ONLY: KIND_R=>KIND_Rdp  !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "05-B-FUN_xyEXPAND_P.f90.hdr"
!!--begin--
INCLUDE "05-B-FUN_xyEXPAND_P.f90.bdy"
!!--end--
END FUNCTION

END MODULE

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


text:
c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
c Copyright (C) INRIA
c
c This file must be used under the terms of the CeCILL.
c This source file is licensed as described in the file COPYING, which
c you should have received as part of this distribution.  The terms
c are also available at
c http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
c
      subroutine zqnbd(indqn,simul,dh,n,binf,bsup,x,f,g,zero,napmax,
     &itmax,indic,izig,nfac,imp,io,epsx,epsf,epsg,x1,x2,g1,dir,df0,
     &ig,in,irel,izag,iact,epsrel,ieps1,izs,rzs,dzs)
c
      implicit double precision (a-h,o-z)
      real rzs(*)
      double precision dzs(*)
      dimension x1(n),x2(n),g1(n),dir(n),epsx(n)
      dimension binf(n),bsup(n),x(n),g(n),dh(*),indic(n),izig(n),
     &izs(*)
      external simul,proj
c
      if(imp.lt.4)go to 3
      write(io,1020)izag,ig,in,irel,iact,epsrel
1020  format(' qnbd :  izag,ig,in,irel,iact,epsrel=',5i3,f11.4)
c
      if(ig.eq.1)write(io,110)
110   format(' test sur gradient pour sortie ib')
      if(in.eq.1)write(io,111)
111   format(' test sur nombre de defactorisations pour sortie ib')
      if(izag.ne.0)write(io,112)izag
112   format(' memorisation de variables izag=',i3)
      if(irel.eq.1)write(io,114)epsrel
114   format(' methode de minimisations incompletes ; epsrel=',d11.4)
      if(iact.eq.1)write(io,116)
116   format(' blocage des variables dans ib')
      if(ieps1.eq.1)write(io,118)
118   format(' parametre eps1 nul')
      if(ieps1.eq.2)write(io,119)
119   format(' parametre eps1 grand')
c
c     cscal1 utilise pour calculer eps(x) = eps1 cf avant 310
      cscal1=1.0d+8
      if(ieps1.eq.2)write(io,120)cscal1
120   format(' parametre eps1=eps(x) calcule avec cscal1=',d11.4)
3     continue
c
      difg0=1.0d+0
      difg1=0.0d+0
c
c     eps0 sert a partitionner les variables
      eps0=0.0d+0
      do 5 i=1,n
      izig(i)=0
5     eps0=eps0+epsx(i)
      eps0=10.*eps0/n
c
c     section 1  mise en forme de dh
c     si indqn=1 on init dh a ident puis scal a it=2
c
      call proj(n,binf,bsup,x)
      ndh=n*(n+1)/2
      if(indqn.eq.1)go to 10
      if(indqn.eq.2)go to 30
c     erreur
      if(imp.gt.0)write(io,105)indqn
105   format(' qnbd  : valeur non admissible de indqn  ',i5)
      indqn=-105
      if(imp.gt.0)write(io,123)indqn
      return
10    continue
c     on initialise dh a l identite puis a l iteration 2
c     on met a l echelle
      nfac=0
      do 11 i=1,n
11    indic(i)=i
      do 12 i=1,ndh
12    dh(i)=0.0d+0
30    continue
c
c     section 2  mise a jour dh
c
c     iter nombre d iterations de descente
      iter=0
      scal=1.0d+0
      nap=1
      indsim=4
      if(indqn.eq.1) call simul(indsim,n,x,f,g,izs,rzs,dzs)
      if(indsim.le.0)then
      indqn=-1
      if(indsim.eq.0)indqn=0
      if(imp.gt.0)write(io,123)indqn
123   format(' qnbd : indqn=',i8)
      return
      endif
      if(indqn.ne.1)go to 200
c     mise a echelle dh
c     df0 decroissance prevue . si mod quad df0=((dh)-1g,g)/2
c     et on cherche dh diag de la forme cst/(dx)**2
c     d ou cst=som((y(i)*(dx))**2))/(2*df0)
      cof1=0.0d+0
      do 80 i=1,n
80    cof1=cof1+(g(i)*epsx(i))**2
      cof1=cof1/(2.0d+0*df0)
      i1=-n
      do 82 i=1,n
      i1=i1+n+2-i
82    dh(i1)=(cof1 + zero)/(epsx(i)**2 + zero)
      iconv=0
200   iter=iter +1
      if(iter.le.itmax)go to 202
      if(imp.gt.0)write(io,1202)
1202  format(' qnbd : maximum d iterations atteint')
      indqn=5
      if(imp.gt.0)write(io,123)indqn
      return
202   if(imp.ge.2)write(io,1210)iter,f
1210  format(/' qnbd : iter=',i3,'  f=',d15.7)
c     x1,g1 valeurs a l iteration precedente
      if(iter.eq.1)go to 300
      cof1=0.0d+0
      do 201 i=1,n
      x1(i)=x(i)-x1(i)
      g1(i)=g(i)-g1(i)
201   cof1=cof1 + x1(i)*g1(i)
      if(cof1.le.zero)go to 250
      if(iter.gt.2.or.indqn.ne.1)go to 250
c     mise a l echelle de dh par methode shanno-phua
c      dh=(y,y)/(y,s)*id
      cof2=0.0d+0
      do 203 i=1,n
203   cof2=cof2 + g1(i)**2
      cof2=cof2/cof1
      if(imp.gt.3)write(io,1203)cof2
1203  format(' qnbd : facteur d echelle=',d11.4)
      dh(1)=cof2
      i1=1
      do 205 i=1,nfac
      i1=i1+n+1-i
205   dh(i1)=cof2
c
c     scal= (y,s)/(y,y)
c     scal sert de coeff a g dans le calcul de dir pour i dans ib
      scal=1.0d+0/cof2
250   continue
c
c     mise a jour dh par methode bfgs (majour) si iter ge 2
c     dh1=dh +y*yt/(y,s) - dh*s*st*dh/(s,dh*s)
c     exprimons ds=x1 et y=g1 dans les nouv variables soit x2 et g1
      do 251 i=1,n
      i1=indic(i)
      x2(i1)=g1(i)
251   dir(i1)=x1(i)
      do 252 i=1,n
252   g1(i)=x2(i)
      do 253 i=1,n
      i1=indic(i)
253   x2(i1)=x1(i)
c     on stocke d abord dh*s dans x2
c     calcul des nfac premieres variables,en deux fois
      continue
      if(nfac.eq.0) go to 2312
      if(nfac.gt.1) go to 2300
      dir(1)=dir(1)*dh(1)
      go to 2312
2300  continue
      np=nfac+1
      ii=1
      n1=nfac-1
      do 2303 i=1,n1
      y=dir(i)
      if(dh(ii).eq.0.0d+0) go to 2302
      ij=ii
      ip=i+1
      do 2301 j=ip,nfac
      ij=ij+1
2301  y=y+dir(j)*dh(ij)
2302  dir(i)=y*dh(ii)
2303  ii=ii+np-i
      dir(nfac)=dir(nfac)*dh(ii)
      do 2311 k=1,n1
      i=nfac-k
      ii=ii-np+i
      if(dir(i).eq.0.0d+0) go to 2311
      ip=i+1
      ij=ii
      y=dir(i)
      do 2310 j=ip,nfac
      ij=ij+1
2310  dir(j)=dir(j)+dh(ij)*dir(i)
2311  continue
2312  continue
      nfac1=nfac+1
      n2fac=(nfac*nfac1)/2
      nnfac=n-nfac
      k=n2fac
      if(nfac.eq.n)go to 268
      do 255 i=nfac1,n
255   dir(i)=0.0d+0
      if(nfac.eq.0)go to 265
      do 260i=1,nfac
      do 260j=nfac1,n
      k=k+1
      if(x2(j).eq.0.)go to 260
      dir(i)= dir(i) + dh(k)*x2(j)
260   continue
c     calcul autres comp de dh*s=d en deux fois
      k=n2fac
      do 264 j=1,nfac
      do 264 i=nfac1,n
      k=k+1
      dir(i)=dir(i) + dh(k)*x2(j)
264   continue
265   continue
      k=n2fac+nfac*nnfac
      do 266 j=nfac1,n
      do 266 i=j,n
      k=k+1
      if(x2(j).eq.0.)go to 266
      dir(i)=dir(i) + dh(k)*x2(j)
266   continue
      if(nfac.eq.n-1)go to 268
      nm1=n-1
      k=n2fac+nfac*nnfac
      do 267 i=nfac1,nm1
      k=k+1
      i1=i+1
      do 267 j=i1,n
      k=k+1
      if(x2(j).eq.0.)go to 267
      dir(i)=dir(i)+dh(k)*x2(j)
267   continue
c     calcul de dh*s fini
c     calcul sig1 pour 2eme mise a jour
268   sig1=0.0d+0
      do 271 i=1,n
271   sig1=sig1+dir(i)*x2(i)
      if(sig1.gt.0.0d+0)go to 272
      if(imp.gt.2)write(io,1272)sig1
1272  format(' qnbd : pb (bs,s) negatif=',d11.4)
c
c     ******************************************************
      indqn=8
      if(iter.eq.1)indqn=-5
      if(imp.gt.0)write(io,123)indqn
      return
272      sig1=-1.0d+0/sig1
c     truc powell si (y,s) negatif
      if(cof1.gt.zero)go to 277
      if(imp.gt.2)write(io,1270)cof1
1270  format(' qnbd : emploi truc powell (y,s)=',d11.4)
      teta=-1.0d+0/sig1
      teta=.8*teta/(teta-cof1)
      teta1=1.0d+0-teta
      do 274 i=1,n
274   g1(i)=teta*g1(i)+teta1*dir(i)
      cof1=-.2/sig1
277   continue
c
c      premiere mise a jour de dh
      sig=1.0d+0/cof1
      zsig1=1.0d+0/sig1
      mk=0
      ir=nfac
      epsmc=1.0d-9
      call calmaj(dh,n,g1,sig,x2,ir,mk,epsmc,nfac)
      if(ir.ne.nfac)go to 280
      call calmaj(dh,n,dir,sig1,x2,ir,mk,epsmc,nfac)
      if(ir.ne.nfac)go to 280
      go to 300
280   if(imp.gt.0)write(io,282)
282   format(' qnbd : pb dans appel majour')
      indqn=8
      if(iter.eq.1)indqn=-5
      if(imp.gt.0)write(io,123)indqn
      return
300   continue
c
c     section 3 determination des variables libres et bloquees
c
c     calcul eps1
c
      scal1=scal
      if(ieps1.eq.1)scal1=0.0d+0
      if(ieps1.eq.2)scal1=scal*cscal1
305   do 310 i=1,n
310   x1(i)=x(i)-scal1*abs(g(i))*g(i)
      call proj(n,binf,bsup,x1)
      eps1=0.0d+0
      do 320 i=1,n
320   eps1=eps1 + abs(x1(i)-x(i))
      eps1=min(eps0,eps1)
      if(ieps1.eq.1)eps1=0.0d+0
      if(ieps1.eq.2)eps1=eps1*1.0d+4
      if(imp.gt.3)write(io,322)eps1
322   format(' qnbd : val de eps1 servant a partitionner les variables'
     &,d11.4)
c     nfac nombre de lignes factorisees (nr pour ajour)
      ifac=0
      idfac=0
      k=0
c
c
      gr=0.0d+0
      if(ig.eq.1)gr=0.2*difg/n
      n3=n
      if(in.eq.1)n3=n/10
c     si irit=1 on peut relacher des variables
      irit=0
      if(difg1.le.epsrel*difg0)irit=1
      if(irel.eq.0.or.iter.eq.1)irit=1
      if(irit*irel.gt.0.and.imp.gt.3)write(io,1320)difg0,epsrel,difg1
1320  format(' qnbd : redemarrage ; difg0,epsrel,difg1=',3d11.4)
c
      tiers=1.0d+0/3.0d+0
      do 340 k=1,n
      izig(k)=izig(k)-1
      if(izig(k).le.0)izig(k)=0
      bi=binf(k)
      bs=bsup(k)
      ic=indic(k)
      d1=x(k)-bi
      d2=bs-x(k)
      dd=(bs-bi)*tiers
      ep=min(eps1,dd)
      if(d1.gt.ep)go to 324
      if(g(k).gt.0.)go to 330
      go to 335
324   if(d2.gt.ep)go to 335
      if(g(k).gt.0.)go to 335
      go to 330
c     on defactorise si necessaire
330   continue
      if(ic.gt.nfac)go to 340
      idfac=idfac+1
      mode=-1
      if(imp.ge.4)write(io,336)k
336   format(' defactorisation de ',i3)
      izig(k)=izig(k) + izag
      call ajour(mode,n,k,nfac,dh,x2,indic)
      if(mode.eq.0) go to 340
      if(imp.gt.0)write(io,333)mode
333   format(' qnbd : pb dans ajour. mode=',i3)
      indqn=8
      if(iter.eq.1)indqn=-5
      if(imp.gt.0)write(io,123)indqn
      return
c     on factorise
335   continue
      if(irit.eq.0)go to 340
      if(ic.le.nfac)go to 340
      if(izig(k).ge.1)go to 340
      mode=1
      if(ifac.ge.n3.and.iter.gt.1)go to 340
      if(abs(g(k)).le.gr)go to 340
      ifac=ifac+1
      if(imp.ge.4)write(io,339)k
339   format(' on factorise l indice ',i3)
      call ajour(mode,n,k,nfac,dh,x2,indic)
      if(mode.eq.0)go to 340
      if(imp.gt.0) write(io,333)mode
      indqn=8
      if(iter.eq.1)indqn=-5
      if(imp.gt.0)write(io,123)indqn
      return
340   continue
      if(imp.ge.2)write(io,350)ifac,idfac,nfac
350   format(' qnbd : nbre fact',i3,' defact',i3,
     &' total var factorisees',i3)
c
c     *********************************************** a voir
      if(iconv.eq.1)return
c
c     section 6 resolution systeme lineaire et expression de dir
c     on inverse le syst correspondant aux nl premieres composantes
c     dans le nouveau syst d indices
      ir=nfac
      do 640 i=1,n
      i1=indic(i)
640   x2(i1)=g(i)
641   continue
      if(ir.lt.nfac) go to 412
      if(nfac.gt.1) go to 400
      x2(1)=x2(1)/dh(1)
      go to 412
400   continue
      do 402 i=2,nfac
      ij=i
      i1=i-1
      v=x2(i)
      do 401 j=1,i1
      v=v-dh(ij)*x2(j)
401   ij=ij+nfac-j
      x2(i)=v
402   x2(i)=v
      x2(nfac)=x2(nfac)/dh(ij)
      np=nfac+1
      do 411 nip=2,nfac
      i=np-nip
      ii=ij-nip
      v=x2(i)/dh(ii)
      ip=i+1
      ij=ii
      do 410 j=ip,nfac
      ii=ii+1
410   v=v-dh(ii)*x2(j)
411   x2(i)=v
412   continue
      if(ir.eq.nfac)go to 660
      if(imp.gt.0)write(io,650)
650   format(' qnbd : pb num dans mult par inverse')
      indqn=7
      if(iter.eq.1)indqn=-6
      if(imp.gt.0)write(io,123)indqn
      return
660   continue
      do 610 i=1,n
      i1=indic(i)
      dir(i)=-g(i)*scal
610   if(i1.le.nfac) dir(i)=-x2(i1)
      continue
c
c     gestion contraintes actives (si iact=1)
      if(iact.ne.1)go to 675
      do 670 i=1,n
      if(izig(i).gt.0)dir(i)=0.
      if(indic(i).gt.nfac)dir(i)=0.0d+0
670   continue
675   continue
c
c     recherche lineaire
c     conservation de x et g . calcul de dir+ et fpn
      do 700 i=1,n
      g1(i)=g(i)
700   x1(i)=x(i)
c     ifp =1 si fpn trop petit. on prend alors d=-g
      ifp=0
      fn=f
709   fpn=0.0d+0
      do 710 i=1,n
      if(x(i)-binf(i).le.epsx(i).and.dir(i).lt.0.)dir(i)=0.0d+0
      if(bsup(i)-x(i).le.epsx(i).and.dir(i).gt.0.)dir(i)=0.0d+0
710   fpn=fpn + g(i)*dir(i)
      if(fpn.gt.0.0d+0) then
         if(ifp.eq.1) then
            if(imp.gt.0)write(io,1705) fpn
1705        format(' qnbd : arret fpn non negatif=',d11.4)
            indqn=6
            if(iter.eq.1)indqn=-3
            if(imp.gt.0)write(io,123)indqn
            return
         else
            ifp=1
            do 707 i=1,n
707         if(izig(i).gt.0)dir(i)=-scal*g(i)
            irit=1
            go to 709
         endif
      endif
c     calcul du t initial suivant une idee de fletcher
      t1=t
      if(iter.eq.1) diff=df0
      t=-2.0d+0*diff/fpn
      if(t.gt.0.30d+0.and.t.lt.3.0d+0)t=1.0d+0
      if(eps1.lt.eps0)t=1.0d+0
      if(indqn.eq.2)t=1.0d+0
      if(iter.gt.1.and.t1.gt.0.010d+0.and.t1.lt.100.0d+0)t=1.0d+0
      tmax=1.0d+10
      t=min(t,tmax)
      t=max(t,10.*zero)
c     amd,amf tests sur h'(t) et diff
      amd=.7
      amf=.1
      napm=15
      napm1=nap + napm
      if(napm1.gt.napmax)napm1=napmax
      call rlbd(indrl,n,simul,x,binf,bsup,fn,fpn,t,tmax,dir,g,
     & tproj,amd,amf,imp,io,zero,nap,napm1,x2,izs,rzs,dzs)
      if(indrl.ge.10)then
         indsim=4
         nap=nap + 1
         call simul(indsim,n,x,f,g,izs,rzs,dzs)
         if(indsim.le.0)then
            indqn=-3
            if(indsim.eq.0)indqn=0
            if(imp.gt.0)write(io,123)indqn
            return
         endif
      endif
      if(indrl.le.0)then
         indqn=10
         if(indrl.eq.0)indqn=0
         if(indrl.eq.-3)indqn=13
         if(indrl.eq.-4)indqn=12
         if(indrl.le.-1000)indqn=11
         if(imp.gt.0)write(io,123)indqn
         return
      endif
c
753   if(imp.lt.6)go to 778
      do 760 i=1,n
760      write(io,777)i,x(i),g(i),dir(i)
777   format(' i=',i2,' xgd ',3f11.4)
c
778   continue
      if(nap.lt.napmax)go to 758
      f=fn
      if(imp.gt.0)write(io,755)napmax
755   format(' qnbd : retour cause max appels simul',i9)
      indqn=4
      if(imp.gt.0)write(io,123)indqn
      return
758   continue
c     section 8 test de convergence
c
      do 805 i=1,n
      if(abs(x(i)-x1(i)).gt.epsx(i))go to 806
805   continue
      f=fn
      if(imp.gt.0)write(io,1805)
1805  format(' qnbd : retour apres convergence de x')
      indqn=3
      if(imp.gt.0)write(io,123)indqn
      return
806   continue
      difg=0.0d+0
      do 810 i=1,n
      aa=g(i)
      if(x(i)-binf(i).le.epsx(i))aa=min(0.0d+0,aa)
      if(bsup(i)-x(i).le.epsx(i))aa=max(0.0d+0,aa)
810   difg=difg + aa**2
      difg1=0.0d+0
      do 820 i=1,n
      if(indic(i).gt.nfac)go to 820
      aa=g(i)
      if(x(i)-binf(i).le.epsx(i))aa=min(0.0d+0,aa)
      if(bsup(i)-x(i).le.epsx(i))aa=max(0.0d+0,aa)
      difg1=difg1 + aa**2
820   continue
      difg1=sqrt(difg1)
      difg=sqrt(difg)
      difg=difg/sqrt(real(n))
      diff=abs(f-fn)
      df0=-diff
      if(irit.eq.1)difg0=difg1
      f=fn
      if(imp.ge.2)write(io,860)epsg,difg,epsf,diff,nap
860   format(' qnbd : epsg,difg=',2d11.4,'  epsf,diff=',2d11.4
     &,'  nap=',i3)
      if(diff.lt.epsf)then
      indqn=2
      if(imp.gt.0)write(io,1865)diff
1865  format(' qnbd : retour cause decroissance f trop petite=',d11.4)
      if(imp.gt.0)write(io,123)indqn
      return
      endif
      if(difg.gt.epsg)go to 200
      indqn=1
      if(imp.gt.0)write(io,1900)difg
1900  format(' qnbd : retour cause gradient projete petit=',d11.4)
      if(imp.gt.0)write(io,123)indqn
      return
      end

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text:
      SUBROUTINE COSFT(Y,N,ISIGN)
      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
      DIMENSION Y(N)
      THETA=3.14159265358979D0/DBLE(N)
      WR=1.0D0
      WI=0.0D0
      WPR=-2.0D0*DSIN(0.5D0*THETA)**2
      WPI=DSIN(THETA)
      SUM=Y(1)
      M=N/2
      DO 11 J=1,M-1
        WTEMP=WR
        WR=WR*WPR-WI*WPI+WR
        WI=WI*WPR+WTEMP*WPI+WI
        Y1=0.5*(Y(J+1)+Y(N-J+1))
        Y2=(Y(J+1)-Y(N-J+1))
        Y(J+1)=Y1-WI*Y2
        Y(N-J+1)=Y1+WI*Y2
        SUM=SUM+WR*Y2
11    CONTINUE
      CALL REALFT(Y,M,+1)
      Y(2)=SUM
      DO 12 J=4,N,2
        SUM=SUM+Y(J)
        Y(J)=SUM
12    CONTINUE
      IF (ISIGN.EQ.-1) THEN
        EVEN=Y(1)
        ODD=Y(2)
        DO 13 I=3,N-1,2
          EVEN=EVEN+Y(I)
          ODD=ODD+Y(I+1)
13      CONTINUE
        ENF0=2.0*(EVEN-ODD)
        SUMO=Y(1)-ENF0
        SUME=(2.0*ODD/FLOAT(N))-SUMO
        Y(1)=0.5*ENF0
        Y(2)=Y(2)-SUME
        DO 14 I=3,N-1,2
          Y(I)=Y(I)-SUMO
          Y(I+1)=Y(I+1)-SUME
14      CONTINUE
      ENDIF
      RETURN
      END

ext: for