text:
      SUBROUTINE MB02JD( JOB, K, L, M, N, P, S, TC, LDTC, TR, LDTR, Q,
     $                   LDQ, R, LDR, DWORK, LDWORK, INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To compute a lower triangular matrix R and a matrix Q with
C     Q^T Q = I such that
C                                    T
C                           T  =  Q R ,
C
C     where T is a K*M-by-L*N block Toeplitz matrix with blocks of size
C     (K,L). The first column of T will be denoted by TC and the first
C     row by TR. It is assumed that the first MIN(M*K, N*L) columns of T
C     have full rank.
C
C     By subsequent calls of this routine the factors Q and R can be
C     computed block column by block column.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     JOB     CHARACTER*1
C             Specifies the output of the routine as follows:
C             = 'Q':  computes Q and R;
C             = 'R':  only computes R.
C
C     Input/Output Parameters
C
C     K       (input)  INTEGER
C             The number of rows in one block of T.  K >= 0.
C
C     L       (input)  INTEGER
C             The number of columns in one block of T.  L >= 0.
C
C     M       (input)  INTEGER
C             The number of blocks in one block column of T.  M >= 0.
C
C     N       (input)  INTEGER
C             The number of blocks in one block row of T.  N >= 0.
C
C     P       (input)  INTEGER
C             The number of previously computed block columns of R.
C             P*L < MIN( M*K,N*L ) + L and P >= 0.
C
C     S       (input)  INTEGER
C             The number of block columns of R to compute.
C             (P+S)*L < MIN( M*K,N*L ) + L and S >= 0.
C
C     TC      (input) DOUBLE PRECISION array, dimension (LDTC, L)
C             On entry, if P = 0, the leading M*K-by-L part of this
C             array must contain the first block column of T.
C
C     LDTC    INTEGER
C             The leading dimension of the array TC.
C             LDTC >= MAX(1,M*K).
C
C     TR      (input)  DOUBLE PRECISION array, dimension (LDTR,(N-1)*L)
C             On entry, if P = 0, the leading K-by-(N-1)*L part of this
C             array must contain the first block row of T without the
C             leading K-by-L block.
C
C     LDTR    INTEGER
C             The leading dimension of the array TR.
C             LDTR >= MAX(1,K).
C
C     Q       (input/output)  DOUBLE PRECISION array, dimension
C                             (LDQ,MIN( S*L, MIN( M*K,N*L )-P*L ))
C             On entry, if JOB = 'Q'  and  P > 0, the leading M*K-by-L
C             part of this array must contain the last block column of Q
C             from a previous call of this routine.
C             On exit, if JOB = 'Q'  and  INFO = 0, the leading
C             M*K-by-MIN( S*L, MIN( M*K,N*L )-P*L ) part of this array
C             contains the P-th to (P+S)-th block columns of the factor
C             Q.
C
C     LDQ     INTEGER
C             The leading dimension of the array Q.
C             LDQ >= MAX(1,M*K), if JOB = 'Q';
C             LDQ >= 1,          if JOB = 'R'.
C
C     R       (input/output)  DOUBLE PRECISION array, dimension
C                             (LDR,MIN( S*L, MIN( M*K,N*L )-P*L ))
C             On entry, if P > 0, the leading (N-P+1)*L-by-L
C             part of this array must contain the nozero part of the
C             last block column of R from a previous call of this
C             routine.
C             One exit, if INFO = 0, the leading
C             MIN( N, N-P+1 )*L-by-MIN( S*L, MIN( M*K,N*L )-P*L )
C             part of this array contains the nonzero parts of the P-th
C             to (P+S)-th block columns of the lower triangular
C             factor R.
C             Note that elements in the strictly upper triangular part
C             will not be referenced.
C
C     LDR     INTEGER
C             The leading dimension of the array R.
C             LDR >= MAX( 1, MIN( N, N-P+1 )*L )
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             On exit, if INFO = -17,  DWORK(1) returns the minimum
C             value of LDWORK.
C             If JOB = 'Q', the first 1 + ( (N-1)*L + M*K )*( 2*K + L )
C             elements of DWORK should be preserved during successive
C             calls of the routine.
C             If JOB = 'R', the first 1 + (N-1)*L*( 2*K + L ) elements
C             of DWORK should be preserved during successive calls of
C             the routine.
C
C     LDWORK  INTEGER
C             The length of the array DWORK.
C             JOB = 'Q':
C                LDWORK >= 1 + ( M*K + ( N - 1 )*L )*( L + 2*K ) + 6*L
C                            + MAX( M*K,( N - MAX( 1,P )*L ) );
C             JOB = 'R':
C                If P = 0,
C                   LDWORK >= MAX( 1 + ( N - 1 )*L*( L + 2*K ) + 6*L
C                                    + (N-1)*L, M*K*( L + 1 ) + L );
C                If P > 0,
C                   LDWORK >= 1 + (N-1)*L*( L + 2*K ) + 6*L + (N-P)*L.
C             For optimum performance LDWORK should be larger.
C
C             If LDWORK = -1, then a workspace query is assumed;
C             the routine only calculates the optimal size of the
C             DWORK array, returns this value as the first entry of
C             the DWORK array, and no error message related to LDWORK
C             is issued by XERBLA.
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             = 1:  the full rank condition for the first MIN(M*K, N*L)
C                   columns of T is (numerically) violated.
C
C     METHOD
C
C     Block Householder transformations and modified hyperbolic
C     rotations are used in the Schur algorithm [1], [2].
C
C     REFERENCES
C
C     [1] Kailath, T. and Sayed, A.
C         Fast Reliable Algorithms for Matrices with Structure.
C         SIAM Publications, Philadelphia, 1999.
C
C     [2] Kressner, D. and Van Dooren, P.
C         Factorizations and linear system solvers for matrices with
C         Toeplitz structure.
C         SLICOT Working Note 2000-2, 2000.
C
C     NUMERICAL ASPECTS
C
C     The implemented method yields a factor R which has comparable
C     accuracy with the Cholesky factor of T^T * T. Q is implicitly
C     computed from the formula Q = T * inv(R^T R) * R, i.e., for ill
C     conditioned problems this factor is of very limited value.
C                                 2
C     The algorithm requires 0(K*L *M*N) floating point operations.
C
C     CONTRIBUTOR
C
C     D. Kressner, Technical Univ. Berlin, Germany, May 2001.
C
C     REVISIONS
C
C     V. Sima, Research Institute for Informatics, Bucharest, June 2001.
C     D. Kressner, Technical Univ. Berlin, Germany, July 2002.
C     V. Sima, Research Institute for Informatics, Bucharest, Mar. 2004,
C     Aug. 2011.
C
C     KEYWORDS
C
C     Elementary matrix operations, Householder transformation, matrix
C     operations, Toeplitz matrix.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, ONE
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
C     .. Scalar Arguments ..
      CHARACTER         JOB
      INTEGER           INFO, K, L, LDQ, LDR, LDTC, LDTR, LDWORK,
     $                  M, N, P, S
C     .. Array Arguments ..
      DOUBLE PRECISION  DWORK(LDWORK), Q(LDQ,*), R(LDR,*), TC(LDTC,*),
     $                  TR(LDTR,*)
C     .. Local Scalars ..
      INTEGER           COLR, I, IERR, KK, LEN, NB, NBMIN, PDQ, PDW,
     $                  PNQ, PNR, PRE, PT, RNK, SHFR, STPS, WRKMIN,
     $                  WRKOPT
      LOGICAL           COMPQ, LQUERY
C     .. Local Arrays ..
      INTEGER           IPVT(1)
C     .. External Functions ..
      LOGICAL           LSAME
      INTEGER           ILAENV
      EXTERNAL          ILAENV, LSAME
C     .. External Subroutines ..
      EXTERNAL          DGELQF, DGEQRF, DLACPY, DLASET, DORGQR, MA02AD,
     $                  MB02CU, MB02CV, MB02KD, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC         DBLE, INT, MAX, MIN
C
C     .. Executable Statements ..
C
C     Decode the scalar input parameters.
C
      INFO  = 0
      COMPQ = LSAME( JOB, 'Q' )
      IF ( COMPQ ) THEN
         WRKMIN = 1 + ( M*K + ( N - 1 )*L )*( L + 2*K ) + 6*L
     $              + MAX( M*K, ( N - MAX( 1, P ) )*L )
      ELSE
         WRKMIN = 1 + ( N - 1 )*L*( L + 2*K ) + 6*L
     $              + ( N - MAX( P, 1 ) )*L
         IF ( P.EQ.0 ) THEN
            WRKMIN = MAX( WRKMIN, M*K*( L + 1 ) + L )
         END IF
      END IF
C
C     Check the scalar input parameters.
C
      IF ( .NOT.( COMPQ .OR. LSAME( JOB, 'R' ) ) ) THEN
         INFO = -1
      ELSE IF ( K.LT.0 ) THEN
         INFO = -2
      ELSE IF ( L.LT.0 ) THEN
         INFO = -3
      ELSE IF ( M.LT.0 ) THEN
         INFO = -4
      ELSE IF ( N.LT.0 ) THEN
         INFO = -5
      ELSE IF ( P*L.GE.MIN( M*K, N*L ) + L .OR. P.LT.0 ) THEN
         INFO = -6
      ELSE IF ( ( P + S )*L.GE.MIN( M*K, N*L ) + L .OR. S.LT.0  ) THEN
         INFO = -7
      ELSE IF ( LDTC.LT.MAX( 1, M*K ) ) THEN
         INFO = -9
      ELSE IF ( LDTR.LT.MAX( 1, K ) ) THEN
         INFO = -11
      ELSE IF ( LDQ.LT.1 .OR. ( COMPQ .AND. LDQ.LT.M*K ) ) THEN
         INFO = -13
      ELSE IF ( LDR.LT.MAX( 1, MIN( N, N - P + 1 )*L ) ) THEN
         INFO = -15
      ELSE
         LQUERY = LDWORK.EQ.-1
         IF ( LQUERY ) THEN
            WRKOPT = 1
            IF ( M*K.LE.L ) THEN
               PDW = M*K*L + M*K
               CALL DGEQRF( M*K, L, DWORK, M*K, DWORK, DWORK, -1, IERR )
               WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) + PDW )
               CALL DORGQR( M*K, M*K, M*K, DWORK, M*K, DWORK, DWORK, -1,
     $                      IERR )
               WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) + PDW )
               IF ( N.GT.1 ) THEN
                  PDW = M*K*M*K
                  CALL MB02KD( 'Row', 'Transpose', K, L, M, N-1, M*K,
     $                         ONE, ZERO, TC, LDTC, TR, LDTR, DWORK,
     $                         M*K, R, LDR, DWORK, -1, IERR )
                  WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) + PDW )
               END IF
C
            ELSE IF ( P.EQ.0 ) THEN
               IF ( COMPQ ) THEN
                  CALL DGEQRF( M*K, L, Q, LDQ, DWORK, DWORK, -1, IERR )
                  WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) + L )
                  CALL DORGQR( M*K, L, L, Q, LDQ, DWORK, DWORK, -1,
     $                        IERR )
                  WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) + L )
                  IF ( N.GT.1 ) THEN
                     CALL MB02KD( 'Row', 'Transpose', K, L, M, N-1, L,
     $                            ONE, ZERO, TC, LDTC, TR, LDTR, Q, LDQ,
     $                            R, LDR, DWORK, LDWORK, IERR )
                     WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) )
                  END IF
               ELSE
                  PDW = M*K*L
                  CALL DGEQRF( M*K, L, DWORK, M*K, DWORK, DWORK, -1,
     $                        IERR )
                  WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) + PDW + L )
                  CALL DORGQR( M*K, L, L, DWORK, M*K, DWORK, DWORK, -1,
     $                         IERR )
                  WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) + PDW + L )
                  IF ( N.GT.1 ) THEN
                     CALL MB02KD( 'Row', 'Transpose', K, L, M, N-1, L,
     $                            ONE, ZERO, TC, LDTC, TR, LDTR, DWORK,
     $                            M*K, R, LDR, DWORK, -1, IERR )
                     WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) + PDW )
                  END IF
               END IF
               PRE = 1
            END IF
C
            IF ( M*K.GT.L .OR. P.GT.0 .OR.
     $         ( P.EQ.0 .AND. N.GT.1 ) ) THEN
               PRE = MAX( 1, P )
               IF ( COMPQ ) THEN
                  PDW = ( 2*K + L )*( ( N - 1 )*L + M*K ) + 1
                  LEN = MAX( L + M*K, ( N - PRE + 1 )*L )
               ELSE
                  PDW = ( 2*K + L )*( N - 1 )*L + 1
                  LEN = ( N - PRE + 1 )*L
               END IF
C
               CALL DGELQF( LEN, L, DWORK, MAX( 1, LEN ), DWORK, DWORK,
     $                      -1, IERR )
               WRKOPT = MAX( WRKOPT, PDW + 6*L + INT( DWORK(1) ) )
            END IF
         ELSE IF ( LDWORK.LT.WRKMIN ) THEN
            DWORK(1) = DBLE( WRKMIN )
            INFO = -17
         END IF
      END IF
C
C     Return if there were illegal values.
C
      IF ( INFO .NE. 0 ) THEN
         CALL XERBLA( 'MB02JD', -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         DWORK(1) = WRKOPT
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF ( MIN( M, N, K*L, S ) .EQ.0 ) THEN
         DWORK(1) = ONE
         RETURN
      END IF
C
C     Catch M*K <= L.
C
      WRKOPT = 1
      IF ( M*K.LE.L ) THEN
         CALL DLACPY( 'All', M*K, L, TC, LDTC, DWORK, M*K )
         PDW = M*K*L + 1
         CALL DGEQRF( M*K, L, DWORK, M*K, DWORK(PDW),
     $                DWORK(PDW+M*K), LDWORK-PDW-M*K+1, IERR )
         WRKOPT = MAX( WRKOPT, INT( DWORK(PDW+M*K) ) + PDW + M*K - 1 )
         CALL MA02AD( 'Upper part', M*K, L, DWORK, M*K, R, LDR )
         CALL DORGQR( M*K, M*K, M*K, DWORK, M*K, DWORK(PDW),
     $                DWORK(PDW+M*K), LDWORK-PDW-M*K+1, IERR )
         WRKOPT = MAX( WRKOPT, INT( DWORK(PDW+M*K) ) + PDW + M*K - 1 )
         IF ( COMPQ ) THEN
            CALL DLACPY( 'All', M*K, M*K, DWORK, M*K, Q, LDQ )
         END IF
         PDW = M*K*M*K + 1
         IF ( N.GT.1 ) THEN
            CALL MB02KD( 'Row', 'Transpose', K, L, M, N-1, M*K, ONE,
     $                   ZERO, TC, LDTC, TR, LDTR, DWORK, M*K, R(L+1,1),
     $                   LDR, DWORK(PDW), LDWORK-PDW+1, IERR )
            WRKOPT = MAX( WRKOPT, INT( DWORK(PDW) ) + PDW - 1 )
         END IF
         DWORK(1) = DBLE( WRKOPT )
         RETURN
      END IF
C
C     Compute the generator if P = 0.
C
      IF ( P.EQ.0 ) THEN
C
C        1st column of the generator.
C
         IF ( COMPQ ) THEN
            CALL DLACPY( 'All', M*K, L, TC, LDTC, Q, LDQ )
            CALL DGEQRF( M*K, L, Q, LDQ, DWORK, DWORK(L+1),
     $                   LDWORK-L, IERR )
            WRKOPT = MAX( WRKOPT, INT( DWORK(L+1) ) + L )
            CALL MA02AD( 'Upper part', L, L, Q, LDQ, R, LDR )
            CALL DORGQR( M*K, L, L, Q, LDQ, DWORK, DWORK(L+1), LDWORK-L,
     $                  IERR )
            WRKOPT = MAX( WRKOPT, INT( DWORK(L+1) ) + L )
            IF ( N.GT.1 ) THEN
               CALL MB02KD( 'Row', 'Transpose', K, L, M, N-1, L, ONE,
     $                      ZERO, TC, LDTC, TR, LDTR, Q, LDQ, R(L+1,1),
     $                      LDR, DWORK, LDWORK, IERR )
               WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) )
            END IF
         ELSE
            PDW = M*K*L + 1
            CALL DLACPY( 'All', M*K, L, TC, LDTC, DWORK, M*K )
            CALL DGEQRF( M*K, L, DWORK, M*K, DWORK(PDW), DWORK(PDW+L),
     $                   LDWORK-PDW-L+1, IERR )
            WRKOPT = MAX( WRKOPT, INT( DWORK(PDW+L) ) + PDW + L - 1 )
            CALL MA02AD( 'Upper part', L, L, DWORK, M*K, R, LDR )
            CALL DORGQR( M*K, L, L, DWORK, M*K, DWORK(PDW),
     $                   DWORK(PDW+L), LDWORK-PDW-L+1, IERR )
            WRKOPT = MAX( WRKOPT, INT( DWORK(PDW+L) ) + PDW + L - 1 )
            IF ( N.GT.1 ) THEN
               CALL MB02KD( 'Row', 'Transpose', K, L, M, N-1, L, ONE,
     $                      ZERO, TC, LDTC, TR, LDTR, DWORK, M*K,
     $                      R(L+1,1), LDR, DWORK(PDW), LDWORK-PDW+1,
     $                      IERR )
               WRKOPT = MAX( WRKOPT, INT( DWORK(PDW) ) + PDW - 1 )
            END IF
         END IF
C
C        Quick return if N = 1.
C
         IF ( N.EQ.1 ) THEN
            DWORK(1) = DBLE( WRKOPT )
            RETURN
         END IF
C
C        2nd column of the generator.
C
         PNR = ( N - 1 )*L*K + 2
         CALL MA02AD( 'All', K, (N-1)*L, TR, LDTR, DWORK(2), (N-1)*L )
C
C        3rd and 4th column of the generator.
C
         CALL DLACPY( 'All', (N-1)*L, L, R(L+1,1), LDR, DWORK(PNR),
     $                (N-1)*L )
         PT  = ( M - 1 )*K + 1
         PDW = PNR + ( N - 1 )*L*L
C
         DO 10  I = 1, MIN( M, N-1 )
            CALL MA02AD( 'All', K, L, TC(PT,1), LDTC, DWORK(PDW),
     $                   (N-1)*L )
            PT  = PT  - K
            PDW = PDW + L
   10    CONTINUE
C
         PT = 1
C
         DO 20  I = M + 1, N - 1
            CALL MA02AD( 'All', K, L, TR(1,PT), LDTR, DWORK(PDW),
     $                   (N-1)*L )
            PT  = PT  + L
            PDW = PDW + L
   20    CONTINUE
C
         IF ( COMPQ ) THEN
            PDQ = ( 2*K + L )*( N - 1 )*L + 2
            PDW = ( 2*K + L )*( ( N - 1 )*L + M*K ) + 2
            PNQ = PDQ + M*K*K
            CALL DLASET( 'All', K, K, ZERO, ONE, DWORK(PDQ), M*K )
            CALL DLASET( 'All', (M-1)*K, K, ZERO, ZERO, DWORK(PDQ+K),
     $                   M*K )
            CALL DLACPY( 'All', M*K, L, Q, LDQ, DWORK(PNQ), M*K )
            CALL DLASET( 'All', M*K, K, ZERO, ZERO, DWORK(PNQ+M*L*K),
     $                   M*K )
         ELSE
            PDW = ( 2*K + L )*( N - 1 )*L + 2
         END IF
         PRE  = 1
         STPS = S - 1
      ELSE
C
C        Set workspace pointers.
C
         PNR = ( N - 1 )*L*K + 2
         IF ( COMPQ ) THEN
            PDQ = ( 2*K + L )*( N - 1 )*L + 2
            PDW = ( 2*K + L )*( ( N - 1 )*L + M*K ) + 2
            PNQ = PDQ + M*K*K
         ELSE
            PDW = ( 2*K + L )*( N - 1 )*L + 2
         END IF
         PRE  = P
         STPS = S
      END IF
C
C     Determine suitable size for the block Housholder reflectors.
C
      IF ( COMPQ ) THEN
         LEN = MAX( L + M*K, ( N - PRE + 1 )*L )
      ELSE
         LEN = ( N - PRE + 1 )*L
      END IF
      NB    = MIN( ( LDWORK - PDW - 6*L + 1 ) / LEN, L )
      NBMIN = MAX( 2, ILAENV( 2, 'DGELQF', ' ', LEN, L, -1, -1 ) )
      IF ( NB.LT.NBMIN )  NB = 0
      COLR = L + 1
C
C     Generator reduction process.
C
      LEN  = ( N - PRE )*L
      SHFR = ( PRE - 1 )*L
      DO 30  I = PRE, PRE + STPS - 1
C
C        IF M*K < N*L the last block might have less than L columns.
C
         KK = MIN( L, M*K - I*L )
         CALL DLACPY( 'Lower', LEN, KK, R(COLR-L,COLR-L), LDR,
     $                R(COLR,COLR), LDR )
         CALL MB02CU( 'Column', KK, KK+K, L+K, NB, R(COLR,COLR), LDR,
     $                 DWORK(SHFR+2), (N-1)*L, DWORK(PNR+SHFR), (N-1)*L,
     $                 RNK, IPVT, DWORK(PDW), ZERO, DWORK(PDW+6*L),
     $                 LDWORK-PDW-6*L+1, IERR )
         IF ( IERR.NE.0 )  THEN
C
C           Error return:  The rank condition is (numerically) not
C                          satisfied.
C
            INFO = 1
            RETURN
         END IF
         IF ( LEN.GT.KK ) THEN
            CALL MB02CV( 'Column', 'NoStructure', KK, LEN-KK, KK+K, L+K,
     $                   NB, -1, R(COLR,COLR), LDR, DWORK(SHFR+2),
     $                   (N-1)*L, DWORK(PNR+SHFR), (N-1)*L,
     $                   R(COLR+KK,COLR), LDR, DWORK(SHFR+KK+2),
     $                   (N-1)*L, DWORK(PNR+SHFR+KK), (N-1)*L,
     $                   DWORK(PDW), DWORK(PDW+6*L), LDWORK-PDW-6*L+1,
     $                   IERR )
         END IF
         IF ( COMPQ ) THEN
            CALL DLASET( 'All', K, KK, ZERO, ZERO, Q(1,COLR), LDQ )
            IF ( M.GT.1 ) THEN
               CALL DLACPY( 'All', (M-1)*K, KK, Q(1,COLR-L), LDQ,
     $                      Q(K+1,COLR), LDQ )
            END IF
            CALL MB02CV( 'Column', 'NoStructure', KK, M*K, KK+K, L+K,
     $                   NB, -1, R(COLR,COLR), LDR, DWORK(SHFR+2),
     $                   (N-1)*L, DWORK(PNR+SHFR), (N-1)*L, Q(1,COLR),
     $                   LDQ, DWORK(PDQ), M*K, DWORK(PNQ), M*K,
     $                   DWORK(PDW), DWORK(PDW+6*L), LDWORK-PDW-6*L+1,
     $                   IERR )
         END IF
         LEN  = LEN  - L
         COLR = COLR + L
         SHFR = SHFR + L
   30 CONTINUE
C
      DWORK(1) = DBLE( WRKOPT )
      RETURN
C
C *** Last line of MB02JD ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


!
! NLTE_Predictor_IO
!
! Module containing routines to read and write Binary format
! NLTE_Predictor data files.
!
!
! CREATION HISTORY:
!       Written by:  Paul van Delst, 15-Mar-2011
!                    paul.vandelst@noaa.gov
!

MODULE NLTE_Predictor_IO

  ! ------------------
  ! Environment set up
  ! ------------------
  ! Module use
  USE Type_Kinds           , ONLY: Long, Double
  USE File_Utility         , ONLY: File_Open, File_Exists
  USE Message_Handler      , ONLY: SUCCESS, FAILURE, WARNING, INFORMATION, Display_Message
  USE Binary_File_Utility  , ONLY: Open_Binary_File
  USE NLTE_Predictor_Define, ONLY: NLTE_Predictor_type        , &
                                   NLTE_Predictor_Destroy     , &
                                   NLTE_Predictor_ValidRelease, &
                                   NLTE_Predictor_Info
  ! Disable implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  PUBLIC :: NLTE_Predictor_InquireFile
  PUBLIC :: NLTE_Predictor_ReadFile
  PUBLIC :: NLTE_Predictor_WriteFile
  PUBLIC :: NLTE_Predictor_IOVersion

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PRIVATE, PARAMETER :: MODULE_VERSION_ID = &
    '$Id$'
  CHARACTER(*), PARAMETER :: WRITE_ERROR_STATUS = 'DELETE'
  ! Default message length
  INTEGER, PARAMETER :: ML = 256

CONTAINS

!################################################################################
!################################################################################
!##                                                                            ##
!##                         ## PUBLIC MODULE ROUTINES ##                       ##
!##                                                                            ##
!################################################################################
!################################################################################

!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       NLTE_Predictor_InquireFile
!
! PURPOSE:
!       Function to inquire a Binary format NLTE_Predictor file.
!
! CALLING SEQUENCE:
!       Error_Status = NLTE_Predictor_InquireFile( &
!                        Filename               , &
!                        n_Profiles = n_Profiles, &
!                        Release    = Release   , &
!                        Version    = Version     )
!
! INPUTS:
!       Filename:      Character string specifying the name of the NLTE
!                      predictor data file to inquire.
!                      UNITS:      N/A
!                      TYPE:       CHARACTER(*)
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(IN)
!
! OPTIONAL OUTPUTS:
!       n_Profiles:    The number of profiles for which there is NLTE
!                      predictor information in the data file.
!                      UNITS:      N/A
!                      TYPE:       INTEGER
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: OPTIONAL, INTENT(OUT)
!
!       Release:       The data/file release number. Used to check
!                      for data/software mismatch.
!                      UNITS:      N/A
!                      TYPE:       INTEGER
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       Version:       The data/file version number. Used for
!                      purposes only in identifying the dataset for
!                      a particular release.
!                      UNITS:      N/A
!                      TYPE:       INTEGER
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(OUT), OPTIONAL
!
! FUNCTION RESULT:
!       Error_Status:  The return value is an integer defining the error
!                      status. The error codes are defined in the
!                      Message_Handler module.
!                      If == SUCCESS the file inquire was successful
!                         == FAILURE an unrecoverable error occurred.
!                      UNITS:      N/A
!                      TYPE:       INTEGER
!                      DIMENSION:  Scalar
!
!:sdoc-:
!------------------------------------------------------------------------------

  FUNCTION NLTE_Predictor_InquireFile( &
    Filename  , &  ! Input
    n_Profiles, &  ! Optional output
    Release   , &  ! Optional Output
    Version   ) &  ! Optional Output
  RESULT( err_stat )
    ! Arguments
    CHARACTER(*),           INTENT(IN)  :: Filename
    INTEGER     , OPTIONAL, INTENT(OUT) :: n_Profiles
    INTEGER     , OPTIONAL, INTENT(OUT) :: Release
    INTEGER     , OPTIONAL, INTENT(OUT) :: Version
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'NLTE_Predictor_InquireFile'
    ! Function variables
    CHARACTER(ML) :: msg
    INTEGER :: io_stat
    INTEGER :: fid
    INTEGER :: rel, ver, m

    ! Setup
    err_stat = SUCCESS
    ! ...Check that the file exists
    IF ( .NOT. File_Exists( Filename ) ) THEN
      msg = 'File '//TRIM(Filename)//' not found.'
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Open the file
    err_stat = Open_Binary_File( Filename, fid )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error opening '//TRIM(Filename)
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Read the release and version
    READ( fid, IOSTAT=io_stat ) rel, ver
    IF ( io_stat /= 0 ) THEN
      WRITE( msg,'("Error reading Release/Version. IOSTAT = ",i0)' ) io_stat
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Read the number of profiles
    READ( fid, IOSTAT=io_stat ) m
    IF ( io_stat /= 0 ) THEN
      WRITE( msg,'("Error reading dimensions from ",a,". IOSTAT = ",i0)' ) TRIM(Filename), io_stat
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Close the file
    CLOSE( fid, IOSTAT=io_stat )
    IF ( io_stat /= 0 ) THEN
      WRITE( msg,'("Error closing ",a,". IOSTAT = ",i0)' ) TRIM(Filename), io_stat
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Assign the return arguments
    IF ( PRESENT(n_Profiles) ) n_Profiles = m
    IF ( PRESENT(Release   ) ) Release    = rel
    IF ( PRESENT(Version   ) ) Version    = ver

  CONTAINS

    SUBROUTINE Inquire_CleanUp()
      ! Close file if necessary
      IF ( File_Open(fid) ) THEN
        CLOSE( fid, IOSTAT=io_stat )
        IF ( io_stat /= 0 ) &
          msg = TRIM(msg)//'; Error closing input file during error cleanup'
      END IF
      ! Set error status and print error message
      err_stat = FAILURE
      CALL Display_Message( ROUTINE_NAME, msg, err_stat )
    END SUBROUTINE Inquire_CleanUp

  END FUNCTION NLTE_Predictor_InquireFile

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       NLTE_Predictor_ReadFile
!
! PURPOSE:
!       Function to read NLTE_Predictor object files in Binary format.
!
! CALLING SEQUENCE:
!       Error_Status = NLTE_Predictor_ReadFile( &
!                        Filename      , &
!                        NLTE_Predictor, &
!                        Quiet      = Quiet     , &
!                        n_Profiles = n_Profiles  )
!
! INPUTS:
!       Filename:       Character string specifying the name of a
!                       NLTE_Predictor format data file to read.
!                       UNITS:      N/A
!                       TYPE:       CHARACTER(*)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       NLTE_Predictor: NLTE_Predictor object containing the NLTE correction
!                       algorithm predictor data.
!                       UNITS:      N/A
!                       TYPE:       NLTE_Predictor_type
!                       DIMENSION:  Rank-1
!                       ATTRIBUTES: INTENT(OUT)
!
! OPTIONAL INPUTS:
!       Quiet:          Set this logical argument to suppress INFORMATION
!                       messages being printed to stdout
!                       If == .FALSE., INFORMATION messages are OUTPUT [DEFAULT].
!                          == .TRUE.,  INFORMATION messages are SUPPRESSED.
!                       If not specified, default is .FALSE.
!                       UNITS:      N/A
!                       TYPE:       LOGICAL
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN), OPTIONAL
!
! OPTIONAL OUTPUTS:
!       n_Profiles:     The number of profiles for which data was read.
!                       UNITS:      N/A
!                       TYPE:       INTEGER
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: OPTIONAL, INTENT(OUT)
!
! FUNCTION RESULT:
!       Error_Status:   The return value is an integer defining the error status.
!                       The error codes are defined in the Message_Handler module.
!                       If == SUCCESS, the file read was successful
!                          == FAILURE, an unrecoverable error occurred.
!                       UNITS:      N/A
!                       TYPE:       INTEGER
!                       DIMENSION:  Scalar
!
!:sdoc-:
!------------------------------------------------------------------------------

  FUNCTION NLTE_Predictor_ReadFile( &
    Filename      , &  ! Input
    NLTE_Predictor, &  ! Output
    Quiet         , &  ! Optional input
    n_Profiles    , &  ! Optional output
    Debug         ) &  ! Optional input (Debug output control)
  RESULT( err_stat )
    ! Arguments
    CHARACTER(*),              INTENT(IN)  :: Filename
    TYPE(NLTE_Predictor_type), INTENT(OUT) :: NLTE_Predictor(:)
    LOGICAL,         OPTIONAL, INTENT(IN)  :: Quiet
    INTEGER,         OPTIONAL, INTENT(OUT) :: n_Profiles
    LOGICAL,         OPTIONAL, INTENT(IN)  :: Debug
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'NLTE_Predictor_ReadFile'
    ! Function variables
    CHARACTER(ML) :: msg
    LOGICAL :: noisy
    INTEGER :: io_stat
    INTEGER :: fid
    INTEGER :: n_file_profiles
    INTEGER :: m, n_input_profiles
    TYPE(NLTE_Predictor_type) :: dummy

    ! Setup
    err_stat = SUCCESS
    ! ...Check Quiet argument
    noisy = .TRUE.
    IF ( PRESENT(Quiet) ) noisy = .NOT. Quiet
    ! ...Override Quiet settings if debug set.
    IF ( PRESENT(Debug) ) THEN
      IF ( Debug ) noisy = .TRUE.
    END IF

    ! Open the file if it exists
    IF ( File_Exists( Filename ) ) THEN
      err_stat = Open_Binary_File( Filename, fid )
      IF ( err_Stat /= SUCCESS ) THEN
        msg = 'Error opening '//TRIM(Filename)
        CALL Read_CleanUp(); RETURN
      END IF
    ELSE
      msg = 'File '//TRIM(Filename)//' not found.'
      CALL Read_CleanUp(); RETURN
    END IF

    ! Read and check the release and version
    READ( fid, IOSTAT=io_stat ) dummy%Release, dummy%Version
    IF ( io_stat /= 0 ) THEN
      WRITE( msg,'("Error reading Release/Version. IOSTAT = ",i0)' ) io_stat
      CALL Read_Cleanup(); RETURN
    END IF
    IF ( .NOT. NLTE_Predictor_ValidRelease( dummy ) ) THEN
      msg = 'NLTE_Predictor Release check failed.'
      CALL Read_Cleanup(); RETURN
    END IF

    ! Read the dimensions
    READ( fid, IOSTAT=io_stat ) n_file_profiles
    IF ( io_stat /= 0 ) THEN
      WRITE( msg,'("Error reading profile dimension from ",a,". IOSTAT = ",i0)' ) &
             TRIM(Filename), io_stat
      CALL Read_Cleanup(); RETURN
    END IF
    ! ...Check if n_Profiles in file is > size of output array
    n_input_profiles = SIZE(NLTE_Predictor)
    IF ( n_file_profiles > n_input_profiles ) THEN
      WRITE( msg,'("Number of profiles, ",i0," > size of the output NLTE_Predictor ", &
                  &" array, ",i0,". Only the first ",i0, &
                  &" profiles will be read.")' ) &
                  n_file_profiles, n_input_profiles, n_input_profiles
      CALL Display_Message( ROUTINE_NAME, msg, WARNING )
    END IF
    n_input_profiles = MIN(n_input_profiles, n_file_profiles)

    ! Loop over all the profiles
    Profile_Loop: DO m = 1, n_input_profiles

      ! Read the NLTE predictor data
      ! ...Read the dimensions
      READ( fid, IOSTAT=io_stat ) &
        NLTE_Predictor(m)%n_Layers    , &
        NLTE_Predictor(m)%n_Predictors
      IF ( io_stat /= 0 ) THEN
        WRITE( msg,'("Error reading data dimensions for profile ",i0, &
                    &". IOSTAT = ",i0)' ) m, io_stat
        CALL Read_Cleanup(); RETURN
      END IF
      ! ...Read the logical indicators
      READ( fid, IOSTAT=io_stat ) &
        NLTE_Predictor(m)%Is_Active , &
        NLTE_Predictor(m)%Compute_Tm
      IF ( io_stat /= 0 ) THEN
        WRITE( msg,'("Error reading logical indicators for profile ",i0, &
                    &". IOSTAT = ",i0)' ) m, io_stat
        CALL Read_Cleanup(); RETURN
      END IF
      ! ...Read the array indices
      READ( fid, IOSTAT=io_stat ) &
        NLTE_Predictor(m)%k1, &
        NLTE_Predictor(m)%k2, &
        NLTE_Predictor(m)%isen, &
        NLTE_Predictor(m)%isol
      IF ( io_stat /= 0 ) THEN
        WRITE( msg,'("Error reading array indices for profile ",i0, &
                    &". IOSTAT = ",i0)' ) m, io_stat
        CALL Read_Cleanup(); RETURN
      END IF
      ! ...Read the predictors and interpolation weights
      READ( fid, IOSTAT=io_stat ) &
        NLTE_Predictor(m)%Tm       , &
        NLTE_Predictor(m)%Predictor, &
        NLTE_Predictor(m)%w
      IF ( io_stat /= 0 ) THEN
        WRITE( msg,'("Error reading predictors and interpolation weights for profile ",i0, &
                    &". IOSTAT = ",i0)' ) m, io_stat
        CALL Read_Cleanup(); RETURN
      END IF

      ! Explicitly assign the version number
      NLTE_Predictor(m)%Version = dummy%Version

    END DO Profile_Loop

    ! Close the file
    CLOSE( fid, IOSTAT=io_stat )
    IF ( io_stat /= 0 ) THEN
      WRITE( msg,'("Error closing ",a,". IOSTAT = ",i0)' ) TRIM(Filename), io_stat
      CALL Read_Cleanup(); RETURN
    END IF

    ! Set the return values
    IF ( PRESENT(n_Profiles) ) n_Profiles = n_input_profiles

    ! Output an info message
    IF ( Noisy ) THEN
      WRITE( msg,'("Number of profiles read from ",a,": ",i0)' ) TRIM(Filename), n_input_profiles
      CALL Display_Message( ROUTINE_NAME, TRIM(msg), INFORMATION )
    END IF

   CONTAINS

     SUBROUTINE Read_CleanUp()
       IF ( File_Open(Filename) ) THEN
         CLOSE( fid, IOSTAT=io_stat )
         IF ( io_stat /= 0 ) &
           msg = TRIM(msg)//'; Error closing input file during error cleanup.'
       END IF
       CALL NLTE_Predictor_Destroy( NLTE_Predictor )
       err_stat = FAILURE
       CALL Display_Message( ROUTINE_NAME, msg, err_stat )
     END SUBROUTINE Read_CleanUp

  END FUNCTION NLTE_Predictor_ReadFile

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       NLTE_Predictor_WriteFile
!
! PURPOSE:
!       Function to write NLTE_Predictor object files in Binary format.
!
! CALLING SEQUENCE:
!       Error_Status = NLTE_Predictor_WriteFile( &
!                        Filename           , &
!                        NLTE_Predictor     , &
!                        Quiet    = Quiet     )
!
! INPUTS:
!       Filename:       Character string specifying the name of a
!                       NLTE_Predictor format data file to read.
!                       UNITS:      N/A
!                       TYPE:       CHARACTER(*)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
!       NLTE_Predictor: NLTE_Predictor object containing the NLTE correction
!                       algorithm predictor data.
!                       UNITS:      N/A
!                       TYPE:       NLTE_Predictor_type
!                       DIMENSION:  Rank-1
!                       ATTRIBUTES: INTENT(OUT)
!
! OPTIONAL INPUTS:
!       Quiet:          Set this logical argument to suppress INFORMATION
!                       messages being printed to stdout
!                       If == .FALSE., INFORMATION messages are OUTPUT [DEFAULT].
!                          == .TRUE.,  INFORMATION messages are SUPPRESSED.
!                       If not specified, default is .FALSE.
!                       UNITS:      N/A
!                       TYPE:       LOGICAL
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN), OPTIONAL
!
! FUNCTION RESULT:
!       Error_Status:   The return value is an integer defining the error status.
!                       The error codes are defined in the Message_Handler module.
!                       If == SUCCESS, the file write was successful
!                          == FAILURE, an unrecoverable error occurred.
!                       UNITS:      N/A
!                       TYPE:       INTEGER
!                       DIMENSION:  Scalar
!
!:sdoc-:
!------------------------------------------------------------------------------

  FUNCTION NLTE_Predictor_WriteFile( &
    Filename      , &  ! Input
    NLTE_Predictor, &  ! Input
    Quiet         , &  ! Optional input
    Debug         ) &  ! Optional input (Debug output control)
  RESULT( err_stat )
    ! Arguments
    CHARACTER(*),              INTENT(IN) :: Filename
    TYPE(NLTE_Predictor_type), INTENT(IN) :: NLTE_Predictor(:)
    LOGICAL,         OPTIONAL, INTENT(IN) :: Quiet
    LOGICAL,         OPTIONAL, INTENT(IN) :: Debug
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'NLTE_Predictor_WriteFile'
    ! Function variables
    CHARACTER(ML) :: msg
    LOGICAL :: noisy
    INTEGER :: io_stat
    INTEGER :: fid
    INTEGER :: m, n_output_profiles

    ! Setup
    err_stat = SUCCESS
    ! ...Check Quiet argument
    noisy = .TRUE.
    IF ( PRESENT(Quiet) ) noisy = .NOT. Quiet
    ! ...Override Quiet settings if debug set.
    IF ( PRESENT(Debug) ) THEN
      IF ( Debug ) noisy = .TRUE.
    END IF

    ! Open the file for output
    err_stat = Open_Binary_File( Filename, fid, For_Output=.TRUE. )
    IF ( err_Stat /= SUCCESS ) THEN
      msg = 'Error opening '//TRIM(Filename)
      CALL Write_CleanUp(); RETURN
    END IF

    ! Write the release and version
    WRITE( fid,IOSTAT=io_stat ) NLTE_Predictor(1)%Release, NLTE_Predictor(1)%Version
    IF ( io_stat /= 0 ) THEN
      WRITE( msg,'("Error writing Release/Version. IOSTAT = ",i0)' ) io_stat
      CALL Write_Cleanup(); RETURN
    END IF

    ! Write the dimensions
    n_output_profiles = SIZE(NLTE_Predictor)
    WRITE( fid, IOSTAT=io_stat ) n_output_profiles
    IF ( io_stat /= 0 ) THEN
      WRITE( msg,'("Error writing profile dimension to ",a,". IOSTAT = ",i0)' ) &
             TRIM(Filename), io_stat
      CALL Write_Cleanup(); RETURN
    END IF

    ! Loop over all the profiles
    Profile_Loop: DO m = 1, n_output_profiles

      ! Write the NLTE predictor data
      ! ...Write the dimensions
      WRITE( fid, IOSTAT=io_stat ) &
        NLTE_Predictor(m)%n_Layers    , &
        NLTE_Predictor(m)%n_Predictors
      IF ( io_stat /= 0 ) THEN
        WRITE( msg,'("Error writing data dimensions for profile ",i0, &
                    &". IOSTAT = ",i0)' ) m, io_stat
        CALL Write_Cleanup(); RETURN
      END IF
      ! ...Write the logical indicators
      WRITE( fid, IOSTAT=io_stat ) &
        NLTE_Predictor(m)%Is_Active , &
        NLTE_Predictor(m)%Compute_Tm
      IF ( io_stat /= 0 ) THEN
        WRITE( msg,'("Error writing logical indicators for profile ",i0, &
                    &". IOSTAT = ",i0)' ) m, io_stat
        CALL Write_Cleanup(); RETURN
      END IF
      ! ...Write the array indices
      WRITE( fid, IOSTAT=io_stat ) &
        NLTE_Predictor(m)%k1, &
        NLTE_Predictor(m)%k2, &
        NLTE_Predictor(m)%isen, &
        NLTE_Predictor(m)%isol
      IF ( io_stat /= 0 ) THEN
        WRITE( msg,'("Error writing array indices for profile ",i0, &
                    &". IOSTAT = ",i0)' ) m, io_stat
        CALL Write_Cleanup(); RETURN
      END IF
      ! ...Write the predictors and interpolation weights
      WRITE( fid, IOSTAT=io_stat ) &
        NLTE_Predictor(m)%Tm       , &
        NLTE_Predictor(m)%Predictor, &
        NLTE_Predictor(m)%w
      IF ( io_stat /= 0 ) THEN
        WRITE( msg,'("Error writing predictors and interpolation weights for profile ",i0, &
                    &". IOSTAT = ",i0)' ) m, io_stat
        CALL Write_Cleanup(); RETURN
      END IF

    END DO Profile_Loop

    ! Close the file
    CLOSE( fid, IOSTAT=io_stat )
    IF ( io_stat /= 0 ) THEN
      WRITE( msg,'("Error closing ",a,". IOSTAT = ",i0)' ) TRIM(Filename), io_stat
      CALL Write_Cleanup(); RETURN
    END IF

    ! Output an info message
    IF ( Noisy ) THEN
      WRITE( msg,'("Number of profiles written to ",a,": ",i0)' ) TRIM(Filename), n_output_profiles
      CALL Display_Message( ROUTINE_NAME, TRIM(msg), INFORMATION )
    END IF

   CONTAINS

     SUBROUTINE Write_Cleanup()
       IF ( File_Open(Filename) ) THEN
         CLOSE( fid, IOSTAT=io_stat, STATUS=WRITE_ERROR_STATUS )
         IF ( io_stat /= 0 ) &
           msg = TRIM(msg)//'; Error closing output file during error cleanup.'
       END IF
       err_stat = FAILURE
       CALL Display_Message( ROUTINE_NAME, msg, err_stat )
     END SUBROUTINE Write_Cleanup

  END FUNCTION NLTE_Predictor_WriteFile

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       NLTE_Predictor_IOVersion
!
! PURPOSE:
!       Subroutine to return the module version information.
!
! CALLING SEQUENCE:
!       CALL NLTE_Predictor_IOVersion( Id )
!
! OUTPUT ARGUMENTS:
!       Id:   Character string containing the version Id information
!             for the module.
!             UNITS:      N/A
!             TYPE:       CHARACTER(*)
!             DIMENSION:  Scalar
!             ATTRIBUTES: INTENT(OUT)
!
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE NLTE_Predictor_IOVersion( Id )
    CHARACTER(*), INTENT(OUT) :: Id
    Id = MODULE_VERSION_ID
  END SUBROUTINE NLTE_Predictor_IOVersion

END MODULE NLTE_Predictor_IO


text:
!
! Copyright (c) 2019, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!
module distance_map_acc_m
   use  dist_cuda_m
   use cudafor

contains
   subroutine distance_map_acc(A, B, C, m, n)
      implicit none
      real, intent(in) :: A(:)
      real, intent(in) :: B(:)
      real, intent(inout) :: C(:,:)
      integer, intent(in) :: m, n

      integer :: i, j

!$acc parallel loop copyin(A, B) copyout(C)
      do j = 1,m
!$acc loop
         do i = 1,n
            C(i,j) = dist_cuda(A(j), B(i) )
         end do
      end do
   end subroutine distance_map_acc
end module distance_map_acc_m

max_stars_repo_name: mozhgan-kch/openacc-training-materials


subroutine acc_loop_2d()
    implicit none
    integer i, j, n
    real, allocatable :: buf(:, :)
    real start, finish

    n = 10000
    allocate(buf(n, n))
    !$acc enter data create(buf)

    print*, "2D access:"

    call cpu_time(start)
    !$acc parallel loop present(buf) collapse(2)
    do j=1, n
    do i=1, n
        buf(i, j) = buf(i, j) + 1
    end do
    end do
    call cpu_time(finish)
    print '("I, J order time = ",f6.3," seconds.")', (finish - start)

    call cpu_time(start)
    !$acc parallel loop present(buf) collapse(2)
    do i=1, n
    do j=1, n
        buf(i, j) = buf(i, j) + 1
    end do
    end do
    call cpu_time(finish)
    print '("J, I order time = ",f6.3," seconds.")', (finish - start)

    !$acc exit data delete(buf)
    deallocate(buf)
end subroutine

subroutine acc_loop_3d()
    implicit none
    integer i, j, k, n
    real, allocatable :: buf(:, :, :)
    real start, finish
    n = 1000
    allocate(buf(n, n, n))
    !$acc enter data create(buf)

    print*, "3D access:"

    call cpu_time(start)
    !$acc parallel loop present(buf) collapse(3)
    do k=1, n
    do j=1, n
    do i=1, n
        buf(i, j, k) = buf(i, j, k) + 1
    end do
    end do
    end do
    call cpu_time(finish)
    print '("I, J, K order time = ",f6.3," seconds.")', (finish - start)

    call cpu_time(start)
    !$acc parallel loop present(buf) collapse(3)
    do j=1, n
    do k=1, n
    do i=1, n
        buf(i, j, k) = buf(i, j, k) + 1
    end do
    end do
    end do
    call cpu_time(finish)
    print '("I, K, J order time = ",f6.3," seconds.")', (finish - start)

    call cpu_time(start)
    !$acc parallel loop present(buf) collapse(3)
    do i=1, n
    do k=1, n
    do j=1, n
          buf(i, j, k) = buf(i, j, k) + 1
    end do
    end do
    end do
    call cpu_time(finish)
    print '("J, K, I order time = ",f6.3," seconds.")', (finish - start)

    call cpu_time(start)
    !$acc parallel loop present(buf) collapse(3)
    do k=1, n
    do i=1, n
    do j=1, n
        buf(i, j, k) = buf(i, j, k) + 1
    end do
    end do
    end do
    call cpu_time(finish)
    print '("J, I, K order time = ",f6.3," seconds.")', (finish - start)

    call cpu_time(start)
    !$acc parallel loop present(buf) collapse(3)
    do i=1, n
    do j=1, n
    do k=1, n
        buf(i, j, k) = buf(i, j, k) + 1
    end do
    end do
    end do
    call cpu_time(finish)
    print '("K, J, I order time = ",f6.3," seconds.")', (finish - start)

    call cpu_time(start)
    !$acc parallel loop present(buf) collapse(3)
    do j=1, n
    do i=1, n
    do k=1, n
        buf(i, j, k) = buf(i, j, k) + 1
    end do
    end do
    end do
    call cpu_time(finish)
    print '("K, I, J order time = ",f6.3," seconds.")', (finish - start)

    !$acc exit data delete(buf)

    deallocate(buf)
end subroutine

program main
    implicit none
    call acc_loop_2d
    call acc_loop_3d
end program main


text:
* --------------- BEGIN gm_td_drvr ---------------------------------
      Program gm_td_drvr

* Obtains input parameters and uses time-domain stochastic model to compute
* a suite of accelerations, and from them velocities, displacements, and
* response spectra, from which average peak motions are obtained.

* At most one set of acceleration, velocity, and displacement time series
* will be saved in this driver; see a_ts_drvr for a program that will save
* more than one times series.

* Notes regarding modification of driver:
*  1. All "include" statements assume all subprograms are in the same folder
*     as this program.
*  2. The logical variable "new_mr" must be explicitly set (= .true. for a
*     new value of magnitude and distance); its value is passed through the
*     block common /misc/ (see SMSIM.FI).
*  3. If values of stress other than those in the input parameter file are
*     to be used, set "stressc" equal to the new value, rather than "stress"
*     (and also set the other parameters if the stress changes with moment).

* Dates: 06/01/95 - Written by D. M. Boore;
*                   Renamed and slightly modified PSV_DRVR
*        06/09/95 - Pass r, amag through common rather than parameter lists
*        06/12/95 - Added optional writing of acc, vel time series to a file
*        08/11/95 - Added frequency column to psv output
*        08/11/95 - Changed places where parameter-file and column-file names
*                   are requested
*        08/18/95 - Added call to Write_Params
*        10/17/95 - Added flag to get_params and write_params that tells
*                   whether are dealing with time domain or rv.
*        11/14/95 - Modified output slightly
*        12/08/95 - Added switch to use 91 standard periods
*        12/14/95 - Print out frequency and duration of excitation
*        12/17/95 - Changed location of asking for col file name and
*                   reorder some output
*        12/28/95 - Used new_mr to correct bug related to use of loop over
*                   amag,r
*        12/31/95 - Modified some of the *.sum output
*        01/05/96 - Added column file stem name to column headings
*        01/22/96 - Added total duration, npts to output
*        02/06/96 - Minor formatting improvements
*        04/12/96 - changed names prv, paa to psv, psa; changed
*                   format of time in accvel file to fixed format.
*        02/28/97 - added output of kappa (because I now allow it to be
*                   magnitude dependent)
*        04/20/97 - changed dimension of per from 91 to 400
*        09/01/98 - added calculation of standard deviations
*        09/11/98 - changed array dimensions from 16400 to 33000
*        01/13/99 - added calculation of Arias intensity
*        01/23/99 - modified output to include Sd, and changed
*                   "prv" and "paa" to "psv" and "psa"
*        02/09/99 - Added computation of displacement
*        02/10/99 - Include option to print pga, pgd, pgd and similarly
*                   pga/w, pga/w*w, etc in a column file to provide peak
*                   motion levels in plots of response spectra.  The file
*                   has so many columns that I decided to write the
*                   information to a separate file rather than include it
*                   in the response spectra file.
*        02/13/99 - Changed "accvel" to "avd"
*        03/05/99 - Added include statements at end and write smsim version
*        03/12/99 - Change from 'a\' to 'a' in format statements.  This is
*                   less pleasing aesthetically, but it allows
*                   for redirection of input/output, which did not
*                   work before (a test program showed that redirection
*                   works as long as the total of the sum of the characters
*                   in the query strings written to the screen do not exceed
*                   some undetermined number). Also, change "call banner(5)"
*                   to "call banner(6)" and remove the "cr to quit" in the
*                   "dist (cr to quit)" query
*        03/12/99 - Added check for nper = 1 in loop over computing log-spaced
*                   periods (this is "fool-proofing"; if only one period is
*                   desired the user should reply "y" to the individual
*                   periods query and the computations with nper-1 in the
*                   denominator should never be reached).
*        02/05/00 - Used trim_c for f_out, f_col and corrected error in
*                   making column heads for f_out stem less than 4 characters
*        06/22/00 - Print out standard error of the mean rather than std
*        01/28/01 - Increased array size from 33000 to 66000 (should use
*                   allocatable arrays)
*        02/12/01 - Use dynamically allocatable arrays; renamed "npts" to "npw2"
*        06/13/02 - New organization of subroutines requires initialization
*                   of seed and determination of time series indices in the
*                   calling program, calling gm_td rather than smsim_td, and
*                   use of include statements to bring in necessary
*                   subroutines at compile time.
*                   Also, changed name from "td_drvr" to "gmtddrvr" to indicate
*                   that ground motion is being computed using time domain
*                   simulations.
*        02/09/03 - Only print stress variable if numsource = 1 or 2
*        08/16/05 - Change name from "gmtddrvr" to "gm_td_drvr"
*        05/11/07 - Removed "\smsim\" from include statements
*        08/19/08 - Suggest name of output files
*        02/20/09 - ntaper was not being passed out of get_npts and into acc_ts.
*                   I caught this only when I compiled the program using the -chk
*                   switch.  This probably means that all of my time-domain runs
*                   up till now used ntaper = 0, even if a nonzero taper was
*                   specified (the standard being 0.05).
*        02/21/09 - Added geometric mean computation.  Let pga, pgd asymptotes
*                   be in terms of geometric mean.
*                 - Add 95%-5% duration
*        02/23/09 - Hardwire 75% rather than 95% in determining duration, as
*                   recommended by Ou and Herrmann (this is also consistent with
*                   some Husid plots I made).
*        02/26/09 - Updated calls to get date and time.
*        07/04/09 - Correct small error in summary output if rs not computed.
*        12/04/09 - Write calculated values of dursource and durpath to summary file.
*        12/04/09 - Change subroutine name Get_Npts to Get_Npts_for_TD_Calcs
*        12/20/09 - Use coefficients relating log10_H_eff to M to modify the distance.
*        12/22/09 - Add m, r, h_eff, rmod, numsource, stress to output
*        02/12/10 - As a result of feedback from John Douglas, I replaced "cr" and "psv"
*                   in the screen prompts with "Enter" and "response spectra", respectively.
!        06/01/10 - Suggest file name for avd output.
!        01/04/11 - Add "gm" to column headers in asymptotes file to make it clear
!                   that geometric means are being used (see 02/21/09 comment above).
!        01/17/11 - Compute 95%-5% duration
!        02/02/11 - Write SD for individual simulations to a file named 'SD_sims.out'
!                   (modification in gm_td.for)
!        04/15/11 - I replaced "h_eff" with "f_ff" (because the finite-fault factor is not
!                   really a pseudo-depth).
!                 - Call rmod_calc to compute rmod (rmod = rmod_calc(r, m, iflag, c1, c2, f_ff) )
!        07/05/11 - Write smsim parameter file name in the first line of the output file.
!        12/02/11 - Add computation of am0 right after specification of amag (previously, it
!                   was contained in spect_scale, but it was not clear from where that routine was
!                   called--it was in get_npts_for_td_calcs, called from the driver program.
!                 - I also corrected an error, in which I called rmod_calc, with amag as an
!                   argument, before specifying amag.  I corrected this by changing the order
!                   in which I ask for amag and r (amag first)
!        01/07/14 - Longer file names
!        10/13/14 - Allow for f_ff being defined by two lines
!        12/21/14 - Allow for f_ff being defined by a transition curve between two lines; this
!                   required a change to the calling arguments of rmod_calc

      character colpsv_head*12, colpsv_sem_head*16, colpsa_head*12
      character colsd_head*12
      character colpsvgm_head*14, colpsvgm_sem_head*18,colpsagm_head*14
      character colsdgm_head*14
      character f_params*200, f_out*200, buf*200, f_pavd*200
      character f_suggest*200
      character f_avd*200, save_avd*10
      character psvcalc*10, indvidper*10, f_col*200, message*80,
     :  standard_periods*10, write_pavd*10
      character date*8, time_start*10, time_stop*10
      logical tdflag, fparam_exist

      real acc_save(:), vel_save(:), dis_save(:)
      allocatable :: acc_save, vel_save, dis_save

      real psvsim(400), psvsim_std(400)
      real psvsimgm(400), psvsimgm_std(400)
      real t_pgd(2), t_pgv(2), t_pga(2)
      REAL PER(400)
      DATA PER/0.040,0.042,0.044,0.046,0.048,
     * 0.050,0.055,0.060,0.065,0.070,
     * 0.075,0.080,0.085,0.090,0.095,0.10,0.11,0.12,0.13,0.14,
     * 0.15,0.16,0.17,0.18,0.19,0.20,0.22,0.24,0.26,0.28,
     * 0.30,0.32,0.34,0.36,0.38,0.40,0.42,0.44,0.46,0.48,
     * 0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,
     * 1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,
     * 2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,
     * 4.0,4.2,4.4,4.6,4.8,5.0,5.5,6.0,6.5,7.0,
     * 7.5,8.0,8.5,9.0,9.5,10.0,11.0,12.0,13.0,14.0,
     * 15.0, 309*0.0/

      include 'smsim.fi'

* pi, twopi now computed in smsim.fi
*      pi = 4.0 * atan(1.0)
*      twopi = 2.0 * pi

      call banner(6)

      write(*, '(a)')
     :    ' Enter a message, if desired: '
      message = ' '
      read(*, '(a)') message

333   continue
      f_params = ' '
      write(*, '(a)')
     : ' Enter name of file with parameters (press "Enter" to quit): '
      read(*, '(a)') f_params
      call trim_c(f_params, nc_f_params)
      if (f_params(1:4) .eq. '    ') stop
      inquire(file=f_params(1:nc_f_params), exist=fparam_exist)
      if (.not. fparam_exist) then
        write(*,'(a)') ' ******* FILE DOES NOT EXIST ******* '
        go to 333
      end if

      do i = nc_f_params, 1, -1
        if (f_params(i:i) .eq. '.') then
           indexp = i
           go to 991
        end if
        indexp = 1
      end do
991   continue

      f_suggest = ' '
      f_suggest = f_params(1:indexp)//'gmtd.sum'
      call trim_c(f_suggest,nc_f_suggest)

      f_out = ' '
      write(*, '(a)')
     :    ' Enter name of summary file ("Enter" = '//
     :     f_suggest(1:nc_f_suggest)//'):'
      read(*, '(a)') f_out
      call trim_c(f_out, nc_f_out)
      if (f_out(1:4) .eq. '    ') then
        f_out = f_suggest(1:nc_f_suggest)
      end if
      call trim_c(f_out,nc_f_out)
      call get_lun(nout)
      open(unit=nout,file=f_out(1:nc_f_out),status='unknown')

      call banner(nout)

      write(nout, '(a)') message
      write(nout, '(2a)') ' summary file: ',
     :   f_out(1:20)

      write(nout, '(a)')
     :    ' *** Results computed using gm_td_drvr ***'

* Standard Fortran 90 intrinsic Subroutine DATE_AND_TIME
      call DATE_AND_TIME( date, time_start )
* Date is returned as 'CCYYMMDD' (character date*8)
* Time is returned as 'hhmmss.sss' (character time*10)
*     character datx*8, time_start*10
      write(nout, *)
      write(nout, '(a)')
     :   ' Date: '//date(1:4)//'/'//date(5:6)//'/'//date(7:8)
      write(nout, '(a)')
     : ' Time Start: '//
     : time_start(1:2)//':'//time_start(3:4)//':'//time_start(5:10)
      write(nout, *)

      write(nout, '(a)') ' file with parameters: '//
     :   f_params(1:nc_f_params)

      write(*, '(a)') ' Compute response spectra? (y/n): '
      psvcalc = ' '
      read(*, '(a)') psvcalc
      call trim_c(psvcalc,nc_psvcalc)
      call upstr(psvcalc)
      if (psvcalc(1:1) .eq. 'Y') then
        write(*, '(a)')
     :  ' Enter individual periods? (y/n): '
        indvidper = ' '
        read(*, '(a)') indvidper
        call trim_c(indvidper, nc_indvidper)
        call upstr(indvidper)
        if (indvidper(1:1) .eq. 'Y') then
          write(*, '(a)')
     :       ' Enter fractional damping, nperiods: '
          read(*, *) damp, nper
          do i = 1, nper
            write(*, '(a)') '    Enter oscillator period: '
            read(*,*) per(i)
          end do
        else
          write(*, '(a)') ' Use standard set of 91 periods? (y/n): '
          standard_periods = ' '
          read(*, '(a)') standard_periods
          call trim_c(standard_periods, nc_standard_periods)
          call upstr(standard_periods)
          nper = 91
          if (standard_periods(1:1) .eq. 'Y') then
            nper = 91
            write(*, '(a)') ' Enter fractional damping: '
            read(*, *) damp
          else
            write(*, '(a)')
     :       ' Enter fractional damping, perlow, perhigh, nper: '
            read(*,*) damp, perlow, perhigh, nper
            do i = 1, nper
              if (nper .eq. 1) then
                per(i) = perlow
              else
                per(i) =
     :            perlow * (perhigh/perlow)**(float(i-1)/float(nper-1))
              end if
            end do
          end if
        end if
      else
        nper = 0
      end if

      tdflag = .true.  ! controls extent of the input file read and written
      call get_params( f_params, tdflag )
      call write_params(nout, tdflag)

100   continue

      write(*, '(a)') ' amag: '
      read(*, *) amag
      am0 = 10.**(1.5*amag + 16.05)

      write(*, '(a)') ' dist: '
      read(*, *) r

!Replaced h_eff with f_ff
      rmod = rmod_calc(r, amag,
     :                  iflag_f_ff, nlines_f_ff,
     :                  c1_log10_f_ff, c2_log10_f_ff, mh_f_ff,
     :                  c3_log10_f_ff, c4_log10_f_ff,
     :                  m1t_f_ff, m2t_f_ff,
     :                  c0t_f_ff, c1t_f_ff, c2t_f_ff, c3t_f_ff,
     :                                                          f_ff)

      new_mr = .true.

      if (psvcalc(1:1) .eq. 'Y') then
        f_suggest = ' '
        f_suggest = f_params(1:indexp-1)//'.mxxxxr00000_rs.td.col'
        write(f_suggest(indexp+2:indexp+5),'(f4.2)') amag
        if (rmod .lt. 10.0) then
          write(f_suggest(indexp+9:indexp+11),'(f3.1)') rmod
        else if (rmod .lt. 100.0) then
          write(f_suggest(indexp+8:indexp+11),'(f4.1)') rmod
        else
          write(f_suggest(indexp+7:indexp+11),'(f5.1)') rmod
        end if
        call trim_c(f_suggest,nc_f_suggest)

        f_col = ' '
        write(*, '(a)')
     :    ' Enter name of column file to contain'//
     :    ' response spectra ("Enter" = ' //
     :     f_suggest(1:nc_f_suggest) // '):'
        read(*, '(a)') f_col
        call trim_c(f_col, nc_f_col)
        if (f_col(1:4) .eq. '    ') then
          f_col = f_suggest(1:nc_f_suggest)
        end if
        call trim_c(f_col,nc_f_col)
        call get_lun(ncol)
        open(unit=ncol, file=f_col(1:nc_f_col), status='unknown')
        write(ncol,'(a)') ' Parameter file: '//f_params(1:nc_f_params)
        nc_f_col_stem = min0(8, index(f_col(1:nc_f_col),'.') - 1)
        colsd_head  =     ' sd:'//f_col(1:nc_f_col_stem)//'________'
        colpsv_head =     'psv:'//f_col(1:nc_f_col_stem)//'________'
        colpsv_sem_head ='sem/psv:'//f_col(1:nc_f_col_stem)//'________'
        colpsa_head =     'psa:'//f_col(1:nc_f_col_stem)//'________'

        colsdgm_head  =     ' sdgm:'//f_col(1:nc_f_col_stem)//'________'
        colpsvgm_head =     'psvgm:'//f_col(1:nc_f_col_stem)//'________'
        colpsvgm_sem_head =
     :                  'sem_psvgm:'//f_col(1:nc_f_col_stem)//'________'
        colpsagm_head =     'psagm:'//f_col(1:nc_f_col_stem)//'________'

        write(ncol, '(
     :    2x,a, 6x,a, 5x,a,
     :    6x,a, 7x,a, 4x,a, 4x,a,
     :    1x,a, 2x,a,
     :    1x,a, 1x,a, 1x,a, 1x,a,
     :    1x,a, 1x,a, 1x,a, 1x,a
     :                                           )')
     :   'damp', 'per', 'freq',
     :   'M', 'R', 'F_FF', 'Rmod',
     :   'numsource', 'stress',
     :    colsd_head, colpsv_head, colpsa_head, colpsv_sem_head,
     :    colsdgm_head, colpsvgm_head, colpsagm_head, colpsvgm_sem_head

        write_pavd = ' '
        if (indvidper(1:1) .ne. 'Y') then
          write(*, '(a)')
     :      ' Write pga, pgv, pgd to a column file? (y/n): '
          read(*, '(a)') write_pavd
          call trim_c(write_pavd,nc_write_pavd)
          call upstr(write_pavd)
          if (write_pavd(1:1) .eq. 'Y') then
            f_suggest= ' '
            f_suggest= f_params(1:indexp-1)//'.mxxxxr00000_pavd.td.col'
            write(f_suggest(indexp+2:indexp+5),'(f4.2)') amag
            if (r .lt. 10.0) then
              write(f_suggest(indexp+9:indexp+11),'(f3.1)') r
            else if (r .lt. 100.0) then
              write(f_suggest(indexp+8:indexp+11),'(f4.1)') r
            else
              write(f_suggest(indexp+7:indexp+11),'(f5.1)') r
            end if
            call trim_c(f_suggest,nc_f_suggest)

            write(*, '(a)')
     :        ' Enter name of column file for to contain' //
     :        ' pga, pgv, pgd ("Enter" = '//
     :                       f_suggest(1:nc_f_suggest) // '):'
            f_pavd = ' '
            read( *, '(a)') f_pavd
            call trim_c(f_pavd, nc_f_pavd)
            if(f_pavd(1:4) .eq. '    ') then
              f_pavd = f_suggest(1:nc_f_suggest)
            end if
            call trim_c(f_pavd,nc_f_pavd)

            call get_lun(npavd)
            open(unit=npavd, file=f_pavd(1:nc_f_pavd), status='unknown')
            tl = per(1)
            th = per(nper)
            tdiv = (th/tl)**(0.25)
            t_pga(1) = tl
            t_pga(2) = tl*tdiv
            t_pgv(1) = sqrt(th*tl/tdiv)
            t_pgv(2) = sqrt(th*tl*tdiv)
            t_pgd(1) = th/tdiv
            t_pgd(2) = th

            print *,' tl, th, tdiv, t_pgd(1,2)=',
     :                tl, th, tdiv, t_pgd(1),t_pgd(2)

            write(npavd,
     :            '(t4,a,    t15,a,    t24,a,    t34,a,
     :             t45,a,    t54,a,    t67,a,    t76,a,
     :             t86,a,    t94,a,   t106,a,   t119,a)')
     :           't_pgd',    'pgdgm',  'w*pgdgm', 'w2*pgdgm',
     :           't_pgv',  'pgvgm/w',    'pgvgm',  'w*pgvgm',
     :           't_pga', 'pgagm/w2',  'pgagm/w',    'pgagm'
          end if
        end if
      end if  ! end loop to compute response spectrum

      write(*, '(a)') ' Save acc & vel time series? (y/n): '
      save_avd = ' '
      read(*, '(a)') save_avd
      call upstr(save_avd)
      call trim_c(save_avd, nc_save_avd)
      if (save_avd(1:1) .eq. 'Y') then
        write(*, '(a,i4,a)')
     : '    Save which simulation (out of ', nsims, ')? '
        read(*, *) nacc_save

        f_suggest = ' '
        f_suggest = f_params(1:indexp-1)//'.mxxxxr00000_avd.td.col'
        write(f_suggest(indexp+2:indexp+5),'(f4.2)') amag
        if (rmod .lt. 10.0) then
          write(f_suggest(indexp+9:indexp+11),'(f3.1)') rmod
        else if (rmod .lt. 100.0) then
          write(f_suggest(indexp+8:indexp+11),'(f4.1)') rmod
        else
          write(f_suggest(indexp+7:indexp+11),'(f5.1)') rmod
        end if
        call trim_c(f_suggest,nc_f_suggest)

        write(*, '(a)')
     : '    Enter file name for time series ("Enter" = ' //
     :                       f_suggest(1:nc_f_suggest) // '):'
        f_avd = ' '
        read(*, '(a)') f_avd
        call trim_c(f_avd, nc_f_avd)
        if (f_avd(1:3) .eq. '   ') then
           f_avd = f_suggest(1:nc_f_suggest)
           call trim_c(f_avd, nc_f_avd)
        end if
      else
        nacc_save = 0
      end if

      call DATE_AND_TIME( date, time_start )

      call Get_Npts_for_TD_Calcs(nstart, nstop, npw2, te, ntaper)

      ndimen = npw2
      allocate (acc_save(ndimen), vel_save(ndimen), dis_save(ndimen))

      call gm_td(per, nper,
     :   psvsim, psvsim_std,
     :   pgasim, pgasim_std,
     :   pgvsim, pgvsim_std,
     :   pgdsim, pgdsim_std,
     :   arias, arias_std,
     :   psvsimgm, psvsimgm_std,
     :   pgasimgm, pgasimgm_std,
     :   pgvsimgm, pgvsimgm_std,
     :   pgdsimgm, pgdsimgm_std,
     :   ariasgm, ariasgm_std,
     :   dur_75_05, dur_75_05gm,
     :   dur_95_05, dur_95_05gm,
     :   psvcalc, nacc_save, acc_save, vel_save, dis_save,
     :   nstart, nstop, npw2, te, ntaper)

      if (save_avd(1:1) .eq. 'Y') then
        call get_lun(nuacc)
        open(unit=nuacc,file=f_avd(1:nc_f_avd),status='unknown')
        write(nuacc, '(t12,a,t24,a,t36,a,t48,a)') 'T', 'A', 'V', 'D'
        do i = 1, npw2
          write(nuacc, '(1x,f11.6,1p3(1x,e11.4))')
     :      float(i-1)*dt, acc_save(i), vel_save(i), dis_save(i)
        end do
      end if

      deallocate (acc_save, vel_save, dis_save)

      write(nout, *)
      write(nout, '(a)') ' *********** NEW R AND M ***********'
      write(nout, '(a, 1p3(1x,e10.3))')
     :      ' r, f_ff, rmod, amag, kappa = ',
     :        r, f_ff, rmod, amag, kappa_f(amag)
      if (psvcalc(1:1) .eq. 'Y') then
        write(nout, '(2x,2a)') 'File with response spectra: ',
     :    f_col(1:30)
      else
        write(nout, '(2x,a)') 'No response spectra computed'
      end if
      write(nout, '(2x,a,1x,1pe10.3)') ' const= ', const
      if (numsource .eq. 1 .or. numsource .eq. 2) then
        write(nout, '(2x,a, f6.3,
     :      1p, 1x,e9.2, 2(1x,e10.3), 3(1x,e9.2))')
     :    ' amag, stress, fa, fb, dursource, durpath, durex= ',
     :      amag, stress, fa, fb, dursource_calc, durpath_calc, durex
      else
        write(nout, '(2x,a, f6.3,
     :      1p, 2(1x,e10.3), 3(1x,e9.2))')
     :    ' amag, fa, fb, dursource, durpath, durex= ',
     :      amag, fa, fb, dursource_calc, durpath_calc, durex
      end if
      write(nout, '(2x,a,1p2(1x,e10.3))') ' am0, am0b_m0fa= ',
     :                                       am0, am0b_m0
      write(nout, '(2x,a, i6, 1x,i6, 1x,f8.5, 2(1x,f6.1))')
     :  ' npw2, ntaper, dt, te, total duration = ',
     :    npw2, ntaper, dt, te, float(npw2)*dt

      write(nout,'(2x,a, 2x,a,
     :             1x,a, 2x,a,
     :             1x,a, 2x,a)')
     :    'pgdam(cm)',    'sem/pgdam',
     :    'pgvam(cm/s)',  'sem/pgvam',
     :    'pgaam(cm/s2)', 'sem/pgaam'
      write(nout,'(1p,
     :         1x,e10.3, 1x,e10.3,
     :         2x,e10.3, 1x,e10.3,
     :         3x,e10.3, 1x,e10.3)')
     :     pgdsim, (pgdsim_std/sqrt(float(nsims)))/pgdsim,
     :     pgvsim, (pgvsim_std/sqrt(float(nsims)))/pgvsim,
     :     pgasim, (pgasim_std/sqrt(float(nsims)))/pgasim

      write(nout,'(2x,a, 2x,a,
     :             1x,a, 2x,a,
     :             1x,a, 2x,a)')
     :    'pgdgm(cm)',    'sem_pgdgm',
     :    'pgvgm(cm/s)',  'sem_pgvgm',
     :    'pgagm(cm/s2)', 'sem_pgagm'
      write(nout,'(1p,
     :         1x,e10.3, 1x,e10.3,
     :         2x,e10.3, 1x,e10.3,
     :         3x,e10.3, 1x,e10.3
     :                                      )')
     :     pgdsimgm,
     :     (pgdsimgm_std**(1.0/sqrt(float(nsims)))),
     :     pgvsimgm,
     :     (pgvsimgm_std**(1.0/sqrt(float(nsims)))),
     :     pgasimgm,
     :     (pgasimgm_std**(1.0/sqrt(float(nsims))))

      write(nout,'( 1x,a, 2x,a)')
     :    'Arias intensity(cm/s)', 'sem/Arias'
      write(nout,'(1p, 10x,e10.3,  10x,e10.3)')
     :     arias, (arias_std/sqrt(float(nsims)))/arias

      write(nout,'( 1x,a, 1x,a)')
     :    'Arias intensitygm(cm/s)', 'sem_Ariasgm'
      write(nout,'(1p, 14x,e10.3,  2x,e10.3)')
     :     ariasgm,
     :     (ariasgm_std**(1.0/sqrt(float(nsims))))

      write(*,*)

      if (psvcalc(1:1) .eq. 'Y') then

        write(nout,'(2x,a,f6.3)')
     :        ' Fractional oscillator damping = ', damp
        write(nout,'(
     :                4x,a, 6x,a,
     :                2x,a,
     :                4x,a,
     :                1x,a, 4x,a,
     :                1x,a,
     :                2x,a,
     :                1x,a, 2x,a
     :                                                        )')
     :             'per(s)', 'freq',
     :             'psv(cm/s)',
     :             'sem/psv',
     :             'psa(cm/s2)', 'psd(cm)',
     :             'psvgm(cm/s)',
     :             'psvgm_sem',
     :             'psagm(cm/s2)', 'psdgm(cm)'

        do i = 1, nper
          write(nout, '(
     :          1x,f9.4,  1x,f9.4, 1p,
     :          1x,e10.3,
     :          1x,e10.3,
     :          1x,e10.3, 1x,e10.3,
     :          2x,e10.3,
     :          1x,e10.3,
     :          2x,e10.3, 1x,e10.3
     :                                       )')
     :      per(i), 1.0/per(i),
     :      psvsim(i),
     :      (psvsim_std(i)/sqrt(float(nsims)))/psvsim(i),
     :      (twopi/per(i))*psvsim(i), psvsim(i)/(twopi/per(i)),
     :      psvsimgm(i),
     :      (psvsimgm_std(i)**(1.0/sqrt(float(nsims)))),
     :      (twopi/per(i))*psvsimgm(i), psvsimgm(i)/(twopi/per(i))

          if (numsource == 1 .or. numsource == 2 .or.
     :        numsource == 11) then
            stress_out = stress
          else
            stress_out = 0.0
          end if

          write(ncol,'(
     :      1x,f5.3, 1x,f8.4, 1x,f8.3,
     :      1x,f6.3, 1x,f7.2, 1x,f7.2, 1x,f7.2,
     :      8x,i2, 1x,f7.2,
     :      1p,
     :              3x,e10.3, 3x,e10.3,
     :              3x,e10.3,
     :              7x,e10.3,
     :              5x,e10.3, 5x,e10.3,
     :              5x,e10.3,
     :              9x,e10.3
     :                                  )')
     :      damp, per(i), 1.0/per(i),
     :      amag, r, f_ff, rmod,
     :      numsource, stress_out,
     :      psvsim(i)/(twopi/per(i)), psvsim(i),
     :      (twopi/per(i))*psvsim(i),
     :      (psvsim_std(i)/sqrt(float(nsims)))/psvsim(i),
     :      psvsimgm(i)/(twopi/per(i)), psvsimgm(i),
     :      (twopi/per(i))*psvsimgm(i),
     :      (psvsimgm_std(i)**(1.0/sqrt(float(nsims))))
        end do

        close(unit=ncol)
        if (write_pavd(1:1) .eq. 'Y') then
          do i = 1, 2
            write(npavd, '(3(0p1x,f7.3,1p3(1x,e10.3)))')
     :       t_pgd(i),pgdsimgm,
     :                pgdsimgm*(twopi/t_pgd(i)),
     :                pgdsimgm*(twopi/t_pgd(i))**2,
     :       t_pgv(i),pgvsimgm/(twopi/t_pgv(i)),
     :                pgvsimgm,
     :                pgvsimgm*(twopi/t_pgv(i)),
     :       t_pga(i),pgasimgm/(twopi/t_pga(i))**2,
     :                pgasimgm/(twopi/t_pga(i)),
     :                pgasimgm
          end do
          close(unit=npavd)
        end if

      end if

      write(nout, *)
      call DATE_AND_TIME( date, time_stop )
      call time_diff(time_start, time_stop, time_elapsed)
      write(nout, '(a,1x,1pe10.3)')
     :     ' Elapsed time (sec): ', time_elapsed

      write(*,*)
      write(*, '(a)')
     : ' Compute results for another r and M? (y/n;"Enter"=quit): '
      read(*, '(a)') buf
      if (buf(1:4) .eq. '    ') go to 999
      if (buf(1:1) .eq. 'n' .or. buf(1:1) .eq. 'N') go to 999

      goto 100

999   continue

      close(unit=nout)
      stop
      end
* --------------- END gm_td_drvr ---------------------------------

!      include 'gm_td.for'
!      include 'acc_ts.for'
!      include 'td_subs.for'
   !NOTE: run make_td_subs.bat to make the file; it now includes rv_td_subs, recipes, etc
!      include 'rv_td_subs.for'
!      include 'recipes.for'
      include 'td_subs.for'
      include 'smsim_util_subs.for'

ext: for

max_stars_repo_name: cliftbar/FlaPyDisaster

max_issues_repo_name: cliftbar/FlaPyDisaster

max_forks_repo_name: cliftbar/FlaPyDisaster


text:
!------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in     !
!  continuous development by various groups and is based on information  !
!  from these groups: Federal Government employees, contractors working  !
!  within a United States Government contract, and non-Federal sources   !
!  including research institutions.  These groups give the Government    !
!  permission to use, prepare derivative works of, and distribute copies !
!  of their work in the CMAQ system to the public and to permit others   !
!  to do so.  The United States Environmental Protection Agency          !
!  therefore grants similar permission to use the CMAQ system software,  !
!  but users are requested to provide copies of derivative works or      !
!  products designed to operate in the CMAQ system to the United States  !
!  Government without restrictions as to use by others.  Software        !
!  that is used with the CMAQ system but distributed under the GNU       !
!  General Public License or the GNU Lesser General Public License is    !
!  subject to their copyright restrictions.                              !
!------------------------------------------------------------------------!

        SUBROUTINE PM3WARN( CALLER, JDATE, JTIME, ERRTXT )
C.....................................................................
C
C  PURPOSE:   Wrapper for M3WARN in parallel environment, to
C             add the processor-id suffix to the name of the
C             caller.
C
C
C  PRECONDITIONS REQUIRED:  Same as for M3WARN.
C
C
C  REVISION  HISTORY:
C       Original version 07/1998 by Al Bourgeois.
C       Modified 12/07/1998 by Al Bourgeois to add EXTERNAL declarations.
C       Modified 02/23/2011 by Shawn Roselle
C          -- Replaced I/O API include files with M3UTILIO; removed
C             deprecated TRIMLEN
C
C
C  ARGUMENT LIST DESCRIPTION:
C  IN:
C     CHARACTER*(*)   CALLER          ! Name of the caller.
C     INTEGER         JDATE           ! Model date for the error.
C     INTEGER         JTIME           ! Model time for the error.
C     CHARACTER*(*)   ERRTXT          ! Error message.
C
C     COMMON BLOCK PIOVARS:
C     INTEGER  MY_PE               !  Local processor id.
C
C
C  SUBROUTINES AND FUNCTIONS CALLED:  M3WARN.
C
C***********************************************************************

      USE M3UTILIO              ! i/o api

      IMPLICIT NONE

C...........   INCLUDES:

      INCLUDE 'PIOVARS.EXT'      ! Parameters for parallel implementation.

C...........   ARGUMENTS and their descriptions:

      CHARACTER*(*)   CALLER          ! Name of the caller.
      INTEGER         JDATE           ! Model date for the error.
      INTEGER         JTIME           ! Model time for the error.
      CHARACTER*(*)   ERRTXT          ! Error message.

C...........   LOCAL VARIABLES

      INTEGER      LENSTR       ! String length of CALLER.
      CHARACTER*3  CMYPE        ! Processor ID string.
      CHARACTER*7  PE_STR       ! String suffix to go with processor ID.
      CHARACTER*16 CALL16       ! First 16 characters of CALLER.
      CHARACTER*26 PCALLER      ! New caller string with PE information.

C.............................................................................
C   begin subroutine PM3WARN

C.......  Create strings to append to CALLER.
      WRITE (PE_STR,'(A7)') ' on PE '
      WRITE(CMYPE,'(I3.3)') MY_PE

C.......  Construct new CALLER string.
      LENSTR = MIN( 16, LEN_TRIM(CALLER) )
      CALL16 = CALLER( 1: LENSTR )
      PCALLER = CALL16(1:LENSTR)//PE_STR//CMYPE

C.......  Pass the new sting to M3WARN.
      CALL M3WARN( PCALLER, JDATE, JTIME, ERRTXT )

      RETURN
      END

ext: f

max_stars_repo_name: Simeng-unique/CMAQ-changed

max_issues_repo_name: Simeng-unique/CMAQ-changed

max_forks_repo_name: Simeng-unique/CMAQ-changed


text:
*DECK XADJ
      SUBROUTINE XADJ (X, IX, IERROR)
C***BEGIN PROLOGUE  XADJ
C***PURPOSE  To provide single-precision floating-point arithmetic
C            with an extended exponent range.
C***LIBRARY   SLATEC
C***CATEGORY  A3D
C***TYPE      SINGLE PRECISION (XADJ-S, DXADJ-D)
C***KEYWORDS  EXTENDED-RANGE SINGLE-PRECISION ARITHMETIC
C***AUTHOR  Lozier, Daniel W., (National Bureau of Standards)
C           Smith, John M., (NBS and George Mason University)
C***DESCRIPTION
C     REAL X
C     INTEGER IX
C
C                  TRANSFORMS (X,IX) SO THAT
C                  RADIX**(-L) .LE. ABS(X) .LT. RADIX**L.
C                  ON MOST COMPUTERS THIS TRANSFORMATION DOES
C                  NOT CHANGE THE MANTISSA OF X PROVIDED RADIX IS
C                  THE NUMBER BASE OF SINGLE-PRECISION ARITHMETIC.
C
C***SEE ALSO  XSET
C***REFERENCES  (NONE)
C***ROUTINES CALLED  XERMSG
C***COMMON BLOCKS    XBLK2
C***REVISION HISTORY  (YYMMDD)
C   820712  DATE WRITTEN
C   881020  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
C   901019  Revisions to prologue.  (DWL and WRB)
C   901106  Changed all specific intrinsics to generic.  (WRB)
C           Corrected order of sections in prologue and added TYPE
C           section.  (WRB)
C           CALLs to XERROR changed to CALLs to XERMSG.  (WRB)
C   920127  Revised PURPOSE section of prologue.  (DWL)
C***END PROLOGUE  XADJ
      REAL X
      INTEGER IX
      REAL RADIX, RADIXL, RAD2L, DLG10R
      INTEGER L, L2, KMAX
      COMMON /XBLK2/ RADIX, RADIXL, RAD2L, DLG10R, L, L2, KMAX
      SAVE /XBLK2/
C
C   THE CONDITION IMPOSED ON L AND KMAX BY THIS SUBROUTINE
C IS
C     2*L .LE. KMAX
C
C THIS CONDITION MUST BE MET BY APPROPRIATE CODING
C IN SUBROUTINE XSET.
C
C***FIRST EXECUTABLE STATEMENT  XADJ
      IERROR=0
      IF (X.EQ.0.0) GO TO 50
      IF (ABS(X).GE.1.0) GO TO 20
      IF (RADIXL*ABS(X).GE.1.0) GO TO 60
      X = X*RAD2L
      IF (IX.LT.0) GO TO 10
      IX = IX - L2
      GO TO 70
   10 IF (IX.LT.-KMAX+L2) GO TO 40
      IX = IX - L2
      GO TO 70
   20 IF (ABS(X).LT.RADIXL) GO TO 60
      X = X/RAD2L
      IF (IX.GT.0) GO TO 30
      IX = IX + L2
      GO TO 70
   30 IF (IX.GT.KMAX-L2) GO TO 40
      IX = IX + L2
      GO TO 70
   40 CALL XERMSG ('SLATEC', 'XADJ', 'overflow in auxiliary index', 107,
     +             1)
      IERROR=107
      RETURN
   50 IX = 0
   60 IF (ABS(IX).GT.KMAX) GO TO 40
   70 RETURN
      END


text:
program test_ewald4

! This test compares lattice energy calculated from the Ewald summation
! against the energy calculated using the Madelung constant for various lattice
! constants L. The agreement is essentially to machine precision.

! This test is testing the example from section 4.1 in [1].
!
! [1] Pask, J. E., Sukumar, N., & Mousavi, S. E. (2012). Linear scaling
! solution of the all-electron Coulomb problem in solids. International Journal
! for Multiscale Computational Engineering, 10(1), 83–99.
! doi:10.1615/IntJMultCompEng.2011002201

use types, only: dp
use constants, only: ang2bohr, kJmol2Ha
use ewald_sums, only: ewald_box, ewald_fft1, fred2fcart
use utils, only: assert, init_random
implicit none

integer :: natom, ntypat
real(dp), parameter :: Llist(4) = [2 / sqrt(3._dp), 1._dp, 2._dp, 3._dp]
real(dp), parameter :: alpha = 0.90980836237716_dp ! Madelung constant
real(dp), parameter :: stress0 = 1.213077816502878_dp

real(dp) :: stress(6)
real(dp), allocatable :: xred(:, :), forces(:, :), q(:)
real(dp) :: L, E_ewald, E_ewald_fft, E_madelung
integer :: i

natom = 2
ntypat = 1
allocate(xred(3, natom), forces(3, natom), q(natom))
xred(:, 1) = [0, 0, 0]
xred(:, 2) = [1, 1, 1] / 2._dp
q = [1, 1]

do i = 1, size(Llist)
    L = Llist(i)
    call ewald_box(L, xred*L, q, E_ewald, forces, stress)
    E_ewald = E_ewald / (natom/ntypat)

    call ewald_fft1(L, xred*L, q, 64, 0.25_dp*L, E_ewald_fft)
    E_ewald_fft = E_ewald_fft / (natom/ntypat)

    E_madelung = -2*alpha/L

    print *, "a =", L/ang2bohr*100, "pm"
    print *, "Madelung: ", E_madelung / kJmol2Ha, "kJ/mol =", E_madelung, "a.u."
    print *, "Ewald:    ", E_ewald / kJmol2Ha, "kJ/mol =", E_ewald, "a.u."
    print *, "Ewald FFT:", E_ewald_fft / kJmol2Ha, "kJ/mol =",E_ewald_fft,"a.u."
    print *, "Ewald error:     ", abs(E_ewald - E_madelung), "a.u."
    print *, "Ewald FFT error: ", abs(E_ewald_fft - E_madelung), "a.u."
    call assert(abs(E_ewald - E_madelung) < 5e-14_dp)
    call assert(abs(E_ewald_fft - E_madelung) < 5e-12_dp)
    call assert(all(abs(forces) < 1e-15_dp))
    call assert(all(abs(stress - (-stress0/L)*[1, 1, 1, 0, 0, 0]) < 1e-15_dp))
end do

end program

max_stars_repo_name: certik/hfsolver

max_issues_repo_name: certik/hfsolver

max_forks_repo_name: certik/hfsolver


text:
\ Calculation of digits of pi without floating-point
\ The algorithm calculate one digit at the time
\ Michel Jean, April 2020

1001 constant nbr-digits                \ number of digit - add one, the last digit is not printed
variable pos                            \ position in the array
variable ret                            \ carry (retenue in french ;-) )
nbr-digits  10 * 3 / constant nbr-cells \ initialisation on the number of cells
variable arr-cells nbr-cells cells allot

: initialisation ( -- ) \ initialisation of the array
      nbr-cells 0 do
        2 arr-cells i cells + !
        loop ;

: algo-base ( -- )
        pos @ 1 = if
        1 cells arr-cells + @ 10 * ret @ +
        dup
        10 mod arr-cells 1 cells + !
        10 / ret !
        else
                pos @ cells arr-cells + @ 10 * ret @ +
                2 pos @ * 1 -
                2dup
                mod arr-cells pos @ cells + !
                / pos @ 1 - * ret !
        then ;

: release-stack                 \ hold the last digit and release other digits of the stack
        depth 1 > if depth 1- roll . recurse else then ;

: release-stack+1               \ hold the last digit and release other digits of the stack + 1 modulo 10
        depth 1 > if
        depth 1- roll 1+ 10 mod .
        recurse else then ;

: set-predigit ( -- n )
        ret @ dup 9 < if release-stack else
        dup 10 = if release-stack+1 drop 0 else
        then then ;

: 1digit ( -- )                 \ calculate 1 digit at the time
        1 nbr-cells do
        i pos ! algo-base -1 +loop set-predigit ;

: run ( -- )
        0 ret ! initialisation
        nbr-digits 50 / 0 do
        50 0 do 1digit loop     \ 50 digits by line
        cr loop
        drop ;                  \ drop the last digit (value not safe)
run


text:
c=======================================================================
c
c     subroutine ALLOCMVTC
c
c     MV allocation srtategy with transaction costs (linear)
c
c     Min[ w*Q*w + alpha*CostTransaction(w)]
c      s.t.
c     C*w <= b                 (linear constraints)
c     Cinf <= w <= Csup        (lower/upper bounds)
c
c        w   : portfolio weights
c        Q   : covariance matrix
c        rho : assets performance
c
c-----------------------------------------------------------------------
      SUBROUTINE allocmvtc ( n, cov, neq, nin, ccst, bcst, cinf, csup,
     &                       cost, wini, alpha,
     &                       iwork, dwork, wopt, info )
c-----------------------------------------------------------------------
c
c     INPUT
c       n      : number of assets                                integer
c       cov    : covariance matrix (n*n)                          double
c       neq    : number of equality constraints                  integer
c       nin    : number of inequality constraints                integer
c       ccst   : matrix of constraints (n*(neq+nin))              double
c       bcst   : vector of constraints (neq+nin)                  double
c       cinf   : lower bound (n)                                  double
c       csup   : upper bound (n)                                  double
c       cost   : cost transaction (n)                             double
c       wini   : initial portfolio (n)                            double
c       alpha  : transaction cost parameter                       double
c
c     WORKSPACE
c       iwork  : 26*n + 2*nin + neq + 1                          integer
c       dwork  : 2*n*(16*n + nin + neq + 44) + 2(2*nin + neq)     double
c
c     OUTPUT
c       wopt   : optimal portfolio (n)                            double
c       info   : diagnostic argument                             integer
c
c     CALL
c       YM      : matrix copy
c       YV      :
c       QP      : quadratic solver
c       CTITCOV :
c       CTITCST :
c       TESTSDP :
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
c
c     arguments i/o
      INTEGER n, info, neq, nin
      DOUBLE PRECISION alpha
      DOUBLE PRECISION cov(*), cinf(*), csup(*), ccst(n,*), bcst(*),
     &                 wopt(*), cost(*), wini(*)
c
c     workspaces
      INTEGER iwork(*)
      DOUBLE PRECISION dwork(*)
c
c     local variables
      INTEGER p, i, neqtot, nintot, piwork, pdwork, piwk,
     &        pigk, pdsp, pdss, pdsa, pdwk, piwq, pdlin, pdwq, pdlagr,
     &        pdbcs, pdccs, pdro, pdwopt, pdcov, pdcinf, pdcsup, test
      DOUBLE PRECISION eps, epsilon, maxv
      PARAMETER ( eps = 1.E-8, epsilon = 1.E-15, maxv = 1.E+15 )
c
c     external subroutines
      EXTERNAL qp, YM, YV, ctitcov, ctitcst, testsdp
c
c-----------------------------------------------------------------------
c
c     initialization
      info = 0
      test = 0
      p = n
c
c     test if cost(i) = 0 for i = 1,...,n
      DO i = 1,n
        IF (cost(i) .GT. epsilon) test = 1
      ENDDO
c
c     case alpha > eps (with trans. costs)
      IF ((alpha .GT. eps).AND.(test .EQ. 1)) p = n + n
c
c     pointers for integer workspace : iwork
c     --------------------------------------
      piwork = 1
      piwk   = piwork
c     piwk  : pointer for workspace of KATAVEP who needs ( 12*p )
      pigk   = piwk + ( 12*p )
c     pigk  : pointer Kato groups, so (p) more
c     the call of KATAVEP needs ( 13*p )
c
c     pointers for QP who uses a part of the same space
      piwq   = piwork
c     piwq  : pointer for QP who needs
c             ( 3*p + 2*nintot + neqtot + 1 ),
c          so ( 3*p + 2*(nin + p) + neq + 1 )
c          so ( 5*p + 2*nin + neq + 1 )
c
c     together KATAVEP and QP need (union of space)
c     ( 13*p + 2*nin + neq + 1 )
c
c     Total size of iwork array = ( 13*p + 2*nin + neq + 1 )
c
c     pointers for double precision workspace  : dwork
c     ------------------------------------------------
      pdwork = 1
      pdwopt  = pdwork
c     pdwopt : pointer for the optimal point, so ( p ) more
      pdcov  = pdwopt + p
c     pdcov  : pointer for quadratic part, so (p*p) more
      pdsp   = pdcov + ( p*p )
c     pdsp   : pointer for vector spectr of KATAVEP, so ( p ) more
      pdss   = pdsp + ( p )
c     pdss   : pointer for vector eigsor of KATAVEP, so ( p ) more
      pdsa   = pdss + ( p )
c     pdsa   : pointer for vector spectr of KATAVEP, so ( p ) more
      pdwk   = pdsa + ( p )
c     pdwk  : pointer for KATAVEP workspace, so p*(4*p + 26) more
      pdro   = pdwk + ( p*(4*p + 26) )
c     pdro   : pointer for the robust quadratic part, so ( p*p ) more
      pdlin  = pdro + ( p*p )
c     pdlin  : pointer for linear part vector, so ( p ) more
      pdwq   = pdlin + ( p )
c     pdwq  : pointer for QP workspace who needs
c                ( p*p + 6*p + 2*nintot)
c             so ( p*p + 6*p + 2*(nin + p) )
c             so ( p*p + 8*p + 2*nin )
      pdlagr  = pdwq + ( p*p + 6*p + 2*nin )
      IF ((alpha .GT. eps).AND.(test .EQ. 1))
     &  pdlagr  = pdwq + ( p*p + 8*p + 2*nin )
c     pdlagr : pointer for Lagrange multipliers vector
c              ( p + nintot + neqtot)
      pdbcs   = pdlagr + ( p + nin + neq )
      IF ((alpha .GT. eps).AND.(test .EQ. 1))
     &   pdbcs   = pdlagr + ( 2*p + nin + neq )
c     pdbcs : pointer for the constraints vector ( nintot + neqtot ),
      pdccs   = pdbcs +  ( nin + neq)
      IF ((alpha .GT. eps).AND.(test .EQ. 1))
     &   pdccs   = pdbcs +  ( p + nin + neq)
c     pdccs : pointer for the constraints matrix
c              p*( nintot + neqtot ),
      pdcinf  = pdccs + ( p*(nin + neq) )
      IF ((alpha .GT. eps).AND.(test .EQ. 1))
     &   pdcinf  = pdccs + ( p*(p + nin + neq) )
c     pdcinf : pointer for the cinf optimal point, so ( p ) more
      pdcsup  = pdcinf + ( p )
c     pdcsup : pointer for the csup optimal point, so ( p ) more
c
c                      = p
c                      + p*p
c                      + p
c                      + p
c                      + p
c                      + P*(4*p + 26)
c                      + p*p
c                      + p
c                      + p*( p + 8 ) + 2*nin
c                      + 2*p + nin + neq
c                      + p + nin + neq
c                      + p*( p + nin + neq )
c                      + p
c                      + p
c                      = p*( 8*p + nin + neq + 44 ) + 4*nin+ + 2*neq
c
c     Total size of dwork array = p*(8*p + nin + neq + 44) + 4*nin + 2*neq
c
c--------------------------------------------------------------------
c
c     linear part and optimal weights initialization
      DO i = 1,n
        dwork(pdlin + i - 1) = epsilon
      ENDDO
c
c     case alpha = 0 (without transaction cost. p=0)
c     quadratic solver (Mean-Variance Optimization)
      IF ((alpha .LE. eps).OR.(test .EQ. 0)) THEN
c
c           copy cov. matrix
            CALL YM ( n, n, cov, dwork(pdro) )
c
c           quadratic term (matrix) SDP test
            CALL testsdp (n, dwork(pdro), iwork(piwk), dwork(pdwk),info)
            IF (info .NE. 0) THEN
                info = -108
                RETURN
            ENDIF
            CALL qp ( n, dwork(pdro), dwork(pdlin),
     &                neq, nin, ccst, bcst, cinf, csup,
     &                iwork(piwq), dwork(pdwq), dwork(pdlagr),
     &                wopt, info )
            RETURN
      ENDIF
c
c     quadratic term (matrix) SDP test
      CALL testsdp (n, cov, iwork(piwk), dwork(pdwk),info)
      IF (info .NE. 0) THEN
        info = -108
        RETURN
      ENDIF
c
c     add transaction cost. in linear part
c     p=[0, ..., 0, alpha*cost(1), ..., alpha*cost(n)]
      DO i = 1,n
        IF ((alpha*cost(i)) .LT. eps) THEN
            dwork(pdlin + n + i - 1) = alpha*eps
        ELSE
            dwork(pdlin + n + i - 1) = alpha*cost(i)
        ENDIF
      ENDDO
c
c     problem formulation: min[x'*Q*x + p'*x]
c
c     quadratic part Gamma := [cov   0             ] (cf. utalloc.f)
c                             [ 0   (alpha*cost)^2 ]
      CALL ctitcov (n, cov, alpha, cost, dwork(pdcov), info)
c
      CALL YM ( p, p, dwork(pdcov), dwork(pdro) )
c
c     quadratic term (matrix) SDP test
      CALL testsdp (p, dwork(pdro), iwork(piwk), dwork(pdwk), info)
      IF (info .NE. 0) THEN
         info = -108
         RETURN
      ENDIF
c
c     constraint modeling (cf. UTALLOC.F)
c     C*w <= b + cost transaction
      CALL ctitcst ( n, cost, neq, nin, ccst, bcst, wini,
     &               neqtot, nintot, dwork(pdccs), dwork(pdbcs), info )
c
c     lower/upper bounds
      DO i = 1,n
        dwork(pdcinf + i - 1)     = cinf(i)
        dwork(pdcinf + n + i - 1) = epsilon
        dwork(pdcsup + i - 1)     = csup(i)
        dwork(pdcsup + n + i - 1) = maxv
      ENDDO
c
c     quadratic solver (QUAPRO)
      CALL qp ( p, dwork(pdro), dwork(pdlin), neqtot, nintot,
     &          dwork(pdccs), dwork(pdbcs), dwork(pdcinf),
     &          dwork(pdcsup),
     &          iwork(piwq), dwork(pdwq), dwork(pdlagr),
     &          dwork(pdwopt), info )
      IF (info .LT. 0) RETURN
c
c     optimal portfolio
      CALL YV (n, dwork(pdwopt), wopt)
      RETURN
      END

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text: i = size (a)

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData