text:
      SUBROUTINE TG01DD( COMPZ, L, N, P, A, LDA, E, LDE, C, LDC, Z, LDZ,
     $                   DWORK, LDWORK, INFO )
C
C     SLICOT RELEASE 5.5.
C
C     Copyright (c) 2002-2012 NICONET e.V.
C
C     PURPOSE
C
C     To reduce the descriptor system pair (C,A-lambda E) to the
C     RQ-coordinate form by computing an orthogonal transformation
C     matrix Z such that the transformed descriptor system pair
C     (C*Z,A*Z-lambda E*Z) has the descriptor matrix E*Z in an upper
C     trapezoidal form.
C     The right orthogonal transformations performed to reduce E can
C     be optionally accumulated.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     COMPZ   CHARACTER*1
C             = 'N':  do not compute Z;
C             = 'I':  Z is initialized to the unit matrix, and the
C                     orthogonal matrix Z is returned;
C             = 'U':  Z must contain an orthogonal matrix Z1 on entry,
C                     and the product Z1*Z is returned.
C
C     Input/Output Parameters
C
C     L       (input) INTEGER
C             The number of rows of matrices A and E.  L >= 0.
C
C     N       (input) INTEGER
C             The number of columns of matrices A, E, and C.  N >= 0.
C
C     P       (input) INTEGER
C             The number of rows of matrix C.  P >= 0.
C
C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
C             On entry, the leading L-by-N part of this array must
C             contain the state dynamics matrix A.
C             On exit, the leading L-by-N part of this array contains
C             the transformed matrix A*Z.
C
C     LDA     INTEGER
C             The leading dimension of array A.  LDA >= MAX(1,L).
C
C     E       (input/output) DOUBLE PRECISION array, dimension (LDE,N)
C             On entry, the leading L-by-N part of this array must
C             contain the descriptor matrix E.
C             On exit, the leading L-by-N part of this array contains
C             the transformed matrix E*Z in upper trapezoidal form,
C             i.e.
C
C                      ( E11 )
C                E*Z = (     ) ,  if L >= N ,
C                      (  R  )
C             or
C
C                E*Z = ( 0  R ),  if L < N ,
C
C             where R is an MIN(L,N)-by-MIN(L,N) upper triangular
C             matrix.
C
C     LDE     INTEGER
C             The leading dimension of array E.  LDE >= MAX(1,L).
C
C     C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
C             On entry, the leading P-by-N part of this array must
C             contain the state/output matrix C.
C             On exit, the leading P-by-N part of this array contains
C             the transformed matrix C*Z.
C
C     LDC     INTEGER
C             The leading dimension of array C.  LDC >= MAX(1,P).
C
C     Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
C             If COMPZ = 'N':  Z is not referenced.
C             If COMPZ = 'I':  on entry, Z need not be set;
C                              on exit, the leading N-by-N part of this
C                              array contains the orthogonal matrix Z,
C                              which is the product of Householder
C                              transformations applied to A, E, and C
C                              on the right.
C             If COMPZ = 'U':  on entry, the leading N-by-N part of this
C                              array must contain an orthogonal matrix
C                              Z1;
C                              on exit, the leading N-by-N part of this
C                              array contains the orthogonal matrix
C                              Z1*Z.
C
C     LDZ     INTEGER
C             The leading dimension of array Z.
C             LDZ >= 1,        if COMPZ = 'N';
C             LDZ >= MAX(1,N), if COMPZ = 'U' or 'I'.
C
C     Workspace
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if INFO = 0, DWORK(1) returns the optimal value
C             of LDWORK.
C
C     LDWORK  INTEGER
C             The length of the array DWORK.
C             LDWORK >= MAX(1, MIN(L,N) + MAX(L,N,P)).
C             For optimum performance
C             LWORK >= MAX(1, MIN(L,N) + MAX(L,N,P)*NB),
C             where NB is the optimal blocksize.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value.
C
C     METHOD
C
C     The routine computes the RQ factorization of E to reduce it
C     the upper trapezoidal form.
C
C     The transformations are also applied to the rest of system
C     matrices
C
C         A <- A * Z,  C <- C * Z.
C
C     NUMERICAL ASPECTS
C
C     The algorithm is numerically backward stable and requires
C     0( L*N*N )  floating point operations.
C
C     CONTRIBUTOR
C
C     C. Oara, University "Politehnica" Bucharest.
C     A. Varga, German Aerospace Center, DLR Oberpfaffenhofen.
C     March 1999. Based on the RASP routine RPDSRQ.
C
C     REVISIONS
C
C     July 1999, V. Sima, Research Institute for Informatics, Bucharest.
C
C     KEYWORDS
C
C     Descriptor system, matrix algebra, matrix operations,
C     orthogonal transformation.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ONE, ZERO
      PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0 )
C     .. Scalar Arguments ..
      CHARACTER          COMPZ
      INTEGER            INFO, L, LDA, LDC, LDE, LDWORK, LDZ, N, P
C     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), C( LDC, * ), DWORK( * ),
     $                   E( LDE, * ), Z( LDZ, * )
C     .. Local Scalars ..
      LOGICAL            ILZ
      INTEGER            ICOMPZ, LN, WRKOPT
C     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
C     .. External Subroutines ..
      EXTERNAL           DGERQF, DLASET, DORMRQ, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC          INT, MAX, MIN
C
C     .. Executable Statements ..
C
C     Decode COMPZ.
C
      IF( LSAME( COMPZ, 'N' ) ) THEN
         ILZ = .FALSE.
         ICOMPZ = 1
      ELSE IF( LSAME( COMPZ, 'U' ) ) THEN
         ILZ = .TRUE.
         ICOMPZ = 2
      ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
         ILZ = .TRUE.
         ICOMPZ = 3
      ELSE
         ICOMPZ = 0
      END IF
C
C     Test the input parameters.
C
      INFO = 0
      WRKOPT = MAX( 1, MIN( L, N ) + MAX( L, N, P ) )
      IF( ICOMPZ.EQ.0 ) THEN
         INFO = -1
      ELSE IF( L.LT.0 ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( P.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDA.LT.MAX( 1, L ) ) THEN
         INFO = -6
      ELSE IF( LDE.LT.MAX( 1, L ) ) THEN
         INFO = -8
      ELSE IF( LDC.LT.MAX( 1, P ) ) THEN
         INFO = -10
      ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
         INFO = -12
      ELSE IF( LDWORK.LT.WRKOPT ) THEN
         INFO = -14
      END IF
      IF( INFO .NE. 0 ) THEN
         CALL XERBLA( 'TG01DD', -INFO )
         RETURN
      END IF
C
C     Initialize Q if necessary.
C
      IF( ICOMPZ.EQ.3 )
     $   CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
C
C     Quick return if possible.
C
      IF( L.EQ.0 .OR. N.EQ.0 ) THEN
         DWORK( 1 ) = ONE
         RETURN
      END IF
C
      LN = MIN( L, N )
C
C     Compute the RQ decomposition of E, E = R*Z.
C
C     Workspace: need   MIN(L,N) + L;
C                prefer MIN(L,N) + L*NB.
C
      CALL DGERQF( L, N, E, LDE, DWORK, DWORK( LN+1 ), LDWORK-LN, INFO )
      WRKOPT = MAX( WRKOPT, INT( DWORK( LN+1 ) ) + LN )
C
C     Apply transformation on the rest of matrices.
C
C     A <--  A * Z'.
C     Workspace: need   MIN(L,N) + L;
C                prefer MIN(L,N) + L*NB.
C
      CALL DORMRQ( 'Right', 'Transpose', L, N, LN, E( L-LN+1,1 ), LDE,
     $              DWORK, A, LDA, DWORK( LN+1 ), LDWORK-LN, INFO )
      WRKOPT = MAX( WRKOPT, INT( DWORK( LN+1 ) ) + LN )
C
C     C <-- C * Z'.
C     Workspace: need   MIN(L,N) + P;
C                prefer MIN(L,N) + P*NB.
C
      CALL DORMRQ( 'Right', 'Transpose', P, N, LN, E( L-LN+1,1 ), LDE,
     $              DWORK, C, LDC, DWORK( LN+1 ), LDWORK-LN, INFO )
      WRKOPT = MAX( WRKOPT, INT( DWORK( LN+1 ) ) + LN )
C
C     Z <-- Z1 * Z'.
C     Workspace: need   MIN(L,N) + N;
C                prefer MIN(L,N) + N*NB.
C
      IF( ILZ ) THEN
         CALL DORMRQ( 'Right', 'Transpose', N, N, LN, E( L-LN+1,1 ),
     $                LDE, DWORK, Z, LDZ, DWORK( LN+1 ), LDWORK-LN,
     $                INFO )
         WRKOPT = MAX( WRKOPT, INT( DWORK( LN+1 ) ) + LN )
      END IF
C
C     Set lower triangle of E to zero.
C
      IF( L.LT.N ) THEN
         CALL DLASET( 'Full', L, N-L, ZERO, ZERO, E, LDE )
         IF( L.GE.2 ) CALL DLASET( 'Lower', L-1, L, ZERO, ZERO,
     $                             E( 2, N-L+1 ), LDE )
      ELSE
         IF( N.GE.2 ) CALL DLASET( 'Lower', N-1, N, ZERO, ZERO,
     $                             E( L-N+2, 1 ), LDE )
      END IF
C
      DWORK(1) = WRKOPT
C
      RETURN
C *** Last line of TG01DD ***
      END

max_stars_repo_path: External/SLICOT/TG01DD - Copy.f

max_stars_repo_name: bgin/MissileSimulation

max_issues_repo_path: External/SLICOT/TG01DD - Copy.f

max_issues_repo_name: bgin/MissileSimulation

max_forks_repo_path: External/SLICOT/TG01DD - Copy.f

max_forks_repo_name: bgin/MissileSimulation


text:
MODULE PAR_DiscreteOrdinates
!!#### PURPOSE
!! Contains parameters for the discrete ordinates method, mostly
!! default values.

INTEGER,PARAMETER :: DOR_LEVEL=1,DOR_PRODUCT=2,DOR_GENERAL=3

INTEGER,PARAMETER :: DEFAULT_QuadratureType = DOR_PRODUCT
CHARACTER(10),PARAMETER :: DEFAULT_pQuadrature = "dGL       "
CHARACTER(10),PARAMETER :: DEFAULT_aQuadrature = "GL        "

CHARACTER(10),PARAMETER :: DEFAULT_Quadrature  = "LQ        "

INTEGER,PARAMETER :: DEFAULT_aOrder = 16
INTEGER,PARAMETER :: DEFAULT_pOrder = 16
LOGICAL,PARAMETER :: DEFAULT_UniformZGeometry = .FALSE.
INTEGER,PARAMETER :: DEFAULT_Order  = 128

ENDMODULE

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


text:
*     MB04FD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2020 NICONET e.V.
*
*     .. Parameters ..
      INTEGER            NIN, NOUT
      PARAMETER          ( NIN = 5, NOUT = 6 )
      INTEGER            NMAX
      PARAMETER          ( NMAX = 50 )
      INTEGER            LDA, LDB, LDDE, LDFG, LDQ, LDWORK
      PARAMETER          ( LDA  = NMAX/2, LDB = NMAX/2, LDDE = NMAX/2,
     $                     LDFG = NMAX/2, LDQ = NMAX, LDWORK =
     $                     3*NMAX*NMAX/4 )
*
*     .. Local Scalars ..
      CHARACTER          COMPQ, JOB
      INTEGER            I, INFO, J, N
*
*     .. Local Arrays ..
      DOUBLE PRECISION   A( LDA, NMAX/2 ),  ALPHAI( NMAX/2 ),
     $                   ALPHAR( NMAX/2 ),  B( LDB, NMAX/2 ),
     $                   BETA( NMAX/2 ),  DE( LDDE, NMAX/2+1 ),
     $                   DWORK( LDWORK ), FG( LDFG, NMAX/2+1 ),
     $                   Q( LDQ, NMAX )
      INTEGER            IWORK( NMAX/2+1 )
*
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*
*     .. External Subroutines ..
      EXTERNAL           MB04FD
*
*     .. Intrinsic Functions ..
      INTRINSIC          MOD
*
*     .. Executable statements ..
*
      WRITE( NOUT, FMT = 99999 )
*
*     Skip first line in data file.
*
      READ( NIN, FMT = * )
      READ( NIN, FMT = * ) JOB, COMPQ, N
      READ( NIN, FMT = * ) ( (  A( I, J ), J = 1, N/2 ),   I = 1, N/2 )
      READ( NIN, FMT = * ) ( ( DE( I, J ), J = 1, N/2+1 ), I = 1, N/2 )
      READ( NIN, FMT = * ) ( (  B( I, J ), J = 1, N/2 ),   I = 1, N/2 )
      READ( NIN, FMT = * ) ( ( FG( I, J ), J = 1, N/2+1 ), I = 1, N/2 )
      IF( LSAME( COMPQ, 'U' ) )
     $   READ( NIN, FMT = * ) ( ( Q( I, J ), J = 1, N ), I = 1, N )
      IF( N.LT.0 .OR. N.GT.NMAX .OR. MOD( N, 2 ).NE.0 ) THEN
         WRITE( NOUT, FMT = 99998 ) N
      ELSE
*
*        Test of MB04FD.
*
         CALL MB04FD( JOB, COMPQ, N, A, LDA, DE, LDDE, B, LDB, FG, LDFG,
     $                Q, LDQ, ALPHAR, ALPHAI, BETA, IWORK, DWORK,
     $                LDWORK, INFO )
         IF( INFO.NE.0 ) THEN
            WRITE( NOUT, FMT = 99997 ) INFO
         ELSE
            IF( LSAME( JOB, 'T' ) ) THEN
               WRITE( NOUT, FMT = 99996 )
               DO 10 I = 1, N/2
                  WRITE( NOUT, FMT = 99995 ) ( A( I, J ), J = 1, N/2 )
   10          CONTINUE
            END IF
            WRITE( NOUT, FMT = 99994 )
            DO 20 I = 1, N/2
               WRITE( NOUT, FMT = 99995 ) ( DE( I, J ), J = 1, N/2+1 )
   20       CONTINUE
            IF( LSAME( JOB, 'T' ) ) THEN
               WRITE( NOUT, FMT = 99993 )
               DO 30 I = 1, N/2
                  WRITE( NOUT, FMT = 99995 ) ( B( I, J ), J = 1, N/2 )
   30          CONTINUE
            END IF
            WRITE( NOUT, FMT = 99992 )
            DO 40 I = 1, N/2
               WRITE( NOUT, FMT = 99995 ) ( FG( I, J ), J = 1, N/2+1 )
   40       CONTINUE
            IF( .NOT.LSAME( COMPQ, 'N' ) ) THEN
               WRITE( NOUT, FMT = 99991 )
               DO 50 I = 1, N
                  WRITE( NOUT, FMT = 99995 ) ( Q( I, J ), J = 1, N )
   50          CONTINUE
            END IF
            WRITE( NOUT, FMT = 99990 )
            WRITE( NOUT, FMT = 99995 ) ( ALPHAR(I), I = 1, N/2 )
            WRITE( NOUT, FMT = 99995 ) ( ALPHAI(I), I = 1, N/2 )
            WRITE( NOUT, FMT = 99995 ) (   BETA(I), I = 1, N/2 )
         END IF
      END IF
      STOP
99999 FORMAT ( 'MB04FD EXAMPLE PROGRAM RESULTS', 1X )
99998 FORMAT ( 'N is out of range.', /, 'N = ', I5 )
99997 FORMAT ( 'INFO on exit from MB04FD = ', I2 )
99996 FORMAT (/' The transformed matrix A is' )
99995 FORMAT ( 51( 1X, F8.4 ) )
99994 FORMAT (/' The transformed matrix DE is' )
99993 FORMAT (/' The transformed matrix B is' )
99992 FORMAT (/' The transformed matrix FG is' )
99991 FORMAT (/' The transformed matrix Q is ' )
99990 FORMAT (/' The real, imaginary, and beta parts of eigenvalues are'
     $       )
      END

ext: f

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: bnavigator/SLICOT-Reference

max_forks_repo_name: bnavigator/SLICOT-Reference


text:
        subroutine idz_housemat(n,vn,scal,h)
c
c       fills h with the Householder matrix
c       identity_matrix - scal * vn * adjoint(vn).
c
c       input:
c       n -- size of vn and h, though the indexing of vn goes
c            from 2 to n
c       vn -- entries 2 to n of the vector vn;
c             vn(1) is assumed to be 1
c       scal -- scalar multiplying vn * adjoint(vn)
c
c       output:
c       h -- identity_matrix - scal * vn * adjoint(vn)
c
        implicit none
        save
        integer n,j,k
        real*8 scal
        complex*16 vn(2:*),h(n,n),factor1,factor2
c
c
c       Fill h with the identity matrix.
c
        do j = 1,n
          do k = 1,n
c
            if(j .eq. k) h(k,j) = 1
            if(j .ne. k) h(k,j) = 0
c
          enddo ! k
        enddo ! j
c
c
c       Subtract from h the matrix scal*vn*adjoint(vn).
c
        do j = 1,n
          do k = 1,n
c
            if(j .eq. 1) factor1 = 1
            if(j .ne. 1) factor1 = vn(j)
c
            if(k .eq. 1) factor2 = 1
            if(k .ne. 1) factor2 = conjg(vn(k))
c
            h(k,j) = h(k,j) - scal*factor1*factor2
c
          enddo ! k
        enddo ! j
c
c
        return
        end

max_issues_repo_name: birichie/catboost

max_forks_repo_name: birichie/catboost


text:
C     ( Last modified on 23 Dec 2000 at 22:01:38 )
      SUBROUTINE CVARTY( N, IVARTY )
      INTEGER N
      INTEGER IVARTY( N )
C
C  Determine the type (continuous, 0-1, integer) of each variable
C
C  Nick Gould, for CGT productions.
C  December 1999.
C
      INTEGER            LIWK  , LWK   , LFUVAL, LLOGIC, LCHARA
C
C  ---------------------------------------------------------------------
C
C  Parameters whose value might be changed by the user:
C
C  The following parameters define the sizes of problem
C  dependent arrays. These may be changed by the user to
C  suit a particular problem or system configuration.
C
C  The TOOLS will issue error messages if any of these sizes
C  is too small, telling which parameter to increase.
C
C  ---------------------------------------------------------------------
C
C#{sizing}
C
      INTEGER                IWK( LIWK    )
      LOGICAL              LOGI ( LLOGIC  )
      CHARACTER * 10       CHA  ( LCHARA  )
CS    REAL                   WK ( LWK     )
CD    DOUBLE PRECISION       WK ( LWK     )
CS    REAL              FUVALS  ( LFUVAL  )
CD    DOUBLE PRECISION  FUVALS  ( LFUVAL  )
C
C  ---------------------------------------------------------------------
C
C  End of parameters which might be changed by the user.
C
C  ---------------------------------------------------------------------
C
C  integer variables from the GLOBAL common block.
C
      INTEGER            NG    , NELNUM, NGEL  , NVARS , NNZA  , NGPVLU
      INTEGER            NEPVLU, NG1   , NEL1  , ISTADG, ISTGP , ISTADA
      INTEGER            ISTAEV, ISTEP , ITYPEG, KNDOFC, ITYPEE
      INTEGER            IELING, IELVAR, ICNA  , ISTADH, INTVAR, IVAR
      INTEGER            ICALCF, ITYPEV, IWRK  , A     , B
      INTEGER            U     , GPVALU, EPVALU
      INTEGER            ESCALE, GSCALE, VSCALE, GVALS , XT    , DGRAD
      INTEGER            Q     , WRK   , INTREP, GXEQX , GNAMES, VNAMES
      INTEGER            LO    , CH    , LIWORK, LWORK , NGNG  , FT
      INTEGER            LA, LB, NOBJGR, LU, LELVAR
      INTEGER            LSTAEV, LSTADH, LNTVAR, LCALCF
      INTEGER            LELING, LINTRE, LFT, LGXEQX, LSTADG, LGVALS
      INTEGER            LICNA , LSTADA, LKNDOF, LGPVLU, LEPVLU
      INTEGER            LGSCAL, LESCAL, LVSCAL, LCALCG
C
C  integer variables from the LOCAL common block.
C
      INTEGER            LFXI  , LGXI  , LHXI  , LGGFX , LDX   , LGRJAC
      INTEGER            LQGRAD, LBREAK, LP    , LXCP  , LX0   , LGX0
      INTEGER            LDELTX, LBND  , LWKSTR, LSPTRS, LSELTS, LINDEX
      INTEGER            LSWKSP, LSTAGV, LSTAJC, LIUSED, LFREEC
      INTEGER            LNNONZ, LNONZ2, LSYMMD, LSYMMH
      INTEGER            LSLGRP, LSVGRP, LGCOLJ, LVALJR, LSEND
      INTEGER            LNPTRS, LNELTS, LNNDEX, LNWKSP, LNSTGV
      INTEGER            LNSTJC, LNIUSE, LNFREC, LNNNON, LNNNO2, LNSYMD
      INTEGER            LNSYMH, LNLGRP, LNVGRP, LNGCLJ, LNVLJR, LNQGRD
      INTEGER            LNBRAK, LNP   , LNBND , LNFXI , LNGXI , LNGUVL
      INTEGER            LNHXI , LNHUVL, LNGGFX, LNDX  , LNGRJC, LIWK2
      INTEGER            LWK2  , MAXSIN, NINVAR, MAXSEL
      INTEGER            NTYPE , NSETS , LSTYPE, LSSWTR, LSSIWT, LSIWTR
      INTEGER            LSWTRA, LNTYPE, LNSWTR, LNSIWT, LNIWTR
      INTEGER            LNWTRA, LSISET, LSSVSE, LNISET, LNSVSE
      LOGICAL            ALTRIV, FIRSTG
      COMMON / GLOBAL /  IWK   , WK    , FUVALS, LOGI  ,
     *                   NG    , NELNUM, NGEL  , NVARS , NNZA  , NGPVLU,
     *                   NEPVLU, NG1   , NEL1  , ISTADG, ISTGP , ISTADA,
     *                   ISTAEV, ISTEP , ITYPEG, KNDOFC, ITYPEE,
     *                   IELING, IELVAR, ICNA  , ISTADH, INTVAR, IVAR  ,
     *                   ICALCF, ITYPEV, IWRK  , A     , B     ,
     *                   U     , GPVALU, EPVALU,
     *                   ESCALE, GSCALE, VSCALE, GVALS , XT    , DGRAD ,
     *                   Q     , WRK   , INTREP, GXEQX , GNAMES, VNAMES,
     *                   LO    , CH    , LIWORK, LWORK , NGNG  , FT    ,
     *                   ALTRIV, FIRSTG,
     *                   LA, LB, NOBJGR, LU, LELVAR,
     *                   LSTAEV, LSTADH, LNTVAR, LCALCF,
     *                   LELING, LINTRE, LFT, LGXEQX, LSTADG, LGVALS,
     *                   LICNA , LSTADA, LKNDOF, LGPVLU, LEPVLU,
     *                   LGSCAL, LESCAL, LVSCAL, LCALCG
      COMMON / CHARA /   CHA
      COMMON / LOCAL /   LFXI  , LGXI  , LHXI  , LGGFX , LDX   , LGRJAC,
     *                   LQGRAD, LBREAK, LP    , LXCP  , LX0   , LGX0  ,
     *                   LDELTX, LBND  , LWKSTR, LSPTRS, LSELTS, LINDEX,
     *                   LSWKSP, LSTAGV, LSTAJC, LIUSED, LFREEC,
     *                   LNNONZ, LNONZ2, LSYMMD, LSYMMH,
     *                   LSLGRP, LSVGRP, LGCOLJ, LVALJR, LSEND ,
     *                   LNPTRS, LNELTS, LNNDEX, LNWKSP, LNSTGV,
     *                   LNSTJC, LNIUSE, LNFREC, LNNNON, LNNNO2, LNSYMD,
     *                   LNSYMH, LNLGRP, LNVGRP, LNGCLJ, LNVLJR, LNQGRD,
     *                   LNBRAK, LNP   , LNBND , LNFXI , LNGXI , LNGUVL,
     *                   LNHXI , LNHUVL, LNGGFX, LNDX  , LNGRJC, LIWK2 ,
     *                   LWK2  , MAXSIN, NINVAR, MAXSEL, NTYPE ,
     *                   NSETS , LSTYPE, LSSWTR, LSSIWT, LSIWTR,
     *                   LSWTRA, LNTYPE, LNSWTR, LNSIWT, LNIWTR,
     *                   LNWTRA, LSISET, LSSVSE, LNISET, LNSVSE
      SAVE             / GLOBAL /, / LOCAL /, / CHARA /
C
C  local variables.
C
      INTEGER            I
C
C  Set the type of each variable.
C
      DO 10 I = 1, N
        IVARTY( I ) = IWK( ITYPEV + I )
   10 CONTINUE
      RETURN
C
C  end of CVARTY.
C
      END

max_forks_repo_name: dynamics-of-stellar-systems/dynamite


text:
program main

   implicit none

   integer :: i, j
   real :: s1

   real, dimension(10, 10) :: c1010, d1010

   s1 = 0.125

   d1010 = 0.0
   d1010(1:9:2, 1:9:2) = reshape([((sin(real(i))*sin(real(j)), i=1, 5), j=1, 5)], [5, 5])
   d1010(2:10:2, 2:10:2) = reshape([((cos(real(i))*cos(real(j)), i=1, 5), j=1, 5)], [5, 5])

   print *, char(10), ' d1010: '
   do i = 1, 10

      print "(*(F9.4))", (d1010(i, j), j=1, 10)
   end do

   c1010 = s1
   c1010(1:9:2, 1) = d1010(1, 1:9:2)

   print *, char(10), 'c1010: '
   print "(*(g0))", ((abs(c1010(i, j)), " ", j=1, 10), new_line("A"), i=1, 10)

end program

max_stars_repo_name: WVUHPC/Modern-Fortran

max_issues_repo_name: WVUHPC/Modern-Fortran

max_forks_repo_name: WVUHPC/Modern-Fortran


text:
MODULE Ace
IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)

! BLOCK DATA
! COMMON /parms/ itape, maxit, nterm, span, alpha, big

!     block data
!     common /parms/ itape,maxit,nterm,span,alpha,big

!------------------------------------------------------------------

! these procedure parameters can be changed in the calling routine
! by defining the above labeled common and resetting the values with
! executable statements.

! itape : fortran file number for printer output.
!         (itape.le.0 => no printer output.)
! maxit : maximum number of iterations.
! nterm : number of consecutive iterations for which
!         rsq must change less than delcor for convergence.
! span, alpha : super smoother parameters (see below).
! big : a large representable floating point number.

!------------------------------------------------------------------

! DATA itape, maxit, nterm, span, alpha, big /6, 20, 3, 0.0, 0.0, 1.0E20/
! END

! COMMON /parms/ itape, maxit, nterm, span, alpha, big
INTEGER, SAVE  :: itape = 6, maxit = 20, nterm = 3
REAL, SAVE     :: span = 0.0, alpha = 0.0, big = 1.E20

! BLOCK DATA
! COMMON /spans/ spans(3) /consts/ big, sml, eps

!---------------------------------------------------------------

! this sets the compile time (default) values for various
! internal parameters :

! spans : span values for the three running linear smoothers.
! spans(1) : tweeter span.
! spans(2) : midrange span.
! spans(3) : woofer span.
! (these span values should be changed only with care.)
! big : a large representable floating point number.
! sml : a small number. should be set so that (sml)**(10.0) does
!       not cause floating point underflow.
! eps : used to numerically stabilize slope calculations for
!       running linear fits.

! these parameter values can be changed by declaring the
! relevant labeled common in the main program and resetting
! them with executable statements.

!-----------------------------------------------------------------

! DATA spans, big, sml, eps /0.05, 0.2, 0.5, 1.0E20, 1.0E-7, 1.0E-3/
! END

REAL, SAVE  :: spans(3) = 0.05, sml = 1.0E-7, eps = 1.0E-3

CONTAINS

SUBROUTINE mace(p, n, x, y, w, l, delrsq, ns, tx, ty, rsq, ierr, m, z)

! Code converted using TO_F90 by Alan Miller
! Date: 2003-07-10  Time: 18:07:25

!------------------------------------------------------------------

! estimate multiple optimal transformations for regression and
! correlation by alternating conditional expectation estimates.

! version 3/28/85.

! Breiman and Friedman `Estimating optimal transformations for multiple
! regression and correlation',
! Journal of the American Statistical Association, Vol.80, pp.580-614 (1985).

! coded  and copywrite (c) 1985 by:

!                        jerome h. friedman
!                     department of statistics
!                               and
!                stanford linear accelerator center
!                        stanford university

! all rights reserved.

! N.B. Jerry Friedman has given permission for this version to be made
!      available on this web site.

! input:

!    n : number of observations.
!    p : number of predictor variables for each observation.
!    x(p,n) : predictor data matrix.
!    y(n) : response values for the observations.
!       missing values are signified by a value (response or
!       predictor) greater than or equal to big.
!       (see below - default, big = 1.0e20)
!    w(n) : weights for the observations.
!    l(p+1) : flag for each variable.
!       l(1) through l(p) : predictor variables.
!       l(p+1) : response variable.
!       l(i)=0 => ith variable not to be used.
!       l(i)=1 => ith variable assumes orderable values.
!       l(i)=2 => ith variable assumes circular (periodic) values
!                 in the range (0.0,1.0) with period 1.0.
!       l(i)=3 => ith variable transformation is to be monotone.
!       l(i)=4 => ith variable transformation is to be linear.
!       l(i)=5 => ith variable assumes categorical (unorderable) values.
!   delrsq : termination threshold. iteration stops when
!       rsq changes less than delrsq in nterm
!       consecutive iterations (see below - default, nterm=3).
!   ns : number of eigensolutions (sets of transformations).

! output:

!   tx(n,p,ns) : predictor transformations.
!      tx(j,i,k) = transformed value of ith predictor for jth obs
!                  for kth eigensolution.
!   ty(n,ns) = response transformations.
!      ty(j,k) = transformed response value for jth observation
!                for kth eigensolution.
!   rsq(ns) = fraction of variance(ty<y>)
!                       p
!         explained by sum tx(i)<x(i)>  for each eigensolution.
!                      i=1
!   ierr : error flag.
!      ierr = 0 : no errors detected.
!      ierr > 0 : error detected - see format statements below.

! scratch:

!    m(n,p+1), z(n,12) : internal working storage.

! note: mace uses an iterative procedure for solving the optimization
!    problem. default starting transformations are ty(j,k)=y(j),
!    tx(j,i,k)=x(i,j) : j=1,n, i=1,p, k=1,ns.  Other starting transformations
!    can be specified (if desired) for either the response and/or any of
!    the predictor variables.  This is signaled by negating the
!    corresponding l(i) value and storing the starting transformed values
!    in the corresponding array (ty(j,k), tx(j,i,k)) before calling mace.

!------------------------------------------------------------------

INTEGER, INTENT(IN)      :: p
INTEGER, INTENT(IN)      :: n
REAL, INTENT(IN)         :: x(p,n)
REAL, INTENT(IN)         :: y(n)
REAL, INTENT(IN)         :: w(n)
INTEGER, INTENT(IN OUT)  :: l(1)
REAL, INTENT(IN OUT)     :: delrsq
INTEGER, INTENT(IN)      :: ns
REAL, INTENT(OUT)        :: tx(n,p,ns)
REAL, INTENT(OUT)        :: ty(n,ns)
REAL, INTENT(OUT)        :: rsq(ns)
INTEGER, INTENT(OUT)     :: ierr
INTEGER, INTENT(OUT)     :: m(n,1)
REAL, INTENT(OUT)        :: z(n,12)

INTEGER   :: i, is, ism1, iter, j, js, k, np, nit, nt, pp1
REAL      :: cmn, cmx, ct(10), rsqi
REAL (dp) :: sm, sv, sw, sw1

ierr = 0
pp1 = p + 1
sm = 0.0
sv = sm
sw = sv
sw1 = sw
DO  i = 1, pp1
  IF (l(i) < -5.OR.l(i) > 5) THEN
    ierr = 6
    IF (itape > 0) WRITE (itape,5800) i, l(i)
  END IF
END DO
IF (ierr /= 0) RETURN
IF (l(pp1) == 0) THEN
  ierr = 4
  IF (itape > 0) WRITE (itape,5600) pp1
  RETURN
END IF
np = 0
DO  i = 1, p
  IF (l(i) /= 0) np = np + 1
END DO
IF (np <= 0) THEN
  ierr = 5
  IF (itape > 0) WRITE (itape,5700) p
  RETURN
END IF
DO  j = 1, n
  sw = sw + w(j)
END DO
IF (sw <= 0.0) THEN
  ierr = 1
  IF (itape > 0) WRITE (itape,5300)
  RETURN
END IF
DO  is = 1, ns
  IF (itape > 0) WRITE (itape,5000) is
  DO  j = 1, n
    IF (l(pp1) > 0) ty(j,is) = y(j)
  END DO
  DO  i = 1, p
    IF (l(i) == 0) THEN
      DO  j = 1, n
        tx(j,i,is) = 0.0
      END DO
    ELSE
      IF (l(i) > 0) THEN
        DO  j = 1, n
          tx(j,i,is) = x(i,j)
        END DO
      END IF
      DO  j = 1, n
        IF (tx(j,i,is) < big) THEN
          sm = sm + w(j) * tx(j,i,is)
          sw1 = sw1 + w(j)
        END IF
      END DO
      IF (sw1 <= 0.0) THEN
        DO  j = 1, n
          tx(j,i,is) = 0.0
        END DO
        sm = 0.0
        sw1 = sm
      ELSE
        sm = sm / sw1
        DO  j = 1, n
          IF (tx(j,i,is) < big) THEN
            tx(j,i,is) = tx(j,i,is) - sm
          ELSE
            tx(j,i,is) = 0.0
          END IF
        END DO
        sm = 0.0
        sw1 = sm
      END IF
    END IF
  END DO
  DO  j = 1, n
    IF (ty(j,is) < big) THEN
      sm = sm + w(j) * ty(j,is)
      sw1 = sw1 + w(j)
    END IF
  END DO
  IF (sw1 <= 0.0) THEN
    ierr = 1
    IF (itape > 0) WRITE (itape,5300)
    RETURN
  END IF
  sm = sm / sw1
  DO  j = 1, n
    IF (ty(j,is) < big) THEN
      ty(j,is) = ty(j,is) - sm
    ELSE
      ty(j,is) = 0.0
    END IF
  END DO
  DO  j = 1, n
    sv = sv + w(j) * ty(j,is) ** 2
  END DO
  sv = sv / sw
  IF (sv > 0.0) THEN
    sv = 1.0 / SQRT(sv)
  ELSE
    IF (l(pp1) > 0) THEN
      ierr = 2
      IF (itape > 0) WRITE (itape,5400)
    ELSE
      ierr = 3
      IF (itape > 0) WRITE (itape,5500) is
    END IF
    RETURN
  END IF
  DO  j = 1, n
    ty(j,is) = ty(j,is) * sv
  END DO
  IF (is == 1) THEN
    DO  j = 1, n
      m(j,pp1) = j
      z(j,2) = y(j)
    END DO
    CALL sort(z(1:,2), m(1:,pp1), 1, n)
    DO  i = 1, p
      IF (l(i) /= 0) THEN
        DO  j = 1, n
          m(j,i) = j
          z(j,2) = x(i,j)
        END DO
        CALL sort(z(1:,2), m(1:,i), 1, n)
      END IF
    END DO
  END IF
  CALL scale(p, n, w, sw, ty(1:,is), tx(1:,1:,is), delrsq, p, z(1:,5), z(1:,6))
  rsq(is) = 0.0
  iter = 0
  nterm = MIN(nterm,10)
  nt = 0
  ct(1:nterm) = 100.0

  190   iter = iter + 1
  nit = 0
  200   rsqi = rsq(is)
  nit = nit + 1
  DO  j = 1, n
    z(j,5) = ty(j,is)
    DO  i = 1, p
      IF (l(i) /= 0) z(j,5) = z(j,5) - tx(j,i,is)
    END DO
  END DO
  DO  i = 1, p
    IF (l(i) /= 0) THEN
      DO  j = 1, n
        k = m(j,i)
        z(j,1) = z(k,5) + tx(k,i,is)
        z(j,2) = x(i,k)
        z(j,4) = w(k)
      END DO
      CALL smothr(ABS(l(i)), n, z(1,2), z, z(1,4), z(1,3), z(1,6))
      sm = 0.0
      DO  j = 1, n
        sm = sm + z(j,4) * z(j,3)
      END DO
      sm = sm / sw
      DO  j = 1, n
        z(j,3) = z(j,3) - sm
      END DO
      sv = 0.0
      DO  j = 1, n
        sv = sv + z(j,4) * (z(j,1)-z(j,3)) ** 2
      END DO
      sv = 1.0 - sv / sw
      IF (sv > rsq(is)) THEN
        rsq(is) = sv
        DO  j = 1, n
          k = m(j,i)
          tx(k,i,is) = z(j,3)
          z(k,5) = z(j,1) - z(j,3)
        END DO
      END IF
    END IF
  END DO
  IF (np /= 1 .AND. rsq(is)-rsqi > delrsq .AND. nit < maxit) THEN
    GO TO 200
  END IF
  DO  j = 1, n
    k = m(j,pp1)
    z(j,2) = y(k)
    z(j,4) = w(k)
    z(j,1) = 0.0
    DO  i = 1, p
      IF (l(i) /= 0) z(j,1) = z(j,1) + tx(k,i,is)
    END DO
  END DO
  CALL smothr(ABS(l(pp1)), n, z(1,2), z, z(1,4), z(1,3), z(1,6))
  IF (is > 1) THEN
    ism1 = is - 1
    DO  js = 1, ism1
      sm = 0.0
      DO  j = 1, n
        k = m(j,pp1)
        sm = sm + w(k) * z(j,3) * ty(k,js)
      END DO
      sm = sm / sw
      DO  j = 1, n
        k = m(j,pp1)
        z(j,3) = z(j,3) - sm * ty(k,js)
      END DO
    END DO
  END IF
  sm = 0.0
  sv = sm
  DO  j = 1, n
    k = m(j,pp1)
    sm = sm + w(k) * z(j,3)
    z(k,2) = z(j,1)
  END DO
  sm = sm / sw
  DO  j = 1, n
    z(j,3) = z(j,3) - sm
    sv = sv + z(j,4) * z(j,3) ** 2
  END DO
  sv = sv / sw
  IF (sv > 0.0) THEN
    sv = 1.0 / SQRT(sv)
  ELSE
    ierr = 3
    IF (itape > 0) WRITE (itape,5500) is
    RETURN
  END IF
  DO  j = 1, n
    k = m(j,pp1)
    ty(k,is) = z(j,3) * sv
  END DO
  sv = 0.0
  DO  j = 1, n
    sv = sv + w(j) * (ty(j,is)-z(j,2)) ** 2
  END DO
  rsq(is) = 1.0 - sv / sw
  IF (itape > 0) WRITE (itape,5200) iter, rsq(is)
  nt = MOD(nt,nterm) + 1
  ct(nt) = rsq(is)
  cmn = 100.0
  cmx = -100.0
  DO  i = 1, nterm
    cmn = MIN(cmn,ct(i))
    cmx = MAX(cmx,ct(i))
  END DO
  IF (cmx-cmn > delrsq .AND. iter < maxit) THEN
    GO TO 190
  END IF
  IF (itape > 0) WRITE (itape,5100) is, rsq(is)
END DO
RETURN
5000 FORMAT ('0eigensolution ',i2,':')
5100 FORMAT (' eigensolution ',i2,'   r**2  =  1 - e**2  =',g12.4)
5200 FORMAT ('     iteration ',i2,'   r**2  =  1 - e**2  =',g12.4)
5300 FORMAT (' ierr=1: sum of weights (w) not positive.')
5400 FORMAT (' ierr=2: y has zero variance.')
5500 FORMAT (' ierr=3: ty(.,',i2,') has zero variance.')
5600 FORMAT (' ierr=4: l(',i2,') must be nonzero.')
5700 FORMAT (' ierr=5: at least one l(1)-l(',i2,') must be nonzero.')
5800 FORMAT (' ierr=6: l(',i2,') =',g12.4, ' must be in the range (-5, 5).')
END SUBROUTINE mace

SUBROUTINE model(p, n, y, w, l, tx, ty, f, t, m, z)

!--------------------------------------------------------------------

! computes response predictive  function f for the model yhat = f(t),
! where
!                                        p
!            f(t) = e(y : t),     t =   sum  tx<i> ( x<i> )
!                                       i=1
! using the x transformations tx constructed by subroutine ace.
! if y is a categorical variable (classification) then
!                                -1
!                       f(t) = ty  (t).
! input:

!    p,n,y,w,l : same input as for subroutine ace.
!    tx,ty,m,z : output from subroutine ace.

! output:

!    f(n),t(n) : input for subroutine acemod.

! note: this subroutine must be called before subroutine acemod.

!-------------------------------------------------------------------

INTEGER, INTENT(IN)      :: p
INTEGER, INTENT(IN)      :: n
REAL, INTENT(IN)         :: y(n)
REAL, INTENT(IN)         :: w(n)
INTEGER, INTENT(IN OUT)  :: l(1)
REAL, INTENT(IN)         :: tx(n,p)
REAL, INTENT(IN)         :: ty(n)
REAL, INTENT(OUT)        :: f(n)
REAL, INTENT(OUT)        :: t(n)
INTEGER, INTENT(OUT)     :: m(n,1)
REAL, INTENT(OUT)        :: z(n,12)

INTEGER  :: i, j, j1, j2, k, pp1
REAL     :: s

pp1 = p + 1
IF (ABS(l(pp1)) == 5) THEN
  DO  j = 1, n
    t(j) = ty(j)
    m(j,pp1) = j
  END DO
ELSE
  DO  j = 1, n
    s = 0.0
    DO  i = 1, p
      s = s + tx(j,i)
    END DO
    t(j) = s
    m(j,pp1) = j
  END DO
END IF
CALL sort(t, m(1,pp1), 1, n)
DO  j = 1, n
  k = m(j,pp1)
  z(j,2) = w(k)
  IF (y(k) < big) THEN
    z(j,1) = y(k)
  ELSE
    j1 = j
    j2 = j1
    40     IF (y(m(j1,pp1)) >= big) THEN
      j1 = j1 - 1
      IF (j1 >= 1) THEN
        GO TO 40
      END IF
    END IF
    50     IF (y(m(j2,pp1)) >= big) THEN
      j2 = j2 + 1
      IF (j2 <= n) THEN
        GO TO 50
      END IF
    END IF
    IF (j1 < 1) THEN
      k = j2
    ELSE
      IF (j2 > n) THEN
        k = j1
      ELSE
        IF (t(j)-t(j1) < t(j2)-t(j)) THEN
          k = j1
        ELSE
          k = j2
        END IF
      END IF
    END IF
    z(j,1) = y(m(k,pp1))
    t(j) = t(k)
  END IF
END DO
IF (ABS(l(pp1)) == 5) THEN
  DO  j = 1, n
    f(j) = z(j,1)
  END DO
ELSE
  CALL smothr(1, n, t, z, z(1,2), f, z(1,6))
END IF
RETURN
END SUBROUTINE model

SUBROUTINE acemod(v, p, n, x, l, tx, f, t, m, yhat)
!--------------------------------------------------------------------

! computes response y estimates from the model

!                yhat =  f ( t( v ) )

! using the x transformations tx constructed by subroutine ace and
! the predictor function (f,t) constructed by subroutine model.

! input:

!       v(p) : vector of predictor values.
!    p,n,x,l : same input as for subroutine ace.
!       tx,m : output from subroutine ace.
!        f,t : output from subroutine model.

! output:

!    yhat : estimated response value for v.

! note: this subroutine must not be called before subroutine model.

!-------------------------------------------------------------------

INTEGER, INTENT(IN)      :: p
REAL, INTENT(IN)         :: v(p)
INTEGER, INTENT(IN)      :: n
REAL, INTENT(IN)         :: x(p,n)
INTEGER, INTENT(IN OUT)  :: l(1)
REAL, INTENT(IN)         :: tx(n,p)
REAL, INTENT(IN)         :: f(n)
REAL, INTENT(IN)         :: t(n)
INTEGER, INTENT(IN)      :: m(n,1)
REAL, INTENT(OUT)        :: yhat

INTEGER  :: i, jh, jl, low, high, place
REAL     :: th, vi, xt

th = 0.0
DO  i = 1, p
  IF (l(i) /= 0) THEN
    vi = v(i)
    IF (vi >= big) THEN
      IF (x(i,m(n,i)) >= big) th = th + tx(m(n,i),i)
    ELSE
      IF (vi <= x(i,m(1,i))) THEN
        place = 1
      ELSE
        IF (vi >= x(i,m(n,i))) THEN
          place = n
        ELSE
          low = 0
          high = n + 1
          10   IF (low+1 < high) THEN
            place = (low+high) / 2
            xt = x(i,m(place,i))
            IF (vi == xt) GO TO 20
            IF (vi < xt) THEN
              high = place
              GO TO 10
            END IF
            low = place
            GO TO 10
          END IF
          IF (ABS(l(i)) == 5) CYCLE
          jl = m(low,i)
          jh = m(high,i)
          IF (x(i,jh) >= big) THEN
            th = th + tx(jl,i)
            CYCLE
          END IF
          th = th + tx(jl,i) + (tx(jh,i)-tx(jl,i)) * (vi-x(i,jl))  &
              / (x(i,jh)-x(i,jl))
          CYCLE
        END IF
      END IF
      20  th = th + tx(m(place,i),i)
    END IF
  END IF
END DO
IF (th <= t(1)) THEN
  yhat = f(1)
  RETURN
END IF
IF (th >= t(n)) THEN
  yhat = f(n)
  RETURN
END IF
low = 0
high = n + 1
40 IF (low+1 < high) THEN
  place = (low+high) / 2
  xt = t(place)
  IF (th == xt) THEN
    yhat = f(place)
    RETURN
  END IF
  IF (th < xt) THEN
    high = place
    GO TO 40
  END IF
  low = place
  GO TO 40
END IF
IF (ABS(l(p+1)) == 5) THEN
  IF (th-t(low) <= t(high)-th) THEN
    yhat = f(low)
    GO TO 50
  END IF
  yhat = f(high)
ELSE
  yhat = f(low) + (f(high)-f(low)) * (th-t(low)) / (t(high)- t(low))
END IF
50 RETURN
END SUBROUTINE acemod

SUBROUTINE smothr(l, nn, x, y, w, smo, scr)

INTEGER, INTENT(IN)  :: l
INTEGER, INTENT(IN)  :: nn
REAL, INTENT(IN)     :: x(nn)
REAL, INTENT(IN)     :: y(nn)
REAL, INTENT(IN)     :: w(nn)
REAL, INTENT(OUT)    :: smo(nn)
REAL, INTENT(OUT)    :: scr(nn,7)

REAL (dp) :: sm, sw, a, b, d
INTEGER   :: i, j, j0, n, np1

n = nn
sm = 0.0
sw = sm
10 IF (x(n) >= big) THEN
  sm = sm + w(n) * y(n)
  sw = sw + w(n)
  n = n - 1
  IF (n >= 1) THEN
    GO TO 10
  END IF
END IF
IF (n < nn) THEN
  np1 = n + 1
  sm = sm / sw
  DO  j = np1, nn
    smo(j) = sm
  END DO
END IF
IF (n < 1) RETURN
IF (l >= 5) THEN
  j = 1
  30   j0 = j
  sm = w(j) * y(j)
  sw = w(j)
  IF (j < n) THEN
    40  IF (x(j+1) <= x(j)) THEN
      j = j + 1
      sm = sm + w(j) * y(j)
      sw = sw + w(j)
      IF (j < n) THEN
        GO TO 40
      END IF
    END IF
  END IF
  sm = sm / sw
  DO  i = j0, j
    smo(i) = sm
  END DO
  j = j + 1
  IF (j > n) GO TO 120
  GO TO 30
END IF
IF (l == 4) THEN
  sm = 0.0
  sw = sm
  b = sw
  d = b
  DO  j = 1, n
    sm = sm + w(j) * x(j) * y(j)
    sw = sw + w(j) * x(j) ** 2
    b = b + w(j) * x(j)
    d = d + w(j)
  END DO
  a = sw - (b**2) / d
  IF (a <= 0.0) THEN
    a = 0.0
  ELSE
    a = sm / a
  END IF
  b = b / d
  DO  j = 1, n
    smo(j) = a * (x(j)-b)
  END DO
ELSE
  CALL supsmu(n, x, y, w, l, span, alpha, smo, scr)
  IF (l == 3) THEN
    DO  j = 1, n
      scr(j,1) = smo(j)
      scr(n-j+1,2) = scr(j,1)
    END DO
    CALL montne(scr, n)
    CALL montne(scr(1,2), n)
    sm = 0.0
    sw = sm
    DO  j = 1, n
      sm = sm + (smo(j)-scr(j,1)) ** 2
      sw = sw + (smo(j)-scr(n-j+1,2)) ** 2
    END DO
    IF (sm < sw) THEN
      DO  j = 1, n
        smo(j) = scr(j,1)
      END DO
    ELSE
      DO  j = 1, n
        smo(j) = scr(n-j+1,2)
      END DO
    END IF
  END IF
END IF
120 RETURN
END SUBROUTINE smothr

SUBROUTINE montne(x, n)

INTEGER, INTENT(IN)   :: n
REAL, INTENT(IN OUT)  :: x(n)

INTEGER  :: bb, eb, br, er, bl, el, i
REAL     :: pmn

bb = 0
eb = bb
10 IF (eb < n) THEN
  bb = eb + 1
  eb = bb
  20   IF (eb < n) THEN
    IF (x(bb) == x(eb+1)) THEN
      eb = eb + 1
      GO TO 20
    END IF
  END IF
  30   IF (eb < n) THEN
    IF (x(eb) > x(eb+1)) THEN
      br = eb + 1
      er = br
      40  IF (er < n) THEN
        IF (x(er+1) == x(br)) THEN
          er = er + 1
          GO TO 40
        END IF
      END IF
      pmn = (x(bb)*(eb-bb+1) + x(br)*(er-br+1)) / (er-bb+1)
      eb = er
      DO  i = bb, eb
        x(i) = pmn
      END DO
    END IF
  END IF
  IF (bb <= 1) GO TO 10
  IF (x(bb-1) <= x(bb)) GO TO 10
  bl = bb - 1
  el = bl
  60   IF (bl > 1) THEN
    IF (x(bl-1) == x(el)) THEN
      bl = bl - 1
      GO TO 60
    END IF
  END IF
  pmn = (x(bb)*(eb-bb+1)+x(bl)*(el-bl+1)) / (eb-bl+1)
  bb = bl
  DO  i = bb, eb
    x(i) = pmn
  END DO
  GO TO 30
END IF
RETURN
END SUBROUTINE montne

SUBROUTINE scale(p, n, w, sw, ty, tx, eps, maxit, r, sc)

INTEGER, INTENT(IN)    :: p
INTEGER, INTENT(IN)    :: n
REAL, INTENT(IN)       :: w(n)
REAL (dp), INTENT(IN)  :: sw
REAL, INTENT(IN OUT)   :: ty(n)
REAL, INTENT(IN OUT)   :: tx(n,p)
REAL, INTENT(IN OUT)   :: eps
INTEGER, INTENT(IN)    :: maxit
REAL, INTENT(OUT)      :: r(n)
REAL, INTENT(OUT)      :: sc(p,5)

REAL (dp) :: s, h, t, u, gama, delta, v
INTEGER   :: i, iter, j, nit

DO  i = 1, p
  sc(i,1) = 0.0
END DO
nit = 0
20 nit = nit + 1
DO  i = 1, p
  sc(i,5) = sc(i,1)
END DO
DO  iter = 1, p
  DO  j = 1, n
    s = 0.0
    DO  i = 1, p
      s = s + sc(i,1) * tx(j,i)
    END DO
    r(j) = (ty(j)-s) * w(j)
  END DO
  DO  i = 1, p
    s = 0.0
    DO  j = 1, n
      s = s + r(j) * tx(j,i)
    END DO
    sc(i,2) = -2.0 * s / sw
  END DO
  s = 0.0
  DO  i = 1, p
    s = s + sc(i,2) ** 2
  END DO
  IF (s <= 0.0) EXIT
  IF (iter == 1) THEN
    DO  i = 1, p
      sc(i,3) = -sc(i,2)
    END DO
    h = s
  ELSE
    gama = s / h
    h = s
    DO  i = 1, p
      sc(i,3) = -sc(i,2) + gama * sc(i,4)
    END DO
  END IF
  s = 0.0
  t = s
  DO  j = 1, n
    u = 0.0
    DO  i = 1, p
      u = u + sc(i,3) * tx(j,i)
    END DO
    s = s + u * r(j)
    t = t + w(j) * u ** 2
  END DO
  delta = s / t
  DO  i = 1, p
    sc(i,1) = sc(i,1) + delta * sc(i,3)
    sc(i,4) = sc(i,3)
  END DO
END DO

v = 0.0
DO  i = 1, p
  v = MAX(v, ABS(sc(i,1)-sc(i,5)))
END DO
IF (v >= eps.AND.nit < maxit) THEN
  GO TO 20
END IF
DO  i = 1, p
  DO  j = 1, n
    tx(j,i) = sc(i,1) * tx(j,i)
  END DO
END DO
RETURN
END SUBROUTINE scale

SUBROUTINE sort(v, a, ii, jj)

!  Puts into a the permutation vector which sorts v into
!  increasing order.  Only elements from ii to jj are considered.
!  arrays iu(k) and il(k) permit sorting up to 2**(k+1)-1 elements

!  This is a modification of CACM algorithm #347 by R. C. Singleton,
!  which is a modified Hoare quicksort.

REAL, INTENT(IN OUT)     :: v(:)
INTEGER, INTENT(IN)      :: jj
INTEGER, INTENT(IN OUT)  :: a(jj)
INTEGER, INTENT(IN)      :: ii

INTEGER  :: iu(20), il(20)
INTEGER  :: i, ij, j, k, l, m, t, tt
REAL     :: vt, vtt

m = 1
i = ii
j = jj
10 IF (i >= j) GO TO 60
20 k = i
ij = (j+i) / 2
t = a(ij)
vt = v(ij)
IF (v(i) > vt) THEN
  a(ij) = a(i)
  a(i) = t
  t = a(ij)
  v(ij) = v(i)
  v(i) = vt
  vt = v(ij)
END IF
l = j
IF (v(j) >= vt) GO TO 40
a(ij) = a(j)
a(j) = t
t = a(ij)
v(ij) = v(j)
v(j) = vt
vt = v(ij)
IF (v(i) <= vt) GO TO 40
a(ij) = a(i)
a(i) = t
t = a(ij)
v(ij) = v(i)
v(i) = vt
vt = v(ij)
GO TO 40

30 a(l) = a(k)
a(k) = tt
v(l) = v(k)
v(k) = vtt

40 l = l - 1
IF (v(l) > vt) GO TO 40
tt = a(l)
vtt = v(l)
50 k = k + 1
IF (v(k) < vt) GO TO 50
IF (k <= l) GO TO 30
IF (l-i > j-k) THEN
  il(m) = i
  iu(m) = l
  i = k
  m = m + 1
  GO TO 70
END IF
il(m) = k
iu(m) = j
j = l
m = m + 1
GO TO 70
60 m = m - 1
IF (m == 0) RETURN
i = il(m)
j = iu(m)
70 IF (j-i > 10) GO TO 20
IF (i == ii) GO TO 10
i = i - 1

80 i = i + 1
IF (i == j) GO TO 60
t = a(i+1)
vt = v(i+1)
IF (v(i) <= vt) GO TO 80
k = i

90 a(k+1) = a(k)
v(k+1) = v(k)
k = k - 1
IF (vt < v(k)) GO TO 90
a(k+1) = t
v(k+1) = vt
GO TO 80
END SUBROUTINE sort

SUBROUTINE supsmu(n, x, y, w, iper, span, alpha, smo, sc)

!------------------------------------------------------------------

! super smoother (Friedman, 1984).

! version 10/10/84

! coded  and copywrite <c> 1984 by:

!                        jerome h. friedman
!                     department of statistics
!                               and
!                stanford linear accelerator center
!                        stanford university

! all rights reserved.

! input:
!    n : number of observations (x,y - pairs).
!    x(n) : ordered abscissa values.
!    y(n) : corresponding ordinate (response) values.
!    w(n) : weight for each (x,y) observation.
!    iper : periodic variable flag.
!       iper=1 => x is ordered interval variable.
!       iper=2 => x is a periodic variable with values
!                 in the range (0.0,1.0) and peroid 1.0.
!    span : smoother span (fraction of observations in window).
!           span=0.0 => automatic (variable) span selection.
!    alpha : controles high frequency (small span) penality
!            used with automatic span selection (bass tone control).
!            (alpha.le.0.0 or alpha.gt.10.0 => no effect.)
! output:
!   smo(n) : smoothed ordinate (response) values.
! scratch:
!   sc(n,7) : internal working storage.

! note:
!    for small samples (n < 40) or if there are substantial serial
!    correlations between obserations close in x - value, then
!    a prespecified fixed span smoother (span > 0) should be
!    used. reasonable span values are 0.3 to 0.5.

!------------------------------------------------------------------

INTEGER, INTENT(IN)   :: n
REAL, INTENT(IN)      :: x(n)
REAL, INTENT(IN)      :: y(n)
REAL, INTENT(IN)      :: w(n)
INTEGER, INTENT(IN)   :: iper
REAL, INTENT(IN OUT)  :: span
REAL, INTENT(IN OUT)  :: alpha
REAL, INTENT(OUT)     :: smo(n)
REAL, INTENT(IN OUT)  :: sc(n,7)

! COMMON /spans/ spans(3) /consts/ big, sml, eps
INTEGER  :: i, j, jper
REAL     :: a, f, h(n), resmin, scale, sw, sy, vsmlsq

IF (x(n) <= x(1)) THEN
  sy = 0.0
  sw = sy
  DO  j = 1, n
    sy = sy + w(j) * y(j)
    sw = sw + w(j)
  END DO
  a = 0.0
  IF (sw > 0.0) a = sy / sw
  DO  j = 1, n
    smo(j) = a
  END DO
  RETURN
END IF
i = n / 4
j = 3 * i
scale = x(j) - x(i)
30 IF (scale <= 0.0) THEN
  IF (j < n) j = j + 1
  IF (i > 1) i = i - 1
  scale = x(j) - x(i)
  GO TO 30
END IF
vsmlsq = (eps*scale) ** 2
jper = iper
IF (iper == 2 .AND. (x(1) < 0.0 .OR. x(n) > 1.0)) jper = 1
IF (jper < 1 .OR. jper > 2) jper = 1
IF (span > 0.0) THEN
  CALL smooth(n, x, y, w, span, jper, vsmlsq, smo, sc)
  RETURN
END IF
DO  i = 1, 3
  CALL smooth(n, x, y, w, spans(i), jper, vsmlsq, sc(1:,2*i-1), sc(1:,7))
  CALL smooth(n, x, sc(1:,7), w, spans(2), -jper, vsmlsq, sc(1:,2*i), h)
END DO
DO  j = 1, n
  resmin = big
  DO  i = 1, 3
    IF (sc(j,2*i) < resmin) THEN
      resmin = sc(j,2*i)
      sc(j,7) = spans(i)
    END IF
  END DO
  IF (alpha > 0.0 .AND. alpha <= 10.0 .AND. resmin < sc(j,6) .AND.   &
      resmin > 0.0) sc(j,7) = sc(j,7) + (spans(3)-sc(j,7)) *  &
      MAX(sml,resmin/sc(j,6)) ** (10.0-alpha)
END DO
CALL smooth(n, x, sc(1,7), w, spans(2), -jper, vsmlsq, sc(1,2), h)
DO  j = 1, n
  IF (sc(j,2) <= spans(1)) sc(j,2) = spans(1)
  IF (sc(j,2) >= spans(3)) sc(j,2) = spans(3)
  f = sc(j,2) - spans(2)
  IF (f < 0.0) THEN
    f = -f / (spans(2)-spans(1))
    sc(j,4) = (1.0-f) * sc(j,3) + f * sc(j,1)
  ELSE
    f = f / (spans(3)-spans(2))
    sc(j,4) = (1.0-f) * sc(j,3) + f * sc(j,5)
  END IF
END DO
CALL smooth(n, x, sc(1,4), w, spans(1), -jper, vsmlsq, smo, h)
RETURN
END SUBROUTINE supsmu

SUBROUTINE smooth(n, x, y, w, span, iper, vsmlsq, smo, acvr)

INTEGER, INTENT(IN)   :: n
REAL, INTENT(IN)      :: x(n)
REAL, INTENT(IN)      :: y(n)
REAL, INTENT(IN)      :: w(n)
REAL, INTENT(IN)      :: span
INTEGER, INTENT(IN)   :: iper
REAL, INTENT(IN OUT)  :: vsmlsq
REAL, INTENT(OUT)     :: smo(n)
REAL, INTENT(OUT)     :: acvr(n)

INTEGER    :: i, ibw, in, it, j, j0, jper, out
REAL (dp)  :: wt, fbo, fbw, xm, ym, tmp, var, cvar, a, h, sy, xti, xto

xm = 0.0
ym = xm
var = ym
cvar = var
fbw = cvar
jper = ABS(iper)
ibw = 0.5 * span * n + 0.5
IF (ibw < 2) ibw = 2
it = 2 * ibw + 1
DO  i = 1, it
  j = i
  IF (jper == 2) j = i - ibw - 1
  xti = x(j)
  IF (j < 1) THEN
    j = n + j
    xti = x(j) - 1.0
  END IF
  wt = w(j)
  fbo = fbw
  fbw = fbw + wt
  IF (fbw > 0.0) xm = (fbo*xm+wt*xti) / fbw
  IF (fbw > 0.0) ym = (fbo*ym+wt*y(j)) / fbw
  tmp = 0.0
  IF (fbo > 0.0) tmp = fbw * wt * (xti-xm) / fbo
  var = var + tmp * (xti-xm)
  cvar = cvar + tmp * (y(j)-ym)
END DO
DO  j = 1, n
  out = j - ibw - 1
  in = j + ibw
  IF (.NOT.((jper /= 2) .AND. (out < 1 .OR. in > n))) THEN
    IF (out < 1) THEN
      out = n + out
      xto = x(out) - 1.0
      xti = x(in)
    ELSE
      IF (in > n) THEN
        in = in - n
        xti = x(in) + 1.0
        xto = x(out)
      ELSE
        xto = x(out)
        xti = x(in)
      END IF
    END IF
    wt = w(out)
    fbo = fbw
    fbw = fbw - wt
    tmp = 0.0
    IF (fbw > 0.0) tmp = fbo * wt * (xto-xm) / fbw
    var = var - tmp * (xto-xm)
    cvar = cvar - tmp * (y(out)-ym)
    IF (fbw > 0.0) xm = (fbo*xm-wt*xto) / fbw
    IF (fbw > 0.0) ym = (fbo*ym-wt*y(out)) / fbw
    wt = w(in)
    fbo = fbw
    fbw = fbw + wt
    IF (fbw > 0.0) xm = (fbo*xm+wt*xti) / fbw
    IF (fbw > 0.0) ym = (fbo*ym+wt*y(in)) / fbw
    tmp = 0.0
    IF (fbo > 0.0) tmp = fbw * wt * (xti-xm) / fbo
    var = var + tmp * (xti-xm)
    cvar = cvar + tmp * (y(in)-ym)
  END IF
  a = 0.0
  IF (var > vsmlsq) a = cvar / var
  smo(j) = a * (x(j)-xm) + ym
  IF (iper > 0) THEN
    h = 0.0
    IF (fbw > 0.0) h = 1.0 / fbw
    IF (var > vsmlsq) h = h + (x(j)-xm) ** 2 / var
    acvr(j) = 0.0
    a = 1.0 - w(j) * h
    IF (a > 0.0) THEN
      acvr(j) = ABS(y(j)-smo(j)) / a
    ELSE
      IF (j > 1) THEN
        acvr(j) = acvr(j-1)
      END IF
    END IF
  END IF
END DO
j = 1
30 j0 = j
sy = smo(j) * w(j)
fbw = w(j)
IF (j < n) THEN
  40   IF (x(j+1) <= x(j)) THEN
    j = j + 1
    sy = sy + w(j) * smo(j)
    fbw = fbw + w(j)
    IF (j < n) THEN
      GO TO 40
    END IF
  END IF
END IF
IF (j > j0) THEN
  a = 0.0
  IF (fbw > 0.0) a = sy / fbw
  DO  i = j0, j
    smo(i) = a
  END DO
END IF
j = j + 1
IF (j <= n) THEN
  GO TO 30
END IF
RETURN
END SUBROUTINE smooth

END MODULE Ace

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


!========================================================================
!
!                   S P E C F E M 2 D  Version 7 . 0
!                   --------------------------------
!
!     Main historical authors: Dimitri Komatitsch and Jeroen Tromp
!                        Princeton University, USA
!                and CNRS / University of Marseille, France
!                 (there are currently many more authors!)
! (c) Princeton University and CNRS / University of Marseille, April 2014
!
! This software is a computer program whose purpose is to solve
! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
! using a spectral-element method (SEM).
!
! This software is governed by the CeCILL license under French law and
! abiding by the rules of distribution of free software. You can use,
! modify and/or redistribute the software under the terms of the CeCILL
! license as circulated by CEA, CNRS and Inria at the following URL
! "http://www.cecill.info".
!
! As a counterpart to the access to the source code and rights to copy,
! modify and redistribute granted by the license, users are provided only
! with a limited warranty and the software's author, the holder of the
! economic rights, and the successive licensors have only limited
! liability.
!
! In this respect, the user's attention is drawn to the risks associated
! with loading, using, modifying and/or developing or reproducing the
! software by the user in light of its specific status of free software,
! that may mean that it is complicated to manipulate, and that also
! therefore means that it is reserved for developers and experienced
! professionals having in-depth computer knowledge. Users are therefore
! encouraged to load and test the software's suitability as regards their
! requirements in conditions enabling the security of their systems and/or
! data to be ensured and, more generally, to use and operate it in the
! same conditions as regards security.
!
! The full text of the license is available in file "LICENSE".
!
!========================================================================

  subroutine compute_arrays_source(ispec_selected_source,xi_source,gamma_source,sourcearray, &
             Mxx,Mzz,Mxz,xix,xiz,gammax,gammaz,xigll,zigll,nspec)

  use specfem_par, only : AXISYM,is_on_the_axis,xiglj

  implicit none

  include "constants.h"

  integer ispec_selected_source
  integer nspec

  double precision xi_source,gamma_source
  double precision Mxx,Mzz,Mxz

  real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz

  real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLZ) :: sourcearray

  double precision xixd,xizd,gammaxd,gammazd

! Gauss-Lobatto-Legendre points of integration and weights
  double precision, dimension(NGLLX) :: xigll
  double precision, dimension(NGLLZ) :: zigll

! source arrays
  double precision, dimension(NGLLX,NGLLZ) :: G11,G13,G31,G33
  double precision, dimension(NGLLX) :: hxis,hpxis
  double precision, dimension(NGLLZ) :: hgammas,hpgammas

  integer k,m
  integer ir,iv

! calculate G_ij for general source location
! the source does not necessarily correspond to a Gauss-Lobatto point
  do m=1,NGLLZ
      do k=1,NGLLX

        xixd    = xix(k,m,ispec_selected_source)
        xizd    = xiz(k,m,ispec_selected_source)
        gammaxd = gammax(k,m,ispec_selected_source)
        gammazd = gammaz(k,m,ispec_selected_source)

        G11(k,m) = Mxx*xixd+Mxz*xizd
        G13(k,m) = Mxx*gammaxd+Mxz*gammazd
        G31(k,m) = Mxz*xixd+Mzz*xizd
        G33(k,m) = Mxz*gammaxd+Mzz*gammazd

!!!!        G21(k,m) = Mxy*xixd+Myz*xizd
!!!!        G23(k,m) = Mxy*gammaxd+Myz*gammazd

      enddo
  enddo

! compute Lagrange polynomials at the source location

  if (AXISYM) then
    if(is_on_the_axis(ispec_selected_source)) then ! TODO verify if we have to add : .and. is_proc_source(i)
      call lagrange_any(xi_source,NGLJ,xiglj,hxis,hpxis)
    else
      call lagrange_any(xi_source,NGLLX,xigll,hxis,hpxis)
    endif
  else
    call lagrange_any(xi_source,NGLLX,xigll,hxis,hpxis)
  endif

  call lagrange_any(gamma_source,NGLLZ,zigll,hgammas,hpgammas)

! calculate source array
  do m=1,NGLLZ
    do k=1,NGLLX

      sourcearray(:,k,m) = ZERO

      do iv=1,NGLLZ
        do ir=1,NGLLX

          sourcearray(1,k,m) = sourcearray(1,k,m) + hxis(ir)*hgammas(iv) &
                                 *(G11(ir,iv)*hpxis(k)*hgammas(m) &
                                 +G13(ir,iv)*hxis(k)*hpgammas(m))

!        sourcearray(2,k,m) = sourcearray(2,k,m) + hxis(ir)*hgammas(iv) &
!                               *(G21(ir,iv)*hpxis(k)*hgammas(m) &
!                               +G23(ir,iv)*hxis(k)*hpgammas(m))

          sourcearray(2,k,m) = sourcearray(2,k,m) + hxis(ir)*hgammas(iv) &
                                 *(G31(ir,iv)*hpxis(k)*hgammas(m) &
                                 +G33(ir,iv)*hxis(k)*hpgammas(m))

        enddo
      enddo

    enddo
  enddo

  end subroutine compute_arrays_source

! ------------------------------------------------------------------------------------------------------

 subroutine compute_arrays_adj_source(xi_rec,gamma_rec)

 use constants
 use specfem_par

 implicit none

  double precision xi_rec, gamma_rec
  logical :: st !YY

  integer icomp, itime, i, k
  double precision :: junk
  character(len=3) :: comp(3)

  adj_sourcearray(:,:,:,:) = 0.

  if(seismotype == 1 .or. seismotype == 2 .or. seismotype == 3) then

  comp = (/"BXX","BXY","BXZ"/)

  do icomp = 1,3

   filename = 'SEM/'//trim(adj_source_file) // '.'// comp(icomp) // '.adj'

   inquire(file=filename,exist=st) !YY
   ! if file exist  !YY
   if(st) then   !YY

    open(unit = IIN, file = trim(filename), iostat = ios)
    if (ios /= 0) call exit_MPI(' file '//trim(filename)//'does not exist')

    do itime = 1, NSTEP
      read(IIN,*) junk, adj_src_s(itime,icomp)
    enddo
    close(IIN)

   else     !YY
     adj_src_s(:,icomp) = 0.0 !YY
   endif !YY

  enddo

  source_adjointe(irec_local,:,2) = adj_src_s(:,3)

  else if (seismotype == 4 ) then

    filename = 'SEM/'//trim(adj_source_file) // '.PRE.adj'
    open(unit = IIN, file = trim(filename), iostat = ios)
    if (ios /= 0) call exit_MPI(' file '//trim(filename)//'does not exist')

    do itime = 1, NSTEP
      read(IIN,*) junk, adj_src_s(itime,1)
    enddo
    close(IIN)

  else if (seismotype == 6 ) then

    filename = 'SEM/'//trim(adj_source_file) // '.POT.adj'
    open(unit = IIN, file = trim(filename), iostat = ios)
    if (ios /= 0) call exit_MPI(' file '//trim(filename)//'does not exist')

    do itime = 1, NSTEP
      read(IIN,*) junk, adj_src_s(itime,1)
    enddo
    close(IIN)

  endif

  source_adjointe(irec_local,:,1) = adj_src_s(:,1)

  if (AXISYM) then
    if(is_on_the_axis(ispec_selected_rec(irec_local))) then ! TODO verify irec_local...
      call lagrange_any(xi_rec,NGLJ,xiglj,hxir,hpxir)
    else
      call lagrange_any(xi_rec,NGLLX,xigll,hxir,hpxir)
    endif
  else
    call lagrange_any(xi_rec,NGLLX,xigll,hxir,hpxir)
  endif

  call lagrange_any(gamma_rec,NGLLZ,zigll,hgammar,hpgammar)

  do k = 1, NGLLZ
      do i = 1, NGLLX
        adj_sourcearray(:,:,i,k) = sngl(hxir(i) * hgammar(k)) * adj_src_s(:,:)
      enddo
  enddo

end subroutine compute_arrays_adj_source


text:
      SUBROUTINE CFE2MY (TPOSE,Y,X,FILE,BUF)
C*******
C     CFE2MY FORMS THE COMPLEX DOUBLE PRECISION MATRIX
C     PRODUCT X = M*Y FOR THE COMPLEX FEER METHOD
C*******
C     DEFINITION OF INPUT AND OUTPUT PARAMETERS
C*******
C     TPOSE    = .FALSE. --- USE MATRIX M
C              = .TRUE.  --- USE MATRIX M-TRANSPOSE
C     Y        = INPUT  VECTOR
C     X        = OUTPUT VECTOR
C     FILE     = INPUT MATRIX CONTROL BLOCK FOR THE
C                REQUIRED MATRIX
C     BUF      = INPUT REQUIRED GINO BUFFER
C*******
      DOUBLE PRECISION  X(1)     ,Y(1)     ,DA
      INTEGER           FILE(7)  ,BUF(1)   ,EOL      ,DIAG
      LOGICAL           TPOSE(1)
      COMMON  /NAMES /  RD       ,RDREW    ,WRT      ,WRTREW
     2                 ,REW      ,NOREW    ,EOFNRW   ,RSP
     3                 ,RDP      ,CSP      ,CDP      ,SQR
     4                 ,RECT     ,DIAG     ,LOWTRI   ,UPRTRI
     5                 ,SYM      ,ROW      ,IDENTY
      COMMON  /ZNTPKX/  DA(2)    ,II       ,EOL
      NCOL2 = FILE(2)+FILE(2)
      IF (FILE(4).EQ.IDENTY) GO TO 50
      CALL GOPEN (FILE(1),BUF(1),RDREW)
      DO 10  I = 1,NCOL2
   10 X(I) = 0.D0
      IF (FILE(4).EQ.DIAG) GO TO 40
      IF (TPOSE(1)) GO TO 31
C*******
C     GENERAL MATRIX*VECTOR PRODUCT
C*******
      DO 30  I = 1,NCOL2,2
      J = I+1
      IF (Y(I).EQ.0.D0.AND.Y(J).EQ.0.D0) GO TO 25
      CALL INTPK(*30,FILE(1),0,CDP,0)
   22 CALL ZNTPKI
      JJ = II+II
      II = JJ-1
      X(II) = X(II) + DA(1)*Y(I) - DA(2)*Y(J)
      X(JJ) = X(JJ) + DA(1)*Y(J) + DA(2)*Y(I)
      IF (EOL) 30,22,30
   25 CALL SKPREC (FILE(1),1)
   30 CONTINUE
      GO TO 80
C*******
C     GENERAL MATRIX-TRANSPOSE*VECTOR PRODUCT
C*******
   31 DO 36  I = 1,NCOL2,2
      J = I+1
      CALL INTPK(*36,FILE(1),0,CDP,0)
   32 CALL ZNTPKI
      JJ = II+II
      II = JJ-1
      X(I) = X(I) + DA(1)*Y(II) - DA(2)*Y(JJ)
      X(J) = X(J) + DA(1)*Y(JJ) + DA(2)*Y(II)
      IF (EOL) 36,32,36
   36 CONTINUE
      GO TO 80
C*******
C     MATRIX IS DIAGONAL
C*******
  40  CALL INTPK(*80,FILE(1),0,CDP,0)
   45 CALL ZNTPKI
      JJ = II+II
      II = JJ-1
      X(II) = DA(1)*Y(II) - DA(2)*Y(JJ)
      X(JJ) = DA(1)*Y(JJ) + DA(2)*Y(II)
      IF (EOL) 80,45,80
C*******
C     MATRIX IS IDENTITY
C*******
   50 DO 55  I = 1,NCOL2
   55 X(I) = Y(I)
      GO TO 90
   80 CALL CLOSE (FILE(1),REW)
   90 RETURN
      END


text:
program main
  use AbstractOperators
  !use SimpleVectors
  use DVectors
  use ExampleOperators

  use KrylovSolver
  implicit none
  ! problem vars
  type(Laplacian) :: laplacian_1d
  type(DiagonalMatrix) :: inverse_diagonal_laplacian
  type(dvector) :: sol
  type(dvector) :: rhs
  integer :: n
  integer :: max_iterations=400
  double precision :: tolerance=1.0d-5,h,sign

  ! local vars
  integer :: ierr,i,m,j
  type(dvector) :: res
  n=100
  h=1.0d0/(n-1)

  !! init matrix
  call laplacian_1d%init(n)

  !! allocate solution and rhs and set the latter
  call sol%init(n)
  call rhs%init(n)
  call res%init(n)

  do i=1,n
     call rhs%set(i,0.d0)
  end do
  call rhs%set(1,1.0d0)
  call rhs%set(n,1.0d0)

  do i=1,n
     call sol%set(i,0.d0)
  end do

  !! init diagonal and set values with inverse of
  !! diagonal of the Laplacian
  call inverse_diagonal_laplacian%init(n)
  inverse_diagonal_laplacian%coefficients=1.0d0/(2.0d0)

  write(*,*) norm_dvector(rhs)
  call pcg(laplacian_1d,rhs,sol,ierr,&
       tolerance=1.0d-3,&
       max_iterations=400)!,&
  !preconditi1.0d0r=inverse_diagonal_laplacian)

  !! check results

  call laplacian_1d%apply(sol,res,ierr)
  call axpby_dvector(-1.0d0,rhs,1.0d0,res)
  write(*,*) 'Res=',norm_dvector(res)

end program main

max_stars_repo_name: enricofacca/KrylovAPI

max_issues_repo_name: enricofacca/KrylovAPI

max_forks_repo_name: enricofacca/KrylovAPI


text:
      SUBROUTINE MC01OD( K, REZ, IMZ, REP, IMP, DWORK, INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To compute the coefficients of a complex polynomial P(x) from its
C     zeros.
C
C     ARGUMENTS
C
C     Input/Output Parameters
C
C     K       (input) INTEGER
C             The number of zeros (and hence the degree) of P(x).
C             K >= 0.
C
C     REZ     (input) DOUBLE PRECISION array, dimension (K)
C     IMZ     (input) DOUBLE PRECISION array, dimension (K)
C             The real and imaginary parts of the i-th zero of P(x)
C             must be stored in REZ(i) and IMZ(i), respectively, where
C             i = 1, 2, ..., K. The zeros may be supplied in any order.
C
C     REP     (output) DOUBLE PRECISION array, dimension (K+1)
C     IMP     (output) DOUBLE PRECISION array, dimension (K+1)
C             These arrays contain the real and imaginary parts,
C             respectively, of the coefficients of P(x) in increasing
C             powers of x. If K = 0, then REP(1) is set to one and
C             IMP(1) is set to zero.
C
C     Workspace
C
C     DWORK   DOUBLE PRECISION array, dimension (2*K+2)
C             If K = 0, this array is not referenced.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value.
C
C     METHOD
C
C     The routine computes the coefficients of the complex K-th degree
C     polynomial P(x) as
C
C        P(x) = (x - r(1)) * (x - r(2)) * ... * (x - r(K))
C
C     where r(i) = (REZ(i),IMZ(i)), using real arithmetic.
C
C     NUMERICAL ASPECTS
C
C     None.
C
C     CONTRIBUTORS
C
C     Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Mar. 1997.
C     Supersedes Release 2.0 routine MC01CD by Alan Brown and
C     A.J. Geurts.
C
C     REVISIONS
C
C     V. Sima, May 2002.
C
C     KEYWORDS
C
C     Elementary polynomial operations, polynomial operations.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, ONE
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
C     .. Scalar Arguments ..
      INTEGER           INFO, K
C     .. Array Arguments ..
      DOUBLE PRECISION  DWORK(*), IMP(*), IMZ(*), REP(*), REZ(*)
C     .. Local Scalars ..
      INTEGER           I, K2
      DOUBLE PRECISION  U, V
C     .. External Subroutines ..
      EXTERNAL          DAXPY, DCOPY, XERBLA
C     .. Executable Statements ..
C
C     Test the input scalar arguments.
C
      IF( K.LT.0 ) THEN
         INFO = -1
C
C        Error return.
C
         CALL XERBLA( 'MC01OD', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      INFO   = 0
      REP(1) = ONE
      IMP(1) = ZERO
      IF ( K.EQ.0 )
     $   RETURN
C
      K2 = K + 2
C
      DO 20 I = 1, K
         U = REZ(I)
         V = IMZ(I)
         DWORK(1)  = ZERO
         DWORK(K2) = ZERO
         CALL DCOPY( I, REP, 1, DWORK(2), 1 )
         CALL DCOPY( I, IMP, 1, DWORK(K2+1), 1 )
C
         IF ( U.NE.ZERO ) THEN
            CALL DAXPY( I, -U, REP, 1, DWORK, 1 )
            CALL DAXPY( I, -U, IMP, 1, DWORK(K2), 1 )
         END IF
C
         IF ( V.NE.ZERO ) THEN
            CALL DAXPY( I,  V, IMP, 1, DWORK, 1 )
            CALL DAXPY( I, -V, REP, 1, DWORK(K2), 1 )
         END IF
C
         CALL DCOPY( I+1, DWORK, 1, REP, 1 )
         CALL DCOPY( I+1, DWORK(K2), 1, IMP, 1 )
   20 CONTINUE
C
      RETURN
C *** Last line of MC01OD ***
      END

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr