text:
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C   FFTPACK 5.0
C   Copyright (C) 1995-2004, Scientific Computing Division,
C   University Corporation for Atmospheric Research
C   Licensed under the GNU General Public License (GPL)
C
C   Authors:  Paul N. Swarztrauber and Richard A. Valent
C
C   $Id: cosqmb.f,v 1.2 2004/06/15 21:14:57 rodney Exp $
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE COSQMB (LOT, JUMP, N, INC, X, LENX, WSAVE, LENSAV,
     1                   WORK, LENWRK, IER)
      INTEGER    LOT, JUMP, N, INC, LENX, LENSAV, LENWRK, IER
      REAL       X(INC,*), WSAVE(LENSAV), WORK(LENWRK)
      LOGICAL    XERCON
C
      IER = 0
C
      IF (LENX .LT. (LOT-1)*JUMP + INC*(N-1) + 1) THEN
        IER = 1
        CALL XERFFT ('COSQMB', 6)
        GO TO 300
      ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))) +4) THEN
        IER = 2
        CALL XERFFT ('COSQMB', 8)
        GO TO 300
      ELSEIF (LENWRK .LT. LOT*N) THEN
        IER = 3
        CALL XERFFT ('COSQMB', 10)
        GO TO 300
      ELSEIF (.NOT. XERCON(INC,JUMP,N,LOT)) THEN
        IER = 4
        CALL XERFFT ('COSQMB', -1)
        GO TO 300
      ENDIF
C
      LJ = (LOT-1)*JUMP+1
      IF (N-2) 101,102,103
 101  DO 201 M=1,LJ,JUMP
      X(M,1) = X(M,1)
 201  CONTINUE
      RETURN
 102  SSQRT2 = 1./SQRT(2.)
      DO 202 M=1,LJ,JUMP
      X1 = X(M,1)+X(M,2)
      X(M,2) = SSQRT2*(X(M,1)-X(M,2))
      X(M,1) = X1
 202  CONTINUE
      RETURN
  103 CALL MCSQB1 (LOT,JUMP,N,INC,X,WSAVE,WORK,IER1)
      IF (IER1 .NE. 0) THEN
        IER = 20
        CALL XERFFT ('COSQMB',-5)
      ENDIF
C
  300 CONTINUE
      RETURN
      END

ext: f

max_stars_repo_name: rpepr/rpEPR

max_issues_repo_name: rpepr/rpEPR

max_forks_repo_name: rpepr/rpEPR


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

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

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

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

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

! End MIT license text.

      SUBROUTINE CRS_NONSYM_TO_CRS_SYM ( NAME_A, NROW_A, NTERM_A, I_A, J_A, A, NAME_B, NTERM_B, I_B, J_B, B )

! Transforms a square symmetric input matrix, A, that is stored in non-symmetric sparse CRS form (i.e. all nonzero terms stored)
! into output matrix B that is stored as symmetric sparse CRS form (i.e. only terms on and above the diagonal stored)

! Input matrix A is stored as:

!      I_A is an array of NROW_A+1 integers that is used to specify the number of nonzero terms in rows of matrix A. That is:
!          I_A(I+1) - I_A(I) are the number of nonzero terms in row I of matrix A

!      J_A is an integer array giving the col numbers of the NTERM_A nonzero terms in matrix A

!        A is a real array of all of the nonzero terms in matrix A.

! Output matrix B is stored as

!      I_B is an array of NROW_A+1 integers that is used to specify the number of nonzero terms in rows of matrix B. That is:
!          I_B(I+1) - I_B(I) are the number of nonzero terms in row I of matrix B

!      J_B is an integer array giving the col numbers of the nonzero terms in matrix B

!        B is a real array of the nonzero terms on and above the diagonal of input matrix A.

! The number of terms in B (i.e. NTERM_B) is related to NTERM_A and NROW_A as: NTERM_B = (NTERM_A - NROW_A)/2 + NROW_A. This assumes
! a square input matrix. The relationship is not checked herein. In addition, symmetry of the input matrix is assumed and only the
! terms on, and above, the diagonal of A are stored in B.

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, F04
      USE SCONTR, ONLY                :  BLNK_SUB_NAM
      USE TIMDAT, ONLY                :  TSEC
      USE CONSTANTS_1, ONLY           :  ZERO
      USE SUBR_BEGEND_LEVELS, ONLY    :  CRS_NONSYM_TO_CRS_SYM_BEGEND

      USE CRS_NONSYM_TO_CRS_SYM_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'CRS_NONSYM_TO_CRS_SYM'
      CHARACTER(LEN=*), INTENT(IN)    :: NAME_A            ! Name of input matrix
      CHARACTER(LEN=*), INTENT(IN)    :: NAME_B            ! Name of output matrix

      INTEGER(LONG), INTENT(IN)       :: NROW_A            ! Number of rows in input matrix, A
      INTEGER(LONG), INTENT(IN)       :: NTERM_A           ! Number of nonzero terms in input  matrix, A
      INTEGER(LONG), INTENT(IN)       :: NTERM_B           ! Number of nonzero terms in output matrix, B
      INTEGER(LONG), INTENT(IN)       :: I_A(NROW_A+1)     ! I_A(I+1) - I_A(I) are the number of nonzeros in A row I
      INTEGER(LONG), INTENT(IN)       :: J_A(NTERM_A)      ! Col numbers for nonzero terms in A
      INTEGER(LONG), INTENT(OUT)      :: I_B(NROW_A+1)     ! I_B(I+1) - I_B(I) are the num of nonzeros in B row I
      INTEGER(LONG), INTENT(OUT)      :: J_B(NTERM_B)      ! Col numbers for nonzero terms in B
      INTEGER(LONG)                   :: I,K               ! DO loop indices or counters
      INTEGER(LONG)                   :: KBEG_A            ! Index into array I_A where a row of matrix A begins
      INTEGER(LONG)                   :: KEND_A            ! Index into array I_A where a row of matrix A ends
      INTEGER(LONG)                   :: KTERM_B           ! Count of number of nonzero terms put into output matrix B
      INTEGER(LONG)                   :: A_NTERM_ROW_I     ! Number of terms in a row of input matrix A
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = CRS_NONSYM_TO_CRS_SYM_BEGEND

      REAL(DOUBLE) , INTENT(IN)       :: A(NTERM_A)        ! Real nonzero values in input  matrix A
      REAL(DOUBLE) , INTENT(OUT)      :: B(NTERM_B)        ! Real nonzero values in output matrix B

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

! **********************************************************************************************************************************
! Initialize outputs

      DO I=1,NROW_A+1
         I_B(I) = 0
      ENDDO

      DO I=1,NTERM_B
         J_B(I) = 0
           B(I) = ZERO
      ENDDO

      I_B(1)  = 1
      KTERM_B = 0
      KBEG_A  = 1
      DO I=1,NROW_A
         I_B(I+1) = I_B(I)
         A_NTERM_ROW_I = I_A(I+1) - I_A(I)
         KEND_A = KBEG_A + A_NTERM_ROW_I - 1               ! KBEG_A to KEND_A is the range of indices of terms in A for row I of A
         DO K=KBEG_A,KEND_A
            IF (J_A(K) >= I) THEN                          ! This is a term from A that is on, or above, the diagonal
               KTERM_B = KTERM_B + 1                       ! Increment the counter for the total number of terms in output matrix B
               IF (KTERM_B > NTERM_B) CALL ARRAY_SIZE_ERROR_1( SUBR_NAME, NTERM_B, NAME_B )
               I_B(I+1) = I_B(I+1) + 1                       ! Increment I_B(I+1) for the number of terms in row I of B
               J_B(KTERM_B) = J_A(K)                       ! Set column number for the term to go into B
                 B(KTERM_B) =   A(K)                       ! Set value for the term to go into B
            ENDIF
         ENDDO
         KBEG_A = KEND_A + 1
      ENDDO

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

      RETURN

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

      END SUBROUTINE CRS_NONSYM_TO_CRS_SYM


text:
      SUBROUTINE VA15AD(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,S,Y,
     *          POINT,W,IFLAG,ITER)
C
C
      DOUBLE PRECISION X(N),G(N),S(M*N),Y(M*N),DIAG(N),W(N+2*M)
      DOUBLE PRECISION FTOL,GTOL,XTOL,STPMIN,STPMAX,STP,F,YS,SQ,
     *              YR,BETA,ONE,ZERO,EPS,XNORM,GNORM,YY,DDOT,STP1
C
      INTEGER BOUND,LP,ITER,NFUN,NFEV,IPRINT(2),POINT,CP,IFLAG
      LOGICAL FINISH,DIAGCO
C      COMMON /VA15DD/MP,LP, GTOL
C MP, LP, GTOL ARE GIVEN HERE TO PREVENT COMMON IN MAIN 7/31/02 Y.HONDA
C      DATA   MP, LP /6, 6/
      COMMON /VA15DD/ MP, LP
      SAVE
      DATA ONE,ZERO/1.0D+0,0.0D+0/
C
C     ------------------------------------------------------------
C     INITIALIZE
C     ------------------------------------------------------------
C
      IF(IFLAG.EQ.0) GO TO 1
      GO TO (72,10) IFLAG
   1  ITER= 0
C MP, LP, GTOL ARE GIVEN HERE TO PREVENT COMMON IN MAIN 7/31/02 Y.HONDA
      CALL VA15CD
C
      IF(N.LE.0.OR.M.LE.0) GO TO 96
      IF(GTOL.LE.1.D-04) THEN
        IF(LP.GT.0) WRITE(LP,145)
        GTOL=1.D-02
      ENDIF
      NFUN= 1
      POINT= 0
      FINISH= .FALSE.
      IF(DIAGCO) THEN
         DO 3 I=1,N
 3       IF (DIAG(I).LE.ZERO) GO TO 95
      ELSE
         DO 4 I=1,N
 4       DIAG(I)= 1.0D0
      ENDIF
      DO 5 I=1,N
 5    S(I)= -G(I)*DIAG(I)
      GNORM= DSQRT(DDOT(N,G,1,G,1))

      STP1= ONE/GNORM
C
C     PARAMETERS FOR LINE SEARCH ROUTINE
C     ----------------------------------
      FTOL= 1.0D-4
      GTOL= 9.0D-1
      XTOL= 1.0D-17
      STPMIN= 1.0D-20
      STPMAX= 1.0D+20
      MAXFEV= 20
C
       CALL VA15BD(IPRINT,ITER,NFUN,
     *                     N,M,X,F,G,STP,FINISH)
C
C    ------------------------------------------------------------
C     MAIN ITERATION LOOP
C    --------------------------------------------------------
C
 8    ITER= ITER+1
      INFO=0
      BOUND=ITER-1
      IF (ITER .GT. M)BOUND=M
      IF(ITER.EQ.1) GO TO 65
C
C     ------------------------------------------------------------
C     COMPUTE -HG AND THE DIAGONAL SCALING MATRIX IN DIAG
C     ------------------------------------------------------------
C
      IF(.NOT.DIAGCO) THEN
CZOU         PRINT*,'form scaling ITER=', ITER
         DO 9 I=1,N
   9     DIAG(I)= YS/YY
      ELSE
         IFLAG=2
         RETURN
      ENDIF
  10  CONTINUE
      DO 11 I=1,N
  11  IF (DIAG(I).LE.ZERO) GO TO 95
C
      CP= POINT
      IF (POINT.EQ.0) CP=M
      W(N+CP)= ONE/YS
      DO 12 I=1,N
  12  W(I)= -G(I)
      CP= POINT
      DO 25 II= 1,BOUND
         CP=CP-1
         IF (CP.EQ. -1)CP=M-1
         SQ= DDOT(N,S(CP*N+1),1,W,1)
         W(N+M+CP+1)= W(N+CP+1)*SQ
         DO 20 K=1,N
  20     W(K)= W(K)-W(N+M+CP+1)*Y(CP*N+K)
  25  CONTINUE
C
      DO 30 I=1,N
  30  W(I)=DIAG(I)*W(I)
      DO 45 II=1,BOUND
         YR= DDOT(N,Y(CP*N+1),1,W,1)
         BETA= W(N+CP+1)*YR
         DO 40 K=1,N
  40     W(K)= W(K)+S(CP*N+K)*(W(N+M+CP+1)-BETA)
         CP=CP+1
         IF (CP.EQ.M)CP=0
  45  CONTINUE
C
C     ------------------------------------------------------------
C     STORE THE NEW DIRECTION IN S
C     ------------------------------------------------------------
C
       DO 60 J=1,N
  60   S(POINT*N+J)= W(J)
C
C     ------------------------------------------------------------
C     OBTAIN THE MINIMIZER OF THE FUNCTION ALONG THE
C     DIRECTION S BY USING THE LINE SEARCH ROUTINE OF VD05AD
C     ------------------------------------------------------------
  65  NFEV=0
      STP=ONE
      IF (ITER.EQ.1) STP=STP1
      DO 70 I=1,N
  70  W(I)=G(I)
  72  CONTINUE
C
      CALL VD05AD(N,X,F,G,S(POINT*N+1),STP,FTOL,GTOL,
     *            XTOL,STPMIN,STPMAX,MAXFEV,INFO,NFEV,DIAG)
C
      IF (INFO .EQ. -1) THEN
        IFLAG=1
        RETURN
      ENDIF
      IF (INFO .NE. 1) GO TO 90
      NFUN= NFUN + NFEV
C
C     ------------------------------------------------------------
C     COMPUTE THE NEW S AND Y
C     ------------------------------------------------------------
C
      NPT=POINT*N
      DO 75 I=1,N
      S(NPT+I)= STP*S(NPT+I)
  75  Y(NPT+I)= G(I)-W(I)
      YS= DDOT(N,Y(NPT+1),1,S(NPT+1),1)
      YY= DDOT(N,Y(NPT+1),1,Y(NPT+1),1)
      POINT=POINT+1
      IF (POINT.EQ.M)POINT=0
C
C     ------------------------------------------------------------
C     CONVERGENCE CHECK
C     ------------------------------------------------------------
C
      GNORM= DDOT(N,G,1,G,1)
      GNORM=DSQRT(GNORM)
      XNORM= DDOT(N,X,1,X,1)
      XNORM=DSQRT(XNORM)
      XNORM= DMAX1(1.0D0,XNORM)
C
CZOU start
C      XIAO=GNORM/XNORM
C      PRINT*,'GNORM/(1, XNORM)=', XIAO
CZOU end
      IF (GNORM/XNORM .LE. EPS) FINISH=.TRUE.
C
      CALL VA15BD(IPRINT,ITER,NFUN,
     *               N,M,X,F,G,STP,FINISH)
      IF (FINISH) THEN
         IFLAG=0
         RETURN
      ENDIF
      GO TO 8
C
C     ------------------------------------------------------------
C     END OF MAIN ITERATION LOOP. ERROR EXITS.
C     ------------------------------------------------------------
C
  90  IF(LP.LE.0) RETURN
      IF (INFO.EQ.0) THEN
           IFLAG= -1
           WRITE(LP,100)IFLAG
      ELSE IF (INFO.EQ.2) THEN
           IFLAG= -2
           WRITE(LP,105)IFLAG
      ELSE IF (INFO.EQ.3) THEN
           IFLAG= -3
           WRITE(LP,110)IFLAG
      ELSE IF (INFO.EQ.4) THEN
           IFLAG= -4
           WRITE(LP,115)IFLAG
      ELSE IF (INFO.EQ.5) THEN
           IFLAG= -5
           WRITE(LP,120)IFLAG
      ELSE IF (INFO.EQ.6) THEN
           IFLAG= -6
           WRITE(LP,125)IFLAG
      ENDIF
      RETURN
C
  95  IFLAG= -7
      IF(LP.GT.0) WRITE(LP,135)IFLAG,I
      RETURN
  96  IFLAG= -8
      IF(LP.GT.0) WRITE(LP,140)IFLAG
C
C     ------------------------------------------------------------
C     FORMATS
C     ------------------------------------------------------------
C
 100  FORMAT(/' IFLAG= ',I2,/' IMPROPER INPUT PARAMETERS DURING',
     .       ' THE LINE SEARCH.')
 105  FORMAT(/' IFLAG= ',I2,/' RELATIVE WIDTH OF THE INTERVAL OF',
     .       ' UNCERTAINTY IN THE LINE SEARCH'/'IS OF THE ORDER OF
     .          MACHINE ROUNDOFF.')
 110  FORMAT(/' IFLAG= ',I2,/' NUMBER OF CALLS TO FUNCTION IN THE',
     .       ' LINE SEARCH HAS REACHED 20.')
 115  FORMAT(/' IFLAG= ',I2,/' THE STEP IN THE LINE SEARCH IS',
     .       ' TOO SMALL.')
 120  FORMAT(/' IFLAG= ',I2,/' THE STEP IN THE LINE SEARCH IS',
     .       ' TOO LARGE.')
 125  FORMAT(/' IFLAG= ',I2,/' ROUNDING ERRORS PREVENT FURTHER',
     .       ' PROGRESS IN THE LINE SEARCH.')
 135  FORMAT(/' IFLAG= ',I2,/' THE',I5,'-TH DIAGONAL ELEMENT OF THE',
     .       ' INVERSE HESSIAN APPROXIMATION IS NOT POSITIVE')
 140  FORMAT(/' IFLAG= ',I2,/' IMPROPER INPUT PARAMETERS (N OR M',
     .       ' ARE NOT POSITIVE)')
 145  FORMAT(/'  GTOL IS LESS THAN OR EQUAL TO 1.D-04',
     .       / 'IT HAS BEEN RESET TO 1.D-02')
      RETURN
      END
C
C
C
C
      SUBROUTINE VA15BD(IPRINT,ITER,NFUN,
     *                     N,M,X,F,G,STP,FINISH)
C
C     ------------------------------------------------------------------
C     THIS ROUTINE PRINTS MONITORING INFORMATION. THE FREQUENCY AND AMOU
C     OF OUTPUT ARE SPECIFIED AS FOLLOWS:
C
C     IPRINT(1) < 0 : NO OUTPUT IS GENERATED
C     IPRINT(1) = 0 : OUTPUT ONLY AT FIRST AND LAST ITERATION
C     IPRINT(1) > 0 : OUTPUT EVERY IPRINT(1) ITERATION
C     IPRINT(2) = 0 : ITERATION COUNT, FUNCTION VALUE, NORM OF THE GRADI
C                     ,NUMBER OF FUNCTION CALLS AND STEP LENGTH
C     IPRINT(2) = 1 : + VECTOR OF VARIABLES AND GRADIENT VECTOR AT THE
C                       INITIAL POINT
C     IPRINT(2) = 2 : + VECTOR OF VARIABLES
C     IPRINT(2) = 3 : + GRADIENT VECTOR
C     ------------------------------------------------------------------
C
      DOUBLE PRECISION X(N),G(N),F,GNORM,STP,FACTOR,DDOT,GTOL
      INTEGER IPRINT(2),ITER,NFUN,PROB,LP
      LOGICAL FINISH
      COMMON /SET/ FACTOR,PROB
      COMMON /VA15DD/MP,LP, GTOL
C
      IF (IPRINT(1).LT.0)RETURN
      GNORM= DDOT(N,G,1,G,1)
      GNORM= DSQRT(GNORM)

      IF (ITER.EQ.0)THEN
           WRITE(MP,10)
           WRITE(MP,20) PROB,N,M
           WRITE(MP,30)F,GNORM
                 IF (IPRINT(2).GE.1)THEN
                     WRITE(MP,40)
                     WRITE(MP,50) (X(I),I=1,N)
                     WRITE(MP,60)
                     WRITE(MP,50) (G(I),I=1,N)
                  ENDIF
           WRITE(MP,10)
           WRITE(MP,70)
      ELSE
          IF ((IPRINT(1).EQ.0).AND.(ITER.NE.1.AND..NOT.FINISH))RETURN
              IF (IPRINT(1).NE.0)THEN
                   IF(MOD(ITER-1,IPRINT(1)).EQ.0.OR.FINISH)THEN
                         WRITE(MP,80)ITER,NFUN,F,GNORM,STP
                   ELSE
                         RETURN
                   ENDIF
              ELSE
                   WRITE(MP,80)ITER,NFUN,F,GNORM,STP
              ENDIF
              IF (IPRINT(2).EQ.2.OR.IPRINT(2).EQ.3)THEN
                    IF (FINISH)THEN
                        WRITE(MP,90)
               WRITE(MP,50) (X(I),I=1,N)
        WRITE(MP,60)
        WRITE(MP,50) (G(I),I=1,N)
                    ENDIF
              ENDIF
            IF (FINISH) WRITE(MP,100)
      ENDIF
C
C10   FORMAT('*************************************************')
 10   FORMAT(/'*************************************************')
C20   FORMAT(' PROB=',I3,'   N=',I5,'   NUMBER OF CORRECTIONS=',I2)
 20   FORMAT(' PROB=',I3,'   N=',I7,'   NUMBER OF CORRECTIONS=',I2)
 30   FORMAT('    ',E15.8,'          ',E15.8)
 40   FORMAT(' VECTOR X= ')
 50   FORMAT(6(2X,1E15.8))
 60   FORMAT(' GRADIENT VECTOR G= ')
C70   FORMAT(/'I    NFN',4X,'FUNC',8X,'GNORM',7X,'STEPLENGTH'/)
 70   FORMAT(/3X,'I',2X,'NFN',11X,'FUNC',12X,'GNORM',9X,'STEPLENGTH'/)
 80   FORMAT(2(I4,1X),3X,3(E15.8,2X))
 90   FORMAT(' FINAL POINT X= ')
 100  FORMAT(/' THE MINIMIZATION TERMINATED WITHOUT DETECTING ERRORS.',
     .       /' IFLAG = 0')
C
      RETURN
      END
C
C   ----------------------------------------------------------
C   DATA BLOCK
C   ----------------------------------------------------------
C
CTT   BLOCK DATA VA15CD
      SUBROUTINE VA15CD
      COMMON /VA15DD/MP,LP, GTOL
      INTEGER LP
      DOUBLE PRECISION GTOL
CTT   DATA MP,LP,GTOL/6,6,9.0D-01/
c     MP   = 6
c     LP   = 6
c     GTOL = 9.0D-01
      MP   = 66
      LP   = 66
      GTOL = 9.0D-01
CTS
      RETURN
CTE
      END
C
C   -------------------------------------------------------------
C
      SUBROUTINE VD05AD(N,X,F,G,S,STP,FTOL,GTOL,XTOL,
     *                  STPMIN,STPMAX,MAXFEV,INFO,NFEV,WA)
      INTEGER N,MAXFEV,INFO,NFEV
      DOUBLE PRECISION F,STP,FTOL,GTOL,XTOL,STPMIN,STPMAX
      DOUBLE PRECISION X(N),G(N),S(N),WA(N)
      SAVE
C     **********
C
C     SUBROUTINE VD05AD
C
C     THE PURPOSE OF VD05AD IS TO FIND A STEP WHICH SATISFIES
C     A SUFFICIENT DECREASE CONDITION AND A CURVATURE CONDITION.
C     THE USER MUST PROVIDE A SUBROUTINE WHICH CALCULATES THE
C     FUNCTION AND THE GRADIENT.
C
C     AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF
C     UNCERTAINTY WITH ENDPOINTS STX AND STY. THE INTERVAL OF
C     UNCERTAINTY IS INITIALLY CHOSEN SO THAT IT CONTAINS A
C     MINIMIZER OF THE MODIFIED FUNCTION
C
C          F(X+STP*S) - F(X) - FTOL*STP*(GRADF(X)'S).
C
C     IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION
C     HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE DERIVATIVE,
C     THEN THE INTERVAL OF UNCERTAINTY IS CHOSEN SO THAT IT
C     CONTAINS A MINIMIZER OF F(X+STP*S).
C
C     THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES
C     THE SUFFICIENT DECREASE CONDITION
C
C           F(X+STP*S) .LE. F(X) + FTOL*STP*(GRADF(X)'S),
C
C     AND THE CURVATURE CONDITION
C
C           ABS(GRADF(X+STP*S)'S)) .LE. GTOL*ABS(GRADF(X)'S).
C
C     IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION
C     IS BOUNDED BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES
C     BOTH CONDITIONS. IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH
C     CONDITIONS, THEN THE ALGORITHM USUALLY STOPS WHEN ROUNDING
C     ERRORS PREVENT FURTHER PROGRESS. IN THIS CASE STP ONLY
C     SATISFIES THE SUFFICIENT DECREASE CONDITION.
C
C     THE SUBROUTINE STATEMENT IS
C
C        SUBROUTINE VD05AD(N,X,F,G,S,STP,FTOL,GTOL,XTOL,
C                          STPMIN,STPMAX,MAXFEV,INFO,NFEV,WA)
C     WHERE
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF VARIABLES.
C
C       X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
C         BASE POINT FOR THE LINE SEARCH. ON OUTPUT IT CONTAINS
C         X + STP*S.
C
C       F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F
C         AT X. ON OUTPUT IT CONTAINS THE VALUE OF F AT X + STP*S.
C
C       G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
C         GRADIENT OF F AT X. ON OUTPUT IT CONTAINS THE GRADIENT
C         OF F AT X + STP*S.
C
C       S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE
C         SEARCH DIRECTION.
C
C       STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN
C         INITIAL ESTIMATE OF A SATISFACTORY STEP. ON OUTPUT
C         STP CONTAINS THE FINAL ESTIMATE.
C
C       FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. TERMINATION
C         OCCURS WHEN THE SUFFICIENT DECREASE CONDITION AND THE
C         DIRECTIONAL DERIVATIVE CONDITION ARE SATISFIED.
C
C       XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS
C         WHEN THE RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
C         IS AT MOST XTOL.
C
C       STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH
C         SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP.
C
C       MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION
C         OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST
C         MAXFEV BY THE END OF AN ITERATION.
C
C       INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
C
C         INFO = 0  IMPROPER INPUT PARAMETERS.
C
C         INFO =-1  A RETURN IS MADE TO COMPUTE THE FUNCTION AND GRADIEN
C
C         INFO = 1  THE SUFFICIENT DECREASE CONDITION AND THE
C                   DIRECTIONAL DERIVATIVE CONDITION HOLD.
C
C         INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
C                   IS AT MOST XTOL.
C
C         INFO = 3  NUMBER OF CALLS TO FCN HAS REACHED MAXFEV.
C
C         INFO = 4  THE STEP IS AT THE LOWER BOUND STPMIN.
C
C         INFO = 5  THE STEP IS AT THE UPPER BOUND STPMAX.
C
C         INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS.
C                   THERE MAY NOT BE A STEP WHICH SATISFIES THE
C                   SUFFICIENT DECREASE AND CURVATURE CONDITIONS.
C                   TOLERANCES MAY BE TOO SMALL.
C
C       NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF
C         CALLS TO FCN.
C
C       WA IS A WORK ARRAY OF LENGTH N.
C
C     SUBPROGRAMS CALLED
C
C       HARWELL-SUPPLIED...VD05BD
C
C       FORTRAN-SUPPLIED...ABS,MAX,MIN
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
C     JORGE J. MORE', DAVID J. THUENTE
C
C     **********
      INTEGER INFOC,J
      LOGICAL BRACKT,STAGE1
      DOUBLE PRECISION DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM,
     *       FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY,
     *       STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,ZERO
      DATA P5,P66,XTRAPF,ZERO /0.5D0,0.66D0,4.0D0,0.0D0/
      IF(INFO.EQ.-1) GO TO 45
      INFOC = 1
C
C     CHECK THE INPUT PARAMETERS FOR ERRORS.
C
      IF (N .LE. 0 .OR. STP .LE. ZERO .OR. FTOL .LT. ZERO .OR.
     *    GTOL .LT. ZERO .OR. XTOL .LT. ZERO .OR. STPMIN .LT. ZERO
     *    .OR. STPMAX .LT. STPMIN .OR. MAXFEV .LE. 0) RETURN
C
C     COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION
C     AND CHECK THAT S IS A DESCENT DIRECTION.
C
      DGINIT = ZERO
      DO 10 J = 1, N
         DGINIT = DGINIT + G(J)*S(J)
   10    CONTINUE
      IF (DGINIT .GE. ZERO) RETURN
C
C     INITIALIZE LOCAL VARIABLES.
C
      BRACKT = .FALSE.
      STAGE1 = .TRUE.
      NFEV = 0
      FINIT = F
      DGTEST = FTOL*DGINIT
      WIDTH = STPMAX - STPMIN
      WIDTH1 = WIDTH/P5
      DO 20 J = 1, N
         WA(J) = X(J)
   20    CONTINUE
C
C     THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP,
C     FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP.
C     THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP,
C     FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF
C     THE INTERVAL OF UNCERTAINTY.
C     THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP,
C     FUNCTION, AND DERIVATIVE AT THE CURRENT STEP.
C
      STX = ZERO
      FX = FINIT
      DGX = DGINIT
      STY = ZERO
      FY = FINIT
      DGY = DGINIT
C
C     START OF ITERATION.
C
   30 CONTINUE
C
C        SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND
C        TO THE PRESENT INTERVAL OF UNCERTAINTY.
C
         IF (BRACKT) THEN
            STMIN = MIN(STX,STY)
            STMAX = MAX(STX,STY)
         ELSE
            STMIN = STX
            STMAX = STP + XTRAPF*(STP - STX)
            END IF
C
C        FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN.
C
         STP = MAX(STP,STPMIN)
         STP = MIN(STP,STPMAX)
C
C        IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET
C        STP BE THE LOWEST POINT OBTAINED SO FAR.
C
         IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX))
     *      .OR. NFEV .GE. MAXFEV-1 .OR. INFOC .EQ. 0
     *      .OR. (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX)) STP = STX
C
C        EVALUATE THE FUNCTION AND GRADIENT AT STP
C        AND COMPUTE THE DIRECTIONAL DERIVATIVE.
C
         DO 40 J = 1, N
            X(J) = WA(J) + STP*S(J)
   40       CONTINUE
         INFO=-1
         RETURN
C
   45    INFO=0
         NFEV = NFEV + 1
         DG = ZERO
         DO 50 J = 1, N
            DG = DG + G(J)*S(J)
   50       CONTINUE
         FTEST1 = FINIT + STP*DGTEST
C
C        TEST FOR CONVERGENCE.
C
         IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX))
     *      .OR. INFOC .EQ. 0) INFO = 6
         IF (STP .EQ. STPMAX .AND.
     *       F .LE. FTEST1 .AND. DG .LE. DGTEST) INFO = 5
         IF (STP .EQ. STPMIN .AND.
     *       (F .GT. FTEST1 .OR. DG .GE. DGTEST)) INFO = 4
         IF (NFEV .GE. MAXFEV) INFO = 3
         IF (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX) INFO = 2
         IF (F .LE. FTEST1 .AND. ABS(DG) .LE. GTOL*(-DGINIT)) INFO = 1
C
C        CHECK FOR TERMINATION.
C
         IF (INFO .NE. 0) RETURN
C
C        IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED
C        FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE.
C
         IF (STAGE1 .AND. F .LE. FTEST1 .AND.
     *       DG .GE. MIN(FTOL,GTOL)*DGINIT) STAGE1 = .FALSE.
C
C        A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF
C        WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED
C        FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE
C        DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN
C        OBTAINED BUT THE DECREASE IS NOT SUFFICIENT.
C
         IF (STAGE1 .AND. F .LE. FX .AND. F .GT. FTEST1) THEN
C
C           DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES.
C
            FM = F - STP*DGTEST
            FXM = FX - STX*DGTEST
            FYM = FY - STY*DGTEST
            DGM = DG - DGTEST
            DGXM = DGX - DGTEST
            DGYM = DGY - DGTEST
C
C           CALL CSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
C           AND TO COMPUTE THE NEW STEP.
C
            CALL VD05BD(STX,FXM,DGXM,STY,FYM,DGYM,STP,FM,DGM,
     *                 BRACKT,STMIN,STMAX,INFOC)
C
C           RESET THE FUNCTION AND GRADIENT VALUES FOR F.
C
            FX = FXM + STX*DGTEST
            FY = FYM + STY*DGTEST
            DGX = DGXM + DGTEST
            DGY = DGYM + DGTEST
         ELSE
C
C           CALL VD05BD TO UPDATE THE INTERVAL OF UNCERTAINTY
C           AND TO COMPUTE THE NEW STEP.
C
            CALL VD05BD(STX,FX,DGX,STY,FY,DGY,STP,F,DG,
     *                 BRACKT,STMIN,STMAX,INFOC)
            END IF
C
C        FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE
C        INTERVAL OF UNCERTAINTY.
C
         IF (BRACKT) THEN
            IF (ABS(STY-STX) .GE. P66*WIDTH1)
     *         STP = STX + P5*(STY - STX)
            WIDTH1 = WIDTH
            WIDTH = ABS(STY-STX)
            END IF
C
C        END OF ITERATION.
C
         GO TO 30
C
C     LAST CARD OF SUBROUTINE VD05AD.
C
      END
      SUBROUTINE VD05BD(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,
     *                 STPMIN,STPMAX,INFO)
      INTEGER INFO
      DOUBLE PRECISION STX,FX,DX,STY,FY,DY,STP,FP,DP,STPMIN,STPMAX
      LOGICAL BRACKT,BOUND
C     **********
C
C     SUBROUTINE VD05BD
C
C     THE PURPOSE OF VD05BD IS TO COMPUTE A SAFEGUARDED STEP FOR
C     A LINESEARCH AND TO UPDATE AN INTERVAL OF UNCERTAINTY FOR
C     A MINIMIZER OF THE FUNCTION.
C
C     THE PARAMETER STX CONTAINS THE STEP WITH THE LEAST FUNCTION
C     VALUE. THE PARAMETER STP CONTAINS THE CURRENT STEP. IT IS
C     ASSUMED THAT THE DERIVATIVE AT STX IS NEGATIVE IN THE
C     DIRECTION OF THE STEP. IF BRACKT IS SET TRUE THEN A
C     MINIMIZER HAS BEEN BRACKETED IN AN INTERVAL OF UNCERTAINTY
C     WITH ENDPOINTS STX AND STY.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE VD05BD(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,
C                        STPMIN,STPMAX,INFO)
C
C     WHERE
C
C       STX, FX, AND DX ARE VARIABLES WHICH SPECIFY THE STEP,
C         THE FUNCTION, AND THE DERIVATIVE AT THE BEST STEP OBTAINED
C         SO FAR. THE DERIVATIVE MUST BE NEGATIVE IN THE DIRECTION
C         OF THE STEP, THAT IS, DX AND STP-STX MUST HAVE OPPOSITE
C         SIGNS. ON OUTPUT THESE PARAMETERS ARE UPDATED APPROPRIATELY.
C
C       STY, FY, AND DY ARE VARIABLES WHICH SPECIFY THE STEP,
C         THE FUNCTION, AND THE DERIVATIVE AT THE OTHER ENDPOINT OF
C         THE INTERVAL OF UNCERTAINTY. ON OUTPUT THESE PARAMETERS ARE
C         UPDATED APPROPRIATELY.
C
C       STP, FP, AND DP ARE VARIABLES WHICH SPECIFY THE STEP,
C         THE FUNCTION, AND THE DERIVATIVE AT THE CURRENT STEP.
C         IF BRACKT IS SET TRUE THEN ON INPUT STP MUST BE
C         BETWEEN STX AND STY. ON OUTPUT STP IS SET TO THE NEW STEP.
C
C       BRACKT IS A LOGICAL VARIABLE WHICH SPECIFIES IF A MINIMIZER
C         HAS BEEN BRACKETED. IF THE MINIMIZER HAS NOT BEEN BRACKETED
C         THEN ON INPUT BRACKT MUST BE SET FALSE. IF THE MINIMIZER
C         IS BRACKETED THEN ON OUTPUT BRACKT IS SET TRUE.
C
C       STPMIN AND STPMAX ARE INPUT VARIABLES WHICH SPECIFY LOWER
C         AND UPPER BOUNDS FOR THE STEP.
C
C       INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
C         IF INFO = 1,2,3,4,5, THEN THE STEP HAS BEEN COMPUTED
C         ACCORDING TO ONE OF THE FIVE CASES BELOW. OTHERWISE
C         INFO = 0, AND THIS INDICATES IMPROPER INPUT PARAMETERS.
C
C     SUBPROGRAMS CALLED
C
C       FORTRAN-SUPPLIED ... ABS,MAX,MIN,SQRT
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
C     JORGE J. MORE', DAVID J. THUENTE
C
C     **********
      DOUBLE PRECISION GAMMA,P,Q,R,S,SGND,STPC,STPF,STPQ,THETA
      INFO = 0
C
C     CHECK THE INPUT PARAMETERS FOR ERRORS.
C
      IF ((BRACKT .AND. (STP .LE. MIN(STX,STY) .OR.
     *    STP .GE. MAX(STX,STY))) .OR.
     *    DX*(STP-STX) .GE. 0.0 .OR. STPMAX .LT. STPMIN) RETURN
C
C     DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN.
C
      SGND = DP*(DX/ABS(DX))
C
C     FIRST CASE. A HIGHER FUNCTION VALUE.
C     THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER
C     TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN,
C     ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN.
C
      IF (FP .GT. FX) THEN
         INFO = 1
         BOUND = .TRUE.
         THETA = 3*(FX - FP)/(STP - STX) + DX + DP
         S = MAX(ABS(THETA),ABS(DX),ABS(DP))
         GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S))
         IF (STP .LT. STX) GAMMA = -GAMMA
         P = (GAMMA - DX) + THETA
         Q = ((GAMMA - DX) + GAMMA) + DP
         R = P/Q
         STPC = STX + R*(STP - STX)
         STPQ = STX + ((DX/((FX-FP)/(STP-STX)+DX))/2)*(STP - STX)
         IF (ABS(STPC-STX) .LT. ABS(STPQ-STX)) THEN
            STPF = STPC
         ELSE
           STPF = STPC + (STPQ - STPC)/2
           END IF
         BRACKT = .TRUE.
C
C     SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF
C     OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC
C     STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP,
C     THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN.
C
      ELSE IF (SGND .LT. 0.0) THEN
         INFO = 2
         BOUND = .FALSE.
         THETA = 3*(FX - FP)/(STP - STX) + DX + DP
         S = MAX(ABS(THETA),ABS(DX),ABS(DP))
         GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S))
         IF (STP .GT. STX) GAMMA = -GAMMA
         P = (GAMMA - DP) + THETA
         Q = ((GAMMA - DP) + GAMMA) + DX
         R = P/Q
         STPC = STP + R*(STX - STP)
         STPQ = STP + (DP/(DP-DX))*(STX - STP)
         IF (ABS(STPC-STP) .GT. ABS(STPQ-STP)) THEN
            STPF = STPC
         ELSE
            STPF = STPQ
            END IF
         BRACKT = .TRUE.
C
C     THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
C     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES.
C     THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY
C     IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC
C     IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE
C     EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO
C     COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP
C     CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN.
C
      ELSE IF (ABS(DP) .LT. ABS(DX)) THEN
         INFO = 3
         BOUND = .TRUE.
         THETA = 3*(FX - FP)/(STP - STX) + DX + DP
         S = MAX(ABS(THETA),ABS(DX),ABS(DP))
C
C        THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND
C        TO INFINITY IN THE DIRECTION OF THE STEP.
C
         GAMMA = S*SQRT(MAX(0.0D0,(THETA/S)**2 - (DX/S)*(DP/S)))
         IF (STP .GT. STX) GAMMA = -GAMMA
         P = (GAMMA - DP) + THETA
         Q = (GAMMA + (DX - DP)) + GAMMA
         R = P/Q
         IF (R .LT. 0.0 .AND. GAMMA .NE. 0.0) THEN
            STPC = STP + R*(STX - STP)
         ELSE IF (STP .GT. STX) THEN
            STPC = STPMAX
         ELSE
            STPC = STPMIN
            END IF
         STPQ = STP + (DP/(DP-DX))*(STX - STP)
         IF (BRACKT) THEN
            IF (ABS(STP-STPC) .LT. ABS(STP-STPQ)) THEN
               STPF = STPC
            ELSE
               STPF = STPQ
               END IF
         ELSE
            IF (ABS(STP-STPC) .GT. ABS(STP-STPQ)) THEN
               STPF = STPC
            ELSE
               STPF = STPQ
               END IF
            END IF
C
C     FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
C     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES
C     NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP
C     IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN.
C
      ELSE
         INFO = 4
         BOUND = .FALSE.
         IF (BRACKT) THEN
            THETA = 3*(FP - FY)/(STY - STP) + DY + DP
            S = MAX(ABS(THETA),ABS(DY),ABS(DP))
            GAMMA = S*SQRT((THETA/S)**2 - (DY/S)*(DP/S))
            IF (STP .GT. STY) GAMMA = -GAMMA
            P = (GAMMA - DP) + THETA
            Q = ((GAMMA - DP) + GAMMA) + DY
            R = P/Q
            STPC = STP + R*(STY - STP)
            STPF = STPC
         ELSE IF (STP .GT. STX) THEN
            STPF = STPMAX
         ELSE
            STPF = STPMIN
            END IF
         END IF
C
C     UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT
C     DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE.
C
      IF (FP .GT. FX) THEN
         STY = STP
         FY = FP
         DY = DP
      ELSE
         IF (SGND .LT. 0.0) THEN
            STY = STX
            FY = FX
            DY = DX
            END IF
         STX = STP
         FX = FP
         DX = DP
         END IF
C
C     COMPUTE THE NEW STEP AND SAFEGUARD IT.
C
      STPF = MIN(STPMAX,STPF)
      STPF = MAX(STPMIN,STPF)
      STP = STPF
      IF (BRACKT .AND. BOUND) THEN
         IF (STY .GT. STX) THEN
            STP = MIN(STX+0.66*(STY-STX),STP)
         ELSE
            STP = MAX(STX+0.66*(STY-STX),STP)
            END IF
         END IF
      RETURN
C
C     LAST CARD OF SUBROUTINE VD05BD.
C
      END

       DOUBLE PRECISION FUNCTION DDOT(N,D,I1,S,I2)
C
C      -------------------------------------------------------
C      THIS FUNCTION COMPUTES THE INNER PRODUCT OF TWO VECTORS
C      -------------------------------------------------------
C
       DOUBLE PRECISION D(N),S(N),PROD
       INTEGER I1,I2
C
        PROD=0.0D0
        DO 10 I=1,N
 10     PROD= PROD+D(I)*S(I)
C
        DDOT= PROD
C
       RETURN
       END

max_stars_repo_name: skotsuki/speedy-lpf

max_issues_repo_name: skotsuki/speedy-lpf

max_forks_repo_name: skotsuki/speedy-lpf


text:
!==============================================================
! Copyright © 2020 Intel Corporation
!
! SPDX-License-Identifier: MIT
! =============================================================
program main
    use omp_lib
    integer, parameter :: N=16
    integer :: i
    integer, allocatable :: x(:)
    logical :: is_cpu = .true.

    !$omp allocate allocator(omp_target_shared_mem_alloc)
    allocate(x(N))

    do i=1,N
        x(i) = i
    end do

    !$omp target map(tofrom: is_cpu) has_device_addr(x)
    !$omp teams distribute parallel do
    do i=1,N
        if ((i==1) .and. (.not.(omp_is_initial_device()))) is_cpu=.false.
        x(i) = x(i) * 2
    end do
    !$omp end target

    if (is_cpu) then
        print *, "Running on CPU"
    else
        print *, "Running on GPU"
    end if

    do i=1,N
        print *, x(i)
    end do

    deallocate(x)
end program main

max_stars_repo_name: yafshar/oneAPI-samples

max_issues_repo_name: yafshar/oneAPI-samples

max_forks_repo_name: yafshar/oneAPI-samples


text:
Travel N Tails Departed Businesses closed in August 2005, and has been replaced by Just KIDding. Previously, this spot held Accents, with the same phone number.

This was a really fun store for an animal lover or anybody who needs to buy a gift for an animal lover. With an excellent selection of Dogs dog related gifts, you could get everything from socks to statues with any given breed of dog. Plenty of dog and cat beds, bowls and other pet things were available.

In a slightly odd twist, they also had a selection of travel related items like luggage tags, travel wallets, small bags, etc.

20041214 22:18:16 nbsp Travel N Tails sells my favoriteLiver biscotti.  Yummmmm....  Users/SheilaTheDog

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
Resident of Davis for 6 years: chose the house I live in because its walking distance to Covell Nugget.

I live here with my lovely, smart U.C.D. grad wife, and two dogs. I am a landscape designer,
local food advocate and wannabe vineyard owner.

I think Davis is one of the top 10 best towns in Cal. to be trapped in during a recession.
(Sonoma, Capitola, Ojai, Pacific Beach are others, but I cant afford them)

The only flaw I can see is a big university that draws international students,
and yet we have no/few good restaurants, wine bars or live music. The Grad
does not count. The Palms left. Tucos is nice, but needs to be 4 times bigger.

City government repression of liquor licenses, townies frowning on live venues?
I cant explain.

20090914 20:52:01 nbsp Hi there. A belated welcome to the wiki. Great new page on landscape architects. I added some language on it and left it open for other LAs to include their names or businesses. Then I created a new business page for your Great Valley Design, with a minimal amount of info that can be expanded upon. Users/robinlaughlin

20090915 14:34:50 nbsp Hi,
are you a landscape architect, or a landscape designer? Id like the page to be accurate and explain the distinction, then perhaps we can add folks like Margot Anderson and others.  Users/DonShor


text:
C -*- Mode: Fortran; -*-
C
C  (C) 2004 by Argonne National Laboratory.
C      See COPYRIGHT in top-level directory.
C
      program main
      implicit none
      include 'mpif.h'
C iooffset.h provides a variable "offset" that is of type MPI_Offset
C (in Fortran 90, kind=MPI_OFFSET_KIND)
      include 'iooffset.h'
C iodisp.h declares disp as an MPI_Offset integer
      include 'iodisp.h'
      integer rank, size
      integer fh, i, group, worldgroup, result
      integer ierr, errs, toterrs
      integer BUFSIZE
      parameter (BUFSIZE=1024)
      integer buf(BUFSIZE)
      character*(50) filename
      character*(MPI_MAX_DATAREP_STRING) datarep
      integer amode
      logical atomicity
      integer newtype, etype, filetype
      integer integer_size, type_size
C
      errs = 0
      call mpi_init( ierr )
      call mpi_comm_rank( MPI_COMM_WORLD, rank, ierr )
      call mpi_comm_size( MPI_COMM_WORLD, size, ierr )
C
C Create a file that we'll then query properties
      filename = "testfile.txt"
      call mpi_file_open( MPI_COMM_WORLD, filename, MPI_MODE_CREATE +
     $     MPI_MODE_RDWR, MPI_INFO_NULL, fh, ierr )
      if (ierr .ne. MPI_SUCCESS) then
         print *, "Unable to create file ", filename
         call mpi_abort( MPI_COMM_WORLD, 1, ierr )
      endif
C
C Fill in some data
      do i=1, BUFSIZE
         buf(i) = i
      enddo
      call mpi_file_write( fh, buf, BUFSIZE, MPI_INTEGER,
     $     MPI_STATUS_IGNORE, ierr )
      call MPI_File_sync( fh, ierr )
C
C Now, query properties of the file and the file handle
      call MPI_File_get_amode(fh, amode, ierr )
      if (amode .ne. MPI_MODE_CREATE + MPI_MODE_RDWR) then
         errs = errs + 1
         print *, " Amode was different than expected"
      endif
C
      call MPI_File_get_atomicity( fh, atomicity, ierr )
      if (atomicity) then
         errs = errs + 1
         print *, " Atomicity was true but should be false"
      endif
C
      call MPI_File_set_atomicity( fh, .true., ierr )
      call MPI_File_get_atomicity( fh, atomicity, ierr )
      if (.not. atomicity) then
         errs = errs + 1
         print *, " Atomicity was set to true but ",
     $        "get_atomicity returned false"
      endif
      call MPI_File_set_atomicity( fh, .false., ierr )
C
C FIXME: original code use 10,10,20, and the following code
C assumed the original
C
C Create a vector type of 10 elements, each of 20 elements, with a stride of
C 30 elements
      call mpi_type_vector( 10, 20, 30, MPI_INTEGER, newtype, ierr )
      call mpi_type_commit( newtype, ierr )
C
C All processes are getting the same view, with a 1000 byte offset
      offset = 1000
      call mpi_file_set_view( fh, offset, MPI_INTEGER, newtype, "native"
     $     , MPI_INFO_NULL, ierr )

      call mpi_file_get_view( fh, offset, etype, filetype, datarep, ierr
     $     )
      if (offset .ne. 1000) then
         print *, " displacement was ", offset, ", expected 1000"
         errs = errs + 1
      endif
      if (datarep .ne. "native") then
         print *, " data representation form was ", datarep,
     $        ", expected native"
         errs = errs + 1
      endif

C Find the byte offset, given an offset of 20 etypes relative to the
C current view (the same as the blockcount of the filetype, which
C places it at the beginning of the next block, hence a stride
C distance away).
      offset = 20
      call mpi_file_get_byte_offset( fh, offset, disp, ierr )
      call mpi_type_size( MPI_INTEGER, integer_size, ierr )
      if (disp .ne. 1000 + 30 * integer_size) then
         errs = errs + 1
         print *, " (offset20)Byte offset = ", disp, ", should be ", 1000+20
     $        *integer_size
      endif
C
C     We should also compare file and etypes.  We just look at the
C     sizes and extents for now

      call mpi_type_size( etype, type_size, ierr )
      if (type_size .ne. integer_size) then
         print *, " Etype has size ", type_size, ", but should be ",
     $        integer_size
         errs = errs + 1
      endif
      call mpi_type_size( filetype, type_size, ierr )
      if (type_size .ne. 10*20*integer_size) then
         print *, " filetype has size ", type_size, ", but should be ",
     $        10*20*integer_size
         errs = errs + 1
      endif
C
C Only free derived type
      call mpi_type_free( filetype, ierr )

      call mpi_file_get_group( fh, group, ierr )
      call mpi_comm_group( MPI_COMM_WORLD, worldgroup, ierr )
      call mpi_group_compare( group, worldgroup, result, ierr )

      if (result .ne. MPI_IDENT) then
         print *, " Group of file does not match group of comm_world"
         errs = errs + 1
      endif
      call mpi_group_free( group, ierr )
      call mpi_group_free( worldgroup, ierr )

      offset = 1000+25*integer_size
      call mpi_file_set_size(fh, offset, ierr )
      call mpi_barrier(MPI_COMM_WORLD, ierr )
      call mpi_file_sync(fh, ierr )
      call mpi_file_get_size( fh, offset, ierr )

      if (offset .ne. 1000+25*integer_size) then
         errs = errs + 1
         print *, " File size is ", offset, ", should be ", 1000 + 25
     $        * integer_size
      endif
C
C File size is 1000+25ints.  Seek to end.  Note that the file size
C places the end of the file into the gap in the view, so seeking
C to the end, which is relative to the view, needs to give the end
C of the first block of 20 ints)
      offset = 0
      call mpi_file_seek( fh, offset, MPI_SEEK_END, ierr )
      call mpi_file_get_position( fh, disp, ierr )
      if (disp .ne. 20) then
         errs = errs + 1
         print *, "File pointer position = ", disp, ", should be 20"
      endif
      call mpi_file_get_byte_offset(fh, disp, offset, ierr )
      if (offset .ne. 1000+30*integer_size) then
         errs = errs + 1
         print *, " (seek)Byte offset = ", offset, ", should be ", 1000
     $        +30*integer_size
      endif

      call mpi_barrier(MPI_COMM_WORLD, ierr )

      offset = -20
      call mpi_file_seek(fh, offset, MPI_SEEK_CUR, ierr )
      call mpi_file_get_position(fh, disp, ierr )
      call mpi_file_get_byte_offset(fh, disp, offset, ierr )
      if (offset .ne. 1000) then
         errs = errs + 1
         print *, " File pointer position in bytes = ", offset,
     $        ", should be 1000"
      endif

      offset = 8192
      call mpi_file_preallocate(fh, offset, ierr )
      offset = 0
      call mpi_file_get_size( fh, offset, ierr )
      if (offset .lt. 8192) then
         errs = errs + 1
         print *, " Size after preallocate is ", offset,
     $        ", should be at least 8192"
      endif
      call mpi_file_close( fh, ierr )

      call mpi_barrier(MPI_COMM_WORLD, ierr )

      if (rank .eq. 0) then
         call MPI_File_delete(filename, MPI_INFO_NULL, ierr )
      endif
C
C Get error summary
      call MPI_Allreduce( errs, toterrs, 1, MPI_INTEGER, MPI_SUM,
     $     MPI_COMM_WORLD, ierr )
      if (rank .eq. 0) then
         if( toterrs .gt. 0) then
            print *, "Found ", toterrs, " errors"
         else
            print *, " No Errors"
         endif
      endif

      call mpi_type_free( newtype, ierr )

      call mpi_finalize( ierr )
      end

max_stars_repo_name: raffenet/mpich-CVS

max_issues_repo_name: raffenet/mpich-CVS

max_forks_repo_name: raffenet/mpich-CVS


text:
subroutine atom(sol,k,eval)
real, intent(in) :: sol
integer, intent(in) :: k(2)
real, intent(out) :: eval(2)
real t1
  t1=sqrt(dble(k(1)**2)-(sol)**2)
  eval(1)=sol**2/sqrt(t1)-sol**2
end subroutine

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      subroutine setvdecay(idv,vcharge)
c--- Set up particle labels and couplings for vector boson
c--- decays according to contents of vdecayid common block
c---
c--- Vector boson is p(idv) with charge vcharge
c---
c--- If vdecayid = FALSE, returns without doing anything
      implicit none
      include 'constants.f'
      include 'ewcharge.f'
      include 'plabel.f'
      include 'vdecayid.f'
      include 'zcouple.f'
      include 'nqcdjets.f'
      integer idv,vcharge
      character*2 vdecay,decay1,decay2
      character*15 decaystring
      double precision decayq,decayl,decayr

c--- return if decay not specified in input file
      if (vdecayid .eqv. .false.) return

      if     (idv .eq. 34) then
        vdecay=v34id
      elseif (idv .eq. 56) then
        vdecay=v56id
      else
        write(6,*) 'Subroutine setvdecay called improperly;'
        write(6,*) 'only idv=34, 56 allowed but idv=',idv
        stop
      endif

c--- allowed string for specifying vector boson decays,
c---  (particle, anti-particle)

c--- Z decays
      if (vcharge .eq. 0) then
      if     ((vdecay .eq. 'el') .or. (vdecay .eq. 'EL')) then
        decaystring='(e-, e+)       '
        decay1='el'
        decay2='ea'
        decayq=-1d0
        decayl=le
        decayr=re
      elseif ((vdecay .eq. 'mu') .or. (vdecay .eq. 'MU')
     &   .or. (vdecay .eq. 'ml') .or. (vdecay .eq. 'ML')) then
        decaystring='(mu-, mu+)     '
        decay1='ml'
        decay2='ma'
        decayq=-1d0
        decayl=le
        decayr=re
      elseif ((vdecay .eq. 'tl') .or. (vdecay .eq. 'TL')) then
        decaystring='(tau-, tau+)   '
        decay1='tl'
        decay2='ta'
        decayq=-1d0
        decayl=le
        decayr=re
      elseif ((vdecay .eq. 'nu') .or. (vdecay .eq. 'NU')
     &   .or. (vdecay .eq. 'nl') .or. (vdecay .eq. 'NL')) then
        decaystring='(nu, nubar) x 3'
        decay1='nl'
        decay2='na'
        decayq=0d0
        decayl=ln*dsqrt(3d0)
        decayr=rn*dsqrt(3d0)
      elseif ((vdecay .eq. 'bq') .or. (vdecay .eq. 'BQ')) then
        decaystring='(b, b-bar)     '
        decay1='bq'
        decay2='ba'
        nqcdjets=nqcdjets+2
        decayq=Q(5)*dsqrt(xn)
        decayl=l(5)*dsqrt(xn)
        decayr=r(5)*dsqrt(xn)
      else
        write(6,*) 'Decay string not recognized, vdecay=',vdecay
        stop
      endif
      endif

c--- W+ decays
      if(vcharge .eq. +1) then
      if     ((vdecay .eq. 'en') .or. (vdecay .eq. 'EN')) then
        decaystring='(ve, e+)       '
        decay1='nl'
        decay2='ea'
      elseif ((vdecay .eq. 'mn') .or. (vdecay .eq. 'MN')) then
        decaystring='(vmu, mu+)     '
        decay1='nm'
        decay2='ma'
      elseif ((vdecay .eq. 'tn') .or. (vdecay .eq. 'TN')) then
        decaystring='(vtau, tau+)   '
        decay1='nt'
        decay2='ta'
      else
        write(6,*) 'Decay string not recognized, vdecay=',vdecay
        stop
      endif
      endif

c--- W- decays
      if(vcharge .eq. -1) then
      if     ((vdecay .eq. 'en') .or. (vdecay .eq. 'EN')) then
        decaystring='(e-, ve~)      '
        decay1='el'
        decay2='na'
      elseif ((vdecay .eq. 'mn') .or. (vdecay .eq. 'MN')) then
        decaystring='(mu-, vmu~)    '
        decay1='ml'
        decay2='bm'
      elseif ((vdecay .eq. 'tn') .or. (vdecay .eq. 'TN')) then
        decaystring='(tau-, vtau~)  '
        decay1='tl'
        decay2='bt'
      else
        write(6,*) 'Decay string not recognized, vdecay=',vdecay
        stop
      endif
      endif

      if (idv .eq. 34) then
        plabel(3)=decay1
        plabel(4)=decay2
        if (vcharge .eq. 0) then
          q1=decayq
          l1=decayl
          r1=decayr
        endif
      endif

      if (idv .eq. 56) then
        plabel(5)=decay1
        plabel(6)=decay2
        if (vcharge .eq. 0) then
          q2=decayq
          l2=decayl
          r2=decayr
        endif
      endif

      write(6,*)
      write(6,*) '****************************************************'
      write(6,99) idv,decaystring
      write(6,*) '****************************************************'

      return

   99 format('*      Vector boson decay ',i2,' -> ',a15,'     *')
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
program atest
   use mytype
   use qsort_c_module
   implicit none
   integer :: numnz,msize
   integer :: k
   integer , parameter :: maxnz=100000
   integer , allocatable:: row_ptr(:)
   type(ent), allocatable:: vals(:)
! Read the number of nonzeros and the size of the matrix
   read(*,*)numnz,msize
! Allocate our vector to store the data triplets.  In
! the end this will also contain col_ind and values
! in the proper order.
! NOTE: numnz <= maxnz
   allocate(vals(maxnz))
! Read our data.  In a real program the data could be
! generated on the fly in any order.
   do k=1,numnz
        read(*,*)vals(k)%i,vals(k)%j,vals(k)%val
   enddo
! On return from QsortC vals%j contains the col_ind
! and vals%val contains the values.  Note we pass in
! only the section of the array that contains data.
   call QsortC(vals(1:numnz))
! By convention row_ptr is one bigger than the
! size of the matrix and the last entry is equal
! to numnz+1.
   allocate(row_ptr(msize+1))
! Fill row_ptr
   call getrow(vals,row_ptr,numnz)
! Write out the row_ptr
   write(*,*)row_ptr
   write(*,*)
! Write out the col_ind
   write(*,*)vals(1:numnz)%j
   write(*,*)
! Write out the values
   write(*,*)vals(1:numnz)%val
 end program

max_stars_repo_name: timkphd/examples

max_issues_repo_name: timkphd/examples

max_forks_repo_name: timkphd/examples


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 09.01.05.01/04                        *
C  *    TEST TITLE : Error 206 handled by user             *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      CALL INITGL ('09.01.05.01/04')
      CALL E05206 (1)
      CALL ENDIT
      END


module stats2
    implicit none
    save
    private
    public :: f, g, h, &
              c, dc, d2c, &
              setf, setg, seth, &
              setc, setdc, setd2c, &
              subproblem, iteration, armijo, norm, angle, &
              initstat, printheader, printstat, setmethod

    interface
        function f_t(x)
            double precision, intent(in) :: x(:)
            double precision             :: f_t
        end function f_t

        subroutine g_t(gx, x)
            double precision, intent(in)  :: x(:)
            double precision, intent(out) :: gx(:)
        end subroutine g_t

        subroutine h_t(hx, x)
            double precision, intent(in)  :: x(:)
            double precision, intent(out) :: hx(:, :)
        end subroutine h_t

        subroutine c_t(cx, x)
            double precision, intent(in)  :: x(:)
            double precision, intent(out) :: cx(:)
        end subroutine c_t

        subroutine dc_t(dcx, x)
            double precision, intent(in)  :: x(:)
            double precision, intent(out) :: dcx(:, :)
        end subroutine dc_t

        subroutine d2c_t(d2cx, x)
            double precision, intent(in)  :: x(:)
            double precision, intent(out) :: d2cx(:, :, :)
        end subroutine d2c_t
    end interface

    procedure   (f_t), pointer ::   f_ptr => null()
    procedure   (g_t), pointer ::   g_ptr => null()
    procedure   (h_t), pointer ::   h_ptr => null()
    procedure   (c_t), pointer ::   c_ptr => null()
    procedure  (dc_t), pointer ::  dc_ptr => null()
    procedure (d2c_t), pointer :: d2c_ptr => null()

    integer          :: n, m
    integer          :: nsub, nit, narm, nnor, nang
    integer          :: nf, ng, nh
    integer          :: nc, ndc, nd2c
    character(len=8) :: method

contains
    subroutine initstat(n_, m_)
        integer, intent(in) :: n_, m_

        n = n_
        m = m_
        nsub = 0
        nit  = 0
        narm = 0
        nnor = 0
        nang = 0
        nf   = 0
        nh   = 0
        ng   = 0

        method = ""
          f_ptr => null()
          g_ptr => null()
          h_ptr => null()
          c_ptr => null()
         dc_ptr => null()
        d2c_ptr => null()
    end subroutine initstat

    subroutine setf(p)
        procedure (f_t) :: p
        f_ptr => p
    end subroutine setf

    subroutine setg(p)
        procedure (g_t) :: p
        g_ptr => p
    end subroutine setg

    subroutine seth(p)
        procedure (h_t) :: p
        h_ptr => p
    end subroutine seth

    subroutine setc(p)
        procedure (c_t) :: p
        c_ptr => p
    end subroutine setc

    subroutine setdc(p)
        procedure (dc_t) :: p
        dc_ptr => p
    end subroutine setdc

    subroutine setd2c(p)
        procedure (d2c_t) :: p
        d2c_ptr => p
    end subroutine setd2c

    function f(x)
        double precision, intent(in) :: x(:)
        double precision             :: f

        nf = nf + 1

        f = f_ptr(x)
    end function f

    subroutine g(gx, x)
        double precision, intent(in)  :: x(:)
        double precision, intent(out) :: gx(:)

        ng = ng + 1

        call g_ptr(gx, x)
    end subroutine g

    subroutine h(hx, x)
        double precision, intent(in)  :: x(:)
        double precision, intent(out) :: hx(:, :)

        nh = nh + 1

        call h_ptr(hx, x)
    end subroutine h

    subroutine c(cx, x)
        double precision, intent(in)  :: x(:)
        double precision, intent(out) :: cx(:)

        nc = nc + 1

        call c_ptr(cx, x)
    end subroutine c

    subroutine dc(dcx, x)
        double precision, intent(in)  :: x(:)
        double precision, intent(out) :: dcx(:, :)

        ndc = ndc + 1

        call dc_ptr(dcx, x)
    end subroutine dc

    subroutine d2c(d2cx, x)
        double precision, intent(in)  :: x(:)
        double precision, intent(out) :: d2cx(:, :, :)

        nd2c = nd2c + 1

        call d2c_ptr(d2cx, x)
    end subroutine d2c

    subroutine subproblem()
        nsub = nsub + 1
        nit = 0
    end subroutine subproblem

    subroutine iteration()
        nit = nit + 1

        if (nit > 1000000) then
            print '(a8, a10)', method, "FC"
            write (57, *), method, " & ", "FC", " \\"
            call exit(1)
        end if
    end subroutine iteration

    subroutine armijo()
        narm = narm + 1
    end subroutine armijo

    subroutine norm()
        nnor = nnor + 1
    end subroutine norm

    subroutine angle()
        nang = nang + 1
    end subroutine angle

    subroutine printheader()
        print '(a22, a16, a10, a10, a12, a13, a10, a10, a11)', "‖c(x)‖", "‖∇L‖", "sub", "f", "∇f", "∇²f", &
                                                           "armijo", "norma", "ângulo"
    end subroutine printheader

    subroutine printstat(x, mu)
        double precision, intent(in) :: x(:)
        double precision, intent(in) :: mu

        double precision, allocatable :: cx(:)
        double precision, allocatable :: dcx(:, :)
        double precision, allocatable :: gx(:)
        double precision, allocatable :: dL(:)

        integer :: i, k

        allocate(cx(m))
        allocate(dcx(m, n))
        allocate(gx(n))
        allocate(dL(n))

        call c_ptr(cx, x)
        call dc_ptr(dcx, x)
        call g_ptr(gx, x)

        dL(:) = gx(:)
        do k = 1, m
            do i = 1, n
                dL(i) = dL(i) + cx(k)*dcx(k, i)/mu
            end do
        end do
        print '(a8, 2e10.2, 7i10)', method, norm2(cx), norm2(dL), nsub, nf, ng, nh, narm, nnor, nang
        write (57,'(2a)',advance='no'), method, " & "
        write (57,'(e10.2, a)',advance='no'), norm2(cx), " & "
        write (57,'(e10.2, a)',advance='no'), norm2(dL), " & "
        write (57,'(i10,a)',advance='no'), nsub, " & "
        write (57,'(i10,a)',advance='no'), nf, " & "
        write (57,'(i10,a)',advance='no'), ng, " & "
        write (57,'(i10,a)',advance='no'), nh, " & "
        write (57,'(i10,a)',advance='no'), narm, " & "
        write (57,'(i10,a)',advance='no'), nnor, " & "
        write (57, *), nang, " \\"

        deallocate(cx)
        deallocate(dcx)
        deallocate(gx)
        deallocate(dL)
    end subroutine printstat

    subroutine setmethod(name)
        character(len=*) :: name

        method = name
    end subroutine setmethod

end module stats2


text:
    subroutine flip_to_cf(xin,xout,nlon,nlat,nz)

        ! Copied from flip_to_cf_and_box_smooth Feb 16, 2014

        implicit none

        integer, parameter :: RKIND = selected_real_kind(12)
        integer, intent(in) :: nlon,nlat,nz
        real (kind=RKIND), intent(in), dimension(nz,nlat,nlon) :: xin
        real (kind=RKIND), intent(out), dimension(nlon,nlat,nz) :: xout

        integer :: i,j,k

        do i=1,nlon
            do j=1,nlat
                do k=1,nz
                    xout(i,j,k) = xin(k,j,i)
                end do
            end do
        end do

    end subroutine flip_to_cf

max_stars_repo_name: NCAR/MPAS-vortex-tracker

max_issues_repo_name: NCAR/MPAS-vortex-tracker

max_forks_repo_name: NCAR/MPAS-vortex-tracker


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2019 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine calcflux(area,vfa,xxna,ipnei,nef,neifa,flux,xxj,
     &  gradpfa,xlet,xle,vel,advfa,ielfa,neiel,ifabou,hfa,nefa,
     &  nefb)
!
!     calculate the mass flow using the newly calculated velocity
!     and based on the Rhie-Chow interpolation
!
!     vfa(1..3,*) is only used for external faces
!
      implicit none
!
      integer i,indexf,ipnei(*),ifa,nef,neifa(*),ielfa(4,*),
     &  neiel(*),iwall,iel,ipointer,ifabou(*),nefa,nefb
!
      real*8 area(*),vfa(0:7,*),xxna(3,*),flux(*),xxj(3,*),gradpfa(3,*),
     &  xlet(*),xle(*),vel(nef,0:7),advfa(*),coef,hfa(3,*)
!
      intent(in) area,vfa,xxna,ipnei,nef,neifa,xxj,
     &  gradpfa,xlet,xle,vel,advfa,ielfa,neiel,ifabou,hfa,nefa,
     &  nefb
!
      intent(inout) flux
!
      do i=nefa,nefb
         do indexf=ipnei(i)+1,ipnei(i+1)
            ifa=neifa(indexf)
            iel=neiel(indexf)
!
            if(iel.eq.0) then
!
!              external face
!
               iwall=0
               ipointer=abs(ielfa(2,ifa))
               if(ipointer.gt.0) then
                  iwall=ifabou(ipointer+5)
               endif
!
               if(iwall.eq.0) then
!
!                 no wall nor sliding conditions
!
                  if(ielfa(3,ifa).gt.0) then
!
!                    inlet: velocity fixed
!
                     flux(indexf)=vfa(5,ifa)*
     &                            (vfa(1,ifa)*xxna(1,indexf)+
     &                             vfa(2,ifa)*xxna(2,indexf)+
     &                             vfa(3,ifa)*xxna(3,indexf))
                  else
!
!                    outlet or 1-layer
!
                     coef=advfa(ifa)*
     &                    ((vfa(4,ifa)-vel(i,4))/xle(indexf)-
     &                     (gradpfa(1,ifa)*xxj(1,indexf)+
     &                      gradpfa(2,ifa)*xxj(2,indexf)+
     &                      gradpfa(3,ifa)*xxj(3,indexf)))
                     flux(indexf)=vfa(5,ifa)*
     &                  ((hfa(1,ifa)-coef*xxj(1,indexf))*xxna(1,indexf)+
     &                   (hfa(2,ifa)-coef*xxj(2,indexf))*xxna(2,indexf)+
     &                   (hfa(3,ifa)-coef*xxj(3,indexf))*xxna(3,indexf))
                  endif
               else
!
!                 wall or sliding: velocity of the wall
!
                  flux(indexf)=0.d0
c                  flux(indexf)=vfa(5,ifa)*
c     &                 (vfa(1,ifa)*xxna(1,indexf)+
c     &                  vfa(2,ifa)*xxna(2,indexf)+
c     &                  vfa(3,ifa)*xxna(3,indexf))
               endif
            else
!
!              internal face
!
               coef=advfa(ifa)*
     &              ((vel(iel,4)-vel(i,4))/xlet(indexf)-
     &              (gradpfa(1,ifa)*xxj(1,indexf)+
     &              gradpfa(2,ifa)*xxj(2,indexf)+
     &              gradpfa(3,ifa)*xxj(3,indexf)))
               flux(indexf)=vfa(5,ifa)*
     &              ((hfa(1,ifa)-coef*xxj(1,indexf))*xxna(1,indexf)+
     &              (hfa(2,ifa)-coef*xxj(2,indexf))*xxna(2,indexf)+
     &              (hfa(3,ifa)-coef*xxj(3,indexf))*xxna(3,indexf))
             endif
         enddo
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
module functions
  implicit none
contains

function binsearch(n,a)
  implicit none
  integer(2) :: binsearch,left,right,mid
  integer(8) :: n,a(:)

  left=1
  right=size(a)
  if ( n >= a(right) ) then
    binsearch=right
    return
  end if
  do while ( .true. )
    mid=(left+right)/2
    if ( left == mid .or. n == a(mid) ) exit
    if ( n > a(mid) ) then
      left=mid
    else
      right=mid
    end if
  end do
  binsearch=mid
  return
end function binsearch

function nchoosek(n,k)
  implicit none
  integer(2) :: k,n,i,ia,ib
  integer(8) :: nchoosek
  real(8) :: A,B

  if ( k < n-k ) then
    ia=n-k+1
    ib=k
  else
    ia=k+1
    ib=n-k
  end if
  A=1;B=1
  do i=ia,n
    A=A*i
  end do
  do i=2,ib
    B=B*i
  end do
  nchoosek=A/B+0.5
  return
end function nchoosek

function complement(a,m)
  implicit none
  integer(2) :: i,j,m
  integer(2) :: a(:),c(m)
  integer(2) :: complement(m)

  c=1
  c(a)=0
  j=0
  complement=0
  do i=1,m
    if ( c(i) /= 0 ) then
      j=j+1
      complement(j)=i
    end if
  end do
  return
end function complement

function width(i_1,i_2,i_4,i_8)
  integer(1) :: width
  integer(1),optional :: i_1
  integer(2),optional :: i_2
  integer(4),optional :: i_4
  integer(8),optional :: i_8
  character(len=20) :: w

  if ( present(i_1) ) write(w,'(I20)') i_1
  if ( present(i_2) ) write(w,'(I20)') i_2
  if ( present(i_4) ) write(w,'(I20)') i_4
  if ( present(i_8) ) write(w,'(I20)') i_8
  width=len_trim(adjustl(w))
end function width

function norm(A)
  implicit none
  real(8) :: norm
  real(8) :: A(3)

  norm=dsqrt(A(1)**2+A(2)**2+A(3)**2)
  return
end function norm

function cross(A,B)
  implicit none
  real(8) :: cross(3)
  real(8) :: A(3),B(3)

  cross(1)=A(2)*B(3)-A(3)*B(2)
  cross(2)=A(3)*B(1)-A(1)*B(3)
  cross(3)=A(1)*B(2)-A(2)*B(1)
  return
end function cross

function trace(A)
  implicit none
  integer(1) :: trace,A(:,:)

  trace=A(1,1)+A(2,2)+A(3,3)
  return
end function trace

function det(A)
  implicit none
  integer(1) :: det
  integer(1) :: A(3,3)

  det=    A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))
  det=det-A(1,2)*(A(2,1)*A(3,3)-A(2,3)*A(3,1))
  det=det+A(1,3)*(A(2,1)*A(3,2)-A(2,2)*A(3,1))
  return
end function det

function inverse(A)
  implicit none
  integer(1) :: i,j,N
  real(8) :: inverse(3,3)
  real(8) :: A(:,:)
  real(8) :: B(3,3),IA(3,3)

  N=3
  forall(i=1:N,j=1:N,i==j) IA(i,j)=1.D0
  forall(i=1:N,j=1:N,i/=j) IA(i,j)=0.D0
  B=A
  call upper(B,IA,N)
  call lower(B,IA,N)
  forall(i=1:N) IA(i,:)=IA(i,:)/B(i,i)
  inverse=IA
  return
end function inverse

subroutine upper(M,S,N)
  implicit none
  integer(1) :: N
  real(8) :: M(N,N)
  real(8) :: S(N,N)
  integer(1) :: I,J
  real(8) :: E

  do I=1,N-1
    do J=I+1,N
      E=M(J,I)/M(I,I)
      M(J,I:N)=M(J,I:N)-M(I,I:N)*E
      S(J,:)=S(J,:)-S(I,:)*E
    end do
  end do
  return
end subroutine upper

subroutine lower(M,S,N)
  implicit none
  integer(1) :: N
  real(8) :: M(N,N)
  real(8) :: S(N,N)
  integer(1) :: I,J
  real(8) :: E

  do I=N,2,-1
    do J=I-1,1,-1
      E=M(J,I)/M(I,I)
      M(J,I:N)=M(J,I:N)-M(I,I:N)*E
      S(J,:)=S(J,:)-S(I,:)*E
    end do
  end do
  return
end subroutine lower

recursive subroutine sort(A,left,right)
  implicit none
  integer(2),intent(in out) :: A(:)
  integer(2),intent(in) :: left,right
  integer(2) :: temp,t
  integer(2) :: i,j

  if ( left > right ) then
    return
  end if
  temp=A(left)
  i=left
  j=right
  do while ( i /= j )
    do while ( A(j) >= temp .and. i < j )
      j=j-1
    end do
    do while ( A(i) <= temp .and. i < j )
      i=i+1
    end do
    if ( i < j ) then
      t=A(i)
      A(i)=A(j)
      A(j)=t
    end if
  end do
  A(left)=A(i)
  A(i)=temp
  call sort(A,left,i-1_2)
  call sort(A,i+1_2,right)
end subroutine sort

end module functions

max_stars_repo_name: Caixinzhang-hnu/disorder

max_issues_repo_name: Caixinzhang-hnu/disorder

max_forks_repo_name: Caixinzhang-hnu/disorder


text:
module data_struct

public

   ! floating point precision
   !integer, parameter :: dp = selected_real_kind(6, 37) ! 32-bit
   !integer, parameter :: dp = selected_real_kind(15, 307) ! 64-bit
   integer, parameter :: dp = selected_real_kind(33, 4931) ! 128-bit

   ! intel 3 vector csr format
   type dcsr3
      integer :: m ! number of rows
      double precision, allocatable :: val(:)
      integer, allocatable :: col(:)
      integer, allocatable :: rowind(:)
   end type

   type fis
      double precision :: val
      double precision :: dir(8)
   end type

   type Mrho
      double precision, allocatable :: eigval(:)
      double precision, allocatable :: eigvec(:,:)
   end type

   type rho
      type( Mrho ), allocatable :: array(:)
      real( kind = dp), allocatable :: weights(:)
   end type

contains

   subroutine EvenBounds( N, M, kmin, kmax )

      implicit none
      integer, intent(in) :: N, M
      integer, intent(out) :: kmin, kmax

      ! local variables
      real, dimension(3) :: btab

      btab = (/ 0.0, real(M)/2.0, -real(M)/2.0 /)
      kmin = ceiling( maxval( btab ) )

      btab = (/ real(N)/2.0, real(N)-real(M)/2.0, real(N)+real(M)/2.0 /)
      kmax = floor( minval( btab ) )

   end subroutine

   subroutine OddBounds( N, M, kmin, kmax )

      implicit none
      integer, intent(in) :: N, M
      integer, intent(out) :: kmin, kmax

      ! local variables
      real btab(3)

      btab = (/ 0.0, -real(M+1)/2.0, real(M-1)/2.0 /)
      kmin = ceiling( maxval( btab ) )

      btab = (/ real(N-1)/2.0, real(N)-real(M+1)/2.0, real(N)+real(M-1)/2.0 /)
      kmax = floor( minval( btab ) )

   end subroutine

end module

max_stars_repo_name: dkajtoch/spinorBEC

max_issues_repo_name: dkajtoch/spinorBEC

max_forks_repo_name: dkajtoch/spinorBEC


      program set_name
! ===============================================
      implicit none

      character*128   ::  buf
      real            ::  lon, lat

      character*1     ::  ew, sn
      character*2     ::  clat
      character*3     ::  clon
      character*7     ::  cname
! ===============================================
      call getarg(1,buf)
       read(buf,*) lon
      call getarg(2,buf)
       read(buf,*) lat

      if( lon<0 )then
        ew='w'
        write(clon,'(i3.3)') int(-lon)
      else
        ew='e'
        write(clon,'(i3.3)')  int(lon)
      endif

      if( lat<0 )then
        sn='s'
        write(clat,'(i2.2)') int(-lat)
      else
        sn='n'
        write(clat,'(i2.2)')  int(lat)
      endif

      cname=sn//clat//ew//clon

      print '(a7)', cname

      end program set_name


program m
    use percolation
    use utilities
    implicit none

    integer :: num_Ls, i, num_samples, fileunit
    integer, dimension(:), allocatable :: Ls
    real(kind=dp) :: x1, x2, tolerance, nu, slope, const
    real(kind=dp), dimension(:,:), allocatable :: invPI
    real(kind=dp), dimension(:), allocatable :: diff, diffapprox
    character(len=:), allocatable :: filename

    Ls = [(2**i, i = 4, 9)]
    num_Ls = size(Ls)

    x1 = 0.8; x2 = 0.3
    num_samples = 2**16; tolerance = 0.0001d0

    allocate(invPI(num_Ls, 2))

    invPI(:,1) = [(spanning_probability_inverse(x1, Ls(i), num_samples/Ls(i), tolerance), &
                   i = 1, num_Ls)]
    invPI(:,2) = [(spanning_probability_inverse(x2, Ls(i), num_samples/Ls(i), tolerance), &
                   i = 1, num_Ls)]

    diff = invPI(:,1) - invPI(:,2)

    call linfit(log(1.0d0*Ls), log(diff), slope, const)
    nu = -1/slope
    diffapprox = exp(const) * Ls**slope

    open(newunit=fileunit, file="tmp/m.dat", status="replace")
    do i = 1, num_Ls
        write(unit=fileunit, fmt=*) Ls(i), diff(i), diffapprox(i)
    end do
    close(unit=fileunit)

    open(newunit=fileunit, file="tmp/m_nu.dat", status="replace")
    write(unit=fileunit, fmt="(f0.3)") nu
    close(unit=fileunit)

end program m


text: #redirect Christian Democratic Party


text:
!
! Copyright (c) 2015, 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.
!

program test_f90_intrinsics
integer, parameter :: n = 4
integer, parameter :: m = 3
integer, parameter :: l = 5
!
integer, parameter :: ntests = 44
logical, dimension(ntests) :: checks
integer ic
!
integer, dimension(n,m,l) :: ia
logical, dimension(n,m,l) :: la
real,    dimension(n,m,l) :: ra
!
! Answers
logical lan1(m,l,2), lan2(n,l,2), lan3(n,m,2)
integer ico1(m,l,5), ico2(n,l,5), ico3(n,m,5)
real    rma1(m,l,4), rma2(n,l,4), rma3(n,m,4)
!
! Checks
logical clan1(m,l), clan2(n,l), clan3(n,m)
logical clal1(m,l), clal2(n,l), clal3(n,m)
integer cico1(m,l), cico2(n,l), cico3(n,m), cic
integer cima1(m,l), cima2(n,l), cima3(n,m), cma
integer cimi1(m,l), cimi2(n,l), cimi3(n,m), cmi
integer cipr1(m,l), cipr2(n,l), cipr3(n,m), cpr
integer ciad1(m,l), ciad2(n,l), ciad3(n,m), cad

integer iexp(ntests)

!
! Input data
data ia /  1,  4,  7, 10,  2,  5,  8, 11,  3,  6,  9, 12, &
           0,  3,  6,  9,  1,  4,  7, 10,  2,  5,  8, 11, &
          -1,  2,  5,  8,  0,  3,  6,  9,  1,  4,  7, 10, &
          -2,  1,  4,  7, -1,  2,  5,  8,  0,  3,  6,  9, &
          -3,  0,  3,  6, -2,  1,  4,  7, -1,  2,  5,  8 /
!
! Result check data
data clan1 / T, T, T, T, T, T, T, T, T, T, T, T, F, T, T /
data clan2 / F, F, T, T, F, F, T, T, F, F, T, T, F, F, F, T, F, F, F, T /
data clan3 / F, F, T, T, F, F, T, T, F, F, T, T /
!
data clal1 / F, F, F, F, F, F, F, F, F, F, F, F, F, F, F /
data clal2 / F, F, T, T, F, F, F, T, F, F, F, T, F, F, F, T, F, F, F, F /
data clal3 / F, F, F, F, F, F, F, T, F, F, F, T /
!
data cic / 20 /
data cico1 / 2, 2, 2, 1, 2, 2, 1, 1, 2, 1, 1, 1, 0, 1, 1 /
data cico2 / 0, 0, 3, 3, 0, 0, 2, 3, 0, 0, 1, 3, 0, 0, 0, 3, 0, 0, 0, 2 /
data cico3 / 0, 0, 1, 4, 0, 0, 2, 5, 0, 0, 3, 5 /
!
data cma / 12 /
data cima1 /10,11,12, 9,10,11, 8, 9,10, 7, 8, 9, 6, 7, 8 /
data cima2 / 3, 6, 9,12, 2, 5, 8,11, 1, 4, 7,10, 0, 3, 6, 9,-1, 2, 5, 8 /
data cima3 / 1, 4, 7,10, 2, 5, 8,11, 3, 6, 9,12 /
!
data cmi / -3 /
data cimi1 / 1, 2, 3, 0, 1, 2,-1, 0, 1,-2,-1, 0,-3,-2,-1 /
data cimi2 / 1, 4, 7,10, 0, 3, 6, 9,-1, 2, 5, 8,-2, 1, 4, 7,-3, 0, 3, 6 /
data cimi3 /-3, 0, 3, 6,-2, 1, 4, 7,-1, 2, 5, 8 /
!
data cpr / 0 /
data cipr1 / 280, 880, 1944, 0, 280, 880, -80, 0, 280, &
             -56, -80, 0, 0, -56, -80 /
data cipr2 / 6, 120, 504, 1320, 0, 60, 336, 990, 0, 24, 210, 720, &
             0, 6, 120, 504, -6, 0, 60, 336 /
data cipr3 / 0, 0, 2520, 30240, 0, 120, 6720, 55440, 0, 720, 15120, 95040 /
!
data cad / 270 /
data ciad1 / 22, 26, 30, 18, 22, 26, 14, 18, 22, 10, 14, 18, 6, 10, 14 /
data ciad2 / 6, 15, 24, 33, 3, 12, 21, 30, 0, 9, 18, 27, &
             -3, 6, 15, 24, -6, 3, 12, 21 /
data ciad3 / -5, 10, 25, 40, 0, 15, 30, 45, 5, 20, 35, 50 /
!
!
!
! Create real and logical arrays
ra = real(ia)
la = ra .gt. 6.0
ic = 1
!
! Do the tests
lan1(:,:,1) = any(la,dim=1)
lan2(:,:,1) = any(la,dim=2)
lan3(:,:,1) = any(la,dim=3)
checks(ic) = .not. any(la); ic = ic + 1
checks(ic) = any(clan1 .neqv. lan1(:,:,1)); ic = ic + 1
checks(ic) = any(clan2 .neqv. lan2(:,:,1)); ic = ic + 1
checks(ic) = any(clan3 .neqv. lan3(:,:,1)); ic = ic + 1
!
lan1(:,:,2) = all(la,dim=1)
lan2(:,:,2) = all(la,dim=2)
lan3(:,:,2) = all(la,dim=3)
checks(ic) = all(la); ic = ic + 1
checks(ic) = any(clal1 .neqv. lan1(:,:,2)); ic = ic + 1
checks(ic) = any(clal2 .neqv. lan2(:,:,2)); ic = ic + 1
checks(ic) = any(clal3 .neqv. lan3(:,:,2)); ic = ic + 1
!
ico1(:,:,1) = count(la,dim=1)
ico2(:,:,1) = count(la,dim=2)
ico3(:,:,1) = count(la,dim=3)
checks(ic) = ((cic - count(la)) .ne. 0); ic = ic + 1
checks(ic) = any((cico1 - ico1(:,:,1)) .ne. 0); ic = ic + 1
checks(ic) = any((cico2 - ico2(:,:,1)) .ne. 0); ic = ic + 1
checks(ic) = any((cico3 - ico3(:,:,1)) .ne. 0); ic = ic + 1
!
ico1(:,:,2) = maxval(ia,dim=1)
ico2(:,:,2) = maxval(ia,dim=2)
ico3(:,:,2) = maxval(ia,dim=3)
checks(ic) = ((cma - maxval(ia)).ne. 0); ic = ic + 1
checks(ic) = any((cima1 - ico1(:,:,2)) .ne. 0); ic = ic + 1
checks(ic) = any((cima2 - ico2(:,:,2)) .ne. 0); ic = ic + 1
checks(ic) = any((cima3 - ico3(:,:,2)) .ne. 0); ic = ic + 1
!
ico1(:,:,3) = minval(ia,dim=1)
ico2(:,:,3) = minval(ia,dim=2)
ico3(:,:,3) = minval(ia,dim=3)
checks(ic) = ((cmi - minval(ia)).ne. 0); ic = ic + 1
checks(ic) = any((cimi1 - ico1(:,:,3)) .ne. 0); ic = ic + 1
checks(ic) = any((cimi2 - ico2(:,:,3)) .ne. 0); ic = ic + 1
checks(ic) = any((cimi3 - ico3(:,:,3)) .ne. 0); ic = ic + 1
!
rma1(:,:,1) = maxval(ra,dim=1)
rma2(:,:,1) = maxval(ra,dim=2)
rma3(:,:,1) = maxval(ra,dim=3)
checks(ic) = ((real(cma) - maxval(ra)).ne. 0); ic = ic + 1
checks(ic) = any((real(cima1) - rma1(:,:,1)) .ne. 0); ic = ic + 1
checks(ic) = any((real(cima2) - rma2(:,:,1)) .ne. 0); ic = ic + 1
checks(ic) = any((real(cima3) - rma3(:,:,1)) .ne. 0); ic = ic + 1
!
rma1(:,:,2) = minval(ra,dim=1)
rma2(:,:,2) = minval(ra,dim=2)
rma3(:,:,2) = minval(ra,dim=3)
checks(ic) = ((real(cmi) - minval(ra)).ne. 0); ic = ic + 1
checks(ic) = any((real(cimi1) - rma1(:,:,2)) .ne. 0); ic = ic + 1
checks(ic) = any((real(cimi2) - rma2(:,:,2)) .ne. 0); ic = ic + 1
checks(ic) = any((real(cimi3) - rma3(:,:,2)) .ne. 0); ic = ic + 1
!
ico1(:,:,4) = product(ia,dim=1)
ico2(:,:,4) = product(ia,dim=2)
ico3(:,:,4) = product(ia,dim=3)
checks(ic) = ((cpr - product(ia)).ne. 0); ic = ic + 1
checks(ic) = any((cipr1 - ico1(:,:,4)) .ne. 0); ic = ic + 1
checks(ic) = any((cipr2 - ico2(:,:,4)) .ne. 0); ic = ic + 1
checks(ic) = any((cipr3 - ico3(:,:,4)) .ne. 0); ic = ic + 1
!
rma1(:,:,3) = product(ra,dim=1)
rma2(:,:,3) = product(ra,dim=2)
rma3(:,:,3) = product(ra,dim=3)
checks(ic) = ((real(cpr) - product(ra)).ne. 0); ic = ic + 1
checks(ic) = any((real(cipr1) - rma1(:,:,3)) .ne. 0); ic = ic + 1
checks(ic) = any((real(cipr2) - rma2(:,:,3)) .ne. 0); ic = ic + 1
checks(ic) = any((real(cipr3) - rma3(:,:,3)) .ne. 0); ic = ic + 1
!
ico1(:,:,5) = sum(ia,dim=1)
ico2(:,:,5) = sum(ia,dim=2)
ico3(:,:,5) = sum(ia,dim=3)
checks(ic) = ((cad - sum(ia)).ne. 0); ic = ic + 1
checks(ic) = any((ciad1 - ico1(:,:,5)) .ne. 0); ic = ic + 1
checks(ic) = any((ciad2 - ico2(:,:,5)) .ne. 0); ic = ic + 1
checks(ic) = any((ciad3 - ico3(:,:,5)) .ne. 0); ic = ic + 1
!
rma1(:,:,4) = sum(ra,dim=1)
rma2(:,:,4) = sum(ra,dim=2)
rma3(:,:,4) = sum(ra,dim=3)
checks(ic) = ((real(cad) - sum(ra)) .ne. 0.0); ic = ic + 1
checks(ic) = any((real(ciad1) - rma1(:,:,4)) .ne. 0.0); ic = ic + 1
checks(ic) = any((real(ciad2) - rma2(:,:,4)) .ne. 0.0); ic = ic + 1
checks(ic) = any((real(ciad3) - rma3(:,:,4)) .ne. 0.0); ic = ic + 1
!
if ((ic-1) .ne. ntests) print *,"Error in number of tests!"
!
do i = 1, ntests
   iexp(i) = 0
end do

call check(checks,iexp,ntests)

!
end

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
known to most as oki
third year student at davis (english creative writing & film studies double major)
is a DJ at KDVS
currently lives east of downtown

http://okijason.blogspot.com Blogs with users/JasonDunne Jason

other profiles:
http://ucdavis.facebook.com/profile.php?id1064520006 Facebook
http://www.myspace.com/okidoki Myspace