text:
!##############################################################################
!# ****************************************************************************
!# <name> meshadapt </name>
!# ****************************************************************************
!#
!# <purpose>
!# This application performs local mesh adaptation based on a given
!# indicator function and prescribed refinement/recoarsening tolerances.
!#
!# Calling example:
!#
!# meshadapt -read2d mymesh -indicator myindicator.dat -refmax 3 \
!#           -reftol 0.8 -crstol 0.2
!#
!# The initial 2D mesh is read from the TRI/PRM file mymesh.tri/prm
!# and refined based on the error element-wise indicator given in file
!# myerror.dat. Cells with error larger than 0.8 are refined and cells
!# with error smalled than 0.2 are re-coarsened. The maximum
!# refinement level is 3, that is, an element of the initial grid can
!# be refined 3 times at most.
!#
!# </purpose>
!##############################################################################

program meshadapt

  use fsystem
  use genoutput
  use signals
  use storage

  use meshadaptbase

  implicit none

  ! Initialise system-wide settings
  call sys_init()

  ! Initialise the output system
  call output_init()

  ! Initialise the FEAT 2.0 storage management
  call storage_init(100, 100)

  if (madapt_signalhandler(SIGUSR1) .eq. 0) then

    ! We are in daemon mode, hence, register signal handler
    call fsignal(SIGUSR2, madapt_signalhandler) ! export mesh
    call fsignal(SIGINT,  madapt_signalhandler) ! single step
    call fsignal(SIGQUIT, madapt_signalhandler) ! finalise
    call fsignal(SIGHUP,  madapt_signalhandler) ! finalise
    call fsignal(SIGTERM, madapt_signalhandler) ! finalise

    daemon: do
      ! Perform mesh adaptation step?
      call fsignal(SIGINT, madapt_signalhandler)

      ! Export mesh to file?
      call fsignal(SIGUSR2, madapt_signalhandler)

      ! Finalise and exit?
      call fsignal(SIGQUIT, madapt_signalhandler)
      call fsignal(SIGTERM, madapt_signalhandler)
      call fsignal(SIGHUP,  madapt_signalhandler)
    end do daemon

  else

    ! We are not in daemon mode, hence, perform a single mesh
    ! adaptation step, export the mesh to file and exit
    if (madapt_signalhandler(SIGINT) .ne. 0) then
      call output_line("An error occured during mesh adaptation!",&
          OU_CLASS_ERROR,OU_MODE_STD,"meshadapt")
      call sys_halt()
    end if

    if (madapt_signalhandler(SIGUSR2) .ne. 0) then
      call output_line("An error occured while exporting mesh to file!",&
          OU_CLASS_ERROR,OU_MODE_STD,"meshadapt")
      call sys_halt()
    end if

    if (madapt_signalhandler(SIGQUIT) .ne. 0) then
      call output_line("An error occured during finalisation!",&
          OU_CLASS_ERROR,OU_MODE_STD,"meshadapt")
      call sys_halt()
    end if
  end if

  ! Clean up the storage management, finish
1 call storage_done()

contains

  subroutine meshadaptQuick(ndim,smesh,nref,nrefmax,dreftol,dcrstol)

    integer, intent(in) :: ndim,nref,nrefmax
    character(len=*), intent(in) :: smesh
    real(DP), intent(in) :: dreftol,dcrstol

    ! local variables
    type(t_meshAdapt) :: rmeshAdapt
    real(DP), dimension(:,:), allocatable :: Dcoords
    real(DP), dimension(:), allocatable :: Dind
    integer, dimension(:,:), allocatable :: IverticesAtElement
    real(DP) :: xc,yc,phi,W
    integer :: nel,nvt,nnve,iref,iel,ive,istatus

    ! Initialisation
    call madapt_init(rmeshAdapt,ndim,smesh)

    ! Get data from adaptation structure
    nel = madapt_getnel(rmeshAdapt)
    nvt = madapt_getnvt(rmeshAdapt)
    nnve = madapt_getnnve(rmeshAdapt)

    ! Get mesh from adaptation structure
    allocate(Dcoords(ndim,nvt), IverticesAtElement(nnve,nel))
    call madapt_getVertexCoords(rmeshAdapt, Dcoords)
    call madapt_getVerticesAtElement(rmeshAdapt, IverticesAtElement)

    W = 0.0_DP
    do iel=1,nel
      do ive=1,nnve
        W = W + sqrt((Dcoords(1,IverticesAtElement(ive,iel))-&
                      Dcoords(1,IverticesAtElement(mod(ive,nnve)+1,iel)))**2&
              +      (Dcoords(2,IverticesAtElement(ive,iel))-&
                      Dcoords(2,IverticesAtElement(mod(ive,nnve)+1,iel)))**2)
      end do
    end do
    W = W/real(nel*nnve)

    ! Deallocate local data
    deallocate(Dcoords,IverticesAtElement)

    ! Perform mesh refinement
    do iref=1,nref

      ! Get mesh from adaptation structure
      allocate(Dcoords(ndim,nvt), IverticesAtElement(nnve,nel),Dind(nel))
      call madapt_getVertexCoords(rmeshAdapt, Dcoords)
      call madapt_getVerticesAtElement(rmeshAdapt, IverticesAtElement)

      do iel=1,nel
        ! Create circular indicator function
        xc = sum(Dcoords(1,IverticesAtElement(:,iel)))/real(nnve,DP)
        yc = sum(Dcoords(2,IverticesAtElement(:,iel)))/real(nnve,DP)
        phi = 0.2_DP - sqrt( (xc-0.5_DP)**2 + (yc-0.5_DP)**2 )

        ! Define indicator array
        if (abs(phi) .le. W) then
          Dind(iel) = 2.0_DP
        else
          Dind(iel) = 1.0_DP
        end if
      end do

      ! Perform one step of mesh adaptation
      istatus = madapt_step(rmeshAdapt,nel,Dind,nrefmax,dreftol,dcrstol)

      ! Deallocate local data
      deallocate(Dcoords,IverticesAtElement,Dind)

      ! Get data from adaptation structure
      nel = madapt_getnel(rmeshAdapt)
      nvt = madapt_getnvt(rmeshAdapt)
      nnve = madapt_getnnve(rmeshAdapt)

      print *, "REFINE:",nel,nvt

      W = W/2.0_DP
    end do

    ! Perform mesh re-coarsening
    do iref=1,nref

      ! Define indicator array
      allocate(Dind(nel)); Dind = 0.1

      ! Perform one step of mesh adaptation
      istatus = madapt_step(rmeshAdapt,nel,Dind,nrefmax,dreftol,dcrstol)

      ! Deallocate local data
      deallocate(Dind)

      ! Get data from adaptation structure
      nel = madapt_getnel(rmeshAdapt)
      nvt = madapt_getnvt(rmeshAdapt)

      print *, "COARSEN:",nel,nvt

      ! Get mesh from adaptation structure
      allocate(Dcoords(ndim,nvt), IverticesAtElement(nnve,nel))

      call madapt_getVertexCoords(rmeshAdapt, Dcoords)
      call madapt_getVerticesAtElement(rmeshAdapt, IverticesAtElement)

      ! Deallocate local data
      deallocate(Dcoords,IverticesAtElement)
    end do

    ! Finalisation
    call madapt_done(rmeshAdapt)

  end subroutine meshadaptQuick

end program meshadapt

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
fl

fswrite btstat.f
\ for the HC05
\ White wire - prop input, hi if connected
[ifndef hcConnect
	16 wconstant hcConnect
]
\ blue wire - prop output, hi to program device
[ifndef hcProgPin
	17 wconstant hcProgPin
]
\ blue wire - prop input, lo when hc05 wants to send data
[ifndef hcRts
	18 wconstant hcRts
]
\ blue wire - prop output, lo when hc05 can send data
[ifndef hcCts
	19 wconstant hcCts
]
\ red wire - prop output, data
[ifndef hcRx
	20 wconstant hcRx
]
\ white wire - prop input, data
[ifndef hcTx
	21 wconstant hcTx
]
[ifndef hcBaud
	d_230400 4/ wconstant hcBaud
\	d_9600 4/ wconstant hcBaud
]
[ifndef hcSerialCog
	0 wconstant hcSerialCog
]

[ifndef hcSend
: hcSend
	cogid hcSerialCog iolink .cstr 1000 delms cogid iounlink cr
;
]

\ hcConnect? ( -- t/f) bluetooth connected
[ifndef hcConnect?
: hcConnect?
	hcConnect dup pinin px?
;
]

\ hcRts? ( -- t/f) bluetooth connected
[ifndef hcRts?
: hcRts?
	hcRts dup pinin px?
;
]

{
: hcQ
	cogid hcSerialCog iolink
	begin
		hcConnect? 0=
		if
			c" NOT CONNECTED~h0D" .concstr
			d_1000 delms
		then
	fkey?
	if
		dup h_1B =
		if
			drop -1
		else
			dup h_30 =
			if
				c" CTS LO~h0D~h0A" hcSerialCog cogx
				hcCts pinlo
			else
				dup h_31 =
				if
					c" CTS hi~h0D~h0A" hcSerialCog cogx
					hcCts pinhi
				then
			then
			drop 0
		then
	else
		drop 0
	then
	until
	cogid iounlink
;
}

hcSerialCog cogreset 100 delms
c" hcRx hcTx hcBaud serial" hcSerialCog cogx 100 delms
1 hcSerialCog sersetflags
hcConnect? .

hcCts pinout hcCts pinlo

hcProgPin pinhi hcProgPin pinout

c" AT+VERSION~h0D~h0A" hcSend
c" AT+ADDR?~h0D~h0A" hcSend
c" AT+NAME?~h0D~h0A" hcSend
c" AT+RNAME?~h0D~h0A" hcSend
c" AT+ROLE?~h0D~h0A" hcSend
c" AT+CLASS?~h0D~h0A" hcSend
c" AT+IAC?~h0D~h0A" hcSend
c" AT+INQM?~h0D~h0A" hcSend
c" AT+PSWD?~h0D~h0A" hcSend
c" AT+UART?~h0D~h0A" hcSend
c" AT+CMODE?~h0D~h0A" hcSend
c" AT+BIND?~h0D~h0A" hcSend
c" AT+POLAR?~h0D~h0A" hcSend
c" AT+IPSCAN?~h0D~h0A" hcSend
c" AT+SNIFF?~h0D~h0A" hcSend
c" AT+SENM?~h0D~h0A" hcSend
c" AT+ADCN?~h0D~h0A" hcSend
c" AT+MRAD?~h0D~h0A" hcSend
c" AT+STATE?~h0D~h0A" hcSend

hcProgPin pinlo

...

fswrite btinit.f

\ for the HC05
\ White wire - prop input, hi if connected
[ifndef hcConnect
	16 wconstant hcConnect
]
\ blue wire - prop output, hi to program device
[ifndef hcProgPin
	17 wconstant hcProgPin
]
\ blue wire - prop input, lo when hc05 wants to send data
[ifndef hcRts
	18 wconstant hcRts
]
\ blue wire - prop output, lo when hc05 can send data
[ifndef hcCts
	19 wconstant hcCts
]
\ red wire - prop output, data
[ifndef hcRx
	20 wconstant hcRx
]
\ white wire - prop input, data
[ifndef hcTx
	21 wconstant hcTx
]

[ifndef hcBaud
	d_230400 4/ wconstant hcBaud
\	d_9600 4/ wconstant hcBaud
]
[ifndef hcSerialCog
	0 wconstant hcSerialCog
]

\ hcConnect? ( -- t/f) bluetooth connected
[ifndef hcConnect?
: hcConnect?
	hcConnect dup pinin px?
;
]

[ifndef hcSend
: hcSend
	cogid hcSerialCog iolink .cstr 1000 delms cogid iounlink cr
;

]

hcSerialCog cogreset 100 delms
c" hcRx hcTx hcBaud serial" hcSerialCog cogx 100 delms
1 hcSerialCog sersetflags
hcConnect? .

hcProgPin pinhi hcProgPin pinout

c" AT+VERSION~h0D~h0A" hcSend
c" AT+ROLE?~h0D~h0A" hcSend
c" AT+UART?~h0D~h0A" hcSend
c" AT+NAME?~h0D~h0A" hcSend
c" AT+PSWD?~h0D~h0A" hcSend
c" AT+NAME=CreativeCompanion~h0D~h0A" hcSend
c" AT+PSWD=password~h0D~h0A" hcSend
c" AT+UART=230400,0,0~h0D~h0A" hcSend
c" AT+VERSION~h0D~h0A" hcSend
c" AT+ROLE?~h0D~h0A" hcSend
c" AT+UART?~h0D~h0A" hcSend
c" AT+NAME?~h0D~h0A" hcSend
c" AT+PSWD?~h0D~h0A" hcSend

hcProgPin pinlo

\ after the hc05 is power cycled, the baud rate will change, not like the hc06

...

fswrite boot.f
hA state orC! version W@ .cstr cr cr cr

: getfile
	." ~h0D........~h0D"
	fsread
	." ~h0D,,,,,,,,~h0D"
;

\ findEETOP ( -- n1 ) the top of the eeprom + 1
: findEETOP
\
\ search 32k block increments until we get a fail
\
	0
	h100000 h8000
	do
		i t0 2 eereadpage
		if
			leave
		else
			i h7FFE + t0 3 eereadpage
			if
				leave
			else
				drop i h8000 +
			then
		then
	h8000 +loop
;

c" boot.f - Finding top of eeprom, " .cstr findEETOP ' fstop 2+ alignl L! forget findEETOP c" Top of eeprom at: " .cstr fstop . cr

c" boot.f - DONE PropForth Loaded~h0D~h0D" .cstr hA state andnC!

\ for the HC05
\ White wire - prop input, hi if connected
[ifndef hcConnect
	16 wconstant hcConnect
]
\ blue wire - prop output, hi to program device
[ifndef hcProgPin
	17 wconstant hcProgPin
]
\ blue wire - prop input, lo when hc05 wants to send data
[ifndef hcRts
	18 wconstant hcRts
]
\ blue wire - prop output, lo when hc05 can send data
[ifndef hcCts
	19 wconstant hcCts
]
\ red wire - prop output, data
[ifndef hcRx
	20 wconstant hcRx
]
\ white wire - prop input, data
[ifndef hcTx
	21 wconstant hcTx
]
[ifndef hcBaud
	d_230400 4/ wconstant hcBaud
\	d_9600 4/ wconstant hcBaud
]
[ifndef hcSerialCog
	0 wconstant hcSerialCog
]

[ifndef xsmDir
	d0 wconstant xsmDir
\ rising edge 2us cycletime max
	d_01 wconstant xsmStep
\ active lo
	d_02 wconstant xsmSlp
\ active lo
	d_03 wconstant xsmRst
	d_04 wconstant xsmMs3
	d_05 wconstant xsmMs2
	d_06 wconstant xsmMs1
\ active lo
	d_07 wconstant xsmEn
	d_24 wconstant xsmTemp
	d_25 wconstant xsmDTemp
]

[ifndef ysmDir
	d_08 wconstant ysmDir
\ rising edge 2us cycletime max
	d_09 wconstant ysmStep
\ active lo
	d_10 wconstant ysmSlp
\ active lo
	d_11 wconstant ysmRst
	d_12 wconstant ysmMs3
	d_13 wconstant ysmMs2
	d_14 wconstant ysmMs1
\ active lo
	d_15 wconstant ysmEn
	d_23 wconstant ysmTemp
	d_22 wconstant ysmDTemp
]

: smInit
	h_F3F3 outa COG! h_FFFF dira COG!
	10 delms h_FBFB outa COG!
	10 delms h_FFFF outa COG!
;

: xsmOn  xsmEn pinlo ;
: xsmOff xsmEn pinhi ;

: ysmOn  ysmEn pinlo ;
: ysmOff ysmEn pinhi ;

: _xStep xsmStep pinhi xsmStep pinlo ;

: xStep xsmDir over 0> if pinhi else pinlo then abs 0 do _xStep loop ;

: _yStep ysmStep pinhi ysmStep pinlo ;

: yStep ysmDir over 0> if pinhi else pinlo then abs 0 do _yStep loop ;

: _temp?
	dup pinlo dup pinout
	1 delms
	cnt COG@ over pinin
	swap >m dup waitpeq
	cnt COG@ swap -
;

: xsmTemp? 0 h10 0 do xsmTemp _temp? + loop h10 u/ ;

: ysmTemp? 0 h10 0 do ysmTemp _temp? + loop h10 u/ ;

: xsmDTemp? 0 h10 0 do xsmDTemp _temp? + loop h10 u/ ;

: ysmDTemp? 0 h10 0 do ysmDTemp _temp? + loop h10 u/ ;

: temp
	begin
		." x: "       xsmTemp? .long space xsmDTemp? .long
		."       y: " ysmTemp? .long space ysmDTemp? .long cr
		fkey?
		if
			h1B =
		else
			drop 0
		then
	until
;

[ifndef _ready_led
	27 wconstant _ready_led
]

: prompt
	c" CreativeCompanion-" prop W!
;

: onreset6
	fkey? and fkey? and or h1B <>
	if
		smInit
		prompt
		$S_con iodis $S_con cogreset 100 delms
		c" _ready_led dup pinout pinhi" $S_con cogx
		c" hcRx hcTx 57600 serial" $S_con cogx 100 delms
		cogid >con
	then
	c" onreset6" (forget)
;

...

ext: f


!*==strevc.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b STREVC
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download STREVC + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE STREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
!                          LDVR, MM, M, WORK, INFO )
!
!       .. Scalar Arguments ..
!       CHARACTER          HOWMNY, SIDE
!       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
!       ..
!       .. Array Arguments ..
!       LOGICAL            SELECT( * )
!       REAL               T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
!      $                   WORK( * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> STREVC computes some or all of the right and/or left eigenvectors of
!> a real upper quasi-triangular matrix T.
!> Matrices of this type are produced by the Schur factorization of
!> a real general matrix:  A = Q*T*Q**T, as computed by SHSEQR.
!>
!> The right eigenvector x and the left eigenvector y of T corresponding
!> to an eigenvalue w are defined by:
!>
!>    T*x = w*x,     (y**H)*T = w*(y**H)
!>
!> where y**H denotes the conjugate transpose of y.
!> The eigenvalues are not input to this routine, but are read directly
!> from the diagonal blocks of T.
!>
!> This routine returns the matrices X and/or Y of right and left
!> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
!> input matrix.  If Q is the orthogonal factor that reduces a matrix
!> A to Schur form T, then Q*X and Q*Y are the matrices of right and
!> left eigenvectors of A.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] SIDE
!> \verbatim
!>          SIDE is CHARACTER*1
!>          = 'R':  compute right eigenvectors only;
!>          = 'L':  compute left eigenvectors only;
!>          = 'B':  compute both right and left eigenvectors.
!> \endverbatim
!>
!> \param[in] HOWMNY
!> \verbatim
!>          HOWMNY is CHARACTER*1
!>          = 'A':  compute all right and/or left eigenvectors;
!>          = 'B':  compute all right and/or left eigenvectors,
!>                  backtransformed by the matrices in VR and/or VL;
!>          = 'S':  compute selected right and/or left eigenvectors,
!>                  as indicated by the logical array SELECT.
!> \endverbatim
!>
!> \param[in,out] SELECT
!> \verbatim
!>          SELECT is LOGICAL array, dimension (N)
!>          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
!>          computed.
!>          If w(j) is a real eigenvalue, the corresponding real
!>          eigenvector is computed if SELECT(j) is .TRUE..
!>          If w(j) and w(j+1) are the real and imaginary parts of a
!>          complex eigenvalue, the corresponding complex eigenvector is
!>          computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
!>          on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
!>          .FALSE..
!>          Not referenced if HOWMNY = 'A' or 'B'.
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the matrix T. N >= 0.
!> \endverbatim
!>
!> \param[in] T
!> \verbatim
!>          T is REAL array, dimension (LDT,N)
!>          The upper quasi-triangular matrix T in Schur canonical form.
!> \endverbatim
!>
!> \param[in] LDT
!> \verbatim
!>          LDT is INTEGER
!>          The leading dimension of the array T. LDT >= max(1,N).
!> \endverbatim
!>
!> \param[in,out] VL
!> \verbatim
!>          VL is REAL array, dimension (LDVL,MM)
!>          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
!>          contain an N-by-N matrix Q (usually the orthogonal matrix Q
!>          of Schur vectors returned by SHSEQR).
!>          On exit, if SIDE = 'L' or 'B', VL contains:
!>          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
!>          if HOWMNY = 'B', the matrix Q*Y;
!>          if HOWMNY = 'S', the left eigenvectors of T specified by
!>                           SELECT, stored consecutively in the columns
!>                           of VL, in the same order as their
!>                           eigenvalues.
!>          A complex eigenvector corresponding to a complex eigenvalue
!>          is stored in two consecutive columns, the first holding the
!>          real part, and the second the imaginary part.
!>          Not referenced if SIDE = 'R'.
!> \endverbatim
!>
!> \param[in] LDVL
!> \verbatim
!>          LDVL is INTEGER
!>          The leading dimension of the array VL.  LDVL >= 1, and if
!>          SIDE = 'L' or 'B', LDVL >= N.
!> \endverbatim
!>
!> \param[in,out] VR
!> \verbatim
!>          VR is REAL array, dimension (LDVR,MM)
!>          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
!>          contain an N-by-N matrix Q (usually the orthogonal matrix Q
!>          of Schur vectors returned by SHSEQR).
!>          On exit, if SIDE = 'R' or 'B', VR contains:
!>          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
!>          if HOWMNY = 'B', the matrix Q*X;
!>          if HOWMNY = 'S', the right eigenvectors of T specified by
!>                           SELECT, stored consecutively in the columns
!>                           of VR, in the same order as their
!>                           eigenvalues.
!>          A complex eigenvector corresponding to a complex eigenvalue
!>          is stored in two consecutive columns, the first holding the
!>          real part and the second the imaginary part.
!>          Not referenced if SIDE = 'L'.
!> \endverbatim
!>
!> \param[in] LDVR
!> \verbatim
!>          LDVR is INTEGER
!>          The leading dimension of the array VR.  LDVR >= 1, and if
!>          SIDE = 'R' or 'B', LDVR >= N.
!> \endverbatim
!>
!> \param[in] MM
!> \verbatim
!>          MM is INTEGER
!>          The number of columns in the arrays VL and/or VR. MM >= M.
!> \endverbatim
!>
!> \param[out] M
!> \verbatim
!>          M is INTEGER
!>          The number of columns in the arrays VL and/or VR actually
!>          used to store the eigenvectors.
!>          If HOWMNY = 'A' or 'B', M is set to N.
!>          Each selected real eigenvector occupies one column and each
!>          selected complex eigenvector occupies two columns.
!> \endverbatim
!>
!> \param[out] WORK
!> \verbatim
!>          WORK is REAL array, dimension (3*N)
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date December 2016
!
!> \ingroup realOTHERcomputational
!
!> \par Further Details:
!  =====================
!>
!> \verbatim
!>
!>  The algorithm used in this program is basically backward (forward)
!>  substitution, with scaling to make the the code robust against
!>  possible overflow.
!>
!>  Each eigenvector is normalized so that the element of largest
!>  magnitude has magnitude 1; here the magnitude of a complex number
!>  (x,y) is taken to be |x| + |y|.
!> \endverbatim
!>
!  =====================================================================
      SUBROUTINE STREVC(Side,Howmny,Select,N,T,Ldt,Vl,Ldvl,Vr,Ldvr,Mm,M,&
     &                  Work,Info)
      IMPLICIT NONE
!*--STREVC226
!
!  -- LAPACK computational routine (version 3.7.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     December 2016
!
!     .. Scalar Arguments ..
      CHARACTER Howmny , Side
      INTEGER Info , Ldt , Ldvl , Ldvr , M , Mm , N
!     ..
!     .. Array Arguments ..
      LOGICAL Select(*)
      REAL T(Ldt,*) , Vl(Ldvl,*) , Vr(Ldvr,*) , Work(*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      REAL ZERO , ONE
      PARAMETER (ZERO=0.0E+0,ONE=1.0E+0)
!     ..
!     .. Local Scalars ..
      LOGICAL allv , bothv , leftv , over , pair , rightv , somev
      INTEGER i , ierr , ii , ip , is , j , j1 , j2 , jnxt , k , ki , n2
      REAL beta , bignum , emax , ovfl , rec , remax , scale , smin ,   &
     &     smlnum , ulp , unfl , vcrit , vmax , wi , wr , xnorm
!     ..
!     .. External Functions ..
      LOGICAL LSAME
      INTEGER ISAMAX
      REAL SDOT , SLAMCH
      EXTERNAL LSAME , ISAMAX , SDOT , SLAMCH
!     ..
!     .. External Subroutines ..
      EXTERNAL SAXPY , SCOPY , SGEMV , SLABAD , SLALN2 , SSCAL , XERBLA
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC ABS , MAX , SQRT
!     ..
!     .. Local Arrays ..
      REAL x(2,2)
!     ..
!     .. Executable Statements ..
!
!     Decode and test the input parameters
!
      bothv = LSAME(Side,'B')
      rightv = LSAME(Side,'R') .OR. bothv
      leftv = LSAME(Side,'L') .OR. bothv
!
      allv = LSAME(Howmny,'A')
      over = LSAME(Howmny,'B')
      somev = LSAME(Howmny,'S')
!
      Info = 0
      IF ( .NOT.rightv .AND. .NOT.leftv ) THEN
         Info = -1
      ELSEIF ( .NOT.allv .AND. .NOT.over .AND. .NOT.somev ) THEN
         Info = -2
      ELSEIF ( N<0 ) THEN
         Info = -4
      ELSEIF ( Ldt<MAX(1,N) ) THEN
         Info = -6
      ELSEIF ( Ldvl<1 .OR. (leftv .AND. Ldvl<N) ) THEN
         Info = -8
      ELSEIF ( Ldvr<1 .OR. (rightv .AND. Ldvr<N) ) THEN
         Info = -10
      ELSE
!
!        Set M to the number of columns required to store the selected
!        eigenvectors, standardize the array SELECT if necessary, and
!        test MM.
!
         IF ( somev ) THEN
            M = 0
            pair = .FALSE.
            DO j = 1 , N
               IF ( pair ) THEN
                  pair = .FALSE.
                  Select(j) = .FALSE.
               ELSEIF ( j<N ) THEN
                  IF ( T(j+1,j)==ZERO ) THEN
                     IF ( Select(j) ) M = M + 1
                  ELSE
                     pair = .TRUE.
                     IF ( Select(j) .OR. Select(j+1) ) THEN
                        Select(j) = .TRUE.
                        M = M + 2
                     ENDIF
                  ENDIF
               ELSE
                  IF ( Select(N) ) M = M + 1
               ENDIF
            ENDDO
         ELSE
            M = N
         ENDIF
!
         IF ( Mm<M ) Info = -11
      ENDIF
      IF ( Info/=0 ) THEN
         CALL XERBLA('STREVC',-Info)
         RETURN
      ENDIF
!
!     Quick return if possible.
!
      IF ( N==0 ) RETURN
!
!     Set the constants to control overflow.
!
      unfl = SLAMCH('Safe minimum')
      ovfl = ONE/unfl
      CALL SLABAD(unfl,ovfl)
      ulp = SLAMCH('Precision')
      smlnum = unfl*(N/ulp)
      bignum = (ONE-ulp)/smlnum
!
!     Compute 1-norm of each column of strictly upper triangular
!     part of T to control overflow in triangular solver.
!
      Work(1) = ZERO
      DO j = 2 , N
         Work(j) = ZERO
         DO i = 1 , j - 1
            Work(j) = Work(j) + ABS(T(i,j))
         ENDDO
      ENDDO
!
!     Index IP is used to specify the real or complex eigenvalue:
!       IP = 0, real eigenvalue,
!            1, first of conjugate complex pair: (wr,wi)
!           -1, second of conjugate complex pair: (wr,wi)
!
      n2 = 2*N
!
      IF ( rightv ) THEN
!
!        Compute right eigenvectors.
!
         ip = 0
         is = M
         DO ki = N , 1 , -1
!
            IF ( ip/=1 ) THEN
               IF ( ki/=1 ) THEN
                  IF ( T(ki,ki-1)/=ZERO ) ip = -1
               ENDIF
!
               IF ( somev ) THEN
                  IF ( ip==0 ) THEN
                     IF ( .NOT.Select(ki) ) GOTO 20
                  ELSEIF ( .NOT.Select(ki-1) ) THEN
                     GOTO 20
                  ENDIF
               ENDIF
!
!           Compute the KI-th eigenvalue (WR,WI).
!
               wr = T(ki,ki)
               wi = ZERO
               IF ( ip/=0 ) wi = SQRT(ABS(T(ki,ki-1)))                  &
     &                           *SQRT(ABS(T(ki-1,ki)))
               smin = MAX(ulp*(ABS(wr)+ABS(wi)),smlnum)
!
               IF ( ip==0 ) THEN
!
!              Real right eigenvector
!
                  Work(ki+N) = ONE
!
!              Form right-hand side
!
                  DO k = 1 , ki - 1
                     Work(k+N) = -T(k,ki)
                  ENDDO
!
!              Solve the upper quasi-triangular system:
!                 (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
!
                  jnxt = ki - 1
                  DO j = ki - 1 , 1 , -1
                     IF ( j<=jnxt ) THEN
                        j1 = j
                        j2 = j
                        jnxt = j - 1
                        IF ( j>1 ) THEN
                           IF ( T(j,j-1)/=ZERO ) THEN
                              j1 = j - 1
                              jnxt = j - 2
                           ENDIF
                        ENDIF
!
                        IF ( j1==j2 ) THEN
!
!                    1-by-1 diagonal block
!
                           CALL SLALN2(.FALSE.,1,1,smin,ONE,T(j,j),Ldt, &
     &                                 ONE,ONE,Work(j+N),N,wr,ZERO,x,2, &
     &                                 scale,xnorm,ierr)
!
!                    Scale X(1,1) to avoid overflow when updating
!                    the right-hand side.
!
                           IF ( xnorm>ONE ) THEN
                              IF ( Work(j)>bignum/xnorm ) THEN
                                 x(1,1) = x(1,1)/xnorm
                                 scale = scale/xnorm
                              ENDIF
                           ENDIF
!
!                    Scale if necessary
!
                           IF ( scale/=ONE )                            &
     &                          CALL SSCAL(ki,scale,Work(1+N),1)
                           Work(j+N) = x(1,1)
!
!                    Update right-hand side
!
                           CALL SAXPY(j-1,-x(1,1),T(1,j),1,Work(1+N),1)
!
                        ELSE
!
!                    2-by-2 diagonal block
!
                           CALL SLALN2(.FALSE.,2,1,smin,ONE,T(j-1,j-1), &
     &                                 Ldt,ONE,ONE,Work(j-1+N),N,wr,    &
     &                                 ZERO,x,2,scale,xnorm,ierr)
!
!                    Scale X(1,1) and X(2,1) to avoid overflow when
!                    updating the right-hand side.
!
                           IF ( xnorm>ONE ) THEN
                              beta = MAX(Work(j-1),Work(j))
                              IF ( beta>bignum/xnorm ) THEN
                                 x(1,1) = x(1,1)/xnorm
                                 x(2,1) = x(2,1)/xnorm
                                 scale = scale/xnorm
                              ENDIF
                           ENDIF
!
!                    Scale if necessary
!
                           IF ( scale/=ONE )                            &
     &                          CALL SSCAL(ki,scale,Work(1+N),1)
                           Work(j-1+N) = x(1,1)
                           Work(j+N) = x(2,1)
!
!                    Update right-hand side
!
                           CALL SAXPY(j-2,-x(1,1),T(1,j-1),1,Work(1+N), &
     &                                1)
                           CALL SAXPY(j-2,-x(2,1),T(1,j),1,Work(1+N),1)
                        ENDIF
                     ENDIF
                  ENDDO
!
!              Copy the vector x or Q*x to VR and normalize.
!
                  IF ( .NOT.over ) THEN
                     CALL SCOPY(ki,Work(1+N),1,Vr(1,is),1)
!
                     ii = ISAMAX(ki,Vr(1,is),1)
                     remax = ONE/ABS(Vr(ii,is))
                     CALL SSCAL(ki,remax,Vr(1,is),1)
!
                     DO k = ki + 1 , N
                        Vr(k,is) = ZERO
                     ENDDO
                  ELSE
                     IF ( ki>1 )                                        &
     &                    CALL SGEMV('N',N,ki-1,ONE,Vr,Ldvr,Work(1+N),1,&
     &                    Work(ki+N),Vr(1,ki),1)
!
                     ii = ISAMAX(N,Vr(1,ki),1)
                     remax = ONE/ABS(Vr(ii,ki))
                     CALL SSCAL(N,remax,Vr(1,ki),1)
                  ENDIF
!
               ELSE
!
!              Complex right eigenvector.
!
!              Initial solve
!                [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
!                [ (T(KI,KI-1)   T(KI,KI)   )               ]
!
                  IF ( ABS(T(ki-1,ki))>=ABS(T(ki,ki-1)) ) THEN
                     Work(ki-1+N) = ONE
                     Work(ki+n2) = wi/T(ki-1,ki)
                  ELSE
                     Work(ki-1+N) = -wi/T(ki,ki-1)
                     Work(ki+n2) = ONE
                  ENDIF
                  Work(ki+N) = ZERO
                  Work(ki-1+n2) = ZERO
!
!              Form right-hand side
!
                  DO k = 1 , ki - 2
                     Work(k+N) = -Work(ki-1+N)*T(k,ki-1)
                     Work(k+n2) = -Work(ki+n2)*T(k,ki)
                  ENDDO
!
!              Solve upper quasi-triangular system:
!              (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
!
                  jnxt = ki - 2
                  DO j = ki - 2 , 1 , -1
                     IF ( j<=jnxt ) THEN
                        j1 = j
                        j2 = j
                        jnxt = j - 1
                        IF ( j>1 ) THEN
                           IF ( T(j,j-1)/=ZERO ) THEN
                              j1 = j - 1
                              jnxt = j - 2
                           ENDIF
                        ENDIF
!
                        IF ( j1==j2 ) THEN
!
!                    1-by-1 diagonal block
!
                           CALL SLALN2(.FALSE.,1,2,smin,ONE,T(j,j),Ldt, &
     &                                 ONE,ONE,Work(j+N),N,wr,wi,x,2,   &
     &                                 scale,xnorm,ierr)
!
!                    Scale X(1,1) and X(1,2) to avoid overflow when
!                    updating the right-hand side.
!
                           IF ( xnorm>ONE ) THEN
                              IF ( Work(j)>bignum/xnorm ) THEN
                                 x(1,1) = x(1,1)/xnorm
                                 x(1,2) = x(1,2)/xnorm
                                 scale = scale/xnorm
                              ENDIF
                           ENDIF
!
!                    Scale if necessary
!
                           IF ( scale/=ONE ) THEN
                              CALL SSCAL(ki,scale,Work(1+N),1)
                              CALL SSCAL(ki,scale,Work(1+n2),1)
                           ENDIF
                           Work(j+N) = x(1,1)
                           Work(j+n2) = x(1,2)
!
!                    Update the right-hand side
!
                           CALL SAXPY(j-1,-x(1,1),T(1,j),1,Work(1+N),1)
                           CALL SAXPY(j-1,-x(1,2),T(1,j),1,Work(1+n2),1)
!
                        ELSE
!
!                    2-by-2 diagonal block
!
                           CALL SLALN2(.FALSE.,2,2,smin,ONE,T(j-1,j-1), &
     &                                 Ldt,ONE,ONE,Work(j-1+N),N,wr,wi, &
     &                                 x,2,scale,xnorm,ierr)
!
!                    Scale X to avoid overflow when updating
!                    the right-hand side.
!
                           IF ( xnorm>ONE ) THEN
                              beta = MAX(Work(j-1),Work(j))
                              IF ( beta>bignum/xnorm ) THEN
                                 rec = ONE/xnorm
                                 x(1,1) = x(1,1)*rec
                                 x(1,2) = x(1,2)*rec
                                 x(2,1) = x(2,1)*rec
                                 x(2,2) = x(2,2)*rec
                                 scale = scale*rec
                              ENDIF
                           ENDIF
!
!                    Scale if necessary
!
                           IF ( scale/=ONE ) THEN
                              CALL SSCAL(ki,scale,Work(1+N),1)
                              CALL SSCAL(ki,scale,Work(1+n2),1)
                           ENDIF
                           Work(j-1+N) = x(1,1)
                           Work(j+N) = x(2,1)
                           Work(j-1+n2) = x(1,2)
                           Work(j+n2) = x(2,2)
!
!                    Update the right-hand side
!
                           CALL SAXPY(j-2,-x(1,1),T(1,j-1),1,Work(1+N), &
     &                                1)
                           CALL SAXPY(j-2,-x(2,1),T(1,j),1,Work(1+N),1)
                           CALL SAXPY(j-2,-x(1,2),T(1,j-1),1,Work(1+n2),&
     &                                1)
                           CALL SAXPY(j-2,-x(2,2),T(1,j),1,Work(1+n2),1)
                        ENDIF
                     ENDIF
                  ENDDO
!
!              Copy the vector x or Q*x to VR and normalize.
!
                  IF ( .NOT.over ) THEN
                     CALL SCOPY(ki,Work(1+N),1,Vr(1,is-1),1)
                     CALL SCOPY(ki,Work(1+n2),1,Vr(1,is),1)
!
                     emax = ZERO
                     DO k = 1 , ki
                        emax = MAX(emax,ABS(Vr(k,is-1))+ABS(Vr(k,is)))
                     ENDDO
!
                     remax = ONE/emax
                     CALL SSCAL(ki,remax,Vr(1,is-1),1)
                     CALL SSCAL(ki,remax,Vr(1,is),1)
!
                     DO k = ki + 1 , N
                        Vr(k,is-1) = ZERO
                        Vr(k,is) = ZERO
                     ENDDO
!
                  ELSE
!
                     IF ( ki>2 ) THEN
                        CALL SGEMV('N',N,ki-2,ONE,Vr,Ldvr,Work(1+N),1,  &
     &                             Work(ki-1+N),Vr(1,ki-1),1)
                        CALL SGEMV('N',N,ki-2,ONE,Vr,Ldvr,Work(1+n2),1, &
     &                             Work(ki+n2),Vr(1,ki),1)
                     ELSE
                        CALL SSCAL(N,Work(ki-1+N),Vr(1,ki-1),1)
                        CALL SSCAL(N,Work(ki+n2),Vr(1,ki),1)
                     ENDIF
!
                     emax = ZERO
                     DO k = 1 , N
                        emax = MAX(emax,ABS(Vr(k,ki-1))+ABS(Vr(k,ki)))
                     ENDDO
                     remax = ONE/emax
                     CALL SSCAL(N,remax,Vr(1,ki-1),1)
                     CALL SSCAL(N,remax,Vr(1,ki),1)
                  ENDIF
               ENDIF
!
               is = is - 1
               IF ( ip/=0 ) is = is - 1
            ENDIF
 20         IF ( ip==1 ) ip = 0
            IF ( ip==-1 ) ip = 1
         ENDDO
      ENDIF
!
      IF ( leftv ) THEN
!
!        Compute left eigenvectors.
!
         ip = 0
         is = 1
         DO ki = 1 , N
!
            IF ( ip/=-1 ) THEN
               IF ( ki/=N ) THEN
                  IF ( T(ki+1,ki)/=ZERO ) ip = 1
               ENDIF
!
               IF ( somev ) THEN
                  IF ( .NOT.Select(ki) ) GOTO 40
               ENDIF
!
!           Compute the KI-th eigenvalue (WR,WI).
!
               wr = T(ki,ki)
               wi = ZERO
               IF ( ip/=0 ) wi = SQRT(ABS(T(ki,ki+1)))                  &
     &                           *SQRT(ABS(T(ki+1,ki)))
               smin = MAX(ulp*(ABS(wr)+ABS(wi)),smlnum)
!
               IF ( ip==0 ) THEN
!
!              Real left eigenvector.
!
                  Work(ki+N) = ONE
!
!              Form right-hand side
!
                  DO k = ki + 1 , N
                     Work(k+N) = -T(ki,k)
                  ENDDO
!
!              Solve the quasi-triangular system:
!                 (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK
!
                  vmax = ONE
                  vcrit = bignum
!
                  jnxt = ki + 1
                  DO j = ki + 1 , N
                     IF ( j>=jnxt ) THEN
                        j1 = j
                        j2 = j
                        jnxt = j + 1
                        IF ( j<N ) THEN
                           IF ( T(j+1,j)/=ZERO ) THEN
                              j2 = j + 1
                              jnxt = j + 2
                           ENDIF
                        ENDIF
!
                        IF ( j1==j2 ) THEN
!
!                    1-by-1 diagonal block
!
!                    Scale if necessary to avoid overflow when forming
!                    the right-hand side.
!
                           IF ( Work(j)>vcrit ) THEN
                              rec = ONE/vmax
                              CALL SSCAL(N-ki+1,rec,Work(ki+N),1)
                              vmax = ONE
                              vcrit = bignum
                           ENDIF
!
                           Work(j+N) = Work(j+N)                        &
     &                                 - SDOT(j-ki-1,T(ki+1,j),1,       &
     &                                 Work(ki+1+N),1)
!
!                    Solve (T(J,J)-WR)**T*X = WORK
!
                           CALL SLALN2(.FALSE.,1,1,smin,ONE,T(j,j),Ldt, &
     &                                 ONE,ONE,Work(j+N),N,wr,ZERO,x,2, &
     &                                 scale,xnorm,ierr)
!
!                    Scale if necessary
!
                           IF ( scale/=ONE )                            &
     &                          CALL SSCAL(N-ki+1,scale,Work(ki+N),1)
                           Work(j+N) = x(1,1)
                           vmax = MAX(ABS(Work(j+N)),vmax)
                           vcrit = bignum/vmax
!
                        ELSE
!
!                    2-by-2 diagonal block
!
!                    Scale if necessary to avoid overflow when forming
!                    the right-hand side.
!
                           beta = MAX(Work(j),Work(j+1))
                           IF ( beta>vcrit ) THEN
                              rec = ONE/vmax
                              CALL SSCAL(N-ki+1,rec,Work(ki+N),1)
                              vmax = ONE
                              vcrit = bignum
                           ENDIF
!
                           Work(j+N) = Work(j+N)                        &
     &                                 - SDOT(j-ki-1,T(ki+1,j),1,       &
     &                                 Work(ki+1+N),1)
!
                           Work(j+1+N) = Work(j+1+N)                    &
     &                        - SDOT(j-ki-1,T(ki+1,j+1),1,Work(ki+1+N), &
     &                        1)
!
!                    Solve
!                      [T(J,J)-WR   T(J,J+1)     ]**T* X = SCALE*( WORK1 )
!                      [T(J+1,J)    T(J+1,J+1)-WR]               ( WORK2 )
!
                           CALL SLALN2(.TRUE.,2,1,smin,ONE,T(j,j),Ldt,  &
     &                                 ONE,ONE,Work(j+N),N,wr,ZERO,x,2, &
     &                                 scale,xnorm,ierr)
!
!                    Scale if necessary
!
                           IF ( scale/=ONE )                            &
     &                          CALL SSCAL(N-ki+1,scale,Work(ki+N),1)
                           Work(j+N) = x(1,1)
                           Work(j+1+N) = x(2,1)
!
                           vmax = MAX(ABS(Work(j+N)),ABS(Work(j+1+N)),  &
     &                            vmax)
                           vcrit = bignum/vmax
!
                        ENDIF
                     ENDIF
                  ENDDO
!
!              Copy the vector x or Q*x to VL and normalize.
!
                  IF ( .NOT.over ) THEN
                     CALL SCOPY(N-ki+1,Work(ki+N),1,Vl(ki,is),1)
!
                     ii = ISAMAX(N-ki+1,Vl(ki,is),1) + ki - 1
                     remax = ONE/ABS(Vl(ii,is))
                     CALL SSCAL(N-ki+1,remax,Vl(ki,is),1)
!
                     DO k = 1 , ki - 1
                        Vl(k,is) = ZERO
                     ENDDO
!
                  ELSE
!
                     IF ( ki<N ) CALL SGEMV('N',N,N-ki,ONE,Vl(1,ki+1),  &
     &                    Ldvl,Work(ki+1+N),1,Work(ki+N),Vl(1,ki),1)
!
                     ii = ISAMAX(N,Vl(1,ki),1)
                     remax = ONE/ABS(Vl(ii,ki))
                     CALL SSCAL(N,remax,Vl(1,ki),1)
!
                  ENDIF
!
               ELSE
!
!              Complex left eigenvector.
!
!               Initial solve:
!                 ((T(KI,KI)    T(KI,KI+1) )**T - (WR - I* WI))*X = 0.
!                 ((T(KI+1,KI) T(KI+1,KI+1))                )
!
                  IF ( ABS(T(ki,ki+1))>=ABS(T(ki+1,ki)) ) THEN
                     Work(ki+N) = wi/T(ki,ki+1)
                     Work(ki+1+n2) = ONE
                  ELSE
                     Work(ki+N) = ONE
                     Work(ki+1+n2) = -wi/T(ki+1,ki)
                  ENDIF
                  Work(ki+1+N) = ZERO
                  Work(ki+n2) = ZERO
!
!              Form right-hand side
!
                  DO k = ki + 2 , N
                     Work(k+N) = -Work(ki+N)*T(ki,k)
                     Work(k+n2) = -Work(ki+1+n2)*T(ki+1,k)
                  ENDDO
!
!              Solve complex quasi-triangular system:
!              ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
!
                  vmax = ONE
                  vcrit = bignum
!
                  jnxt = ki + 2
                  DO j = ki + 2 , N
                     IF ( j>=jnxt ) THEN
                        j1 = j
                        j2 = j
                        jnxt = j + 1
                        IF ( j<N ) THEN
                           IF ( T(j+1,j)/=ZERO ) THEN
                              j2 = j + 1
                              jnxt = j + 2
                           ENDIF
                        ENDIF
!
                        IF ( j1==j2 ) THEN
!
!                    1-by-1 diagonal block
!
!                    Scale if necessary to avoid overflow when
!                    forming the right-hand side elements.
!
                           IF ( Work(j)>vcrit ) THEN
                              rec = ONE/vmax
                              CALL SSCAL(N-ki+1,rec,Work(ki+N),1)
                              CALL SSCAL(N-ki+1,rec,Work(ki+n2),1)
                              vmax = ONE
                              vcrit = bignum
                           ENDIF
!
                           Work(j+N) = Work(j+N)                        &
     &                                 - SDOT(j-ki-2,T(ki+2,j),1,       &
     &                                 Work(ki+2+N),1)
                           Work(j+n2) = Work(j+n2)                      &
     &                                  - SDOT(j-ki-2,T(ki+2,j),1,      &
     &                                  Work(ki+2+n2),1)
!
!                    Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
!
                           CALL SLALN2(.FALSE.,1,2,smin,ONE,T(j,j),Ldt, &
     &                                 ONE,ONE,Work(j+N),N,wr,-wi,x,2,  &
     &                                 scale,xnorm,ierr)
!
!                    Scale if necessary
!
                           IF ( scale/=ONE ) THEN
                              CALL SSCAL(N-ki+1,scale,Work(ki+N),1)
                              CALL SSCAL(N-ki+1,scale,Work(ki+n2),1)
                           ENDIF
                           Work(j+N) = x(1,1)
                           Work(j+n2) = x(1,2)
                           vmax = MAX(ABS(Work(j+N)),ABS(Work(j+n2)),   &
     &                            vmax)
                           vcrit = bignum/vmax
!
                        ELSE
!
!                    2-by-2 diagonal block
!
!                    Scale if necessary to avoid overflow when forming
!                    the right-hand side elements.
!
                           beta = MAX(Work(j),Work(j+1))
                           IF ( beta>vcrit ) THEN
                              rec = ONE/vmax
                              CALL SSCAL(N-ki+1,rec,Work(ki+N),1)
                              CALL SSCAL(N-ki+1,rec,Work(ki+n2),1)
                              vmax = ONE
                              vcrit = bignum
                           ENDIF
!
                           Work(j+N) = Work(j+N)                        &
     &                                 - SDOT(j-ki-2,T(ki+2,j),1,       &
     &                                 Work(ki+2+N),1)
!
                           Work(j+n2) = Work(j+n2)                      &
     &                                  - SDOT(j-ki-2,T(ki+2,j),1,      &
     &                                  Work(ki+2+n2),1)
!
                           Work(j+1+N) = Work(j+1+N)                    &
     &                        - SDOT(j-ki-2,T(ki+2,j+1),1,Work(ki+2+N), &
     &                        1)
!
                           Work(j+1+n2) = Work(j+1+n2)                  &
     &                        - SDOT(j-ki-2,T(ki+2,j+1),1,Work(ki+2+n2),&
     &                        1)
!
!                    Solve 2-by-2 complex linear equation
!                      ([T(j,j)   T(j,j+1)  ]**T-(wr-i*wi)*I)*X = SCALE*B
!                      ([T(j+1,j) T(j+1,j+1)]               )
!
                           CALL SLALN2(.TRUE.,2,2,smin,ONE,T(j,j),Ldt,  &
     &                                 ONE,ONE,Work(j+N),N,wr,-wi,x,2,  &
     &                                 scale,xnorm,ierr)
!
!                    Scale if necessary
!
                           IF ( scale/=ONE ) THEN
                              CALL SSCAL(N-ki+1,scale,Work(ki+N),1)
                              CALL SSCAL(N-ki+1,scale,Work(ki+n2),1)
                           ENDIF
                           Work(j+N) = x(1,1)
                           Work(j+n2) = x(1,2)
                           Work(j+1+N) = x(2,1)
                           Work(j+1+n2) = x(2,2)
                           vmax = MAX(ABS(x(1,1)),ABS(x(1,2)),          &
     &                            ABS(x(2,1)),ABS(x(2,2)),vmax)
                           vcrit = bignum/vmax
!
                        ENDIF
                     ENDIF
                  ENDDO
!
!              Copy the vector x or Q*x to VL and normalize.
!
                  IF ( .NOT.over ) THEN
                     CALL SCOPY(N-ki+1,Work(ki+N),1,Vl(ki,is),1)
                     CALL SCOPY(N-ki+1,Work(ki+n2),1,Vl(ki,is+1),1)
!
                     emax = ZERO
                     DO k = ki , N
                        emax = MAX(emax,ABS(Vl(k,is))+ABS(Vl(k,is+1)))
                     ENDDO
                     remax = ONE/emax
                     CALL SSCAL(N-ki+1,remax,Vl(ki,is),1)
                     CALL SSCAL(N-ki+1,remax,Vl(ki,is+1),1)
!
                     DO k = 1 , ki - 1
                        Vl(k,is) = ZERO
                        Vl(k,is+1) = ZERO
                     ENDDO
                  ELSE
                     IF ( ki<N-1 ) THEN
                        CALL SGEMV('N',N,N-ki-1,ONE,Vl(1,ki+2),Ldvl,    &
     &                             Work(ki+2+N),1,Work(ki+N),Vl(1,ki),1)
                        CALL SGEMV('N',N,N-ki-1,ONE,Vl(1,ki+2),Ldvl,    &
     &                             Work(ki+2+n2),1,Work(ki+1+n2),       &
     &                             Vl(1,ki+1),1)
                     ELSE
                        CALL SSCAL(N,Work(ki+N),Vl(1,ki),1)
                        CALL SSCAL(N,Work(ki+1+n2),Vl(1,ki+1),1)
                     ENDIF
!
                     emax = ZERO
                     DO k = 1 , N
                        emax = MAX(emax,ABS(Vl(k,ki))+ABS(Vl(k,ki+1)))
                     ENDDO
                     remax = ONE/emax
                     CALL SSCAL(N,remax,Vl(1,ki),1)
                     CALL SSCAL(N,remax,Vl(1,ki+1),1)
!
                  ENDIF
!
               ENDIF
!
               is = is + 1
               IF ( ip/=0 ) is = is + 1
            ENDIF
 40         IF ( ip==-1 ) ip = 0
            IF ( ip==1 ) ip = -1
!
         ENDDO
!
      ENDIF
!
!
!     End of STREVC
!
      END SUBROUTINE STREVC


program simpson_test_wf
      implicit none
      integer, parameter :: n = 1001, dp = selected_real_kind(14)
      integer :: i
      real(kind=dp), dimension(n) :: x, y
      real(kind=dp) :: sum, high, low, dx, x_sq
      real(kind=dp), parameter :: pi = real(4, kind=dp)*atan(real(1, kind=dp))

      interface
        function simpson(x,y,n)
                integer, parameter :: dp = selected_real_kind(14)
                integer, intent(in) :: n
                real(kind=dp) :: simpson
                real(kind=dp), dimension(n), intent(in) :: x, y
        end function simpson
      end interface

      sum=0.0d0
      high = 5.0d0; low = -high;
      dx = (high - low)/(n-1)
      do i=1, n
        x(i) = low + dfloat(i)*dx
        x_sq = x(i)*x(i)
        y(i) = (x_sq - 1.0d0)*exp(-x_sq)
        !write(*, '(i, 2f10.3)') i, x(i), y(i)
      end do

      sum = simpson(x,y,n)
      write(*,*) sum, sqrt(pi)/-2.0d0

end program simpson_test_wf


text:
!----- 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: d3dpublish.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/io_dol_f/src/d3dpublish.f90 $
!
! Utility module for DOL / Delft OnLine
!
! It contains some global data and some functions

module D3DPublish
use D3DOnline

implicit none

private    !All viriables and routines are private, unless declared public

public publishGDP
public publishUtils
public setRunningFlag
public setEndTimeStep
public setEndFlag
public timeStepPrint
public printArray

logical :: first       =  .FALSE.
integer :: runningFlag =  0       !0=Flow not running (waiting)
                                  !1=Flow is iterating time steps
integer :: currentStep =  0
integer :: endTimeStep =  0
integer :: endFlag     =  0       !0 = Flow simulation has not finished yet
                                  !1 = Flow has finished
integer :: timeStepInt = -1
!=============================================================================
contains

    subroutine publishGDP(olv_handle, gdp, runid, zmodel)

        use precision
        use globaldata
        use d3d_olv_class

        implicit none

        !-----  Global data  ------!

        type(globdat), target :: gdp ! Declares the gdp and has a USE statement
        include 'fsm.i'
        include 'tri-dyn.igd'
        include 'delftonline.i'

        !-----  Parameters  ------!

        type(olvhandle) :: olv_handle
        logical         :: zmodel
        character(*)    :: runid

        !-----  Local data  ------!

        type(olv_handle_t), pointer  :: olv
        integer                      :: dimens(3)     ! 3D array to contain various D3D array dimensions
        integer                      :: dimens_k(1)   ! 1D array to contain kmax array dimension
        integer                      :: j
        integer                      :: k
        integer(pntrsize)            :: indx
        integer                      :: realDouble
        integer                      :: size_const
        integer                      :: r_offset
        integer                      :: ch_offset
        character(20)                :: name          ! character length MUST be 20 according namcon!!!
        character(120)               :: cident

        if (isNull_olv(olv_handle)) return

        !-----  Executable statements  ------!
        olv => olv_handle%fields

        !----- Set the precision for reals

        if (fp /= prec) then
            return
            !
            ! This causes the Main Development Stream version to abort on single precision mode
            ! ToDo: What DO you want?
            !
            !write (*,*) 'Delft3D Online Abort: Mixed precision between FLOW and BODSED/DPS'
            !call d3stop (1, gdp)    ! ToDo: A proper way to terminate
        endif

        if (fp == hp) then
            realDouble = DOL_DOUBLE
        else
            realDouble = DOL_REAL
        endif

        !-----  Set the run title

        call getfullversionstring_flow2d3d(cident)
        call FLOWOL_SetDescription (cident, "Delft3D-FLOW")

        !----- Size of the array dimensions
        !
        dimens(1)   = gdp%d%nub - gdp%d%nlb + 1
        dimens(2)   = gdp%d%mub - gdp%d%mlb + 1
        dimens(3)   = gdp%d%kmax
        dimens_k(1) = gdp%d%kmax

        call FLOWOL_ArrayShape ("nm"  , 2, dimens)
        call FLOWOL_ArrayShape ("nmk" , 3, dimens)
        call FLOWOL_ArrayShape ("kmax", 1, dimens_k)

        ! For some arrays dimens(3)=kmax+1; So, the third dimension varies, quantity is defined at interfaces.
        !
        dimens(3) = gdp%d%kmax+1
        call FLOWOL_ArrayShape ("nmk1", 3, dimens)
        !
        ! Do not use dimens(3) below (is it kmax or kmax+1?) but use gdp%d%kmax
        ! For safety:
        dimens(3) = -999

        !----- Description parameters D3DOnline_Publish (full description and declaration in esm_alloc_real.f90):
        ! A      internal name
        ! B      user name (description)
        ! C      units as string
        ! D      code "definedon"
        ! E   array shape, defined earlier or empty, which means a scalair
        ! F   type of the variable or array, an integer (enumeration), See file "delftonline.i"
        ! G   the scalair or the array (first array position)
        ! H   integer (enumeration). For OLV it is only DOL_OUT
        !
        ! parameter codes     A               B                        C    D   E       F            G                      H
        !
        !The following parameters, always need to be published (descriptions can be found in esm_alloc_real.f90)!
        call FLOWOL_Publish( 'nlb'  , ''                             , '' , '', ''    , DOL_INTEGER, gdp%d%nlb            , DOL_OUT)
        call FLOWOL_Publish( 'nub'  , ''                             , '' , '', ''    , DOL_INTEGER, gdp%d%nub            , DOL_OUT)
        call FLOWOL_Publish( 'mlb'  , ''                             , '' , '', ''    , DOL_INTEGER, gdp%d%mlb            , DOL_OUT)
        call FLOWOL_Publish( 'mub'  , ''                             , '' , '', ''    , DOL_INTEGER, gdp%d%mub            , DOL_OUT)
        call FLOWOL_Publish( 'nmax' , ''                             , '' , '', ''    , DOL_INTEGER, gdp%d%nmax           , DOL_OUT)
        call FLOWOL_Publish( 'mmax' , ''                             , '' , '', ''    , DOL_INTEGER, gdp%d%mmax           , DOL_OUT)
        call FLOWOL_Publish( 'kmax' , ''                             , '' , '', ''    , DOL_INTEGER, gdp%d%kmax           , DOL_OUT)
        call FLOWOL_Publish( 'thick', 'Rel. layer thickness'         , '' , '', 'kmax', realDouble , r(gdp%gdr_i_ch%thick), DOL_OUT)
        call FLOWOL_Publish( 'xcor' , 'Grid X coordinates'           , 'm', '', 'nm'  , realDouble , r(gdp%gdr_i_ch%xcor) , DOL_OUT)
        call FLOWOL_Publish( 'ycor' , 'Grid Y coordinates'           , 'm', '', 'nm'  , realDouble , r(gdp%gdr_i_ch%ycor) , DOL_OUT)
        call FLOWOL_Publish( 'gvz'  , 'Zeta-Grid cell distance (eta)', 'm', '', 'nm'  , realDouble , r(gdp%gdr_i_ch%gvz)  , DOL_OUT)
        call FLOWOL_Publish( 'guz'  , 'Zeta-Grid cell distance (ksi)', 'm', '', 'nm'  , realDouble , r(gdp%gdr_i_ch%guz)  , DOL_OUT)
        call FLOWOL_Publish( 'xz'   , 'Zeta-Grid X coord.'           , 'm', '', 'nm'  , realDouble , r(gdp%gdr_i_ch%xz)   , DOL_OUT)
        call FLOWOL_Publish( 'yz'   , 'Zeta-Grid Y coord.'           , 'm', '', 'nm'  , realDouble , r(gdp%gdr_i_ch%yz)   , DOL_OUT)
        call FLOWOL_Publish( 'kcs'  , 'KCS mask array'               , '' , '', 'nm'  , DOL_INTEGER, i(gdp%gdr_i_ch%kcs)  , DOL_OUT)
        call FLOWOL_Publish( 'kfs'  , 'KFS mask array'               , '' , '', 'nm'  , DOL_INTEGER, i(gdp%gdr_i_ch%kfs)  , DOL_OUT)
        call FLOWOL_Publish( 'kfu'  , 'KFU mask array'               , '' , '', 'nm'  , DOL_INTEGER, i(gdp%gdr_i_ch%kfu)  , DOL_OUT)
        call FLOWOL_Publish( 'kfv'  , 'KFV mask array'               , '' , '', 'nm'  , DOL_INTEGER, i(gdp%gdr_i_ch%kfv)  , DOL_OUT)
        call FLOWOL_Publish( 'alfas', 'Cell orientations'            , '' , '', 'nm'  , realDouble , r(gdp%gdr_i_ch%alfas), DOL_OUT)

        !The following variables/parameters are optional
        call FLOWOL_Publish( 's1'   , 'Water level'                  , '' , '', 'nm'  , realDouble , r(gdp%gdr_i_ch%s1)   , DOL_OUT)
        call FLOWOL_Publish( 'dp'   , 'Depth'                        , '' , '', 'nm'  , realDouble , r(gdp%gdr_i_ch%dp)   , DOL_OUT)
        call FLOWOL_Publish( 'dps'  , 'Bedlevel'                     , '' , '', 'nm'  , DOL_DOUBLE , d(gdp%gdr_i_ch%dps)  , DOL_OUT)
        call FLOWOL_Publish( 'u1'   , 'Velocity ksi component'       , '' , '', 'nmk' , realDouble , r(gdp%gdr_i_ch%u1)   , DOL_OUT)
        call FLOWOL_Publish( 'v1'   , 'Velocity eta component'       , '' , '', 'nmk' , realDouble , r(gdp%gdr_i_ch%v1)   , DOL_OUT)

        !call FLOWOL_Publish( 'rho',  'Density',              '',    '', 'nm',   realDouble,  r(rho),   DOL_INOUT)

        ! The following variables are only necessary when the z-model is used. It has to be done in this manner,
        ! otherwise things will go wrong in the client (test if zmodel exists).
        if ( zmodel ) then
            call FLOWOL_Publish('zbot'  , ''               , '', '', ''   , realDouble , gdp%gdzmodel%zbot    , DOL_OUT)
            call FLOWOL_Publish('ztop'  , ''               , '', '', ''   , realDouble , gdp%gdzmodel%ztop    , DOL_OUT)
            call FLOWOL_Publish('zmodel', ''               , '', '', ''   , DOL_LOGICAL, zmodel               , DOL_OUT)
            call FLOWOL_Publish('kfsz1' , 'KFSZ mask array', '', '', 'nmk', DOL_INTEGER, i(gdp%gdr_i_ch%kfsz1), DOL_OUT)
        endif

        ! Publish so called "transportable quantities"
        ! The array character(20), dimension(lstsci + ltur) :: namcom
        ! contains the names of these "transportable quantities" AND the turbulence quantities.
        ! Variable  lstsci  is the number of quanties and including  Sal., temp. en spiraalstroming,
        ! which has the names 'Salinity', 'Temperature', 'Secondary flow', 'Turbulent energy' and
        ! 'Energy dissipation '
        ! Variable ltur is the number of turbulence quantities.
        ! lmax = lstsci + ltur
        !
        ! The "transportable quantities" are in array r(r1) as first postion and
        ! the turbulence quantities are in array r(rtur1) as first position.
        !

        ! Character array namcon is part of the large memory pool array "ch" and indicated
        ! by namcon.
        ! Character array namcon is defined as character(len=20).
        ! So here you must calculate the right offset in the large array "ch"
        ! To pass the right character length to routine FLOWOL_Publish the
        ! local character variable name is used.
        !
        ! Publish the "transportable quantities" by calculating their relative position:
        !
        size_const = dimens(1) * dimens(2) * gdp%d%kmax   !number of gridpoints per quantity

        ! Do an allocate of the same array size as ch from esm/fsm
        ! character length is 20 for both arrays
        !

        do j = 1, gdp%d%lstsci
            ch_offset = (j-1) * 20
            r_offset  = (j-1) * size_const
            indx = gdp%gdr_i_ch%namcon+ch_offset
            do k = 1, 20
                name(k:k) = ch(indx)(k:k)
            enddo
            call FLOWOL_Publish(name, name, '', '', 'nmk', realDouble, r(gdp%gdr_i_ch%r1+r_offset), DOL_OUT)
        enddo

        ! Publish the "turbulence quantities" by calculating their relative position:
        !
        size_const = dimens(1) * dimens(2) * (gdp%d%kmax+1)     !number of gridpoints per quantity
        do j = 1, gdp%d%ltur
            ch_offset = (j-1 + gdp%d%lstsci) * 20
            r_offset  = (j-1)                * size_const
            indx = gdp%gdr_i_ch%namcon+ch_offset
            do k = 1, 20
                name(k:k) = ch(indx)(k:k)
            enddo
            call FLOWOL_Publish(name, name, '', '', 'nm', realDouble, r(gdp%gdr_i_ch%rtur1+r_offset), DOL_OUT)
        enddo

    end subroutine publishGDP
!
!---------------------------------------------------------------------------------------

    subroutine publishUtils(olv_handle)
        use d3d_olv_class

        include 'delftonline.i'
        type(olvhandle)             :: olv_handle
        type(olv_handle_t), pointer :: olv

        if (isNull_olv(olv_handle)) return

        olv => olv_handle%fields

        call FLOWOL_Publish( 'runningFlag', '', '', '', '', DOL_INTEGER, olv%runningFlag,  DOL_OUT)
        call FLOWOL_Publish( 'currentStep', '', '', '', '', DOL_INTEGER, olv%currentStep,  DOL_OUT)
        call FLOWOL_Publish( 'endTimeStep', '', '', '', '', DOL_INTEGER, olv%endTimeStep,  DOL_OUT)
        call FLOWOL_Publish( 'endFlag'    , '', '', '', '', DOL_INTEGER, olv%endFlag    ,  DOL_OUT)
        call FLOWOL_Publish( 'timeStepInt', '', '', '', '', DOL_INTEGER, olv%timeStepInt,  DOL_INOUT)

    end subroutine publishUtils
!---------------------------------------------------------------------------------------

    subroutine setRunningFlag(olv_handle, value, nst)
        use d3d_olv_class

        integer :: value
        integer :: nst

        type(olvhandle) :: olv_handle
        type(olv_handle_t), pointer :: olv

        if (isNull_olv(olv_handle)) return

        olv => olv_handle%fields

        olv%runningFlag = value
        olv%currentStep = nst

    end subroutine setRunningFlag
!---------------------------------------------------------------------------------------

    subroutine setEndTimeStep(olv_handle, itstop)
        use d3d_olv_class

        integer :: itstop

        type(olvhandle) :: olv_handle
        type(olv_handle_t), pointer :: olv

        if (isNull_olv(olv_handle)) return

        olv => olv_handle%fields

        olv%endTimeStep = itstop

    end subroutine setEndTimeStep
!---------------------------------------------------------------------------------------

    subroutine setEndFlag(olv_handle, value)
        use d3d_olv_class
        integer :: value

        type(olvhandle) :: olv_handle
        type(olv_handle_t), pointer :: olv

        if (isNull_olv(olv_handle)) return

        olv => olv_handle%fields

        olv%endFlag = value

    end subroutine

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

    subroutine timeStepPrint(olv_handle)
        use d3d_olv_class

        type(olvhandle) :: olv_handle
        type(olv_handle_t), pointer :: olv

        olv => olv_handle%fields

        write(*,*) 'Current Timestep =' , olv%timeStepInt + 1

    end subroutine

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

    subroutine printArray(nlb, nub, mlb, mub, value, name)

        integer :: nlb, nub, mlb, mub
        integer, dimension( nlb:nub, mlb:mub ) :: value
        character(len=*) :: name

        integer :: m, n

        write(*,*) name
        do m = mlb, mub
            write(*, *) 'm=' , m
            do n = nlb, nub
                write(*,'(i4)') value(n,m)
            enddo
        enddo

    end subroutine printArray
!---------------------------------------------------------------------------------------

end module D3DPublish

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
module main

   use iso_fortran_env, only: real64
   implicit none
   private

  type, public :: test_t
      real(real64) :: a, b
  contains
      procedure, private :: init_sub
      generic :: init => init_sub
  end type test_t

   type :: node
      private
      type(node), pointer :: next => null()
      class(*), allocatable :: item

      contains
         final :: node_finalizer

   end type node

  contains

   subroutine init_sub(this, a, b)

      class( test_t ) :: this
      real(real64),intent(in) :: a, b

      this%a = a
      this%b = b

   end subroutine init_sub

   subroutine node_finalizer(a)
      type(node) :: a

   end subroutine node_finalizer

end module main

max_stars_repo_name: quanzhuo/vscode-fortran-support

max_issues_repo_name: quanzhuo/vscode-fortran-support

max_forks_repo_name: quanzhuo/vscode-fortran-support