! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! 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.

!>    @brief modify coefficents for the temperature equation
!>    @param[in] ismpl local sample index
!>    @details
!> modify coefficents for the temperature equation according to the prescribed sources and sinks\n
!> rhs stored in w.\n
      SUBROUTINE set_tq(ismpl)
        use arrays
        use mod_genrl
        use mod_time
        IMPLICIT NONE
        integer :: ismpl
        integer :: i, j, k

        DOUBLE PRECISION deltat, deltt, qt
        EXTERNAL deltat, qt
! rhs: sources
        IF (transient .AND. tr_switch(ismpl)) THEN
          deltt = deltat(simtime(ismpl),ismpl)
!$OMP     do schedule(static) collapse(3)
          DO k = 1, k0
            DO j = 1, j0
              DO i = 1, i0
                w(i,j,k,ismpl) = w(i,j,k,ismpl) - qt(i,j,k,ismpl)
              END DO
            END DO
          END DO
!$OMP     end do nowait
        ELSE
!$OMP     do schedule(static) collapse(3)
          DO k = 1, k0
            DO j = 1, j0
              DO i = 1, i0
                w(i,j,k,ismpl) = w(i,j,k,ismpl) - qt(i,j,k,ismpl)
              END DO
            END DO
          END DO
!$OMP     end do nowait
        END IF

! Heat sources associated with flow dirichlet nodes - - - - - - -

!      do ib=first_flow,last_flow
!          i=ibc_data(ib,cbc_i)
!          j=ibc_data(ib,cbc_j)
!          k=ibc_data(ib,cbc_k)
!          bctype=ibc_data(ib,cbc_bt)
!          if (bctype.eq.bt_diri) then
!             w(i,j,k)=w(i,j,k)-qheadbcd(i,j,k,ismpl) ??? Zeit-Abhaengikeit ???
!             write(99,'(3i6,4x,e15.5)') i,j,k,qheadbcd(i,j,k,ismpl)
!          endif
!      end do

        RETURN
      END


text:
c
c     Adaptation hpf de onde24, Fabien Coelho et Ronan Keryell, 1993
c
c     - taille NP passees en parame`tres et verifiee aux IO~;
c     - idem NT~;
c     - les surdeclarations inutiles ont ete corrigees (NPMAX+1->NP)
c     - tableaux B et UB coupes en 3~: E, S, W~;
c     - tableau temporaire UINT duplique en 3~: E, S, W~;
c     - reductions remplacees par des appels de fonctions~;
c     - inclusion dans le code de l ouverture des fichiers~;
c     - directives de placement proposees par Ronan, et modifiees~;
c     - boucles paralle`les independantes,
c     - avec les variables privees~;
c     - mise en commentaire des mesures de performances~;
c     - & -> $ car marche mieux avec emacs~;
c     - fonctions REDMIN et REDMAX codees en 2D~;
c     - CONTINUE -> ENDDO
c     - boucle 71 distribuee
c
c     + hint to the compiler with independent do scalar loops
c
c     Il faudrait encore
c
c     - revoir les I/O...

      PROGRAM hpftest48

C  ------------------------------------------------------
C        MODELISATION PAR EQUATION D ONDES D UN PHENOMENE
C              DE PROPAGATION EN DIMENSION 2.
C
C  HYPOTHESES  :   - ON CONSIDERE UN DOMAINE CARRE .
C  ----------      - LA SURFACE EST LIBRE.
C                  - ON CONSIDERE DES CONDITIONS ABSORBANTES DU SECOND
C                    ORDRE SUR LES BORDS  (cf. these CHALINDAR).
C                  - LE MILIEU EST HETEROGENE.
C
C
C  UTILISATION DUN SCHEMA AUX DIFFERENCES FINIES DORDRE (2,4)
C  (ORDRE 2 EN TEMPS ET 4 EN ESPACE) COMME METHODE DE RESOLUTION :
C
C         - SCHEMA DORDRE (2,4) POUR LE DOMAINE INTERIEUR.
C         - SCHEMA DORDRE (2,2) POUR LE BORD INTERIEUR DU DOMAINE.
C         - CONDITIONS PARAXIALES DORDRE 2 POUR LE BORD DU DOMAINE.
C
C  -----------------------------------------------------

C  DECLARATION DE CONSTANTES
C  -------------------------
C  NPMAX      : NOMBRE DE POINTS MAXIMUM DU MAILLAGE SPATIAL
C  NTMAX      : NOMBRE MAXIMUM DE PAS DE TEMPS
C  NBTRAC   : NOMBRE DE TRACES SISMIQUES
C  TMAX        : TEMPS DE PRISE EN COMPTE DE LA FONCTION SOURCE (S)

      INTEGER    NPMAX
      INTEGER    NTMAX
      INTEGER    NBTRAC
      INTEGER    TMAX
      REAL*8     PI

      PARAMETER      (NPMAX     = 809)
      PARAMETER      (NTMAX     = 4000)
      PARAMETER      (NBTRAC  = 50)
      PARAMETER      (TMAX       = 1.0)
      PARAMETER      (PI         = 3.141592653589)

C  DECLARATION DE TABLEAUX
C  -----------------------
C  U      :    SOLUTION DE LA MODELISATION PAR EQUATION DONDE DUN
C              PHENOMENE DE PROPAGATION. UTILISATION DUN SCHEMA AUX
C              DIFFERENCES FINIES COMME METHODE DE RESOLUTION.
C                U(.,.,KM) REPRESENTE U(t)
C              U(.,.,KP) REPRESENTE U(t-1) et U(t+1)
C
C  V      :    CHAMP DE VITESSE CARACTERISTIQUE DU MILIEU.
C              - LE TABLEAU MEMORISE LA CONSTANTE LIEE A LA VITESSE
C                     1/12*((V(i,j)*DELTAT/H)**2)
C                POUR LES POINTS INTERIEURS DU DOMAINE.
C              - LE TABLEAU MEMORISE LA CONSTANTE LIEE A LA VITESSE
C                     (V(i,j)*DELTAT/H)**2
C                POUR LES POINTS CONSTITUANT LE BORD INTERNE DU DOMAINE
C              - DANS LE CAS DES POINTS CONSTITUANT LES FRONTIERES
C                E-W-S, LE TABLEAU CONTIENT LA MOYENNE HARMONIQUE
C                DES VITESSES DES POINTS FRONTIERES ET DES POINTS
C                CONNEXES A LA FRONTIERE.
C              - DANS LE CAS DES POINTS CONSTITUANT LES COINS INF.
C                DU DOMAINE, ET LES DEUX POINTS FRONTIERES CONNEXES
C                AUX COINS,  LE TABLEAU CONTIENT LA MOYENNE HARMONIQUE
C                DES VITESSES DES 4 POINTS CONNEXES AU POINT CONCERNE.
C              - LES COINS SUPERIEURS SONT CONSIDERES COMME PARTIE
C                DE LA SURFACE LIBRE.
C
C  B      :    TABLEAUX DE 3 VECTEURS CONTENANT LES CONSTANTES LIEES
C              AUX EQUATIONS ASSOCIEES AUX BORDS.
C              B1 = 1 - (1/(V`*DELTAT + H)) * (2*H + (V`*DELTAT)**2/H)
C              B2 = 2*H / (V`*DELTAT + H)
C              B3 = (V`*DELTAT)*(V`*DELTAT) / 2*H*(V`*DELTAT + H)
C              AVEC V`: VITESSE HARMONIQUE SUR DEUX POINTS.
C
C              DANS LE CAS PARTICULIER DES 3 POINTS COMPOSANT
C              CHAQUE COIN E-S ET S-W, LES TABLEAUX B1 ET B2
C              CONTIENNENT RESPECTIVEMENT LES CONSTANTES :
C                  B1   = 2*V`*DELTAT / ( 4*V`*DELTAT + 3*H)
C                  B2   = 3*H / (4*V`*DELTAT + 3*H)
C              AVEC V`: VITESSE HARMONIQUE SUR QUATRE POINTS.
C
C  UB     :    TABLEAU DE VECTEURS INTERMEDIAIRES UTILISES POUR LE
C              CALCUL DES BORDS E-S-W ; LE CODAGE EST LE MEME QUE POUR
C              LES VECTEURS B.
C
C  UINT   :    VECTEUR INTERMEDIAIRE UTILISE LORS DU CALCUL DE
C              LA SOLUTION AUX BORDS
C
C

c
c     nouvelles declarations avec NP et NT
c
c     pour BE, BW, UINT, (NP-1) est suffisant
c
      INTEGER   NP
      INTEGER   NT
      PARAMETER (NP = 800)
      PARAMETER (NT = 4000)

      external redmin, redmax
      real*8 redmin, redmax

      REAL*8    U     (NP,NP,2)
      REAL*8    V     (NP,NP)

      REAL*8    BE   (NP,3)
      REAL*8    BS   (NP,3)
      REAL*8    BW   (NP,3)

      REAL*8    UBE  (NP)
      REAL*8    UBS  (NP)
      REAL*8    UBW  (NP)

      REAL*8    UINTE(NP)
      REAL*8    UINTS(NP)
      REAL*8    UINTW(NP)

      REAL*8    SISMO (NBTRAC,NT)

c
c     directives de placement HPF
c

CHPF$ TEMPLATE DOMAIN(NP,NP)
CHPF$ PROCESSORS PE(4,4)
CHPF$ DISTRIBUTE DOMAIN(BLOCK,BLOCK) ONTO PE

CHPF$ ALIGN V(I,J), U(I,J,*) WITH DOMAIN(I,J)

CHPF$ ALIGN UBE(I) WITH DOMAIN(I,799)
CHPF$ ALIGN UBS(I) WITH DOMAIN(799,I)
CHPF$ ALIGN UBW(I) WITH DOMAIN(I,2)

CHPF$ ALIGN BE(I,*), UINTE(I) WITH DOMAIN(I,800)
CHPF$ ALIGN BS(I,*), UINTS(I) WITH DOMAIN(800,I)
CHPF$ ALIGN BW(I,*), UINTW(I) WITH DOMAIN(I,1)

CHPF$ ALIGN SISMO(I,*) WITH DOMAIN(1,750+I)

C  DECLARATION DE SCALAIRES
C  ------------------------
C  NP       : NOMBRE DE POINT SELON LA DIMENSION X DU MAILLAGE SPATIAL
C  NT       : NOMBRE DE PAS DE TEMPS DE LA SIMULATION
C  DELTAT   : INTERVALLE DE TEMPS               (0.0005 S)
C  H        : INTERVALLE DESPACE               (1.5    M)
C  F        : FREQUENCE DE LA SOURCE            (30     HERTZ)
C  IS,JS    : COORDONNEES DE LA SOURCE
C
C  KM,KP    : POINTEURS PERMETTANT DE PERMUTER LES PLANS ESPACE.
C  I,J,K,N  : INDICES DE BOUCLES
C  Wi       : VARIABLES DE TRAVAIL
C  G0       : INTENSITE DE LA SOURCE A LINSTANT t
C  TINIT,
C  TFIN,
C  TMESUR   : COMPTEURS DE GESTION DU TEMPS
C  TiCPU    : COMPTEURS DE GESTION DU TEMPS CPU
C  TiELAP   : COMPTEURS DE GESTION DU TEMPS ELAPSED
C  NBOPi    : COMPTEURS DU NOMBRE DOPERATIONS FLOTTANTES
C  MFLOPS   : NOMBRE DE MFLOPS MESURES
C  UMIN     : VALEUR MIN DE LA SOLUTION U AU DERNIER PAS DE TEMPS
C  UMAX     : VALEUR MAX DE LA SOLUTION U AU DERNIER PAS DE TEMPS
C  VMIN     : VALEUR MIN DE LA SOLUTION DU CHAMP DE VITESSE
C  VMAX     : VALEUR MAX DE LA SOLUTION DU CHAMP DE VITESSE
C  SISMIN   : VALEUR MIN DE LA SOLUTION DU SISMOGRAMME
C  SISMAX   : VALEUR MAX DE LA SOLUTION DU SISMOGRAMME
C
      REAL*8    DELTAT
      REAL*8    H
      REAL*8    F
      INTEGER   IS,JS
      INTEGER   KM,KP
      INTEGER   I,J,K,N
      REAL*8    W
      REAL*8    C1,C2,C3

      REAL*8    G0

c      REAL*8
c     $     TMESUR, TINIT, TFIN, T1CPU, T2CPU, T3CPU,
c     $     NBOP1, NBOP2, MFLOPS

      REAL*8    UMIN, UMAX
      REAL*8    VMIN, VMAX
      REAL*8    SISMIN, SISMAX

      INTEGER   NPTMP, NTTMP
C
C             CORPS DU PROGRAMME
C             ------------------

c
c     ouverture des fichiers
c

      OPEN(UNIT=4,FILE='MODELES/param2.dat')
      OPEN(UNIT=3,FILE='MODELES/modele2.dat')
      OPEN(UNIT=1,FILE='SORTIES/stat2')
      OPEN(UNIT=2,FILE='SORTIES/sismo2')

C
C ACQUISITION DES PARAMETRES DE LA MODELISATION
C -----------------------------------------

C WRITE (6,'('' ENTRER LE NB DE POINTS DU MAILLAGE CARRE (NP) ? '',$)')
      READ (4,*) NPTMP
      IF (NPTMP.NE.NP) STOP
C WRITE (6,'('' ENTRER LE NB DE PAS DE TEMPS TOTAL (NT) ? '',$)')
      READ (4,*) NTTMP
      IF (NTTMP.NE.NT) STOP
C WRITE (6,'('' ENTRER LE PAS DE DISCR. EN ESPACE (H) ? '',$)')
      READ (4,*) H
C WRITE (6,'('' ENTRER LE PAS DE DISCR. EN TEMPS (DELTAT) ? '',$)')
      READ (4,*) DELTAT
C WRITE (6,'('' ENTRER LA POSITION EN X DE LA SOURCE (JS) ? '',$)')
      READ (4,*) JS
C WRITE (6,'('' ENTRER LA POSITION EN Z DE LA SOURCE (IS) ? '',$)')
      READ (4,*) IS
C WRITE (6,'('' ENTRER LA FREQUENCE DE LA SOURCE (F) ? '',$)')
      READ (4,*) F

      WRITE(1,1100)
      WRITE(1,1110)
      WRITE(1,1120)
      WRITE(1,1130)NP,NP
      WRITE(1,1135)NT
      WRITE(1,1140)IS,JS
      WRITE(1,1171)H
      WRITE(1,1170)DELTAT
      WRITE(1,1180)F

C
C INITIALISATION DES VARIABLES SCALAIRES
C ---------------------------------
      KM = 1
      KP = 2
      G0 = 0

C INITIALISATION DU CHAMP VITESSE
C V(I,J)     = READ (Vi) (CAS DUN MILIEU HETEROGENE)
C ----------------------------------------------

      READ (3,*) NPTMP
      IF (NPTMP.NE.NP) STOP
      READ (3,*) NPTMP
      IF (NPTMP.NE.NP) STOP
      DO I = 1,NP
         DO J = 1,NP
            READ (3,*) V(I,J)
         ENDDO
      ENDDO

c
c     reductions MIN et MAX
c
      VMIN = REDMIN(V(1,1), 1, NP, 1, NP)
      VMAX = REDMAX(V(1,1), 1, NP, 1, NP)

      WRITE (1,1051) VMIN,VMAX

C     INITIALISATION DE LA SOLUTION U
C     U(.,.,t=0) = U(.,.,t=1) = 0
C     ET DES ENREGISTREMENTS
C     SISMO (T=1) = 0.0
C     ------------------------------

chpf$ independent(k,j,i)
      DO K = 1,2
         DO J = 1,NP
            DO I = 1,NP
               U(I,J,K) = 0.0
            ENDDO
         ENDDO
      ENDDO

chpf$ independent(i)
      DO I = 1,NBTRAC
         SISMO(I,1) = 0.0
      ENDDO

C     INITIALISATION DES CONSTANTES LIEES AUX COINS INFERIEURS
C            V(I,J)  <--  moyenne harmonique sur 4 points
C            B(I,.,.)<--  cf declarations
C     ET AUX BORDS E-S-W
C            V(I,J)  <--  moyenne harmonique sur 2 points
C            B(I,.,.)<--  cf declarations
C     -------------------------------------------

C     CAS DES TROIS POINTS DU COIN E-S
C     --------------------------------

c
c how to help the compiler...
c
chpf$ independent(i), new(C1,W)
      do i=0,0
         C1          = 4/((1/V(NP-1,NP))+(1/V(NP-1,NP-1)) +
     $        (1/V(NP-2,NP))+(1/V(NP-2,NP-1)) )
         W           = C1*DELTAT
         BE(NP-1,1) = 2*W/(4*W+3*H)
         BE(NP-1,2) = 3*H/(4*W+3*H)
         V(NP-1,NP)  = C1
      enddo

chpf$ independent(i), new(C2,W)
      do i=0,0
         C2          = 4/((1/V(NP,NP))+(1/V(NP,NP-1)) +
     $        (1/V(NP-1,NP))+(1/V(NP-1,NP-1)) )
         W           = C2*DELTAT
         BS(NP,1) = 2*W/(4*W+3*H)
         BS(NP,2) = 3*H/(4*W+3*H)
         V(NP,NP)    = C2
      enddo

chpf$ independent(i), new(C3,W)
      do i=0,0
         C3          = 4/((1/V(NP,NP))+(1/V(NP,NP-1)) +
     $        (1/V(NP-1,NP))+(1/V(NP-1,NP-1)) )
         W           = C3*DELTAT
         BS(NP-1,1) = 2*W/(4*W+3*H)
         BS(NP-1,2) = 3*H/(4*W+3*H)
         V(NP,NP-1)  = C3
      enddo

C     CAS DES 3 POINTS DU COIN W-S
C     ----------------------------

chpf$ independent(i), new(C1,W)
      do i=0,0
         C1          = 4/((1/V(NP,2))+(1/V(NP,3)) +
     $        (1/V(NP-1,2))+(1/V(NP-1,3)) )
         W           = C1*DELTAT
         BS(2,1) = 2*W/(4*W+3*H)
         BS(2,2) = 3*H/(4*W+3*H)
         V(NP,2)    = C1
      enddo

chpf$ independent(i), new(C2,W)
      do i=0,0
         C2          = 4/((1/V(NP,1))+(1/V(NP,2)) +
     $        (1/V(NP-1,1))+(1/V(NP-1,2)) )
         W           = C2*DELTAT
         BS(1,1) = 2*W/(4*W+3*H)
         BS(1,2) = 3*H/(4*W+3*H)
         V(NP,1)    = C2
      enddo

chpf$ independent(i), new(C3,W)
      do i=0,0
         C3          = 4/((1/V(NP-1,1))+(1/V(NP-1,2)) +
     $        (1/V(NP-2,1))+(1/V(NP-2,2)) )
         W           = C3*DELTAT
         BW(NP-1,1) = 2*W/(4*W+3*H)
         BW(NP-1,2) = 3*H/(4*W+3*H)
         V(NP-1,1)  = C3
      enddo

C     CAS DES BORDS   E - S - W
C     -------------------------

c
c     les variables reparties ecrites sont bien alignees
c

chpf$ independent(i), new(W)
      DO I = 2,NP-2
         V(I,NP)   = 2/((1/V(I,NP))+(1/V(I,NP-1)))
         W        = V(I,NP)*DELTAT
         BE(I,1) = 1 - (1/(W+H)) * ( 2*H + W*W/H)
         BE(I,2) = 2*H/(W+H)
         BE(I,3) = W*W/(2*H*(W+H))
      ENDDO

chpf$ independent(i), new(W)
      DO I = 3,NP-2
         V(NP,I)   = 2/((1/V(NP,I))+(1/V(NP-1,I)))
         W        = V(NP,I)*DELTAT
         BS(I,1) = 1 - (1/(W+H)) * ( 2*H + W*W/H)
         BS(I,2) = 2*H/(W+H)
         BS(I,3) = W*W/(2*H*(W+H))
      ENDDO

chpf$ independent(i), new(W)
      DO I = 2,NP-2
         V(I,1)   = 2/((1/V(I,1))+(1/V(I,2)))
         W        = V(I,1)*DELTAT
         BW(I,1) = 1 - (1/(W+H)) * ( 2*H + W*W/H)
         BW(I,2) = 2*H/(W+H)
         BW(I,3) = W*W/(2*H*(W+H))
      ENDDO

C
C     INITIALISATION DES CONSTANTES LIEES AUX BORDS INTERNES
C     ------------------------------------------------------

chpf$ independent(i)
      DO I = 3,NP-1
         V(I,2) = (V(I,2)*DELTAT/H)**2
      ENDDO

chpf$ independent(i)
      DO I = 3,NP-1
         V(I,NP-1) = (V(I,NP-1)*DELTAT/H)**2
      ENDDO

chpf$ independent(i)
      DO I = 3,NP-2
         V(NP-1,I) = (V(NP-1,I)*DELTAT/H)**2
      ENDDO

C     INITIALISATION DES CONSTANTES LIEES AU DOMAINE INTERIEUR
C            V(I,J) <--  1/12*((V(I,J)*DELTAT/H)**2)
C     -------------------------------------------

chpf$ independent(j,i)
      DO J = 3,NP-2
         DO I = 3,NP-2
            V(I,J) = ((V(I,J)*DELTAT/H)**2)/12
         ENDDO
      ENDDO

C     -----------------------------------------
C     !     MODELISATION DE LA PROPAGATION DES ONDES
C     !               PAR EQUATION DONDE.
C     !        RESOLUTION PAR DIFFERENCES FINIES.
C     -------------------------------------------

C     ----------------------------
C     !     BOUCLE EN TEMPS      !
C     ----------------------------

      DO N = 2,NT

C PRISE EN COMPTE DE LA FONCTION SOURCE PENDANT UN TEMPS DE TMAX SEC.
C Un = Un + SOURCE
C ---------------------------------------------------
         W = F * (N * DELTAT - 1.0/F)
         IF (W**2.LT.TMAX) THEN
            W = (W*PI)**2
            G0 = (1.0 - 2.0*W) * EXP (-W)
            U(IS,JS,KM) = U(IS,JS,KM) + G0 * 12*V(IS,JS) * H**2
         ENDIF

C     --------------------------------------------------------
C     SAUVEGARDE DE LA FRANGE CONNEXE AUX BORDS POUR LE CALCUL
C     ULTERIEUR DE LA SOLUTION AUX BORDS
C     --------------------------------------------------------
C     boucle 71

chpf$ independent(i)
         DO I = 2,NP-1
            UBE(I) = U(I,NP-1,KP)
         ENDDO

chpf$ independent(i)
         DO I = 2,NP-1
            UBS(I) = U(NP-1,I,KP)
         ENDDO

chpf$ independent(i)
         DO I = 2,NP-1
            UBW(I) = U(I,2,KP)
         ENDDO

C
C     ------------------------------------------------
C     SOLUTION SUR LE DOMAINE INTERIEUR
C     SCHEMA DORDRE 2 EN TEMPS ET DORDRE 4 EN ESPACE
C     ------------------------------------------------

chpf$ independent(j,i)
         DO J = 3,NP-2
            DO I = 3, NP-2

               U(I,J,KP) =
     $              2 * U(I,J,KM) - U(I,J,KP)
     $              - V(I,J) * ( 60 * U(I,J,KM)
     $              - 16 * ( U(I+1,J,KM) + U(I-1,J,KM)
     $                     + U(I,J-1,KM) + U(I,J+1,KM))
     $              +   U(I+2,J,KM) + U(I-2,J,KM)
     $              + U(I,J-2,KM) + U(I,J+2,KM))

            ENDDO
         ENDDO

C
C     ------------------------------------------------
C     SOLUTION SUR LES BORDS INTERIEURS
C     SCHEMA DORDRE 2 EN TEMPS ET DORDRE 2 EN ESPACE
C     ------------------------------------------------

chpf$ independent(i)
         DO I = 3,NP-1
            U(I,2,KP) =
     $           2 * U(I,2,KM) - U(I,2,KP) - V(I,2) *
     $           ( 4 * U(I,2,KM) - ( U(I+1,2,KM) + U(I-1,2,KM)
     $           + U(I,1,KM) + U(I,3,KM)))
         ENDDO

chpf$ independent(i)
         DO I = 3,NP-1
            U(I,NP-1,KP) =
     $           2 * U(I,NP-1,KM) - U(I,NP-1,KP) - V(I,NP-1) *
     $           ( 4 * U(I,NP-1,KM) - ( U(I+1,NP-1,KM) + U(I-1,NP-1,KM)
     $           + U(I,NP-2,KM) + U(I,NP,KM)))
         ENDDO

chpf$ independent(i)
         DO I = 3,NP-2
            U(NP-1,I,KP) =
     $           2 * U(NP-1,I,KM) - U(NP-1,I,KP) - V(NP-1,I) *
     $           ( 4 * U(NP-1,I,KM) - ( U(NP,I,KM) + U(NP-2,I,KM)
     $           + U(NP-1,I-1,KM) + U(NP-1,I+1,KM)))
         ENDDO

C
C     --------------------------------------------------------
C     SOLUTION AU BORD EAST. CONDITIONS ABSORBANTES DORDRE 2.
C     --------------------------------------------------------

chpf$ independent(i)
         DO I = 3,NP-2
            UINTE(I) =
     $           BE(I,1)*( U(I,NP,KP) + U(I,NP-1,KP) ) - UBE(I)
     $           + BE(I,2)*( U(I,NP,KM) + U(I,NP-1,KM) )
     $           + BE(I,3)*( U(I+1,NP-1,KP) + U(I+1,NP,KP)
     $           + U(I-1,NP-1,KP) + U(I-1,NP,KP) )
         ENDDO

chpf$ independent(i)
         DO I = 3,NP-2
            U(I,NP,KP) = UINTE(I)
         ENDDO

C
C     --------------------------------------------------------
C     SOLUTION AU BORD SOUTH. CONDITIONS ABSORBANTES DORDRE 2.
C     --------------------------------------------------------

chpf$ independent(i)
         DO I = 3,NP-2
            UINTS(I) =
     $           BS(I,1)*( U(NP,I,KP) + U(NP-1,I,KP) ) - UBS(I)
     $           + BS(I,2)*( U(NP,I,KM) + U(NP-1,I,KM) )
     $           + BS(I,3)*( U(NP-1,I+1,KP) + U(NP,I+1,KP)
     $           + U(NP-1,I-1,KP) + U(NP,I-1,KP) )
         ENDDO

chpf$ independent(i)
         DO I = 3,NP-2
            U(NP,I,KP) = UINTS(I)
         ENDDO

C
C     --------------------------------------------------------
C     SOLUTION AU BORD WEST. CONDITIONS ABSORBANTES DORDRE 2.
C     --------------------------------------------------------

chpf$ independent(i)
         DO I = 3,NP-2
            UINTW(I) =
     $           BW(I,1)*( U(I,1,KP) + U(I,2,KP)) - UBW(I)
     $           + BW(I,2)*( U(I,1,KM) + U(I,2,KM) )
     $           + BW(I,3)*( U(I+1,2,KP) + U(I+1,1,KP)
     $           + U(I-1,2,KP) + U(I-1,1,KP) )
         ENDDO

chpf$ independent(i)
         DO I = 3,NP-2
            U(I,1,KP) = UINTW(I)
         ENDDO

C     -----------------------------------------------
C     CALCUL DU BORD DU COIN
C     U(NP-1,NP,KP) ,  U(NP,NP-1,KP) ET U(NP-1,NP,KP)
C     -----------------------------------------------

chpf$ independent(i)
         do i=0,0
            U(NP-1,NP,KP) = BE(NP-1,1)*(U(NP-1,NP-1,KP)+U(NP-2,NP,KP))
     $           + BE(NP-1,2)*U(NP-1,NP,KM)
         enddo

chpf$ independent(i)
         do i=0,0
            U(NP,NP-1,KP) = BS(NP-1,1)*(U(NP,NP-2,KP)+U(NP-1,NP-1,KP))
     $           + BS(NP-1,2)*U(NP,NP-1,KM)
         enddo

chpf$ independent(i)
         do i=0,0
            U(NP,NP,KP)   = BS(NP,1)*(U(NP,NP-1,KP)+U(NP-1,NP,KP))
     $           + BS(NP,2)*U(NP,NP,KM)
         enddo

C     -----------------------------------------------
C     CALCUL DU BORD DU COIN
C     U(NP,2,KP) ,  U(NP-1,1,KP) ET U(NP,1,KP)
C     -----------------------------------------------

chpf$ independent(i)
         do i=0,0
            U(NP-1,1,KP)  = BW(NP-1,1)*(U(NP-1,2,KP)+U(NP-2,1,KP))
     $           + BW(NP-1,2)*U(NP-1,1,KM)
         enddo

chpf$ independent(i)
         do i=0,0
            U(NP,2,KP)    = BS(2,1)*(U(NP,3,KP)+U(NP-1,2,KP))
     $           + BS(2,2)*U(NP,2,KM)
         enddo

chpf$ independent(i)
         do i=0,0
            U(NP,1,KP)    = BS(1,1)*(U(NP-1,1,KP)+U(NP,2,KP))
     $           + BS(1,2)*U(NP,1,KM)
         enddo

C
C     ENREGISTREMENT DU SISMOGRAMME
C     -----------------------------

chpf$ independent(i)
         DO I = 1,NBTRAC
            SISMO(I,N) = U(IS, NP - NBTRAC + I, KP)
         ENDDO

C     PERMUTATION A CHAQUE ITERATION DES PLANS ESPACE
C     -----------------------------------------------
C
         KM = KP
         KP = 3 - KP

      ENDDO

C
C     SORTIE FICHIER DU SISMOGRAMME
C     -----------------------------

      WRITE (2,*) NBTRAC
      WRITE (2,*) NT
      DO K = 1,NBTRAC
         DO J = 1,NT
            WRITE (2,*) SISMO(K,J)
         ENDDO
      ENDDO

C
C     CALCUL DE STATISTIQUES SUR LE RESULTAT
C     --------------------------------------

c
c     reductions
c

      UMIN = REDMIN(U(1,1,KM), 1, NP, 1, NP)
      UMAX = REDMAX(U(1,1,KM), 1, NP, 1, NP)

      WRITE (1,1050) UMIN,UMAX

      SISMIN = REDMIN(SISMO(1,1), 1, NBTRAC, 1, NT)
      SISMAX = REDMAX(SISMO(1,1), 1, NBTRAC, 1, NT)

      WRITE (1,1055) SISMIN,SISMAX

      WRITE(1,1200)

      CLOSE(UNIT=1)
      CLOSE(UNIT=2)
      CLOSE(UNIT=3)
      CLOSE(UNIT=4)

      STOP

C     DECLARATIONS DE FORMATS DE SORTIE
C     ---------------------------------

1100  FORMAT (1X,'DEBUT DU PROGRAMME'//)

1120  FORMAT (1X,'MODELISATION PAR EQUATION D"ONDES',/,
     $        1X,'SCHEMA D"ORDRE 4 EN ESPACE ET 2 EN TEMPS'//)
1130  FORMAT (1X,'TAILLE DU MAILLAGE NPXNP       : ',I4,'  X',I4)
1135  FORMAT (1X,'NOMBRE DE PAS DE TEMPS         : ',I4)
1140  FORMAT (1X,'POSITION DE LA SOURCE          : (',I3,',',I3,')')
1171  FORMAT (1X,'PAS DE DISCR. EN ESPACE (H)    : ',E12.5)
1170  FORMAT (1X,'PAS DE DISCR. EN TEMPS (DELTAT): ',E12.5)
1180  FORMAT (1X,'FREQUENCE DE LA SOURCE (F)     : ',E12.5)
1051  FORMAT (1X,'MIN DE V:',E11.5,3X,'MAX DE V:  ',E11.5//)
1200  FORMAT (1X,/,1X,'FIN DU PROGRAMME')

!1000  FORMAT (1X,'TEMPS CPU DE LA RESOLUTION     (RESOL) : ',E15.5)
!1010  FORMAT (1X,'TEMPS CPU DES INITIALISATIONS   (INIT) : ',E15.5)
!1020  FORMAT (1X,'TEMPS CPU DES ENTREES/SORTIES    (E/S) : ',E15.5)
!1030  FORMAT (1X,'TEMPS CPU TOTAL       (RESOL+INIT+E/S) : ',E15.5//)
!1015  FORMAT (1X,'TEMPS CPU DE L"APPEL AU TIMER          : ',E15.5)
1050  FORMAT (1X,'MIN DE U     :',E15.5,3X,'MAX DE U     : ',E15.5)
1055  FORMAT (1X,'MIN DE SISMO :',E15.5,3X,'MAX DE SISMO : ',E15.5//)
!1060  FORMAT (1X,'MFLOPS CPU PARTIE VECTORIELLE  (RESOL) : ',E15.5)
!1065  FORMAT (1X,'MFLOPS CPU RESOLUTION     (RESOL+INIT) : ',E15.5)
!1070  FORMAT (1X,'MFLOPS CPU GLOBAL     (RESOL+INIT+E/S) : ',E15.5//)

C/////////////////////////////////////////////////////////////////////
C
C   PARTIE SPECIFIQUE A CHAQUE MACHINE ET A CHAQUE MODELISATION
C
C            - FORMATS DES COMMENTAIRES
C            - OUVERTURE DE FICHIERS
C            - FONCTIONS DE GESTION DES TIMERS
C
C/////////////////////////////////////////////////////////////////////

1110  FORMAT (1X,'*********************************************',/,
     $        1X,'       RESOLUTION SUR SUN 4                  ',/,
     $        1X,'     DOUBLE PRECISION (64 bits)              ',/,
     $        1X,'         VERSION STANDART                    ',/,
     $        1X,'                                             ',/,
     $        1X,'*********************************************'//)

      END

c
c     calcul du min
c
      REAL*8 FUNCTION REDMIN(A, LO1, UP1, LO2, UP2)
      INTEGER LO1, UP1, LO2, UP2
      REAL*8 A(LO1:UP1,LO2:UP2), MINA
      INTEGER I,J
      MINA = A(LO1, LO2)
      DO J=LO2, UP2
         DO I=LO1, UP1
            IF (A(I,J).LT.MINA) MINA = A(I,J)
         ENDDO
      ENDDO
      REDMIN = MINA
      RETURN
      END
c
c     calcul du max
c
      REAL*8 FUNCTION REDMAX(A, LO1, UP1, LO2, UP2)
      INTEGER LO1, UP1, LO2, UP2
      REAL*8 A(LO1:UP1,LO2:UP2), MAXA
      INTEGER I,J
      MAXA = A(LO1, LO2)
      DO J=LO2, UP2
         DO I=LO1, UP1
            IF (A(I,J).GT.MAXA) MAXA = A(I,J)
         ENDDO
      ENDDO
      REDMAX = MAXA
      RETURN
      END

ext: f


text:
C Copyright (c) 2007-2017 National Technology & Engineering Solutions of
C Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C
C     * Neither the name of NTESS nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C

C************************************************************************
      subroutine rwpval1(ntxt, ndb, flag, nump, nval, id, names,
     &  values, *)
C************************************************************************
C     This subroutine with read properties from an input text file and
C     write properties to an output exodusII file
C     ntxt   - IN - file id to read from
C     ndb    - IN - file id to write to
C     flag   - IN - flag indicating what type of properties to read/write
C     nump   - IN - number of properties
C     nval   - IN - number of values
C     names  - IN - names of the properties
C     values - IN - values of the properties

      include 'exodusII.inc'

      integer ntxt, ndb, flag, nump, nval
      integer id(*)
      character*(mxstln) names(*)
      integer values(*)

      if (nump .eq. 0) return

      do 100 i = 1, nump
C ... Skip comment line
        read (ntxt, *, end=130, err=130)
        read (ntxt, '(A)', end=120, err=120) names(i)
C ... Skip comment line
        read (ntxt, *, end=130, err=130)
        read (ntxt, *, end=130, err=130) (values(j), j = 1, nval)

C ... Write the property (unless it is 'ID')
        if (names(i) .ne. 'ID                              ') then
          do 90 j = 1, nval
            call expp(ndb, flag, id(j), names(i), values(j), ierr)
 90       continue
        end if
 100  continue
      return

 120  continue
 130  continue
      CALL PRTERR('FATAL', 'Problem reading property')
      return 1
      end

max_stars_repo_name: brian-kelley/seacas

max_issues_repo_name: brian-kelley/seacas

max_forks_repo_name: brian-kelley/seacas


text:
*###[ aaxcx :
	subroutine aaxcx(caxi,cbxi,ccxi,d0,xmm,xpi,level,ier)
***#[*comment:***********************************************************
*									*
*	Calculation of formfactors resulting from decvert.sub		*
*	or decvert.frm (up to third rank)				*
*	21-dec-1993: switched to ffxc1 for C1, added numerical checks.	*
*	Definitions:							*
*       C0 = 1/(i pi^2)*\int d^4 Q					*
*		1/((Q^2-m_1^2)((Q+p1)^2-m2^2)((Q-p3)^2-m3^2))		*
*	C1 = 1/(i pi^2)*\int d^n Q Q(mu)/(...)				*
*	   = C11*p1 + C12*p2						*
*	C2 = C21*p1*p1 + C22*p2*p2 + C23*(p1*p2+p2*p1) + C24*g		*
*	C3 = C31*p1*p1*p1 + C32*p2*p2*p2 + C33*(p1*p1*p2 + p1*p2*p1 +	*
*	     p2*p1*p1) + C34*(p1*p2*p2 + p2*p1*p2 + p1*p2*p2) + C35*	*
*	     (p1*g + g*p1 + 'g*p1*g') + C36*(p2*g + g*p2 + 'g*p2*g')	*
*									*
*	Input:	xpi	  the same as in Geert Jan's routines		*
*		level	  rank of tensor(integral)			*
*	Output: caxi(3)  : ca0i 		     i=1,2,3		*
*		cbxi(12) : (cb0i,cb11i,cb21i,cb22i)  i=1,2,3		*
*		ccxi(13) : cc0,cc1(2),cc2(4),cc3(6)			*
*									*
***#]*comment:***********************************************************
*  #[ declarations :
	implicit none
*
*	arguments
*
	integer ier,level
	DOUBLE PRECISION xpi(6),d0,xmm
	DOUBLE COMPLEX caxi(3),cbxi(12),ccxi(13)
*
*	local variables
*
	DOUBLE PRECISION maxi(3),mbxi(12),mcxi(13)
*
*	common blocks
*
	include 'ff.h'
	include 'aa.h'
*
*  #] declarations :
*  #[ call ffxcx:
*
	call ffxcx(caxi,maxi,cbxi,mbxi,ccxi,mcxi,d0,xmm,xpi,level,ier)
*
*  #] call ffxcx:
*###] aaxcx :
	end
*###[ ffxcx:
	subroutine ffxcx(caxi,maxi,cbxi,mbxi,ccxi,mcxi,d0,xmm,xpi,level,
     +		ier)
***#[*comment:***********************************************************
*									*
*	Calculation of three point form factors with more accurate	*
*	error estimates.  Calls ffxc1, the rest is still here.		*
*	Definitions:							*
*       C0 = 1/(i pi^2)*\int d^4 Q					*
*		1/((Q^2-m_1^2)((Q+p1)^2-m2^2)((Q-p3)^2-m3^2))		*
*	C1 = 1/(i pi^2)*\int d^n Q Q(mu)/(...)				*
*	   = C11*p1 + C12*p2						*
*	C2 = C21*p1*p1 + C22*p2*p2 + C23*(p1*p2+p2*p1) + C24*g		*
*	C3 = C31*p1*p1*p1 + C32*p2*p2*p2 + C33*(p1*p1*p2 + p1*p2*p1 +	*
*	     p2*p1*p1) + C34*(p1*p2*p2 + p2*p1*p2 + p1*p2*p2) + C35*	*
*	     (p1*g + g*p1 + 'g*p1*g') + C36*(p2*g + g*p2 + 'g*p2*g')	*
*									*
*	Input:	xpi(6)	     real	m_i^2 (1:3), p_{i-3}^2 (4:6)	*
*		d0,xmu	     real	renormalisation constants	*
*		level	     integer	rank of tensor (integral)	*
*	Output: caxi(3)	     complex	A0(m_i^2)			*
*		maxi(3)	     real	max term in sum to caxi()	*
*		cbxi(12)     complex	3x(B0,B11,B21,B22)(p_i^2)	*
*		mbxi(12)     real	max term in sum to cbxi()	*
*		ccxi(13)     complex	C0,C1(2),C2(4),C3(6)		*
*		mcxi(13)     real	max term in sum to ccxi()	*
*	Note that if level<3 some of these are not defined.		*
*									*
***#]*comment:***********************************************************
*  #[ declarations :
	implicit none
*
*	arguments
*
	integer ier,level
	DOUBLE PRECISION maxi(3),mbxi(12),mcxi(13),xpi(6),d0,xmm
	DOUBLE COMPLEX caxi(3),cbxi(12),ccxi(13)
*
*	local variables
*
	integer i,bl,ier0,ier1
	logical adesav
	DOUBLE PRECISION absc,ma0i(6),mabxi(2),big
	DOUBLE COMPLEX acbxi(2),ca0i(6),cc,cc0
*
*	common blocks
*
	include 'ff.h'
	include 'aa.h'
*
*	statement function
*
	absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
*
*  #] declarations :
*  #[ initialisations:
*
*	initialization to nonsense
*
	big = 1/(1d20*xclogm)
	if ( ltest ) then
	    do 10 i=1,3
		caxi(i) = big
   10	    continue
	    do 20 i=1,12
		cbxi(i) = big
   20	    continue
	    do 30 i=1,13
		ccxi(i) = big
   30	    continue
	endif
*
*  #] initialisations:
*  #[ get C0:
*
*	C0-function
*
	ldot=.TRUE.
	ier1 = ier
	call ffxc0(ccxi(1),xpi,ier1)
	if ( ier1.gt.10 ) then
	    if ( ltest ) then
		print *,'ffxcx: id = ',id,', nevent = ',nevent
		print *,'ffxcx: lost ',ier1,' digits in C0 with isgnal '
     +			,isgnal,', trying other roots, isgnal ',-isgnal
		print *,'       if OK (no further messages) adding this'
     +			,' to your code will improve speed'
	    endif
	    isgnal = -isgnal
	    ier0 = ier
	    call ffxc0(cc0,xpi,ier0)
	    isgnal = -isgnal
	    if ( ier0.lt.ier1 ) then
	    	ier1 = ier0
	    	ccxi(1) = cc0
	    endif
	endif
	if ( ier1 .gt. 10 ) then
	    print *,'ffxcx: id = ',id,', nevent = ',nevent
	    print *,'ffxcx: error: C0 not stable, lost ',ier1,' digits'
	    print *,'       please contact author (t19@nikhef.nl)'
	    print *,'xpi = ',xpi
	endif
*	note that we may have lost another factor xloss**3 or so
	mcxi(1) = absc(ccxi(1))*DBLE(10)**mod(ier1,50)
	if ( awrite ) then
*	 #[ for debugging: imported stuff from ff
	    print *,' '
	    print *,'ffxcx : level 0 '
	    print *,'C0 =',ccxi(1),mcxi(1),ier1
	    print *,'used:',( xpi(i),i=1,3 )
	    print *,'     ',( xpi(i),i=4,6 )
	    print *,'imported stuff via ff.h:'
	    print *,'kin det = ',fdel2
	    print *,'dotpr1,1= ',fpij3(4,4)
	    print *,'dotpr2,2= ',fpij3(5,5)
	    print *,'dotpr1,2= ',fpij3(4,5)
*	 #] for debugging:
	endif

	if ( level.eq.0 ) goto 990
*
*  #] get C0:
*  #[ need B-functions till b-level=(level-1):
	bl=level-1
	if ( awrite ) then
	    print '(a,i1)','  ##[ B-function output: up to level ',bl
	endif
	adesav = aderiv
	aderiv = .FALSE.
	ier0 = ier
	call ffxbx( ca0i(1),ma0i(1),cbxi(1),mbxi(1),acbxi,mabxi,
     +		d0,xmm,xpi(5),xpi(2),xpi(3),bl,ier0)
	ier1 = max(ier1,ier0)
	ier0 = ier
	call ffxbx( ca0i(3),ma0i(3),cbxi(5),mbxi(5),acbxi,mabxi,
     +		d0,xmm,xpi(6),xpi(1),xpi(3),bl,ier0)
	ier1 = max(ier1,ier0)
	ier0 = ier
	call ffxbx( ca0i(5),ma0i(5),cbxi(9),mbxi(9),acbxi,mabxi,
     +		d0,xmm,xpi(4),xpi(1),xpi(2),bl,ier0)
	ier1 = max(ier1,ier0)
	aderiv = adesav
	if ( awrite ) then
	    print '(a)','  ##] B-function output:'
	endif
*		symmetry in A0(i,j)
	caxi(1)=ca0i(1)
	caxi(2)=ca0i(2)
	caxi(3)=ca0i(3)
	maxi(1)=ma0i(1)
	maxi(2)=ma0i(2)
	maxi(3)=ma0i(3)
	if ( lwarn .and. atest ) then
	    if ((ca0i(4)-ca0i(2)) .ne. 0. .or.
     +		(ca0i(5)-ca0i(3)) .ne. 0. .or.
     +		(ca0i(6)-ca0i(1)) .ne. 0. ) then
		print *,'error in A0-calculations in aaxbx.for'
	    endif
	endif
*  #] need B-functions till b-level=(level-1):
*  #[ break to let ffzcz tie in:
	call ffxcxp(caxi,maxi,cbxi,mbxi,ccxi,mcxi,xpi,level,ier1)
*  #] break to let ffzcz tie in:
  990	ier = ier1
	end
	subroutine ffxcxp(caxi,maxi,cbxi,mbxi,ccxi,mcxi,xpi,level,ier)
*  #[ declarations :
	implicit none
*
*	arguments
*
	integer ier,level
	DOUBLE PRECISION maxi(3),mbxi(12),mcxi(13),xpi(6)
	DOUBLE COMPLEX caxi(3),cbxi(12),ccxi(13)
*
*	local variables
*
	integer i,j,ier1,ier2
	DOUBLE PRECISION xi3(3),f1,f2,absc,xmax,R1m,R2m,R3m,R4m,
     +		R5m,R6m,R11m,R12m,R13m,R14m,R15m,R16m,R17m,R18m
	DOUBLE PRECISION mb0i(3),mb11i(3),mxy(2),mb21i(3),mb22i(3)
	DOUBLE COMPLEX R1,R2,R3,R4,R5,R6,R11,R12,R13,R14,R15,R16,R17,
     +		R18,cb0i(3),cb11i(3),cb21i(3),cb22i(3),cc,cxy(2)
*
*	common blocks
*
	include 'ff.h'
	include 'aa.h'
*
*	statement function
*
	absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
*
*  #] declarations :
*  #[ kinematical quantities for 3pv-red :
*	       inverse kinematical matrix xi3  (2X2)
	ier2 = ier
	call aaxi3(xi3,xpi,ier2)
	ier2 = ier2 - ier
*
*	       f-functions:
	f1 = 2*fpij3(1,4)
	f2 = 2*fpij3(1,5)
*  #] kinematical quantities for 3pv-red :
*  #[ level 1 : C11,C12,B0(I)
*		need 3 diff B0(I)-functions,I=1,2,3
	cb0i(1)=cbxi(1)
	cb0i(2)=cbxi(5)
	cb0i(3)=cbxi(9)
	mb0i(1)=mbxi(1)
	mb0i(2)=mbxi(5)
	mb0i(3)=mbxi(9)
	call ffxc1a(ccxi(2),mcxi(2),ccxi(1),mcxi(1),cb0i,mb0i,
     +		xpi,fpij3,fdel2,ier)
	if ( awrite ) then
	    print *,'GEERT JANs-scheme:'
	    print *,'C11=',ccxi(2),mcxi(2),ier
	    print *,'C12=',ccxi(3),mcxi(3),ier
	    print *,' '
	endif
	if ( lwarn .and. atest ) then
*		PV-reduction
	    R1=( f1*ccxi(1)+cb0i(2)-cb0i(1) )/2
	    R2=( f2*ccxi(1)+cb0i(3)-cb0i(2) )/2
	    R1m=max(abs(f1)*mcxi(1),mb0i(2),mb0i(1))/2
	    R2m=max(abs(f2)*mcxi(1),mb0i(3),mb0i(2))/2
	    cxy(1)=xi3(1)*R1+xi3(3)*R2
	    cxy(2)=xi3(3)*R1+xi3(2)*R2
	    mxy(1)=max(abs(xi3(1))*R1m,abs(xi3(3))*R2m)
	    mxy(2)=max(abs(xi3(3))*R1m,abs(xi3(2))*R2m)
	    if ( xloss*absc(ccxi(2)-cxy(1)) .gt. precc*
     +		max(mcxi(2),mxy(1)) )
     +		print *,'ffxcxp: error: FF C11 disagrees with PV: ',
     +		ccxi(2),cxy(1),ccxi(2)-cxy(1),ier
	    if ( xloss*absc(ccxi(3)-cxy(2)) .gt. precc*
     +		max(mcxi(3),mxy(2)) )
     +		print *,'ffxcxp: error: FF C12 disagrees with PV: ',
     +		ccxi(3),cxy(2),ccxi(3)-cxy(2),ier
	    if (awrite) then
		print *,' '
		print *,'ffxcxp : level 1: id,nevent ',id,nevent
		print *,'C11=',ccxi(2)
		print *,'C12=',ccxi(3)
	    endif
	endif
*
	if ( level.eq.1 ) then
	    if ( lwarn ) then
		xmax = 0
		do i=1,3
		    if ( absc(ccxi(i)).ne.0 ) then
			xmax = max(xmax,mcxi(i)/absc(ccxi(i)))
		    elseif ( mcxi(i).ne.0 ) then
			xmax = max(xmax,1/precc)
		    endif
		enddo
		ier1 = int(log10(xmax))
		if ( awrite ) print *,'ier = ',ier1
	    else
		ier1 = 0
	    endif
	    goto 990
	endif
*
*  #] level 1 :
*  #[ level 2 : C21,C22,C23,C24,B11(I),A0(I,J)
*		need 3 diff B1-functions and 3 diff A0-fuctions
	do 12 i=1,3
	    j=(i+1)+(i-1)*3
	    cb11i(i) = cbxi(j)
	    mb11i(i) = mbxi(j)
   12	continue
*		PV-reduction
	ccxi(7)=1/4.d0 + 1/2.d0*xpi(1)*ccxi(1) -
     +		1/4.d0*( f1*ccxi(2)+f2*ccxi(3)-cb0i(1)  )
	ier1 = ier
	if ( lwarn ) then
***	    c7max = max(x1,2*xpi(1)*absc(ccxi(1)),absc(f1*ccxi(2)),
***     +		absc(f2*ccxi(3)),absc(cb0i(1)))/4
***	    if ( absc(ccxi(7)) .lt. xloss*c7max ) then
***		call ffwarn(293,ier1,absc(ccxi(7)),c7max)
***	    endif
	    mcxi(7) = max(x1,2*xpi(1)*mcxi(1),abs(f1)*mcxi(2),
     +		abs(f2)*mcxi(3),mb0i(1))/4
	endif
	R3=( f1*ccxi(2) + cb11i(2) + cb0i(1)  )/2 - ccxi(7)
	R4=( f2*ccxi(2) + cb11i(3) - cb11i(2) )/2
	R5=( f1*ccxi(3) + cb11i(2) - cb11i(1) )/2
	R6=( f2*ccxi(3)            - cb11i(2) )/2 - ccxi(7)
	ccxi(4)=xi3(1)*R3 + xi3(3)*R4
	ccxi(5)=xi3(3)*R5 + xi3(2)*R6
	ccxi(6)=xi3(3)*R3 + xi3(2)*R4
	if ( lwarn ) then
***	    R3m = max(absc(f1*ccxi(2)),absc(cb11i(2)),absc(cb0i(1)),
***     +		2*c7max)/2
***	    R4m = max(absc(f2*ccxi(2)),absc(cb11i(3)),absc(cb11i(2)))/2
***	    R5m = max(absc(f1*ccxi(3)),absc(cb11i(2)),absc(cb11i(1)))/2
***	    R6m = max(absc(f2*ccxi(3)),absc(cb11i(2)),2*c7max)/2
***	    xmax = max(abs(xi3(1))*R3m,abs(xi3(3))*R4m)
***	    if ( absc(ccxi(4)).lt.xloss*xmax ) then
***		ier0 = ier
***		call ffwarn(292,ier0,absc(ccxi(4)),xmax)
***		ier1 = max(ier1,ier0)
***	    endif
***	    xmax = max(abs(xi3(3))*R5m,abs(xi3(2))*R6m)
***	    if ( absc(ccxi(5)).lt.xloss*xmax ) then
***		ier0 = ier
***		call ffwarn(291,ier0,absc(ccxi(5)),xmax)
***		ier1 = max(ier1,ier0)
***	    endif
***	    xmax = max(abs(xi3(3))*R3m,abs(xi3(2))*R4m)
***	    if ( absc(ccxi(6)).lt.xloss*xmax ) then
***		ier0 = ier
***		call ffwarn(290,ier0,absc(ccxi(6)),xmax)
***		ier1 = max(ier1,ier0)
***	    endif
*
*	    for the whole chain
*
	    R3m = max(abs(f1)*mcxi(2),mb11i(2),mb0i(1),2*mcxi(7))/2
	    R4m = max(abs(f2)*mcxi(2),mb11i(3),mb11i(2))/2
	    R5m = max(abs(f1)*mcxi(3),mb11i(2),mb11i(1))/2
	    R6m = max(abs(f2)*mcxi(3),mb11i(2),2*mcxi(7))/2
	    mcxi(4) = max(abs(xi3(1))*R3m,abs(xi3(3))*R4m)
	    mcxi(5) = max(abs(xi3(3))*R5m,abs(xi3(2))*R6m)
	    mcxi(6) = max(abs(xi3(3))*R3m,abs(xi3(2))*R4m)
	endif
	if ( lwarn .and. atest ) then
	    cxy(1) = xi3(1)*R5 + xi3(3)*R6
	    mxy(1) = abs(xi3(1))*R5m + abs(xi3(3))*R6m
	    if ( xloss*absc(cxy(1)-ccxi(6)).gt.precc*max(mcxi(6),mxy(1))
     +			) then
		print *,'ffxcxp: error: id/nevent ',id,'/',nevent
		print *,'redundancy check at level 2 failed: '
		print *,cxy(1),mxy(1)
		print *,ccxi(6),mcxi(6)
		print *,absc(cxy(1)-ccxi(6))
	    endif
	endif
	if ( awrite ) then
	    print *,' '
	    print *,'ffxcxp : level 2: id,nevent ',id,nevent
	    print *,'C21=',ccxi(4),mcxi(4)
	    print *,'C22=',ccxi(5),mcxi(5)
	    print *,'C23=',ccxi(6),mcxi(6)
	    print *,'    ',cxy(1),mxy(1)
	    print *,'C24=',ccxi(7),mcxi(7)
	endif

	if ( level.eq.2 ) then
	    if ( lwarn ) then
		xmax = 0
		do i=1,7
		    if ( absc(ccxi(i)).ne.0 ) then
			xmax = max(xmax,mcxi(i)/absc(ccxi(i)))
		    elseif ( mcxi(i).ne.0 ) then
			xmax = max(xmax,1/precc)
		    endif
		enddo
		ier1 = int(log10(xmax))
		if ( awrite ) print *,'ier = ',ier1
	    else
		ier1 = 0
	    endif
	    goto 990
	endif

*  #] level 2 :
*  #[ level 3 : C31,C32,C33,C34,C35,C36,B21(I),B22(I)
	do 13 i=1,3
	    j = (i+1)+(i-1)*3
	    cb21i(i)=cbxi(j+1)
	    cb22i(i)=cbxi(j+2)
	    mb21i(i)=mbxi(j+1)
	    mb22i(i)=mbxi(j+2)
   13	continue
*		PV-reduction
	R17=( f1*ccxi(7)+cb22i(2)-cb22i(1) )/2
	R18=( f2*ccxi(7)+cb22i(3)-cb22i(2) )/2
	ccxi(12)=xi3(1)*R17+xi3(3)*R18
	ccxi(13)=xi3(3)*R17+xi3(2)*R18
	if ( lwarn ) then
***	    R17m = max(abs(f1)*c7max,absc(cb22i(2)),absc(cb22i(1)))/2
***	    R18m = max(abs(f2)*c7max,absc(cb22i(3)),absc(cb22i(2)))/2
***	    c12max = max(abs(xi3(1))*R17m,abs(xi3(3))*R18m)
***	    if ( absc(ccxi(12)).lt.xloss*c12max ) then
***		ier0 = ier
***		call ffwarn(289,ier0,absc(ccxi(12)),c12max)
***		ier1 = max(ier1,ier0)
***	    endif
***	    c13max = max(abs(xi3(3))*R17m,abs(xi3(2))*R18m)
***	    if ( absc(ccxi(13)).lt.xloss*c13max ) then
***		ier0 = ier
***		call ffwarn(288,ier0,absc(ccxi(13)),c13max)
***		ier1 = max(ier1,ier0)
***	    endif
	    R17m = max(abs(f1)*mcxi(7),mb22i(2),mb22i(1))/2
	    R18m = max(abs(f2)*mcxi(7),mb22i(3),mb22i(2))/2
	    mcxi(12) = max(abs(xi3(1))*R17m,abs(xi3(3))*R18m)
	    mcxi(13) = max(abs(xi3(3))*R17m,abs(xi3(2))*R18m)
	endif
	R11=( f1*ccxi(4)+cb21i(2)-cb0i(1)  )/2 - 2*ccxi(12)
	R12=( f2*ccxi(4)+cb21i(3)-cb21i(2) )/2
	R13=( f1*ccxi(5)+cb21i(2)-cb21i(1) )/2
	R14=( f2*ccxi(5)         -cb21i(2) )/2 - 2*ccxi(13)
	R15=( f1*ccxi(6)+cb21i(2)+cb11i(1) )/2 - ccxi(13)
	R16=( f2*ccxi(6)         -cb21i(2) )/2 - ccxi(12)
	ccxi(8) =xi3(1)*R11 + xi3(3)*R12
	ccxi(9) =xi3(3)*R13 + xi3(2)*R14
	ccxi(10)=xi3(3)*R11 + xi3(2)*R12
	ccxi(11)=xi3(1)*R13 + xi3(3)*R14
	if ( lwarn ) then
***	    R11m = max(absc(f1*ccxi(4)),absc(cb21i(2)),absc(cb0i(1)),
***     +		2*c12max)/2
***	    R12m = max(absc(f2*ccxi(4)),absc(cb21i(3)),absc(cb21i(2)))/2
***	    R13m = max(absc(f1*ccxi(5)),absc(cb21i(2)),absc(cb21i(1)))/2
***	    R14m = max(absc(f2*ccxi(5)),absc(cb21i(2)),4*c13max)/2
***	    R15m = max(absc(f1*ccxi(6)),absc(cb21i(2)),absc(cb11i(1)),
***     +		2*c13max)/2
***	    R16m = max(absc(f2*ccxi(6)),absc(cb21i(2)),2*c12max)/2
***	    xmax = max(abs(xi3(1))*R11m,abs(xi3(3))*R12m)
***	    if ( absc(ccxi(8)).lt.xloss*xmax ) then
***		ier0 = ier
***		call ffwarn(287,ier0,absc(ccxi(8)),xmax)
***		ier1 = max(ier1,ier0)
***	    endif
***	    xmax = max(abs(xi3(3))*R13m,abs(xi3(2))*R14m)
***	    if ( absc(ccxi(9)).lt.xloss*xmax ) then
***		ier0 = ier
***		call ffwarn(286,ier0,absc(ccxi(9)),xmax)
***		ier1 = max(ier1,ier0)
***	    endif
***	    xmax = max(abs(xi3(3))*R11m,abs(xi3(2))*R12m)
***	    if ( absc(ccxi(10)).lt.xloss*xmax ) then
***		ier0 = ier
***		call ffwarn(285,ier0,absc(ccxi(10)),xmax)
***		ier1 = max(ier1,ier0)
***	    endif
***	    xmax = max(abs(xi3(1))*R13m,abs(xi3(3))*R14m)
***	    if ( absc(ccxi(11)).lt.xloss*xmax ) then
***		ier0 = ier
***		call ffwarn(284,ier0,absc(ccxi(11)),xmax)
***		ier1 = max(ier1,ier0)
***	    endif
*
*	    for the whole chain
*
	    R11m = max(abs(f1)*mcxi(4),mb21i(2),mb0i(1),2*mcxi(12))/2
	    R12m = max(abs(f2)*mcxi(4),mb21i(3),mb21i(2))/2
	    R13m = max(abs(f1)*mcxi(5),mb21i(2),mb21i(1))/2
	    R14m = max(abs(f2)*mcxi(5),mb21i(2),4*mcxi(13))/2
	    R15m = max(abs(f1)*mcxi(6),mb21i(2),mb11i(1),2*mcxi(13))/2
	    R16m = max(abs(f2)*mcxi(6),mb21i(2),2*mcxi(12))/2
	    mcxi(8) = max(abs(xi3(1))*R11m,abs(xi3(3))*R12m)
	    mcxi(9) = max(abs(xi3(3))*R13m,abs(xi3(2))*R14m)
	    mcxi(10)= max(abs(xi3(3))*R11m,abs(xi3(2))*R12m)
	    mcxi(11)= max(abs(xi3(1))*R13m,abs(xi3(3))*R14m)
	endif
*	redundancy check
	if ( lwarn .and. atest ) then
	    cxy(1) = xi3(1)*R15 + xi3(3)*R16
	    cxy(2) = xi3(3)*R15 + xi3(2)*R16
	    mxy(1) = abs(xi3(1))*R15m + abs(xi3(3))*R16m
	    mxy(2) = abs(xi3(3))*R15m + abs(xi3(2))*R16m
	    if ( xloss*absc(cxy(1)-ccxi(10)).gt.precc*max(mxy(1),
     +			mcxi(10))
     +	    .or. xloss*absc(cxy(2)-ccxi(11)).gt.precc*max(mxy(2),
     +			mcxi(11)) ) then
		print *,'ffxcxp: error: id/nevent ',id,'/',nevent
		print *,'redundancy check at level 3 failed: '
		print *,cxy(1),mxy(1)
		print *,ccxi(10),mcxi(10)
		print *,absc(cxy(1)-ccxi(10))
		print *,cxy(2),mxy(2)
		print *,ccxi(11),mcxi(11)
		print *,absc(cxy(1)-ccxi(11))
	    endif
	endif
	if ( awrite ) then
	    print *,' '
	    print *,'ffxcxp: level 3: id,nevent ',id,nevent
	    print *,'C31=',ccxi(8),mcxi(8)
	    print *,'C32=',ccxi(9),mcxi(9)
	    print *,'C33=',ccxi(10),mcxi(10)
	    print *,'    ',cxy(1),mxy(1)
	    print *,'C34=',ccxi(11),mcxi(11)
	    print *,'    ',cxy(2),mxy(2)
	    print *,'C35=',ccxi(12),mcxi(12)
	    print *,'C36=',ccxi(13),mcxi(13)
	endif

	if ( level.eq.3 ) then
	    if ( lwarn ) then
		xmax = 0
		do i=1,13
		    if ( absc(ccxi(i)).ne.0 ) then
			xmax = max(xmax,mcxi(i)/absc(ccxi(i)))
		    elseif ( mcxi(i).ne.0 ) then
			xmax = max(xmax,1/precc)
		    endif
		enddo
		ier1 = int(log10(xmax))
		if ( awrite ) print *,'ier = ',ier1
	    else
		ier1 = 0
	    endif
	    goto 990
	endif

*  #] level 3 :
*  #[ end:
	print *,'ffxcxp: level ',level,' not supported.'
	stop
  990	continue
	ier = ier1 + ier2
*  #] end:
*###] ffxcx:
	end


text:
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc                                                                    cc
cc              mnsuev : surev test program                           cc
cc                                                                    cc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real fac
      integer i,idim,ier,j,m,mu,mv,m0,m1,m2,m3,nc,nu4,nv4,nu,nv,l
      real tu(11),tv(10),c(126),u(6),v(6),f(108),wrk(48)
      integer iwrk(12)
c  we set up the grid points for evaluating the spline surface.
      mu = 6
      mv = 6
      do 10 i=1,6
      u(i) = (i-1)*0.2
      v(i) = u(i)
  10  continue
c  the interior knots with respect to the u-variable.
      tu(5) = 0.4
      tu(6) = 0.7
      tu(7) = 0.9
      nu = 11
c  the interior knots with respect to the v-variable.
      tv(5) = 0.3
      tv(6) = 0.8
      nv = 10
c  the boundary knots
      do 20 i=1,4
        tu(i) = 0.
        tv(i) = 0.
        tu(i+7) = 1.
        tv(i+6) = 1.
  20  continue
c  we generate the b-spline coefficients for the test surface
c        x = u*v    y = v**2    z = u+v     0 <= u,v <= 1
c  the dimension of the surface
      idim = 3
c  the number of b-spline coefficients for each co-ordinate
      nu4 = nu-4
      nv4 = nv-4
      nc = nu4*nv4
c  the coefficients for x = u*v
      do 30 i=1,nv4
        c(i) = 0.
  30  continue
      do 40 i=2,nu4
        c((i-1)*nv4+1) = 0.
  40  continue
      m0 = 1
      do 60 i=2,nu4
        m1 = m0+nv4
        fac = (tu(i+3)-tu(i))/9.
        do 50 j=2,nv4
          m2 = m0+1
          m3 = m1+1
          c(m3) = c(m1)+c(m2)-c(m0)+fac*(tv(j+3)-tv(j))
          m0 = m0+1
          m1 = m1+1
  50    continue
        m0 = m0+1
  60  continue
c  the coefficients for y = v**2.
      l = nc
      m0 = l+1
      m1 = m0+1
      c(m0) = 0.
      c(m1) = 0.
      do 70 i=3,nv4
        c(m1+1) = c(m1)+(tv(i+3)-tv(i))*((c(m1)-c(m0))/(tv(i+2)-tv(i-1))
     *    +(tv(i+2)-tv(i))/3.)
        m0 = m1
        m1 = m0+1
  70  continue
      do 80 i=1,nv4
        m0 = l+i
        fac = c(m0)
        do 80 j=1,nu4
          m0 = m0+nv4
          c(m0) = fac
  80  continue
c  the coefficients for z = u+v
      l = l+nc
      m0 = l+1
      c(m0) = 0.
      do 90 i=2,nv4
        m1 = m0+1
        c(m1) = c(m0)+(tv(i+3)-tv(i))/3.
        m0 = m1
  90  continue
      do 100 i=1,nv4
        m0 = l+i
        do 100 j=2,nu4
          m1 = m0+nv4
          c(m1) = c(m0)+(tu(j+3)-tu(j))/3.
          m0 = m1
 100  continue
c  evaluation of the spline surface
      call surev(idim,tu,nu,tv,nv,c,u,mu,v,mv,f,108,wrk,48,iwrk,12,ier)
c  printing of the results
      write(6,900)
      write(6,910)
      write(6,920) (tu(i),i=1,nu)
      write(6,930)
      write(6,920) (tv(i),i=1,nv)
      write(6,940)
      m1 = 0
      do 110 l=1,idim
        m0 = m1+1
        m1 = m1+nc
        write(6,950) (c(j),j=m0,m1)
 110  continue
      write(6,960)
      write(6,970) (v(i),i=1,mv)
      write(6,980)
      m = mu*mv
      m0 = 0
      do 130 i=1,mu
        write(6,990) u(i)
        m1 = m0
        do 120 l=1,idim
          m2 = m1+1
          m3 = m1+mv
          write(6,995) (f(j),j=m2,m3)
          m1 = m1+m
 120    continue
        m0 = m0+mv
 130  continue
      stop
c  format statements.
 900  format(23h0bicubic spline surface)
 910  format(1x,40hposition of the knots in the u-direction)
 920  format(1x,15f5.1)
 930  format(1x,40hposition of the knots in the v-direction)
 940  format(23h b-spline coefficients )
 950  format(5x,8f9.4)
 960  format(1h0,37hspline values at selected grid points)
 970  format(1h0,2x,1hv,6(3x,f4.1))
 980  format(1h ,1x,1hu)
 990  format(1h ,f4.1)
 995  format(5x,6f7.3)
      end

max_issues_repo_name: tbeu/ApproxSpline


text:
      subroutine chol1d ( p, v, qty, npoint, ncol, u, qu )
c  from  * a practical guide to splines *  by c. de boor
c  from  * a practical guide to splines *  by c. de boor
c to be called in  s m o o t h
constructs the upper three diags. in v(i,j), i=2,,npoint-1, j=1,3, of
c  the matrix  6*(1-p)*q-transp.*(d**2)*q + p*r, then computes its
c  l*l-transp. decomposition and stores it also in v, then applies
c  forward and backsubstitution to the right side q-transp.*y in  qty
c  to obtain the solution in  u .
      integer ncol,npoint,   i,npm1,npm2
      real p,qty(npoint),qu(npoint),u(npoint),v(npoint,7),   prev,ratio
     *	  ,six1mp,twop
      npm1 = npoint - 1
c     construct 6*(1-p)*q-transp.*(d**2)*q  +  p*r
      six1mp = 6.*(1.-p)
      twop = 2.*p
      do 2 i=2,npm1
	 v(i,1) = six1mp*v(i,5) + twop*(v(i-1,4)+v(i,4))
	 v(i,2) = six1mp*v(i,6) + p*v(i,4)
    2	 v(i,3) = six1mp*v(i,7)
      npm2 = npoint - 2
      if (npm2 .ge. 2)			go to 10
      u(1) = 0.
      u(2) = qty(2)/v(2,1)
      u(3) = 0.
					go to 41
c  factorization
   10 do 20 i=2,npm2
	 ratio = v(i,2)/v(i,1)
	 v(i+1,1) = v(i+1,1) - ratio*v(i,2)
	 v(i+1,2) = v(i+1,2) - ratio*v(i,3)
	 v(i,2) = ratio
	 ratio = v(i,3)/v(i,1)
	 v(i+2,1) = v(i+2,1) - ratio*v(i,3)
   20	 v(i,3) = ratio
c
c  forward substitution
      u(1) = 0.
      v(1,3) = 0.
      u(2) = qty(2)
      do 30 i=2,npm2
   30	 u(i+1) = qty(i+1) - v(i,2)*u(i) - v(i-1,3)*u(i-1)
c  back substitution
      u(npoint) = 0.
      u(npm1) = u(npm1)/v(npm1,1)
      i = npm2
   40	 u(i) = u(i)/v(i,1)-u(i+1)*v(i,2)-u(i+2)*v(i,3)
	 i = i - 1
	 if (i .gt. 1)			go to 40
c  construct q*u
   41 prev = 0.
      do 50 i=2,npoint
	 qu(i) = (u(i) - u(i-1))/v(i-1,4)
	 qu(i-1) = qu(i) - prev
   50	 prev = qu(i)
      qu(npoint) = -qu(npoint)
					return
      end

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
      INTEGER FUNCTION PIPARMQ( ICTXT, ISPEC, NAME, OPTS, N, ILO, IHI,
     $                          LWORKNB )
*
*     Contribution from the Department of Computing Science and HPC2N,
*     Umea University, Sweden
*
*  -- ScaLAPACK auxiliary routine (version 2.0.1) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     Univ. of Colorado Denver and University of California, Berkeley.
*     January, 2012
*
      IMPLICIT NONE
*
*     .. Scalar Arguments ..
      INTEGER            ICTXT, IHI, ILO, ISPEC, LWORKNB, N
      CHARACTER          NAME*( * ), OPTS*( * )
*
*  Purpose
*  =======
*
*       This program sets problem and machine dependent parameters
*       useful for PxHSEQR and its subroutines. It is called whenever
*       PILAENVX is called with 12 <= ISPEC <= 16
*
*  Arguments
*  =========
*
*       ICTXT  (local input) INTEGER
*              On entry,  ICTXT  specifies the BLACS context handle,
*              indicating the global  context of the operation. The
*              context itself is global, but the value of ICTXT is
*              local.
*
*       ISPEC  (global input) INTEGER
*              ISPEC specifies which tunable parameter PIPARMQ should
*              return.
*
*              ISPEC=12: (INMIN)  Matrices of order nmin or less
*                        are sent directly to PxLAHQR, the implicit
*                        double shift QR algorithm.  NMIN must be
*                        at least 11.
*
*              ISPEC=13: (INWIN)  Size of the deflation window.
*                        This is best set greater than or equal to
*                        the number of simultaneous shifts NS.
*                        Larger matrices benefit from larger deflation
*                        windows.
*
*              ISPEC=14: (INIBL) Determines when to stop nibbling and
*                        invest in an (expensive) multi-shift QR sweep.
*                        If the aggressive early deflation subroutine
*                        finds LD converged eigenvalues from an order
*                        NW deflation window and LD.GT.(NW*NIBBLE)/100,
*                        then the next QR sweep is skipped and early
*                        deflation is applied immediately to the
*                        remaining active diagonal block.  Setting
*                        PIPARMQ(ISPEC=14) = 0 causes PxLAQR0 to skip a
*                        multi-shift QR sweep whenever early deflation
*                        finds a converged eigenvalue.  Setting
*                        PIPARMQ(ISPEC=14) greater than or equal to 100
*                        prevents PxLAQR0 from skipping a multi-shift
*                        QR sweep.
*
*              ISPEC=15: (NSHFTS) The number of simultaneous shifts in
*                        a multi-shift QR iteration.
*
*              ISPEC=16: (IACC22) PIPARMQ is set to 1 or 2 with the
*                        following meanings.
*                        1:  During the multi-shift QR sweep,
*                            PxLAQR5 and/or xLAQR6 accumulates reflections
*                            and uses matrix-matrix multiply to update
*                            the far-from-diagonal matrix entries.
*                        2:  During the multi-shift QR sweep.
*                            PxLAQR5 accumulates reflections and takes
*                            advantage of 2-by-2 block structure during
*                            matrix-matrix multiplies.
*
*                        ( IACC22=0 is valid in LAPACK but not here.
*                        Householder reflections are always accumulated
*                        for the performance consideration.
*                        If xTRMM is slower than xGEMM or NB is small,
*                        PIPARMQ(ISPEC=16)=1 may be more efficient than
*                        PIPARMQ(ISPEC=16)=2 despite the greater level of
*                        arithmetic work implied by the latter choice. )
*
*       NAME    (global input) character string
*               Name of the calling subroutine
*
*       OPTS    (global input) character string
*               This is a concatenation of the string arguments to
*               TTQRE.
*
*       N       (global input) integer scalar
*               N is the order of the Hessenberg matrix H.
*
*       ILO     (global input) INTEGER
*       IHI     (global input) INTEGER
*               It is assumed that H is already upper triangular
*               in rows and columns 1:ILO-1 and IHI+1:N.
*
*       LWORKNB   (global input) INTEGER
*               The amount of workspace available or the blockfactor.
*
*  Further Details
*  ===============
*
*       Little is known about how best to choose these parameters.
*       It is possible to use different values of the parameters
*       for each of PCHSEQR, PDHSEQR, PSHSEQR and PZHSEQR.
*
*       It is probably best to choose different parameters for
*       different matrices and different parameters at different
*       times during the iteration, but this has not been fully
*       implemented --- yet.
*
*
*       The best choices of most of the parameters depend
*       in an ill-understood way on the relative execution
*       rate of PxLAQR3 and PxLAQR5 and on the nature of each
*       particular eigenvalue problem.  Experiment may be the
*       only practical way to determine which choices are most
*       effective.
*
*       Following is a list of default values supplied by PIPARMQ.
*       These defaults may be adjusted in order to attain better
*       performance in any particular computational environment.
*
*       PIPARMQ(ISPEC=12) The PxLAQR1 vs PxLAQR0 crossover point.
*                         Default: 220. (Must be at least 11.)
*
*       PIPARMQ(ISPEC=13) Recommended deflation window size.
*                         This depends on ILO, IHI and NS, the
*                         number of simultaneous shifts returned
*                         by PIPARMQ(ISPEC=15).  The default for
*                         (IHI-ILO+1).LE.500 is NS.  The default
*                         for (IHI-ILO+1).GT.500 is 3*NS/2.
*
*       PIPARMQ(ISPEC=14) Nibble crossover point.
*                         The default for the serial case is 14.
*                         The default for the parallel case is
*                         335 * N**(-0.44) * NPROCS.
*
*       PIPARMQ(ISPEC=15) Number of simultaneous shifts, NS.
*                         a multi-shift QR iteration.
*
*                         If IHI-ILO+1 is ...
*
*                         greater than      ...but less    ... the
*                         or equal to ...      than        default is
*
*                                 0               30       NS =    2+
*                                30               60       NS =    4+
*                                60              150       NS =   10
*                               150              590       NS =   **
*                               590             3000       NS =   64
*                              3000             6000       NS =  128
*                              6000            12000       NS =  256
*                             12000            24000       NS =  512
*                             24000            48000       NS = 1024
*                             48000            96000       NS = 2048
*                             96000         INFINITY       NS = 4096
*
*                     (+)  By default matrices of this order are
*                          passed to the implicit double shift routine
*                          PxLAQR1.  See PIPARMQ(ISPEC=12) above. These
*                          values of NS are used only in case of a rare
*                          PxLAQR1 failure.
*
*                     (**) The asterisks (**) indicate an ad-hoc
*                          function increasing from 10 to 64.
*
*       PIPARMQ(ISPEC=16) Select structured matrix multiply.
*                         (See ISPEC=16 above for details.)
*                         Default: 3.
*
*     ================================================================
*     .. Parameters ..
      INTEGER            INMIN, INWIN, INIBL, ISHFTS, IACC22
      PARAMETER          ( INMIN = 12, INWIN = 13, INIBL = 14,
     $                   ISHFTS = 15, IACC22 = 16 )
      INTEGER            NMIN, NMIN2, K22MIN, KACMIN, NIBBLE, KNWSWP
      PARAMETER          ( NMIN = 220, K22MIN = 14, KACMIN = 14,
     $                   NIBBLE = 14, KNWSWP = 500, NMIN2 = 770 )
      REAL               TWO
      PARAMETER          ( TWO = 2.0 )
*     ..
*     .. Local Scalars ..
      INTEGER            NH, NS, MYROW, MYCOL, NPROW, NPCOL, NP
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          LOG, MAX, MOD, NINT, REAL
*     ..
*     .. External functions ..
      INTEGER            ICEIL
      EXTERNAL           ICEIL
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO
*     ..
*     .. Executable Statements ..
      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
*
      IF( ( ISPEC.EQ.ISHFTS ) .OR. ( ISPEC.EQ.INWIN ) .OR.
     $    ( ISPEC.EQ.IACC22 ) ) THEN
*
*        ==== Set the number simultaneous shifts ====
*
         NH = IHI - ILO + 1
         NS = 2
         IF( NH.GE.30 )
     $        NS = 4
         IF( NH.GE.60 )
     $        NS = 10
         IF( NH.GE.150 )
     $        NS = MAX( 10, NH / NINT( LOG( REAL( NH ) ) / LOG( TWO ) ))
         IF( NH.GE.590 )
     $        NS = 64
         IF( NH.GE.3000 )
     $        NS = 128
         IF( NH.GE.6000 )
     $        NS = 256
         IF( NH.GE.12000 )
     $        NS = 512
         IF( NH.GE.24000 )
     $        NS = 1024
         IF( NH.GE.48000 )
     $        NS = 2048
         IF( NH.GE.96000 )
     $        NS = 4096
         IF( NH.GE.192000 )
     $        NS = 8192
         IF( NH.GE.384000 )
     $        NS = 16384
         IF( NH.GE.768000 )
     $        NS = 32768
         IF( NH.GE.1000000 )
     $        NS = ICEIL( NH, 25 )
         NS = MAX( NS, 2*MIN(NPROW,NPCOL) )
         NS = MAX( 2, NS-MOD( NS, 2 ) )
      END IF
*
      IF( ISPEC.EQ.INMIN ) THEN
*
*
*        ===== Submatrices of order smaller than NMIN*min(P_r,P_c)
*        .     get sent to PxLAHQR, the classic ScaLAPACK algorithm.
*        .     This must be at least 11. ====
*
         PIPARMQ = NMIN * MIN( NPROW, NPCOL )
*
      ELSE IF( ISPEC.EQ.INIBL ) THEN
*
*        ==== INIBL: skip a multi-shift QR iteration and
*        .    whenever aggressive early deflation finds
*        .    at least (NIBBLE*(window size)/100) deflations. ====
*
         NP = MIN( NPROW, NPCOL )
         IF( NP.EQ.1 ) THEN
            PIPARMQ = NIBBLE
         ELSE
            NH = IHI - ILO + 1
            PIPARMQ = MIN( 100,
     $           CEILING( 335.0D+0 * NH**(-0.44D+0) * NP ) )
         END IF
*
      ELSE IF( ISPEC.EQ.ISHFTS ) THEN
*
*        ==== NSHFTS: The number of simultaneous shifts =====
*
         PIPARMQ = NS
*
      ELSE IF( ISPEC.EQ.INWIN ) THEN
*
*        ==== NW: deflation window size.  ====
*
         IF( NH.LE.KNWSWP ) THEN
            PIPARMQ = NS
         ELSE
            PIPARMQ = 3*NS / 2
         END IF
*
      ELSE IF( ISPEC.EQ.IACC22 ) THEN
*
*        ==== IACC22: Whether to use 2-by-2 block structure while
*        .     doing it.  A small amount of work could be saved
*        .     by making this choice dependent also upon the
*        .     NH=IHI-ILO+1.
*
         PIPARMQ = 1
c         PIPARMQ = 0
c         IF( NS.GE.KACMIN )
c     $      PIPARMQ = 1
         IF( NS.GE.K22MIN )
     $      PIPARMQ = 2
*
      ELSE
*        ===== invalid value of ispec =====
         PIPARMQ = -1
*
      END IF
*
*     ==== End of PIPARMQ ====
*
      END

max_stars_repo_name: jschueller/scalapack

max_issues_repo_name: jschueller/scalapack

max_forks_repo_name: jschueller/scalapack


        !COMPILER-GENERATED INTERFACE MODULE: Fri Mar  6 15:34:54 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE GRIDPOINTA__genmod
          INTERFACE
            SUBROUTINE GRIDPOINTA
            END SUBROUTINE GRIDPOINTA
          END INTERFACE
        END MODULE GRIDPOINTA__genmod


text:
!
! read_files_from_CMD.f90
! Copyright (C) 2018 Pi-Yueh Chuang <pychuang@gwu.edu>
!
! Distributed under terms of the MIT license.
!

program main
use, intrinsic:: iso_c_binding, only: c_associated, c_null_ptr
use gdal
implicit none

type(GDALDatasetH):: hDS
type(CSLConstList):: list

type(OGRLayerH):: layer
type(CString):: layer_name_c

type(OGRSpatialReferenceH):: srs
type(CString):: srs_wkt
type(OGRErr):: err

type(OGRFeatureDefnH):: layer_defn

type(OGRFieldDefnH):: field_defn
type(CString):: field_name

type(OGRFeatureH):: feature

type(OGRGeometryH):: geometry
type(CString):: geometry_name
type(CString):: geometry_wkt

integer(kind=c_int):: n_layers
integer(kind=c_int):: n_files
integer(kind=c_int):: n_fields

integer(kind=c_int):: i, j, k

character(len=255):: arg
integer(kind=4):: arg_len, arg_stat

! initialize GDAL
call GDALAllRegister()

! get file name/path from CMD arg
call get_command_argument(1, arg, arg_len, arg_stat)
if (arg_stat .gt. 0) stop "Error: failed to get the CMD argument."
if (arg_stat .eq. -1) stop "Error: CMD argument is too long."

! open files
hDS = GDALOpenEx(&
    trim(arg)//char(0), GDAL_OF_VECTOR+GDAL_OF_VERBOSE_ERROR, &
    c_null_ptr, c_null_ptr, c_null_ptr)
if (.not. c_associated(hDS%ptr)) stop

! get number of layers
n_layers = GDALDatasetGetLayerCount(hDS)
print *, "Number of Layers: ", n_layers
print *,

! get file list
list = GDALGetFileList(hDS)
n_files = CSLCount(list)
print *, "Number of Files: ", n_files
print *,
print *, "List of Files:"
n_files = CSLPrint(list, c_null_ptr)
print *,
call CSLDestroy(list)

do i = 0, n_layers-1
    print *, "========================================================================"
    print *, "Layer", i
    print *, "------------------------------------------------------------------------"

    layer = GDALDatasetGetLayer(hDS, i)

    if (.not. c_associated(layer%ptr)) then
        print *, "ERROR: layer", i, "is expected to be valid but not."
        stop
    endif

    layer_name_c = OGR_L_GetName(layer)
    print *, "Name: ", cstr2fstr(layer_name_c)
    print *,

    srs = OGR_L_GetSpatialRef(layer)
    err = OSRExportToPrettyWkt(srs, srs_wkt, 0)
    print *, "SRC (WKT): ", cstr2fstr(srs_wkt)
    print *,
    call CPLFree(srs_wkt)

    layer_defn = OGR_L_GetLayerDefn(layer)
    n_fields = OGR_FD_GetFieldCount(layer_defn)
    print *, "Number of fields: ", n_fields
    print *,

    print *, "List of fields: "

    do j = 0, n_fields-1
        field_defn = OGR_FD_GetFieldDefn(layer_defn, j)
        field_name = OGR_Fld_GetNameRef(field_defn)
        print *, j, ": ", cstr2fstr(field_name)
    enddo
    print *,

    j = 0
    feature%ptr = c_null_ptr
    call OGR_L_ResetReading(layer)
    do while(.true.)
        if(c_associated(feature%ptr)) call OGR_F_Destroy(feature)

        feature = OGR_L_GetNextFeature(layer)
        if(.not. c_associated(feature%ptr)) exit

        print *, "Feature", j, ":"

        geometry = OGR_F_GetGeometryRef(feature)
        write(*, "(A)", advance="no") "     Geometry type: "

        if (c_associated(geometry%ptr)) then
            geometry_name = OGR_G_GetGeometryName(geometry)
            print *, cstr2fstr(geometry_name)

            print *, "    Spatial DIM: ", OGR_G_GetDimension(geometry)
            print *, "    Coordinate DIM: ", OGR_G_CoordinateDimension(geometry)
            print *, "    Point Count: ", OGR_G_GetPointCount(geometry)

            ! get the Wkt of the geometry, but not printing it because it's messy
            err = OGR_G_ExportToWkt(geometry, geometry_wkt)
            call CPLFree(geometry_wkt)
        else
            print *, "N/A"
        endif

        write(*, "(A)", advance="no") "     Field values: "
        do k = 0, n_fields-1
            field_name = OGR_F_GetFieldAsString(feature, k)
            write(*, "(A, A)", advance="no") cstr2fstr(field_name), ", "
        enddo
        print *,
         j = j + 1
    enddo
    call OGR_F_Destroy(feature)
    print *, "========================================================================"
enddo
print *,

call GDALClose(hDS)

end program main

max_stars_repo_name: piyueh/YAFBGDAL

max_issues_repo_name: piyueh/YAFBGDAL

max_forks_repo_name: piyueh/YAFBGDAL


! { dg-do run }
! PR60810 Bogus end-of-file
program readstrlist
  character(len=80), dimension(2) :: ver
  integer :: a, b, c
  a = 1
  b = 2
  c = 3
  ver(1) = '285 383'
  ver(2) = '985'
  read( ver, *) a, b, c
  if (a /= 285 .or. b /= 383 .or. c /= 985) STOP 1
  !write ( *, *) a, b, c
end