text:
! Computed goto
      INTEGER NEXT

      GO TO (100,200) NEXT
10    PRINT *,'Control transfers here if NEXT does not equal 1 or 2'

100   PRINT *,'Control transfers here if NEXT = 1'

200   PRINT *,'Control transfers here if NEXT = 2'

    ! This second use of the GOTO fails in the new version of OFP/ROSE
      GO TO (100,200) NEXT

      END

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      SUBROUTINE HPGL_DONE

c***********************************************************************

c  Function:
c     Terminate plot to HPGL device.

c  Parameters passed:

c  Parameters returned:

c  Common blocks referenced:
c     graphics

c  Functions and subroutines referenced:
c     close
c     write

c  Entry:

c  References:

c  Change History:

c*******************!***************************************************

      include      'graphics.cmn'

c     Put pen away.
      write(graphics_LogicalUnit,'('' sp;'')')

c     Close output file
      close(graphics_LogicalUnit)

      RETURN
      END      ! HPGL_DONE

ext: for

max_stars_repo_name: spinorkit/LWPC

max_issues_repo_name: spinorkit/LWPC

max_forks_repo_name: spinorkit/LWPC


text:
\ -----------------------------------------------------------------------------
\  OpenGl for Windows
\ -----------------------------------------------------------------------------

S" lib\include\float2.f" INCLUDED

MODULE: HIDDEN  \   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

WINAPI: CreateWindowExA		USER32.DLL
WINAPI: GetSystemMetrics	USER32.DLL
WINAPI: GetDC			USER32.DLL
WINAPI: ReleaseDC		USER32.DLL
WINAPI: ShowCursor		USER32.DLL
WINAPI: GetAsyncKeyState	USER32.DLL

WINAPI: ChoosePixelFormat	GDI32.DLL
WINAPI: SetPixelFormat		GDI32.DLL
WINAPI: SwapBuffers		GDI32.DLL

WINAPI: wglCreateContext	OPENGL32.DLL
WINAPI: wglMakeCurrent		OPENGL32.DLL
WINAPI: glHint			OPENGL32.DLL
WINAPI: glMatrixMode		OPENGL32.DLL
WINAPI: glClear			OPENGL32.DLL
WINAPI: glLoadIdentity		OPENGL32.DLL
WINAPI: glTranslated		OPENGL32.DLL
WINAPI: glPointSize		OPENGL32.DLL
WINAPI: glBegin			OPENGL32.DLL
WINAPI:	glVertex2d		OPENGL32.DLL
WINAPI:	glEnd			OPENGL32.DLL
WINAPI: glColor3b		OPENGL32.DLL
WINAPI: glRotated		OPENGL32.DLL
WINAPI: glLineWidth		OPENGL32.DLL
WINAPI: glPolygonMode		OPENGL32.DLL
WINAPI: glEnable		OPENGL32.DLL
WINAPI: glDisable		OPENGL32.DLL

WINAPI: gluPerspective		GLU32.DLL

0
2 -- nSize
2 -- nVersion
CELL -- dwFlags
1 -- iPixelType
1 -- cColorBits
1 -- cRedBits
1 -- cRedShift
1 -- cGreenBits
1 -- cGreenShift
1 -- cBlueBits
1 -- cBlueShift
1 -- cAlphaBits
1 -- cAlphaShift
1 -- cAccumBits
1 -- cAccumRedBits
1 -- cAccumGreenBits
1 -- cAccumBlueBits
1 -- cAccumAlphaBits
1 -- cDepthBits
1 -- cStencilBits
1 -- cAuxBuffers
1 -- iLayerType
1 -- bReserved
CELL -- dwLayerMask
CELL -- dwVisibleMask
CELL -- dwDamageMask
CONSTANT PIXELFORMATDESCRIPTOR
0 VALUE pfd		\ structure for opengl
PIXELFORMATDESCRIPTOR ALLOCATE THROW TO pfd
0x25 pfd dwFlags !
32 pfd cColorBits C!
pfd FREE THROW

0 VALUE glhandle	\ handle window
0 VALUE glhdc		\ context handle

EXPORT

0.0E FVALUE X		\ x coordinate
0.0E FVALUE Y		\ Y coordinate
0.0E FVALUE X1		\ x1
0.0E FVALUE Y1		\ y1
0.0E FVALUE X2		\ x2
0.0E FVALUE Y2		\ y2
0.0E FVALUE X3		\ x3
0.0E FVALUE Y3		\ y3

\ The height of the screen image in pixels
: VScreen ( -> n )
	1 GetSystemMetrics ;

\ Screen width in pixels
: HScreen ( -> n )
	0 GetSystemMetrics ;
\ Cls
\
\ Number of float from float stack to data stack
: F>FL  ( -> f ) ( F: f -> )
	[
	0x8D C, 0x6D C, 0xFC C,
	0xD9 C, 0x5D C, 0x00 C,
	0x87 C, 0x45 C, 0x00 C,
	0xC3 C, ] ;

\ Number of float from float stack to data stack
: F>DL ( -> d ) ( F: f -> )
	FLOAT>DATA SWAP ;

\ Convert number on stack to float
: S>FL ( n -> f )
	DS>F F>FL ;

\ Convert number on stack to double
: S>DL ( n -> d )
	DS>F F>DL ;

\ Aspect ratio of width and height
: AScreen ( -> d )
	HScreen DS>F VScreen S>D D>F F/ F>DL ;

\ Output the specified number of lines on the screen (from point to 4)
: LineLoop ( n -> )
	DUP 0 > IF Y F>DL X F>DL glVertex2d DROP THEN
	DUP 1 > IF Y1 F>DL X2 F>DL glVertex2d DROP THEN
	DUP 2 > IF Y2 F>DL X2 F>DL glVertex2d DROP THEN
	DUP 3 > IF Y3 F>DL X3 F>DL glVertex2d DROP THEN
	DROP ;

\ Show cursor on screen
: ShowCursore ( -> )
	1 ShowCursor DROP ;

\ Hide screen cursor
: HideCursore ( -> )
	0 ShowCursor DROP ;

\ Output
: glClose ( -> )
	glhdc glhandle ReleaseDC 0 ExitProcess ;

: glClose2 ( -> )
	glhdc glhandle ReleaseDC  ;  \  0 ExitProcess ;

\ Checks if a key is pressed by a given code
: key ( n -- flag )
	GetAsyncKeyState ;

\ &&&&&&&&&&&&&&&&&&&&&       SCREEN  SIZE   &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
\ &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
\ Initialize gl window
: glOpen ( -> )
  600 800 0 0  VScreen  2 /   HScreen   2 /  140 140 \ border
  0x90000000 S" edit" DROP DUP 8 CreateWindowExA
   \  WILL OPEN AT 1/4 SCREEN TOP LEFT
  \ 0 0 800 660 VScreen HScreen 0 0 0x90000000 S" edit" DROP DUP 8 CreateWindowExA

  \  TO MAKE  FULL SCREEN  USE THIS LINE :
\ 0  0  0 0  VScreen  ( 2 /  )  HScreen  ( 2 / )  0 0 0x90000000 S" edit" DROP DUP 8 CreateWindowExA
\ &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
\ &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

 DUP TO glhandle GetDC TO glhdc
 pfd DUP glhdc ChoosePixelFormat glhdc SetPixelFormat DROP
 glhdc wglCreateContext glhdc wglMakeCurrent DROP
 0x1102 0xC50 glHint DROP
 0x1701 glMatrixMode DROP
 100 S>DL 1 DS>F 10 DS>F F/ F>DL AScreen 90 S>DL gluPerspective DROP
 0x1700 glMatrixMode DROP ;

\ Clear image buffer
: Cls ( -> )
	0x4000 glClear DROP ;

 \ Show image buffer (where we drew)
: View ( -> )
	glhdc SwapBuffers DROP ;

\ The unit matrix
: SingleMatrix ( -> )
	glLoadIdentity DROP ;

\ Shifts the current matrix by a vector (x, y, z).
: ShiftMatrix ( F: f f f -> )
	F>DL F>DL F>DL glTranslated DROP ;

\ Point size
: PointSize ( n -> )
	S>FL glPointSize DROP ;

\  DUP 0 > IF Y F>DL X F>DL glVertex2d DROP THEN

\  Display 2D point on the screen (X, Y)
: Point ( -> )
	0 glBegin DROP 1 LineLoop glEnd DROP ;

\ SET COLOR  (B G R)
: Color ( n n n -> )
	glColor3b DROP ;

\  Rotates the current matrix at a given angle around the specified vector.
\ z y x angle
: RotatedMatrix ( F: f f f f -> )
	F>DL F>DL F>DL F>DL glRotated DROP ;

\ \ Output 2D line (X, Y, X1, Y1)
: Line ( -> )
	1 glBegin DROP
	Y F>DL X F>DL glVertex2d DROP
	Y1 F>DL X1 F>DL glVertex2d DROP
	glEnd DROP ;

\ Draw line width
: LineSize ( n -> )
	S>FL glLineWidth DROP ;

\ Print 2D triangle (X,Y,X1,Y1,X2,Y2)
: Triangle ( -> )
	4 glBegin DROP
	Y F>DL X F>DL glVertex2d DROP
	Y1 F>DL X1 F>DL glVertex2d DROP
	Y2 F>DL X2 F>DL glVertex2d DROP
	glEnd DROP ;

\ Drawing shapes with a wire style
: GlLine ( -> )
	0x1B01 0x408 glPolygonMode DROP ;

\ Drawing shapes painted style
: GlFill ( -> )
	0x1B02 0x408 glPolygonMode DROP ;

\ 2D rectangle (X,Y,X1,Y1)
: Rectangle ( -> )
	7 glBegin DROP
	Y F>DL X F>DL glVertex2d DROP
	Y F>DL X1 F>DL glVertex2d DROP
	Y1 F>DL X1 F>DL glVertex2d DROP
	Y1 F>DL X F>DL glVertex2d DROP
	glEnd DROP ;

\ Smoothing
: Smoothing ( -> )
	0xB10 glEnable DROP ;

\ Remove smoothing
: NoSmoothing ( -> )
	0xB10 glDisable DROP ;

\ 2D quadrilateral (X,Y,X1,Y1,X2,Y2,X3,Y3)
: Tetragon ( -> )
	7 glBegin DROP
	Y F>DL X F>DL glVertex2d DROP
	Y1 F>DL X1 F>DL glVertex2d DROP
	Y2 F>DL X2 F>DL glVertex2d DROP
	Y3 F>DL X3 F>DL glVertex2d DROP
	glEnd DROP ;

\ 2D angle (X,Y,X1,Y1,X2,Y2)
: Corner ( -> )
	3 glBegin DROP
	Y F>DL X F>DL glVertex2d DROP
	Y1 F>DL X1 F>DL glVertex2d DROP
	Y2 F>DL X2 F>DL glVertex2d DROP
	glEnd DROP ;

\ basic retro colors

: green ( -- ) \ light-green
	0 100 0 Color ;

: dark-green ( -- )
	0 300 0 Color ;

: red ( -- )
	0 0 100 Color ;

: blue ( -- )
	100 0 0 Color ;

: cyan ( -- )
	100 0 0 Color ;

: yellow ( -- )
	0 100 100 Color ;

: white ( -- )
	100 100 100 Color ;

: black ( -- )
	0 0 0 Color ;

\ STARTLOG

;MODULE  \ >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

0 VALUE msec					\ for synchronization
0 VALUE msei
0.0E FVALUE theta				\ matrix rotation angle

: S>F    S>D D>F ;         \ *******  this is very important  !!!
\ ----------------------------------------------------------------

: Dessin   \ start the routine

\ --------------------------------------: main ( -- )
       glOpen
          GlLine
              Smoothing
                   HideCursore
                         SingleMatrix
           0.0E 0.0E -5.0E ShiftMatrix	\ grandeur 5
              theta 0.0E 0.0E 1.0E RotatedMatrix
                  TIMER@ DROP TO msei  msei msec <> IF   msei TO msec
      Cls
            VIEW
\ -----------------------ATTENTION HERE  -------------------------
;

ext: f


text:
      SUBROUTINE STREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
     $                   INFO )
*
*  -- LAPACK routine (version 3.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2006
*
*     .. Scalar Arguments ..
      CHARACTER          COMPQ
      INTEGER            IFST, ILST, INFO, LDQ, LDT, N
*     ..
*     .. Array Arguments ..
      REAL               Q( LDQ, * ), T( LDT, * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  STREXC reorders the real Schur factorization of a real matrix
*  A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
*  moved to row ILST.
*
*  The real Schur form T is reordered by an orthogonal similarity
*  transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
*  is updated by postmultiplying it with Z.
*
*  T must be in Schur canonical form (as returned by SHSEQR), that is,
*  block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
*  2-by-2 diagonal block has its diagonal elements equal and its
*  off-diagonal elements of opposite sign.
*
*  Arguments
*  =========
*
*  COMPQ   (input) CHARACTER*1
*          = 'V':  update the matrix Q of Schur vectors;
*          = 'N':  do not update Q.
*
*  N       (input) INTEGER
*          The order of the matrix T. N >= 0.
*
*  T       (input/output) REAL array, dimension (LDT,N)
*          On entry, the upper quasi-triangular matrix T, in Schur
*          Schur canonical form.
*          On exit, the reordered upper quasi-triangular matrix, again
*          in Schur canonical form.
*
*  LDT     (input) INTEGER
*          The leading dimension of the array T. LDT >= max(1,N).
*
*  Q       (input/output) REAL array, dimension (LDQ,N)
*          On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
*          On exit, if COMPQ = 'V', Q has been postmultiplied by the
*          orthogonal transformation matrix Z which reorders T.
*          If COMPQ = 'N', Q is not referenced.
*
*  LDQ     (input) INTEGER
*          The leading dimension of the array Q.  LDQ >= max(1,N).
*
*  IFST    (input/output) INTEGER
*  ILST    (input/output) INTEGER
*          Specify the reordering of the diagonal blocks of T.
*          The block with row index IFST is moved to row ILST, by a
*          sequence of transpositions between adjacent blocks.
*          On exit, if IFST pointed on entry to the second row of a
*          2-by-2 block, it is changed to point to the first row; ILST
*          always points to the first row of the block in its final
*          position (which may differ from its input value by +1 or -1).
*          1 <= IFST <= N; 1 <= ILST <= N.
*
*  WORK    (workspace) REAL array, dimension (N)
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*          = 1:  two adjacent blocks were too close to swap (the problem
*                is very ill-conditioned); T may have been partially
*                reordered, and ILST points to the first row of the
*                current position of the block being moved.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO
      PARAMETER          ( ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            WANTQ
      INTEGER            HERE, NBF, NBL, NBNEXT
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     .. External Subroutines ..
      EXTERNAL           SLAEXC, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX
*     ..
*     .. Executable Statements ..
*
*     Decode and test the input arguments.
*
      INFO = 0
      WANTQ = LSAME( COMPQ, 'V' )
      IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
         INFO = -4
      ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN
         INFO = -6
      ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
         INFO = -7
      ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
         INFO = -8
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'STREXC', -INFO )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( N.LE.1 )
     $   RETURN
*
*     Determine the first row of specified block
*     and find out it is 1 by 1 or 2 by 2.
*
      IF( IFST.GT.1 ) THEN
         IF( T( IFST, IFST-1 ).NE.ZERO )
     $      IFST = IFST - 1
      END IF
      NBF = 1
      IF( IFST.LT.N ) THEN
         IF( T( IFST+1, IFST ).NE.ZERO )
     $      NBF = 2
      END IF
*
*     Determine the first row of the final block
*     and find out it is 1 by 1 or 2 by 2.
*
      IF( ILST.GT.1 ) THEN
         IF( T( ILST, ILST-1 ).NE.ZERO )
     $      ILST = ILST - 1
      END IF
      NBL = 1
      IF( ILST.LT.N ) THEN
         IF( T( ILST+1, ILST ).NE.ZERO )
     $      NBL = 2
      END IF
*
      IF( IFST.EQ.ILST )
     $   RETURN
*
      IF( IFST.LT.ILST ) THEN
*
*        Update ILST
*
         IF( NBF.EQ.2 .AND. NBL.EQ.1 )
     $      ILST = ILST - 1
         IF( NBF.EQ.1 .AND. NBL.EQ.2 )
     $      ILST = ILST + 1
*
         HERE = IFST
*
   10    CONTINUE
*
*        Swap block with next one below
*
         IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
*
*           Current block either 1 by 1 or 2 by 2
*
            NBNEXT = 1
            IF( HERE+NBF+1.LE.N ) THEN
               IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
     $            NBNEXT = 2
            END IF
            CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT,
     $                   WORK, INFO )
            IF( INFO.NE.0 ) THEN
               ILST = HERE
               RETURN
            END IF
            HERE = HERE + NBNEXT
*
*           Test if 2 by 2 block breaks into two 1 by 1 blocks
*
            IF( NBF.EQ.2 ) THEN
               IF( T( HERE+1, HERE ).EQ.ZERO )
     $            NBF = 3
            END IF
*
         ELSE
*
*           Current block consists of two 1 by 1 blocks each of which
*           must be swapped individually
*
            NBNEXT = 1
            IF( HERE+3.LE.N ) THEN
               IF( T( HERE+3, HERE+2 ).NE.ZERO )
     $            NBNEXT = 2
            END IF
            CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT,
     $                   WORK, INFO )
            IF( INFO.NE.0 ) THEN
               ILST = HERE
               RETURN
            END IF
            IF( NBNEXT.EQ.1 ) THEN
*
*              Swap two 1 by 1 blocks, no problems possible
*
               CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT,
     $                      WORK, INFO )
               HERE = HERE + 1
            ELSE
*
*              Recompute NBNEXT in case 2 by 2 split
*
               IF( T( HERE+2, HERE+1 ).EQ.ZERO )
     $            NBNEXT = 1
               IF( NBNEXT.EQ.2 ) THEN
*
*                 2 by 2 Block did not split
*
                  CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1,
     $                         NBNEXT, WORK, INFO )
                  IF( INFO.NE.0 ) THEN
                     ILST = HERE
                     RETURN
                  END IF
                  HERE = HERE + 2
               ELSE
*
*                 2 by 2 Block did split
*
                  CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
     $                         WORK, INFO )
                  CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1,
     $                         WORK, INFO )
                  HERE = HERE + 2
               END IF
            END IF
         END IF
         IF( HERE.LT.ILST )
     $      GO TO 10
*
      ELSE
*
         HERE = IFST
   20    CONTINUE
*
*        Swap block with next one above
*
         IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
*
*           Current block either 1 by 1 or 2 by 2
*
            NBNEXT = 1
            IF( HERE.GE.3 ) THEN
               IF( T( HERE-1, HERE-2 ).NE.ZERO )
     $            NBNEXT = 2
            END IF
            CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
     $                   NBF, WORK, INFO )
            IF( INFO.NE.0 ) THEN
               ILST = HERE
               RETURN
            END IF
            HERE = HERE - NBNEXT
*
*           Test if 2 by 2 block breaks into two 1 by 1 blocks
*
            IF( NBF.EQ.2 ) THEN
               IF( T( HERE+1, HERE ).EQ.ZERO )
     $            NBF = 3
            END IF
*
         ELSE
*
*           Current block consists of two 1 by 1 blocks each of which
*           must be swapped individually
*
            NBNEXT = 1
            IF( HERE.GE.3 ) THEN
               IF( T( HERE-1, HERE-2 ).NE.ZERO )
     $            NBNEXT = 2
            END IF
            CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
     $                   1, WORK, INFO )
            IF( INFO.NE.0 ) THEN
               ILST = HERE
               RETURN
            END IF
            IF( NBNEXT.EQ.1 ) THEN
*
*              Swap two 1 by 1 blocks, no problems possible
*
               CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1,
     $                      WORK, INFO )
               HERE = HERE - 1
            ELSE
*
*              Recompute NBNEXT in case 2 by 2 split
*
               IF( T( HERE, HERE-1 ).EQ.ZERO )
     $            NBNEXT = 1
               IF( NBNEXT.EQ.2 ) THEN
*
*                 2 by 2 Block did not split
*
                  CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1,
     $                         WORK, INFO )
                  IF( INFO.NE.0 ) THEN
                     ILST = HERE
                     RETURN
                  END IF
                  HERE = HERE - 2
               ELSE
*
*                 2 by 2 Block did split
*
                  CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
     $                         WORK, INFO )
                  CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1,
     $                         WORK, INFO )
                  HERE = HERE - 2
               END IF
            END IF
         END IF
         IF( HERE.GT.ILST )
     $      GO TO 20
      END IF
      ILST = HERE
*
      RETURN
*
*     End of STREXC
*
      END

max_issues_repo_name: cpatrick/calatk

max_forks_repo_name: cpatrick/calatk


text:
subroutine turwav(nmmax     ,kmax      ,kfs       , &
                & vicuv     ,dis       ,dps       ,s0        ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  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 3.
!
!  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, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: turwav.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/kernel/src/compute_roller/turwav.f90 $
!!--description-----------------------------------------------------------------
!
! Computes eddy viscosity as function of roller energy,
! phase velocity and wave height
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
   !
    use globaldata
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    real(fp)               , pointer :: rhow
    real(fp)               , pointer :: thr
    logical                , pointer :: wave
    logical                , pointer :: htur2d
!
! Global variables
!
    integer                                       , intent(in) :: kmax  !  Description and declaration in esm_alloc_int.f90
    integer                                       , intent(in) :: nmmax !  Description and declaration in dimens.igs
    integer   , dimension(gdp%d%nmlb:gdp%d%nmub)  , intent(in) :: kfs   !  Description and declaration in esm_alloc_int.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub,4), intent(in) :: dis   !  Description and declaration in esm_alloc_real.f90
    real(prec), dimension(gdp%d%nmlb:gdp%d%nmub)  , intent(in) :: dps   !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)  , intent(in) :: s0    !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub, kmax + 2)     :: vicuv !  Description and declaration in esm_alloc_real.f90
!
! Local variables
!
    integer  :: khtur   ! denotes the k-index of vicuv/dicuv containing the Horizontal eddy values
    integer  :: nm
    real(fp) :: h
!
!! executable statements -------------------------------------------------------
!
    thr        => gdp%gdbetaro%thr
    rhow       => gdp%gdphysco%rhow
    wave       => gdp%gdprocs%wave
    htur2d     => gdp%gdprocs%htur2d
    !
    ! Initialization
    !
    khtur = kmax + 2
    !
    ! HORIZONTAL EDDY VISCOSITIES
    !
    do nm = 1, nmmax
       if (kfs(nm) == 1) then
          h = max(0.01_fp, s0(nm) + real(dps(nm),fp))
          !
          ! check on HLES
          ! in that case wave generated visc. should be added
          !
          if (htur2d) then
            vicuv(nm, khtur) = vicuv(nm, khtur)+(dis(nm,2)/rhow)**(1.0/3.0)*h
          else
            vicuv(nm, khtur) = (dis(nm,2)/rhow)**(1.0/3.0)*h
          endif
          vicuv(nm, khtur) = max(vicuv(nm, khtur), thr)
       else
          vicuv(nm, khtur) = thr
       endif
    enddo
end subroutine turwav

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


MODULE m_statist

! ----------------------------------------------------------------------
! MODULE: m_statist.f03
! ----------------------------------------------------------------------
! Purpose:
!  Module for computing of representative statistical quantities of a data series
! ----------------------------------------------------------------------
! Author :	Dr. Thomas Papanikolaou, Cooperative Research Centre for Spatial Information, Australia
! Created:	13 October 2017
! ----------------------------------------------------------------------

      IMPLICIT NONE
      !SAVE

Contains

SUBROUTINE statist (dx, RMSdx, Sigmadx, MEANdx, MINdx, MAXdx)

! ----------------------------------------------------------------------
! SUBROUTINE: statist.f03
! ----------------------------------------------------------------------
! Purpose:
!  Compute the statistical quantities (RMS, Min, Max) of a series matrix of numerical differences
! ----------------------------------------------------------------------
! Input arguments:
! - dx:			Matrix with the x series differences (dynamic allocatable array)
!
! Output arguments:
! - RMSdx: 		Root Mean Square of dx
! - Sigmadx: 	Standard Deviation of dx
! - MEANdx: 	Mean value of dx
! - MINdx: 		Minimum value of dx
! - MAXdx: 		Maximum value of dx
! ----------------------------------------------------------------------
! Author :	Dr. Thomas Papanikolaou, Cooperative Research Centre for Spatial Information, Australia
! Created:	13 October 2017
! ----------------------------------------------------------------------

      USE mdl_precision
      USE mdl_num
      IMPLICIT NONE

! ----------------------------------------------------------------------
! Dummy arguments declaration
! ----------------------------------------------------------------------
! IN
      REAL (KIND = prec_d), INTENT(IN), ALLOCATABLE, DIMENSION(:) :: dx
! ----------------------------------------------------------------------
! OUT
      REAL (KIND = prec_d), INTENT(OUT) :: RMSdx, Sigmadx, MEANdx, MINdx, MAXdx
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Local variables declaration
! ----------------------------------------------------------------------
      REAL (KIND = prec_d) :: sum_dx, sum_dx2, sum_dxx
      INTEGER (KIND = prec_int8) :: sz1, sz2, size_dx, i
      INTEGER (KIND = prec_int2) :: AllocateStatus, DeAllocateStatus
! ----------------------------------------------------------------------

sz1 = size(dx, DIM = 1)
!sz2 = size(dx, DIM = 2)

size_dx = sz1
sum_dx = 0.0D0
sum_dx2 = 0.0D0
MINdx = dx(1)
MAXdx = dx(1)
Do i = 1 , size_dx
   sum_dx2 = sum_dx2 + dx(i)**2
   sum_dx = sum_dx + dx(i)
   If (dx(i) < MINdx) Then
      MINdx = dx(i)
   End IF
   If (dx(i) > MAXdx) Then
      MAXdx = dx(i)
   End IF

End Do

! RMS
RMSdx = sqrt(sum_dx2 / size_dx)

! Mean
MEANdx = sum_dx / size_dx

! Standard Deviation
sum_dxx = 0.0D0
DO i = 1 , size_dx
   sum_dxx = sum_dxx + (dx(i) - MEANdx)**2
END DO
Sigmadx = sqrt(sum_dxx / size_dx)
Sigmadx = sqrt(sum_dxx / (size_dx - 1) )

END SUBROUTINE

End Module


text:
program FFTcf
use omp_lib
use myomp
use p2grid
use head
implicit none

integer(4)::LL(3)
real(4),allocatable::pos(:,:)
real(4),allocatable::ran(:,:)
real(8),allocatable::gden(:,:,:)
real(8),allocatable::rden(:,:,:)
complex(4),allocatable::vc1(:,:,:)
complex(4),allocatable::vc2(:,:,:)

real(4)::boxinfo(2,3)

integer(4)::i,j,k,dir

integer(kind(ngal))::pid

real(8),allocatable::cf8(:,:)
real(8),allocatable::cf2d(:,:,:)
real(4)::fs
real(4)::npeff
logical(4)::debug=.false.

call omp_set_num_threads(NUM_THREADS)
LL=[L,L,L]
boxinfo(1:2,1)=[0.,box]
boxinfo(1:2,2)=[0.,box]
boxinfo(1:2,3)=[0.,box]

call memo(0)
call readdata

fs=L/box
call ompshift(pos,3,ngal,[0.,0.,0.],boxinfo)
if (flag_r) call ompshift(ran,3,nran,[0.,0.,0.],boxinfo)

call ompscale(pos,3,ngal,[fs,fs,fs])
if (flag_r) call ompscale(ran,3,nran,[fs,fs,fs])

call massassign(ngal,pos,L,gden,painter)  ! get delta_g
if (flag_r) call massassign(nran,ran,L,rden,painter)  ! get delta_r

if (flag_r) gden=gden-rden

vc1=cmplx(gden,0.0)
call fftc2c_inplace(vc1,LL)  ! get delta(k)

vc1=vc1*conjg(vc1)/float(L)**3  ! get P(\vec{k})

call ifftc2c_inplace(vc1,LL)  ! get \xi(\vec{r})

if (debug) then
write(*,*) vc1(1,1,1)
do i=2,L/2
  write(*,*) vc1(i,1,1),vc1(L+2-i,1,1)
enddo
write(*,*) vc1(L/2+1,1,1)
endif

call calcxiell(vc1,LL,cf8,rbin,outname)

if (flag_xi2d) call calcxirmu(vc1,LL,cf2d,rbin,ubin,outname2)

call memo(1)

contains

subroutine memo(command)
implicit none
integer(4)::command
if (command.eq.0) then
  write(*,*) 'memo for pos:',(float(ngal)*3)*4/1024.**3,'G'
  allocate(pos(3,ngal))
  write(*,*) 'memo for gden:',(float(L)**3)*8/1024.**3,'G'
  allocate(gden(L,L,L))
  if (flag_r) then
    write(*,*) 'memo for ran:',(float(nran)*3)*4/1024.**3,'G'
    allocate(ran(3,nran))
    write(*,*) 'memo for rden:',(float(L)**3)*8/1024.**3,'G'
    allocate(rden(L,L,L))
  endif
  write(*,*) 'memo for vc1:',(float(L)**3)*8/1024.**3,'G'
  allocate(vc1(L,L,L))
  write(*,*) 'memo for vc2:',(float(L)**3)*8/1024.**3,'G'
  allocate(vc2(L,L,L))

  allocate(cf8(9,rbin))
  if (flag_xi2d) allocate(cf2d(7,rbin,ubin))
elseif(command.eq.1) then
  deallocate(pos)
  deallocate(gden)
  if (flag_r) then
    deallocate(ran)
    deallocate(rden)
  endif
  deallocate(vc1)
  deallocate(vc2)

  deallocate(cf8)
  if (flag_xi2d) deallocate(cf2d)
endif
endsubroutine memo

subroutine readdata
implicit none
integer(4)::colmax
real(4),allocatable::line(:)
! read in galaxy
colmax=maxval(gxyz)
allocate(line(colmax))
write(*,*) 'reading: ',trim(infg)
open(31,file=infg,status='old',form='formatted')
do pid=1,nghead
  read(31,*)
enddo
do pid=1,ngal
  read(31,*) line
  pos(1:3,pid)=line(gxyz(1:3))
enddo
deallocate(line)

! read in random
if (flag_r) then
colmax=maxval(rxyz)
allocate(line(colmax))
write(*,*) 'reading: ',trim(infr)
open(31,file=infr,status='old',form='formatted')
do pid=1,nrhead
  read(31,*)
enddo
do pid=1,nran
  read(31,*) line
  ran(1:3,pid)=line(rxyz(1:3))
enddo
deallocate(line)
endif
endsubroutine readdata

subroutine calcxiell(vc,LL,cf8,rbin,filename)
implicit none
integer(4)::LL(3)
complex(4)::vc(LL(1),LL(2),LL(3))
integer(4)::rbin
character(*)::filename
real(8)::cf8(9,rbin)
real(4)::ri,rj,rk,rmode,mu
real(4)::pl0,pl2,pl4
integer(4)::bin
real(4)::dr
real(4)::rbasic
real(4)::factor

dr=(rmax-rmin)/rbin
rbasic=box/L
cf8=0.
write(*,*) 'calculating correlation function'
write(*,*) 'rmin,rmax=',rmin,rmax
write(*,*) 'rbin=',rbin
write(*,*) 'dr=',dr
!$omp parallel do default(private) shared(vc,LL,rbasic,rmin,dr,rbin) &
!$omp reduction(+:cf8) schedule(guided)
do k=1,LL(3)
  if (k.le.LL(3)/2+1) then
    rk=k-1
  else
    rk=k-1-LL(3)
  endif
  rk=rbasic*rk
  do j=1,LL(2)
    if (j.le.LL(2)/2+1) then
      rj=j-1
    else
      rj=j-1-LL(2)
    endif
    rj=rbasic*rj
    do i=1,LL(1)
      if (i.le.LL(1)/2+1) then
        ri=i-1
      else
        ri=i-1-LL(1)
      endif
      ri=rbasic*ri
      if (i.eq.1.and.j.eq.1.and.k.eq.1) cycle
      if (i.eq.LL(1)/2+1.or.j.eq.LL(2)/2+1.or.k.eq.LL(3)/2+1) then
        factor=2
      else
        factor=1
      endif
      rmode=sqrt(ri*ri+rj*rj+rk*rk)
      mu=rk/rmode
      bin=floor((rmode-rmin)/dr)+1
      if (bin.lt.1.or.bin.gt.rbin) cycle
      pl0=1.*vc(i,j,k)
      pl2=0.5*(3*mu**2-1)*vc(i,j,k)
      pl4=0.125*(35*mu**4-30*mu**2+3)*vc(i,j,k)
      cf8(1,bin)=cf8(1,bin)+rmode*factor
      cf8(2,bin)=cf8(2,bin)+rmode**2*factor
      cf8(3,bin)=cf8(3,bin)+pl0*factor
      cf8(4,bin)=cf8(4,bin)+pl0**2*factor
      cf8(5,bin)=cf8(5,bin)+pl2*factor
      cf8(6,bin)=cf8(6,bin)+pl2**2*factor
      cf8(7,bin)=cf8(7,bin)+pl4*factor
      cf8(8,bin)=cf8(8,bin)+pl4**2*factor
      cf8(9,bin)=cf8(9,bin)+1.*factor
    enddo
  enddo
enddo
!$omp end parallel do

cf8(1,:)=cf8(1,:)/cf8(9,:)
cf8(2,:)=cf8(2,:)/cf8(9,:)-cf8(1,:)**2
cf8(3,:)=cf8(3,:)/cf8(9,:)
cf8(4,:)=cf8(4,:)/cf8(9,:)-cf8(3,:)**2
cf8(5,:)=cf8(5,:)/cf8(9,:)
cf8(6,:)=cf8(6,:)/cf8(9,:)-cf8(5,:)**2
cf8(7,:)=cf8(7,:)/cf8(9,:)
cf8(8,:)=cf8(8,:)/cf8(9,:)-cf8(7,:)**2

if (usecenter) then
  do i=1,rbin
    cf8(1,i)=rmin+(i-0.5)*dr
  enddo
  cf8(2,:)=0.
endif

write(*,*) 'writing: ',trim(filename)
open(32,file=filename,form='formatted',status='replace')
do i=1,rbin
  write(32,'(9e14.6)') cf8(1:9:2,i)
  if (debug) write(*,*) real(cf8(1:9:2,i))
enddo
close(32)
endsubroutine calcxiell

subroutine calcxirmu(vc,LL,cf2d,rbin,ubin,filename)
  implicit none
  integer(4)::LL(3)
  complex(4)::vc(LL(1),LL(2),LL(3))
  integer(4)::rbin
  integer(4)::ubin
  character(*)::filename
  real(8)::cf2d(7,rbin,ubin)
  real(4)::ri,rj,rk,rmode,mu
  integer(4)::bin1,bin2
  real(4)::dr,du
  real(4)::rbasic
  real(4)::factor
  real(4)::rvalue(rbin)
  real(4)::uvalue(ubin)

  dr=(rmax-rmin)/rbin
  du=(umax-umin)/ubin
  rbasic=box/L
  cf2d=0.
  write(*,*) 'calculating 2D correlation function'
  write(*,*) 'rmin,rmax=',rmin,rmax
  write(*,*) 'rbin=',rbin
  write(*,*) 'dr=',dr
  write(*,*) 'umin,rmax=',umin,umax
  write(*,*) 'ubin=',ubin
  write(*,*) 'du=',du
  !$omp parallel do default(private) shared(vc,LL,rbasic,rmin,umin,dr,du,rbin,ubin) &
  !$omp reduction(+:cf2d) schedule(guided)
  do k=1,LL(3)
    if (k.le.LL(3)/2+1) then
      rk=k-1
    else
      rk=k-1-LL(3)
    endif
    rk=rbasic*rk
    do j=1,LL(2)
      if (j.le.LL(2)/2+1) then
        rj=j-1
      else
        rj=j-1-LL(2)
      endif
      rj=rbasic*rj
      do i=1,LL(1)
        if (i.le.LL(1)/2+1) then
          ri=i-1
        else
          ri=i-1-LL(1)
        endif
        ri=rbasic*ri
        if (i.eq.1.and.j.eq.1.and.k.eq.1) cycle
        if (i.eq.LL(1)/2+1.or.j.eq.LL(2)/2+1.or.k.eq.LL(3)/2+1) then
          factor=2
        else
          factor=1
        endif
        rmode=sqrt(ri*ri+rj*rj+rk*rk)
        mu=rk/rmode
        bin1=floor((rmode-rmin)/dr)+1
        bin2=floor((mu-umin)/du)+1
        if (bin1.lt.1.or.bin1.gt.rbin) cycle
        if (bin2.lt.1.or.bin2.gt.ubin) cycle
        cf2d(1,bin1,bin2)=cf2d(1,bin1,bin2)+rmode*factor
        cf2d(2,bin1,bin2)=cf2d(2,bin1,bin2)+rmode**2*factor
        cf2d(3,bin1,bin2)=cf2d(3,bin1,bin2)+mu*factor
        cf2d(4,bin1,bin2)=cf2d(4,bin1,bin2)+mu**2*factor
        cf2d(5,bin1,bin2)=cf2d(5,bin1,bin2)+vc(i,j,k)*factor
        cf2d(6,bin1,bin2)=cf2d(6,bin1,bin2)+vc(i,j,k)**2*factor
        cf2d(7,bin1,bin2)=cf2d(7,bin1,bin2)+1.*factor
      enddo
    enddo
  enddo
  !$omp end parallel do

  cf2d(1,:,:)=cf2d(1,:,:)/cf2d(7,:,:)
  cf2d(2,:,:)=cf2d(2,:,:)/cf2d(7,:,:)-cf2d(1,:,:)**2
  cf2d(3,:,:)=cf2d(3,:,:)/cf2d(7,:,:)
  cf2d(4,:,:)=cf2d(4,:,:)/cf2d(7,:,:)-cf2d(3,:,:)**2
  cf2d(5,:,:)=cf2d(5,:,:)/cf2d(7,:,:)
  cf2d(6,:,:)=cf2d(6,:,:)/cf2d(7,:,:)-cf2d(5,:,:)**2

  do i=1,rbin
    rvalue(i)=rmin+(i-0.5)*dr
  enddo
  do i=1,ubin
    uvalue(i)=umin+(i-0.5)*du
  enddo

  write(*,*) 'writing: ',trim(filename)
  open(32,file=filename,form='formatted',status='replace')
  do i=1,rbin
    do j=1,ubin
      write(32,'(4e14.6)') rvalue(i),uvalue(j),cf2d(5,i,j),cf2d(7,i,j)
      if (debug) write(*,*) rvalue(i),uvalue(j),real(cf2d(5,i,j)),cf2d(7,i,j)
    enddo
  enddo
  close(32)

endsubroutine calcxirmu

endprogram FFTcf

include 'mkl_dfti.f90'
include 'intel_3d_fft.f90'

max_stars_repo_name: ColdThunder/FFTcf

max_issues_repo_name: ColdThunder/FFTcf

max_forks_repo_name: ColdThunder/FFTcf


text:
      SUBROUTINE DT_NUC2CM

C***********************************************************************
C Lorentz-transformation of all wounded nucleons from Lab. to nucl.-   *
C nucl. cms. (This subroutine replaces NUCMOM.)                        *
C This version dated 15.01.95 is written by S. Roesler                 *
C***********************************************************************

      IMPLICIT NONE
      INTEGER i , id , ist , mode , nend
      DOUBLE PRECISION pe , px , py , pz , TINY3 , ZERO
      SAVE

      INCLUDE 'inc/dtflka'

      PARAMETER (ZERO=0.0D0,TINY3=1.0D-3)

C event history

      INCLUDE 'inc/dtevt1'
C extended event history
      INCLUDE 'inc/dtevt2'
C statistics
      INCLUDE 'inc/dtsta1'
C properties of photon/lepton projectiles
      INCLUDE 'inc/dtgpro'
C particle properties (BAMJET index convention)
      INCLUDE 'inc/dtpart'
C Glauber formalism: collision properties
      INCLUDE 'inc/dtglcp'
C*temporary
C statistics: Glauber-formalism
      INCLUDE 'inc/dtsta3'
C*

      ICWp = 0
      ICWt = 0
      NWTacc = 0
      NWAacc = 0
      NWBacc = 0

      NPOint(1) = NHKk + 1
      nend = NHKk
      DO i = 1 , nend
         IF ( (ISThkk(i).EQ.11) .OR. (ISThkk(i).EQ.12) ) THEN
            IF ( ISThkk(i).EQ.11 ) NWAacc = NWAacc + 1
            IF ( ISThkk(i).EQ.12 ) NWBacc = NWBacc + 1
            mode = ISThkk(i) - 9
C            IF (IDHKK(I).EQ.22) THEN
C* VDM assumption
C               PEIN = AAM(33)
C               IDB  = 33
C            ELSE
C               PEIN = PHKK(4,I)
C               IDB  = IDBAM(I)
C            ENDIF
C            CALL DT_LTRANS(PHKK(1,I),PHKK(2,I),PHKK(3,I),PEIN,
C     &           PX,PY,PZ,PE,IDB,MODE)
            IF ( PHKk(5,i).GT.ZERO ) THEN
               CALL DT_LTRANS(PHKk(1,i),PHKk(2,i),PHKk(3,i),PHKk(4,i),
     &                        px,py,pz,pe,IDBam(i),mode)
            ELSE
               px = PGAmm(1)
               py = PGAmm(2)
               pz = PGAmm(3)
               pe = PGAmm(4)
            END IF
            ist = ISThkk(i) - 2
            id = IDHkk(i)
C* VDM assumption
C            IF (ID.EQ.22) ID = 113
            CALL DT_EVTPUT(ist,id,i,0,px,py,pz,pe,0,0,0)
            IF ( ISThkk(i).EQ.11 ) ICWp = ICWp + 1
            IF ( ISThkk(i).EQ.12 ) ICWt = ICWt + 1
         END IF
      END DO

      NWTacc = MAX(NWAacc,NWBacc)
      ICDpr = ICDpr + ICWp
      ICDta = ICDta + ICWt
C*temporary
      IF ( (ICWp.EQ.0) .OR. (ICWt.EQ.0) ) THEN
         CALL DT_EVTOUT(4)
         STOP
      END IF

      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
      SUBROUTINE ZTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S,
     $                   SEP, WORK, LWORK, INFO )
*
*  -- LAPACK routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
*
*     .. Scalar Arguments ..
      CHARACTER          COMPQ, JOB
      INTEGER            INFO, LDQ, LDT, LWORK, M, N
      DOUBLE PRECISION   S, SEP
*     ..
*     .. Array Arguments ..
      LOGICAL            SELECT( * )
      COMPLEX*16         Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  ZTRSEN reorders the Schur factorization of a complex matrix
*  A = Q*T*Q**H, so that a selected cluster of eigenvalues appears in
*  the leading positions on the diagonal of the upper triangular matrix
*  T, and the leading columns of Q form an orthonormal basis of the
*  corresponding right invariant subspace.
*
*  Optionally the routine computes the reciprocal condition numbers of
*  the cluster of eigenvalues and/or the invariant subspace.
*
*  Arguments
*  =========
*
*  JOB     (input) CHARACTER*1
*          Specifies whether condition numbers are required for the
*          cluster of eigenvalues (S) or the invariant subspace (SEP):
*          = 'N': none;
*          = 'E': for eigenvalues only (S);
*          = 'V': for invariant subspace only (SEP);
*          = 'B': for both eigenvalues and invariant subspace (S and
*                 SEP).
*
*  COMPQ   (input) CHARACTER*1
*          = 'V': update the matrix Q of Schur vectors;
*          = 'N': do not update Q.
*
*  SELECT  (input) LOGICAL array, dimension (N)
*          SELECT specifies the eigenvalues in the selected cluster. To
*          select the j-th eigenvalue, SELECT(j) must be set to .TRUE..
*
*  N       (input) INTEGER
*          The order of the matrix T. N >= 0.
*
*  T       (input/output) COMPLEX*16 array, dimension (LDT,N)
*          On entry, the upper triangular matrix T.
*          On exit, T is overwritten by the reordered matrix T, with the
*          selected eigenvalues as the leading diagonal elements.
*
*  LDT     (input) INTEGER
*          The leading dimension of the array T. LDT >= max(1,N).
*
*  Q       (input/output) COMPLEX*16 array, dimension (LDQ,N)
*          On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
*          On exit, if COMPQ = 'V', Q has been postmultiplied by the
*          unitary transformation matrix which reorders T; the leading M
*          columns of Q form an orthonormal basis for the specified
*          invariant subspace.
*          If COMPQ = 'N', Q is not referenced.
*
*  LDQ     (input) INTEGER
*          The leading dimension of the array Q.
*          LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
*
*  W       (output) COMPLEX*16 array, dimension (N)
*          The reordered eigenvalues of T, in the same order as they
*          appear on the diagonal of T.
*
*  M       (output) INTEGER
*          The dimension of the specified invariant subspace.
*          0 <= M <= N.
*
*  S       (output) DOUBLE PRECISION
*          If JOB = 'E' or 'B', S is a lower bound on the reciprocal
*          condition number for the selected cluster of eigenvalues.
*          S cannot underestimate the true reciprocal condition number
*          by more than a factor of sqrt(N). If M = 0 or N, S = 1.
*          If JOB = 'N' or 'V', S is not referenced.
*
*  SEP     (output) DOUBLE PRECISION
*          If JOB = 'V' or 'B', SEP is the estimated reciprocal
*          condition number of the specified invariant subspace. If
*          M = 0 or N, SEP = norm(T).
*          If JOB = 'N' or 'E', SEP is not referenced.
*
*  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
*  LWORK   (input) INTEGER
*          The dimension of the array WORK.
*          If JOB = 'N', LWORK >= 1;
*          if JOB = 'E', LWORK = max(1,M*(N-M));
*          if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).
*
*          If LWORK = -1, then a workspace query is assumed; the routine
*          only calculates the optimal size of the WORK array, returns
*          this value as the first entry of the WORK array, and no error
*          message related to LWORK is issued by XERBLA.
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*
*  Further Details
*  ===============
*
*  ZTRSEN first collects the selected eigenvalues by computing a unitary
*  transformation Z to move them to the top left corner of T. In other
*  words, the selected eigenvalues are the eigenvalues of T11 in:
*
*                Z'*T*Z = ( T11 T12 ) n1
*                         (  0  T22 ) n2
*                            n1  n2
*
*  where N = n1+n2 and Z' means the conjugate transpose of Z. The first
*  n1 columns of Z span the specified invariant subspace of T.
*
*  If T has been obtained from the Schur factorization of a matrix
*  A = Q*T*Q', then the reordered Schur factorization of A is given by
*  A = (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span the
*  corresponding invariant subspace of A.
*
*  The reciprocal condition number of the average of the eigenvalues of
*  T11 may be returned in S. S lies between 0 (very badly conditioned)
*  and 1 (very well conditioned). It is computed as follows. First we
*  compute R so that
*
*                         P = ( I  R ) n1
*                             ( 0  0 ) n2
*                               n1 n2
*
*  is the projector on the invariant subspace associated with T11.
*  R is the solution of the Sylvester equation:
*
*                        T11*R - R*T22 = T12.
*
*  Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
*  the two-norm of M. Then S is computed as the lower bound
*
*                      (1 + F-norm(R)**2)**(-1/2)
*
*  on the reciprocal of 2-norm(P), the true reciprocal condition number.
*  S cannot underestimate 1 / 2-norm(P) by more than a factor of
*  sqrt(N).
*
*  An approximate error bound for the computed average of the
*  eigenvalues of T11 is
*
*                         EPS * norm(T) / S
*
*  where EPS is the machine precision.
*
*  The reciprocal condition number of the right invariant subspace
*  spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
*  SEP is defined as the separation of T11 and T22:
*
*                     sep( T11, T22 ) = sigma-min( C )
*
*  where sigma-min(C) is the smallest singular value of the
*  n1*n2-by-n1*n2 matrix
*
*     C  = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
*
*  I(m) is an m by m identity matrix, and kprod denotes the Kronecker
*  product. We estimate sigma-min(C) by the reciprocal of an estimate of
*  the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
*  cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
*
*  When SEP is small, small changes in T can cause large changes in
*  the invariant subspace. An approximate bound on the maximum angular
*  error in the computed right invariant subspace is
*
*                      EPS * norm(T) / SEP
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY, WANTBH, WANTQ, WANTS, WANTSP
      INTEGER            IERR, K, KASE, KS, LWMIN, N1, N2, NN
      DOUBLE PRECISION   EST, RNORM, SCALE
*     ..
*     .. Local Arrays ..
      INTEGER            ISAVE( 3 )
      DOUBLE PRECISION   RWORK( 1 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   ZLANGE
      EXTERNAL           LSAME, ZLANGE
*     ..
*     .. External Subroutines ..
      EXTERNAL           XERBLA, ZLACN2, ZLACPY, ZTREXC, ZTRSYL
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, SQRT
*     ..
*     .. Executable Statements ..
*
*     Decode and test the input parameters.
*
      WANTBH = LSAME( JOB, 'B' )
      WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
      WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
      WANTQ = LSAME( COMPQ, 'V' )
*
*     Set M to the number of selected eigenvalues.
*
      M = 0
      DO 10 K = 1, N
         IF( SELECT( K ) )
     $      M = M + 1
   10 CONTINUE
*
      N1 = M
      N2 = N - M
      NN = N1*N2
*
      INFO = 0
      LQUERY = ( LWORK.EQ.-1 )
*
      IF( WANTSP ) THEN
         LWMIN = MAX( 1, 2*NN )
      ELSE IF( LSAME( JOB, 'N' ) ) THEN
         LWMIN = 1
      ELSE IF( LSAME( JOB, 'E' ) ) THEN
         LWMIN = MAX( 1, NN )
      END IF
*
      IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.WANTS .AND. .NOT.WANTSP )
     $     THEN
         INFO = -1
      ELSE IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
         INFO = -8
      ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
         INFO = -14
      END IF
*
      IF( INFO.EQ.0 ) THEN
         WORK( 1 ) = LWMIN
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'ZTRSEN', -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( M.EQ.N .OR. M.EQ.0 ) THEN
         IF( WANTS )
     $      S = ONE
         IF( WANTSP )
     $      SEP = ZLANGE( '1', N, N, T, LDT, RWORK )
         GO TO 40
      END IF
*
*     Collect the selected eigenvalues at the top left corner of T.
*
      KS = 0
      DO 20 K = 1, N
         IF( SELECT( K ) ) THEN
            KS = KS + 1
*
*           Swap the K-th eigenvalue to position KS.
*
            IF( K.NE.KS )
     $         CALL ZTREXC( COMPQ, N, T, LDT, Q, LDQ, K, KS, IERR )
         END IF
   20 CONTINUE
*
      IF( WANTS ) THEN
*
*        Solve the Sylvester equation for R:
*
*           T11*R - R*T22 = scale*T12
*
         CALL ZLACPY( 'F', N1, N2, T( 1, N1+1 ), LDT, WORK, N1 )
         CALL ZTRSYL( 'N', 'N', -1, N1, N2, T, LDT, T( N1+1, N1+1 ),
     $                LDT, WORK, N1, SCALE, IERR )
*
*        Estimate the reciprocal of the condition number of the cluster
*        of eigenvalues.
*
         RNORM = ZLANGE( 'F', N1, N2, WORK, N1, RWORK )
         IF( RNORM.EQ.ZERO ) THEN
            S = ONE
         ELSE
            S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )*
     $          SQRT( RNORM ) )
         END IF
      END IF
*
      IF( WANTSP ) THEN
*
*        Estimate sep(T11,T22).
*
         EST = ZERO
         KASE = 0
   30    CONTINUE
         CALL ZLACN2( NN, WORK( NN+1 ), WORK, EST, KASE, ISAVE )
         IF( KASE.NE.0 ) THEN
            IF( KASE.EQ.1 ) THEN
*
*              Solve T11*R - R*T22 = scale*X.
*
               CALL ZTRSYL( 'N', 'N', -1, N1, N2, T, LDT,
     $                      T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
     $                      IERR )
            ELSE
*
*              Solve T11'*R - R*T22' = scale*X.
*
               CALL ZTRSYL( 'C', 'C', -1, N1, N2, T, LDT,
     $                      T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
     $                      IERR )
            END IF
            GO TO 30
         END IF
*
         SEP = SCALE / EST
      END IF
*
   40 CONTINUE
*
*     Copy reordered eigenvalues to W.
*
      DO 50 K = 1, N
         W( K ) = T( K, K )
   50 CONTINUE
*
      WORK( 1 ) = LWMIN
*
      RETURN
*
*     End of ZTRSEN
*
      END

max_stars_repo_name: pkeuchel/oomph-lib

max_issues_repo_name: pkeuchel/oomph-lib

max_forks_repo_name: pkeuchel/oomph-lib


!
!  Copyright 2017 SALMON developers
!
!  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.
!
SUBROUTINE simple_mixing(c1,c2)
use scf_data
implicit none

integer :: ix,iy,iz
real(8),intent(IN) :: c1,c2

!rho = c1*rho + c2*matmul( psi**2, occ )
if(ilsda == 0)then
  do iz=mg_sta(3),mg_end(3)
  do iy=mg_sta(2),mg_end(2)
  do ix=mg_sta(1),mg_end(1)
    rho(ix,iy,iz) = c1*rho_in(ix,iy,iz,num_rho_stock) + c2*rho_out(ix,iy,iz,num_rho_stock)
  end do
  end do
  end do
else if(ilsda == 1)then
  do iz=mg_sta(3),mg_end(3)
  do iy=mg_sta(2),mg_end(2)
  do ix=mg_sta(1),mg_end(1)
    rho_s(ix,iy,iz,1) = c1*rho_s_in(ix,iy,iz,1,num_rho_stock) + c2*rho_s_out(ix,iy,iz,1,num_rho_stock)
    rho_s(ix,iy,iz,2) = c1*rho_s_in(ix,iy,iz,2,num_rho_stock) + c2*rho_s_out(ix,iy,iz,2,num_rho_stock)
    rho(ix,iy,iz) = rho_s(ix,iy,iz,1)+rho_s(ix,iy,iz,2)
  end do
  end do
  end do
end if

return

END SUBROUTINE simple_mixing


text:
      real*8 function isr_tilde(x,beta)
      implicit none
      real*8 x,beta
      real*8 res
      real*8 PI
      PI=4.D0*DATAN(1.D0)
      res=0d0
      if (x .lt. 0.9999999d0) then
      res=-beta*(1d0+x)/2d0
      res=res-beta*beta*(
     c   (1d0+3d0*x*x)/(1d0-x)*dlog(x)
     c      +4d0*(1d0+x)*dlog(1d0-x)+5d0+x)
     c    /8d0
      res=res*(1d0-x)**(1-beta)
      endif
      res=res+
     c (1d0+3d0/4d0*beta+(27d0-8d0*PI*PI)/96d0*beta*beta)
     c  *beta
      isr_tilde=res
      end