text:
	SUBROUTINE svdcmp_sp(a,w,v)
	USE nrtype
 USE nrutil, ONLY : assert_eq,nrerror,outerprod
	USE Debug_module
 USE   nr, ONLY : pythag
	IMPLICIT NONE
	REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a
	REAL(SP), DIMENSION(:), INTENT(OUT) :: w
	REAL(SP), DIMENSION(:,:), INTENT(OUT) :: v
	INTEGER(I4B) :: i,its,j,k,l,m,n,nm
	REAL(SP) :: anorm,c,f,g,h,s,scale,x,y,z
	REAL(SP), DIMENSION(size(a,1)) :: tempm
	REAL(SP), DIMENSION(size(a,2)) :: rv1,tempn
	m=size(a,1)
	n=assert_eq(size(a,2),size(v,1),size(v,2),size(w),'svdcmp_sp')
	g=0.0
	scale=0.0
	do i=1,n
		l=i+1
		rv1(i)=scale*g
		g=0.0
		scale=0.0
		if (i <= m) then
			scale=sum(abs(a(i:m,i)))
			if (scale /= 0.0) then
				a(i:m,i)=a(i:m,i)/scale
				s=dot_product(a(i:m,i),a(i:m,i))
				f=a(i,i)
				g=-sign(sqrt(s),f)
				h=f*g-s
				a(i,i)=f-g
				tempn(l:n)=matmul(a(i:m,i),a(i:m,l:n))/h
				a(i:m,l:n)=a(i:m,l:n)+outerprod(a(i:m,i),tempn(l:n))
				a(i:m,i)=scale*a(i:m,i)
			end if
		end if
		w(i)=scale*g
		g=0.0
		scale=0.0
		if ((i <= m) .and. (i /= n)) then
			scale=sum(abs(a(i,l:n)))
			if (scale /= 0.0) then
				a(i,l:n)=a(i,l:n)/scale
				s=dot_product(a(i,l:n),a(i,l:n))
				f=a(i,l)
				g=-sign(sqrt(s),f)
				h=f*g-s
				a(i,l)=f-g
				rv1(l:n)=a(i,l:n)/h
				tempm(l:m)=matmul(a(l:m,l:n),a(i,l:n))
				a(l:m,l:n)=a(l:m,l:n)+outerprod(tempm(l:m),rv1(l:n))
				a(i,l:n)=scale*a(i,l:n)
			end if
		end if
	end do
	anorm=maxval(abs(w)+abs(rv1))
	do i=n,1,-1
		if (i < n) then
			if (g /= 0.0) then
				v(l:n,i)=(a(i,l:n)/a(i,l))/g
				tempn(l:n)=matmul(a(i,l:n),v(l:n,l:n))
				v(l:n,l:n)=v(l:n,l:n)+outerprod(v(l:n,i),tempn(l:n))
			end if
			v(i,l:n)=0.0
			v(l:n,i)=0.0
		end if
		v(i,i)=1.0
		g=rv1(i)
		l=i
	end do
	do i=min(m,n),1,-1
		l=i+1
		g=w(i)
		a(i,l:n)=0.0
		if (g /= 0.0) then
			g=1.0_sp/g
			tempn(l:n)=(matmul(a(l:m,i),a(l:m,l:n))/a(i,i))*g
			a(i:m,l:n)=a(i:m,l:n)+outerprod(a(i:m,i),tempn(l:n))
			a(i:m,i)=a(i:m,i)*g
		else
			a(i:m,i)=0.0
		end if
		a(i,i)=a(i,i)+1.0_sp
	end do
	do k=n,1,-1
		do its=1,30
			do l=k,1,-1
				nm=l-1
				if ((abs(rv1(l))+anorm) == anorm) exit
				if ((abs(w(nm))+anorm) == anorm) then
					c=0.0
					s=1.0
					do i=l,k
						f=s*rv1(i)
						rv1(i)=c*rv1(i)
						if ((abs(f)+anorm) == anorm) exit
						g=w(i)
						h=pythag(f,g)
						w(i)=h
						h=1.0_sp/h
						c= (g*h)
						s=-(f*h)
						tempm(1:m)=a(1:m,nm)
						a(1:m,nm)=a(1:m,nm)*c+a(1:m,i)*s
						a(1:m,i)=-tempm(1:m)*s+a(1:m,i)*c
					end do
					exit
				end if
			end do
			z=w(k)
			if (l == k) then
				if (z < 0.0) then
					w(k)=-z
					v(1:n,k)=-v(1:n,k)
				end if
				exit
			end if
			if (its == 30) then
 Debug_errorflags(flag_svdcmp)=1
   return
   end if
			x=w(l)
			nm=k-1
			y=w(nm)
			g=rv1(nm)
			h=rv1(k)
			f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0_sp*h*y)
			g=pythag(f,1.0_sp)
			f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x
			c=1.0
			s=1.0
			do j=l,nm
				i=j+1
				g=rv1(i)
				y=w(i)
				h=s*g
				g=c*g
				z=pythag(f,h)
				rv1(j)=z
				c=f/z
				s=h/z
				f= (x*c)+(g*s)
				g=-(x*s)+(g*c)
				h=y*s
				y=y*c
				tempn(1:n)=v(1:n,j)
				v(1:n,j)=v(1:n,j)*c+v(1:n,i)*s
				v(1:n,i)=-tempn(1:n)*s+v(1:n,i)*c
				z=pythag(f,h)
				w(j)=z
				if (z /= 0.0) then
					z=1.0_sp/z
					c=f*z
					s=h*z
				end if
				f= (c*g)+(s*y)
				x=-(s*g)+(c*y)
				tempm(1:m)=a(1:m,j)
				a(1:m,j)=a(1:m,j)*c+a(1:m,i)*s
				a(1:m,i)=-tempm(1:m)*s+a(1:m,i)*c
			end do
			rv1(l)=0.0
			rv1(k)=f
			w(k)=x
		end do
	end do
	END SUBROUTINE svdcmp_sp

	SUBROUTINE svdcmp_dp(a,w,v)
	USE nrtype
 USE nrutil, ONLY : assert_eq,nrerror,outerprod
	USE Debug_module
 USE   nr, ONLY : pythag
	IMPLICIT NONE
	REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: a
	REAL(DP), DIMENSION(:), INTENT(OUT) :: w
	REAL(DP), DIMENSION(:,:), INTENT(OUT) :: v
	INTEGER(I4B) :: i,its,j,k,l,m,n,nm
	REAL(DP) :: anorm,c,f,g,h,s,scale,x,y,z
	REAL(DP), DIMENSION(size(a,1)) :: tempm
	REAL(DP), DIMENSION(size(a,2)) :: rv1,tempn
	m=size(a,1)
	n=assert_eq(size(a,2),size(v,1),size(v,2),size(w),'svdcmp_dp')
	g=0.0
	scale=0.0
	do i=1,n
		l=i+1
		rv1(i)=scale*g
		g=0.0
		scale=0.0
		if (i <= m) then
			scale=sum(abs(a(i:m,i)))
			if (scale /= 0.0) then
				a(i:m,i)=a(i:m,i)/scale
				s=dot_product(a(i:m,i),a(i:m,i))
				f=a(i,i)
				g=-sign(sqrt(s),f)
				h=f*g-s
				a(i,i)=f-g
				tempn(l:n)=matmul(a(i:m,i),a(i:m,l:n))/h
				a(i:m,l:n)=a(i:m,l:n)+outerprod(a(i:m,i),tempn(l:n))
				a(i:m,i)=scale*a(i:m,i)
			end if
		end if
		w(i)=scale*g
		g=0.0
		scale=0.0
		if ((i <= m) .and. (i /= n)) then
			scale=sum(abs(a(i,l:n)))
			if (scale /= 0.0) then
				a(i,l:n)=a(i,l:n)/scale
				s=dot_product(a(i,l:n),a(i,l:n))
				f=a(i,l)
				g=-sign(sqrt(s),f)
				h=f*g-s
				a(i,l)=f-g
				rv1(l:n)=a(i,l:n)/h
				tempm(l:m)=matmul(a(l:m,l:n),a(i,l:n))
				a(l:m,l:n)=a(l:m,l:n)+outerprod(tempm(l:m),rv1(l:n))
				a(i,l:n)=scale*a(i,l:n)
			end if
		end if
	end do
	anorm=maxval(abs(w)+abs(rv1))
	do i=n,1,-1
		if (i < n) then
			if (g /= 0.0) then
				v(l:n,i)=(a(i,l:n)/a(i,l))/g
				tempn(l:n)=matmul(a(i,l:n),v(l:n,l:n))
				v(l:n,l:n)=v(l:n,l:n)+outerprod(v(l:n,i),tempn(l:n))
			end if
			v(i,l:n)=0.0
			v(l:n,i)=0.0
		end if
		v(i,i)=1.0
		g=rv1(i)
		l=i
	end do
	do i=min(m,n),1,-1
		l=i+1
		g=w(i)
		a(i,l:n)=0.0
		if (g /= 0.0) then
			g=1.0_dp/g
			tempn(l:n)=(matmul(a(l:m,i),a(l:m,l:n))/a(i,i))*g
			a(i:m,l:n)=a(i:m,l:n)+outerprod(a(i:m,i),tempn(l:n))
			a(i:m,i)=a(i:m,i)*g
		else
			a(i:m,i)=0.0
		end if
		a(i,i)=a(i,i)+1.0_dp
	end do
	do k=n,1,-1
		do its=1,30
			do l=k,1,-1
				nm=l-1
				if ((abs(rv1(l))+anorm) == anorm) exit
				if ((abs(w(nm))+anorm) == anorm) then
					c=0.0
					s=1.0
					do i=l,k
						f=s*rv1(i)
						rv1(i)=c*rv1(i)
						if ((abs(f)+anorm) == anorm) exit
						g=w(i)
						h=pythag(f,g)
						w(i)=h
						h=1.0_dp/h
						c= (g*h)
						s=-(f*h)
						tempm(1:m)=a(1:m,nm)
						a(1:m,nm)=a(1:m,nm)*c+a(1:m,i)*s
						a(1:m,i)=-tempm(1:m)*s+a(1:m,i)*c
					end do
					exit
				end if
			end do
			z=w(k)
			if (l == k) then
				if (z < 0.0) then
					w(k)=-z
					v(1:n,k)=-v(1:n,k)
				end if
				exit
			end if
			if (its == 30) then
 Debug_errorflags(flag_svdcmp)=1
   return
   end if
			x=w(l)
			nm=k-1
			y=w(nm)
			g=rv1(nm)
			h=rv1(k)
			f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0_dp*h*y)
			g=pythag(f,1.0_dp)
			f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x
			c=1.0
			s=1.0
			do j=l,nm
				i=j+1
				g=rv1(i)
				y=w(i)
				h=s*g
				g=c*g
				z=pythag(f,h)
				rv1(j)=z
				c=f/z
				s=h/z
				f= (x*c)+(g*s)
				g=-(x*s)+(g*c)
				h=y*s
				y=y*c
				tempn(1:n)=v(1:n,j)
				v(1:n,j)=v(1:n,j)*c+v(1:n,i)*s
				v(1:n,i)=-tempn(1:n)*s+v(1:n,i)*c
				z=pythag(f,h)
				w(j)=z
				if (z /= 0.0) then
					z=1.0_dp/z
					c=f*z
					s=h*z
				end if
				f= (c*g)+(s*y)
				x=-(s*g)+(c*y)
				tempm(1:m)=a(1:m,j)
				a(1:m,j)=a(1:m,j)*c+a(1:m,i)*s
				a(1:m,i)=-tempm(1:m)*s+a(1:m,i)*c
			end do
			rv1(l)=0.0
			rv1(k)=f
			w(k)=x
		end do
	end do
	END SUBROUTINE svdcmp_dp

max_stars_repo_name: aasensio/DeepLearning

max_issues_repo_name: aasensio/DeepLearning

max_forks_repo_name: aasensio/DeepLearning


text:
      SUBROUTINE DCHKQ3( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL,
     $                   THRESH, A, COPYA, S, COPYS, TAU, WORK, IWORK,
     $                   NOUT )
*
*  -- LAPACK test routine (version 3.1.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     January 2007
*
*     .. Scalar Arguments ..
      INTEGER            NM, NN, NNB, NOUT
      DOUBLE PRECISION   THRESH
*     ..
*     .. Array Arguments ..
      LOGICAL            DOTYPE( * )
      INTEGER            IWORK( * ), MVAL( * ), NBVAL( * ), NVAL( * ),
     $                   NXVAL( * )
      DOUBLE PRECISION   A( * ), COPYA( * ), COPYS( * ), S( * ),
     $                   TAU( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  DCHKQ3 tests DGEQP3.
*
*  Arguments
*  =========
*
*  DOTYPE  (input) LOGICAL array, dimension (NTYPES)
*          The matrix types to be used for testing.  Matrices of type j
*          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
*          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
*
*  NM      (input) INTEGER
*          The number of values of M contained in the vector MVAL.
*
*  MVAL    (input) INTEGER array, dimension (NM)
*          The values of the matrix row dimension M.
*
*  NN      (input) INTEGER
*          The number of values of N contained in the vector NVAL.
*
*  NVAL    (input) INTEGER array, dimension (NN)
*          The values of the matrix column dimension N.
*
*  NNB     (input) INTEGER
*          The number of values of NB and NX contained in the
*          vectors NBVAL and NXVAL.  The blocking parameters are used
*          in pairs (NB,NX).
*
*  NBVAL   (input) INTEGER array, dimension (NNB)
*          The values of the blocksize NB.
*
*  NXVAL   (input) INTEGER array, dimension (NNB)
*          The values of the crossover point NX.
*
*  THRESH  (input) DOUBLE PRECISION
*          The threshold value for the test ratios.  A result is
*          included in the output file if RESULT >= THRESH.  To have
*          every test ratio printed, use THRESH = 0.
*
*  A       (workspace) DOUBLE PRECISION array, dimension (MMAX*NMAX)
*          where MMAX is the maximum value of M in MVAL and NMAX is the
*          maximum value of N in NVAL.
*
*  COPYA   (workspace) DOUBLE PRECISION array, dimension (MMAX*NMAX)
*
*  S       (workspace) DOUBLE PRECISION array, dimension
*                      (min(MMAX,NMAX))
*
*  COPYS   (workspace) DOUBLE PRECISION array, dimension
*                      (min(MMAX,NMAX))
*
*  TAU     (workspace) DOUBLE PRECISION array, dimension (MMAX)
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension
*                      (MMAX*NMAX + 4*NMAX + MMAX)
*
*  IWORK   (workspace) INTEGER array, dimension (2*NMAX)
*
*  NOUT    (input) INTEGER
*          The unit number for output.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            NTYPES
      PARAMETER          ( NTYPES = 6 )
      INTEGER            NTESTS
      PARAMETER          ( NTESTS = 3 )
      DOUBLE PRECISION   ONE, ZERO
      PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0 )
*     ..
*     .. Local Scalars ..
      CHARACTER*3        PATH
      INTEGER            I, IHIGH, ILOW, IM, IMODE, IN, INB, INFO,
     $                   ISTEP, K, LDA, LW, LWORK, M, MNMIN, MODE, N,
     $                   NB, NERRS, NFAIL, NRUN, NX
      DOUBLE PRECISION   EPS
*     ..
*     .. Local Arrays ..
      INTEGER            ISEED( 4 ), ISEEDY( 4 )
      DOUBLE PRECISION   RESULT( NTESTS )
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH, DQPT01, DQRT11, DQRT12
      EXTERNAL           DLAMCH, DQPT01, DQRT11, DQRT12
*     ..
*     .. External Subroutines ..
      EXTERNAL           ALAHD, ALASUM, DGEQP3, DLACPY, DLAORD, DLASET,
     $                   DLATMS, ICOPY, XLAENV
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN
*     ..
*     .. Scalars in Common ..
      LOGICAL            LERR, OK
      CHARACTER*32       SRNAMT
      INTEGER            INFOT, IOUNIT
*     ..
*     .. Common blocks ..
      COMMON             / INFOC / INFOT, IOUNIT, OK, LERR
      COMMON             / SRNAMC / SRNAMT
*     ..
*     .. Data statements ..
      DATA               ISEEDY / 1988, 1989, 1990, 1991 /
*     ..
*     .. Executable Statements ..
*
*     Initialize constants and the random number seed.
*
      PATH( 1: 1 ) = 'Double precision'
      PATH( 2: 3 ) = 'Q3'
      NRUN = 0
      NFAIL = 0
      NERRS = 0
      DO 10 I = 1, 4
         ISEED( I ) = ISEEDY( I )
   10 CONTINUE
      EPS = DLAMCH( 'Epsilon' )
      INFOT = 0
*
      DO 90 IM = 1, NM
*
*        Do for each value of M in MVAL.
*
         M = MVAL( IM )
         LDA = MAX( 1, M )
*
         DO 80 IN = 1, NN
*
*           Do for each value of N in NVAL.
*
            N = NVAL( IN )
            MNMIN = MIN( M, N )
            LWORK = MAX( 1, M*MAX( M, N )+4*MNMIN+MAX( M, N ),
     $                   M*N + 2*MNMIN + 4*N )
*
            DO 70 IMODE = 1, NTYPES
               IF( .NOT.DOTYPE( IMODE ) )
     $            GO TO 70
*
*              Do for each type of matrix
*                 1:  zero matrix
*                 2:  one small singular value
*                 3:  geometric distribution of singular values
*                 4:  first n/2 columns fixed
*                 5:  last n/2 columns fixed
*                 6:  every second column fixed
*
               MODE = IMODE
               IF( IMODE.GT.3 )
     $            MODE = 1
*
*              Generate test matrix of size m by n using
*              singular value distribution indicated by `mode'.
*
               DO 20 I = 1, N
                  IWORK( I ) = 0
   20          CONTINUE
               IF( IMODE.EQ.1 ) THEN
                  CALL DLASET( 'Full', M, N, ZERO, ZERO, COPYA, LDA )
                  DO 30 I = 1, MNMIN
                     COPYS( I ) = ZERO
   30             CONTINUE
               ELSE
                  CALL DLATMS( M, N, 'Uniform', ISEED, 'Nonsymm', COPYS,
     $                         MODE, ONE / EPS, ONE, M, N, 'No packing',
     $                         COPYA, LDA, WORK, INFO )
                  IF( IMODE.GE.4 ) THEN
                     IF( IMODE.EQ.4 ) THEN
                        ILOW = 1
                        ISTEP = 1
                        IHIGH = MAX( 1, N / 2 )
                     ELSE IF( IMODE.EQ.5 ) THEN
                        ILOW = MAX( 1, N / 2 )
                        ISTEP = 1
                        IHIGH = N
                     ELSE IF( IMODE.EQ.6 ) THEN
                        ILOW = 1
                        ISTEP = 2
                        IHIGH = N
                     END IF
                     DO 40 I = ILOW, IHIGH, ISTEP
                        IWORK( I ) = 1
   40                CONTINUE
                  END IF
                  CALL DLAORD( 'Decreasing', MNMIN, COPYS, 1 )
               END IF
*
               DO 60 INB = 1, NNB
*
*                 Do for each pair of values (NB,NX) in NBVAL and NXVAL.
*
                  NB = NBVAL( INB )
                  CALL XLAENV( 1, NB )
                  NX = NXVAL( INB )
                  CALL XLAENV( 3, NX )
*
*                 Get a working copy of COPYA into A and a copy of
*                 vector IWORK.
*
                  CALL DLACPY( 'All', M, N, COPYA, LDA, A, LDA )
                  CALL ICOPY( N, IWORK( 1 ), 1, IWORK( N+1 ), 1 )
*
*                 Compute the QR factorization with pivoting of A
*
                  LW = MAX( 1, 2*N+NB*( N+1 ) )
*
*                 Compute the QP3 factorization of A
*
                  SRNAMT = 'DGEQP3'
                  CALL DGEQP3( M, N, A, LDA, IWORK( N+1 ), TAU, WORK,
     $                         LW, INFO )
*
*                 Compute norm(svd(a) - svd(r))
*
                  RESULT( 1 ) = DQRT12( M, N, A, LDA, COPYS, WORK,
     $                          LWORK )
*
*                 Compute norm( A*P - Q*R )
*
                  RESULT( 2 ) = DQPT01( M, N, MNMIN, COPYA, A, LDA, TAU,
     $                          IWORK( N+1 ), WORK, LWORK )
*
*                 Compute Q'*Q
*
                  RESULT( 3 ) = DQRT11( M, MNMIN, A, LDA, TAU, WORK,
     $                          LWORK )
*
*                 Print information about the tests that did not pass
*                 the threshold.
*
                  DO 50 K = 1, NTESTS
                     IF( RESULT( K ).GE.THRESH ) THEN
                        IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                     CALL ALAHD( NOUT, PATH )
                        WRITE( NOUT, FMT = 9999 )'DGEQP3', M, N, NB,
     $                     IMODE, K, RESULT( K )
                        NFAIL = NFAIL + 1
                     END IF
   50             CONTINUE
                  NRUN = NRUN + NTESTS
*
   60          CONTINUE
   70       CONTINUE
   80    CONTINUE
   90 CONTINUE
*
*     Print a summary of the results.
*
      CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
*
 9999 FORMAT( 1X, A, ' M =', I5, ', N =', I5, ', NB =', I4, ', type ',
     $      I2, ', test ', I2, ', ratio =', G12.5 )
*
*     End of DCHKQ3
*
      END

ext: f


text:
! From: "James Van Buskirk" <torsop@ix.netcom.com>
! Newsgroups: comp.lang.fortran
! Subject: Re: ? Strassen algorithm for multiplying matrix
! Date: Thu, 14 Sep 2000 06:27:20 -0600

! Latest revision - 16 September 2000

module mykinds
   implicit none
   integer, parameter :: wp = selected_real_kind(12, 60)
   integer            :: threshold
end module mykinds

module matrix_multiply
   use mykinds
   implicit none

contains

recursive subroutine strassen(a, b, c, N)
   integer, intent(in)   :: N
   real(wp), intent(in)  :: a(N,N), b(N,N)
   real(wp), intent(out) :: c(N,N)

   real(wp), dimension(N/2,N/2) :: a11,a21,a12,a22,b11,b21,b12,b22
   real(wp), dimension(N/2,N/2) :: q1,q2,q3,q4,q5,q6,q7

   if(iand(N,1) /= 0 .OR. N < threshold) then
      c = matmul(a,b)
   else
      a11 = a(:N/2,:N/2)
      a21 = a(N/2+1:,:N/2)
      a12 = a(:N/2,N/2+1:)
      a22 = a(N/2+1:,N/2+1:)
      b11 = b(:N/2,:N/2)
      b21 = b(N/2+1:,:N/2)
      b12 = b(:N/2,N/2+1:)
      b22 = b(N/2+1:,N/2+1:)
      call strassen(a11+a22, b11+b22, q1, N/2)
      call strassen(a21+a22, b11, q2, N/2)
      call strassen(a11, b12-b22, q3, N/2)
      call strassen(a22, -b11+b21, q4, N/2)
      call strassen(a11+a12, b22, q5, N/2)
      call strassen(-a11+a21, b11+b12, q6, N/2)
      call strassen(a12-a22, b21+b22, q7, N/2)
      c(:N/2,:N/2) = q1+q4-q5+q7
      c(N/2+1:,:N/2) = q2+q4
      c(:N/2,N/2+1:) = q3+q5
      c(N/2+1:,N/2+1:) = q1+q3-q2+q6
   end if

   return
end subroutine strassen

recursive subroutine simple(a, b, c, N)
   use mykinds
   implicit none
   integer, intent(in)   :: N
   real(wp), intent(in)  :: a(N,N), b(N,N)
   real(wp), intent(out) :: c(N,N)
   real(wp), dimension(N/2,N/2) :: a11,a21,a12,a22,b11,b21,b12,b22,q1,q2

   if(iand(N,1) /= 0 .OR. N < threshold) then
      c = matmul(a,b)
   else
      a11 = a(:N/2,:N/2)
      a21 = a(N/2+1:,:N/2)
      a12 = a(:N/2,N/2+1:)
      a22 = a(N/2+1:,N/2+1:)
      b11 = b(:N/2,:N/2)
      b21 = b(N/2+1:,:N/2)
      b12 = b(:N/2,N/2+1:)
      b22 = b(N/2+1:,N/2+1:)
      call simple(a11, b11, q1, N/2)
      call simple(a12, b21, q2, N/2)
      c(:N/2,:N/2) = q1+q2
      call simple(a21, b11, q1, N/2)
      call simple(a22, b21, q2, N/2)
      c(N/2+1:,:N/2) = q1+q2
      call simple(a11, b12, q1, N/2)
      call simple(a12, b22, q2, N/2)
      c(:N/2,N/2+1:) = q1+q2
      call simple(a21, b12, q1, N/2)
      call simple(a22, b22, q2, N/2)
      c(N/2+1:,N/2+1:) = q1+q2
   end if

   return
end subroutine simple

end module matrix_multiply

program strassen_test
   use mykinds
   use matrix_multiply
   implicit none
   real(wp), allocatable :: a(:,:), b(:,:), c(:,:), d(:,:), e(:,:)
   integer  :: N
   real     :: t0, t1

   write(*, '(a)', advance='no') ' Enter the order of the matrices:> '
   read(*,*) N
   write(*, '(a)', advance='no') ' Enter the threshold:> '
   read(*,*) threshold
   allocate(a(N,N), b(N,N), c(N,N), d(N,N), e(N,N))
   call random_seed()
   call random_number(a)
   call random_number(b)
   call cpu_time(t0)
   c = matmul(a,b)
   call cpu_time(t1)
   write(*,'(a, f10.3)') ' Time for matmul        = ', t1-t0
   call cpu_time(t0)
   call strassen(a, b, d, N)
   call cpu_time(t1)
   write(*,'(a, f10.3)') ' Time for strassen      = ', t1-t0
   write(*,'(a, g12.4)') ' RMS error for strassen = ', sqrt(sum((c-d)**2))/N
   call cpu_time(t0)
   call simple(a,b,e,N)
   call cpu_time(t1)
   write(*,'(a, f10.3)') ' Time for simple        = ', t1-t0
   write(*,'(a, g12.4)') ' RMS error for simple   = ', sqrt(sum((c-e)**2))/N
   stop
end program strassen_test

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%
!%  Description:
!%      +   Run the Monte Carlo sampler of the ParaMonte library given the input log-target density function `getLogFunc()`.
!%  Output:
!%      +   The simulation output files.
!%  Author:
!%      +   Computational Data Science Lab, Monday 9:03 AM, May 16 2016, ICES, The University of Texas at Austin
!%  Visit:
!%      +   https://www.cdslab.org/paramonte
!%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

!   Keep in mind that Fortran is case-insensitive, except for character values and string values.
!   So, feel free to use upper-case or lower-case for the Fortran syntax and entity names.
!   The ParaMonte library uses camelCase convention for variable naming.

#if IS_COMPATIBLE_COMPILER

    program main

        ! This is the Object-Oriented-Programming (OOP) style interface to the ParaMonte routines.
        ! This is more flexible but less portable, as its compilation requires the same compiler
        ! brand and version with which the ParaMonte library has been built.

        use paramonte, only: IK
        use LogFunc_mod, only: getLogFunc, NDIM

        implicit none

        ! Rune the ParaDRAM Adaptive MCMC sampler.

        block
            use paramonte, only: ParaDRAM
            type(ParaDRAM) :: sim
            call sim%runSampler ( ndim = NDIM &
                                , getLogFunc = getLogFunc &
                                , inputFile = "./paramonte.in" &    ! this is optional argument
                                ! You can also specify simulation specifications as input arguments, like
                                ! the following. This is possible only from the OOP interface to ParaDRAM.
                                , proposalAdaptationCountGreedy = 0_IK &    ! this is optional argument
                                , description = "an example run" &  ! this is optional argument
                                ! More optional arguments can appear here.
                                ! See the ParaDRAM routine's list of input arguments.
                                )
            if (sim%Err%occurred) error stop
        end block

        ! Rune the ParaNest Nested sampler.

        block
            use paramonte, only: ParaNest
            type(ParaNest) :: sim
            call sim%runSampler ( ndim = NDIM &
                                , getLogFunc = getLogFunc &
                                , inputFile = "./paramonte.in" &    ! this is optional argument
                                ! You can also specify simulation specifications as input arguments, like
                                ! the following. This is possible only from the OOP interface to ParaNest.
                                , description = "an example run" &  ! this is optional argument
                                ! More optional arguments can appear here.
                                ! See the ParaNest routine's list of input arguments.
                                )
            if (sim%Err%occurred) error stop
        end block

    end program main

#else

    ! This is the default simple procedural interface to the ParaMonte routines.
    ! The first two arguments `ndim` and `getLogFunc` to the sampler routines are mandatory.
    ! The last argument `inputFile` is optional.

    program main

        use paramonte, only: runParaDRAM, runParaNest, IK
        use LogFunc_mod, only: getLogFunc, NDIM

        implicit none

        integer(IK) :: ierr

        ierr = runParaDRAM  ( NDIM &
                            , getLogFunc &
                            , "./paramonte.in" & ! this is optional argument
                            )
        if (ierr /= 0_IK) error stop

        ierr = runParaNest  ( NDIM &
                            , getLogFunc &
                            , "./paramonte.in" & ! this is optional argument
                            )
        if (ierr /= 0_IK) error stop

    end program main

#endif

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

max_stars_repo_name: cdslaborg/paramonte-api-kernel

max_issues_repo_name: cdslaborg/paramonte-api-kernel

max_forks_repo_name: cdslaborg/paramonte-api-kernel


text:
      subroutine stmget(mxint)
c******************************************************************
c                                                                 *
c  This subroutine is called from SR CONTIN to read in storm data *
c  on a storm-by-storm basis. It calls SR SUMRNF.                 *
c                                                                 *
c******************************************************************
c******************************************************************
c                                                                 *
c  Arguments                                                      *
c    ibrkpt - flag for breakpoint rainfall inputs 0 - no brkpts   *
c    mxint - maximum rainfall intensity based on average          *
c            intensity value and Ip value                         *
c                                                                 *
c******************************************************************
c                                                                 *
c  Parameters                                                     *
c      mxplan : maximum number of over land flow elements         *
c      mxtlsq : maximum number of tillage sequences               *
c      mxtill : maximum number of tillage operations per tillage  *
c                                                                 *
c******************************************************************
c
      include 'pmxelm.inc'
      include 'pmxhil.inc'
      include 'pmxpln.inc'
      include 'pmxtls.inc'
      include 'pmxtil.inc'
c
c******************************************************************
c
c  Common Blocks                                                  *
c                                                                 *
c******************************************************************
c                                                                 *
      include 'cclig.inc'
c        read: iclig
c
      include 'cclim.inc'
c******************************************************************
c
c  clim variables updated                                         *
c     tave,tmnavg,tmxavg                                          *
c                                                                 *
c******************************************************************
c
      include 'cdiss1.inc'
      include 'cdiss3.inc'
c
      include 'chydrol.inc'
c
c******************************************************************
c                                                                 *
c  hydrol variables updated                                       *
c    stmdur,avrint                                                *
c                                                                 *
c******************************************************************
c
      include 'cstmflg.inc'
c
c******************************************************************
c                                                                 *
c   stmflg variables updated                                      *
c      norain                                                     *
c                                                                 *
c******************************************************************
c
      include 'cupdate.inc'
c
      include 'cstmgt.inc'
c
      include 'cstruc.inc'
c
c******************************************************************
c                                                                 *
c   local variables                                               *
c   mxint                                                         *
c                                                                 *
c******************************************************************
c
c     real mxint,totrn
      real mxint
      integer xday(10), xmon(10), xyear(10), nbrkpt, i, jkjk, jkl, kjl
      real xrain(10), xstmdr(10), xtimep(10), xip(10), xtmax(10),
     1    xtmin(10), xradly(10), xvwind(10), xwind(10), xtdpt(10)
c    1    xnint(10)
c
c XXX NOTE - use of global save does not follow WEPP coding
c            conventions.  Only local variables whose values need
c            to be retained should be saved.   dcf  5/20/94
      save
c      data totrn /0.0/
c
c  Variables needed to "lump" disk reads.
c
      mxint = 0.
      stmstr = 0.
c
c     Added if-else to prevent numeric overflow of integer variable
c     IDTOT in very long simulations (model bombed in year 90). dcf 2/94
c
      if(jyear.le.1 .and. sdate.le.1)then
        idtot = 0
      elseif(idtot.lt.6)then
        idtot = idtot + 1
      endif
c
      do 10 kjl =1,nplane
        norain(kjl) = 1
 10   continue
c
      if (ibrkpt.eq.0) then
c       read(13,*) day,mon,year,rain,stmdur,timep,ip,
c       1             tmax,tmin,radly,vwind,wind,tdpt
c
c       Try to speed up disk I/O by doing it in larger chunks.
        nrec = nrec + 1
        if (nrec.ge.11) then
          nrec = 1
          read (13,*,end=20) (xday(i),xmon(i),xyear(i),xrain(i),
     1        xstmdr(i),xtimep(i),xip(i),xtmax(i),xtmin(i),xradly(i),
     1        xvwind(i),xwind(i),xtdpt(i),i = 1,10)
        end if
c
   20   continue
c
c       Most of the variables below are never used in STMGET (except RAIN).
c       They are just stuffed into common.
c
        day = xday(nrec)
        mon = xmon(nrec)
        year = xyear(nrec)
c
c       If a leap year check to see if 29 days in Feb
c
        if (mod(year,4).eq.0) then
          if (mon.eq.3.and.day.eq.1.and.sdate.lt.61) then
            write (6,1000)
            write (6,1100) year
          end if
        end if
c
        do 30 jkjk=1,nplane
          rain(jkjk) = xrain(nrec)
 30     continue
        prcp = xrain(nrec)
        stmdur = xstmdr(nrec)
c
c       CLIGEN versions 2.3 to 3.1 correction factor for storm duration
c       from Risse  7/1/94  dcf
c
c       STILL NEED THIS IF USER IS RUNNING OLD (pre-CLIGEN 4) FILES
        if (imodel.eq.1 .and. iclig.eq.2) stmdur = stmdur * 2.06
c
        if (stmdur .gt. 23.999) stmdur = 23.999
c
        timep = xtimep(nrec)
        ip = xip(nrec)
c
        if (imodel.eq.1 .and. iclig.eq.2)then
c
c         correction factor for CLIGEN versions 2.3 to 3.1
c         From Risse and Nearing.   7/1/94.   dcf
c
          ip=ip*1.44
c
        elseif(iclig.eq.1)then
c
c         correction factor for CLIGEN version 4.0+ to
c         account for effects of steady-state erosion equation
c         assumption.   mn/dcf  3/95
c
          ip=ip*0.70
c
        endif
c
        tmax = xtmax(nrec)
        tmin = xtmin(nrec)
        radly = xradly(nrec)
        vwind = xvwind(nrec)
        wind = xwind(nrec)
        tdpt = xtdpt(nrec)
c
c       ...rain is changed from  mm to meters
c
c       CONVERT RAIN FROM MM TO METERS AND STMDUR FROM HR TO SEC
c
        if (stmdur.gt.0.0) then
c
          do 40 jkjk=1,nplane
            rain(jkjk) = rain(jkjk) * 0.001
 40       continue
c
          prcp = prcp * 0.001
          stmdur = stmdur * 3600.0
          avrint = rain(1) / stmdur
          mxint = ip * avrint
        else
          if(prcp.gt.0.0)then
            write (6,1300)mon,day,year
            prcp = 0.0
            do 45 jkjk=1,nplane
               rain(jkjk) = 0.0
 45         continue
          end if
        end if
c
      else
c       READ IN BREAKPOINT DATA
        read (13,*) day, mon, year, nbrkpt, tmax, tmin, radly, vwind,
     1              wind, tdpt
        if (mod(year,4).eq.0) then
          if (mon.eq.3.and.day.eq.1.and.sdate.lt.61) then
            write (6,1000)
            write (6,1200) year
          end if
        end if
c
c       NINT IS THE NUMBER OF BREAKPOINTS OCCURRING ON THE GIVEN DAY.
c       READ IN BREAKPOINT DATA IF RAINFALL OCCURS
c
        if (nbrkpt.gt.0) then
          call brkpt(nbrkpt,mxint)
c
          do 50 jkjk=1,nplane
            rain(jkjk) = prcp
            ninten(jkjk) = nbrkpt
 50       continue
c
        else
c
          do 60 jkjk=1,nplane
            rain(jkjk) = 0.0
            ninten(jkjk) = 0
 60       continue
c
          prcp = 0.0
        end if
      end if
c
c     ...if continuous simulation then get five day antecedent
c     moisture for plant growth
c******************************************************************
c
      if (imodel.eq.1) then
        if (idtot.lt.6) then
          tmp(sdate) = tmin
          tmpx(sdate) = tmax
          tmnavg = (tmnavg*(sdate-1)+tmp(sdate)) / float(sdate)
          tmxavg = (tmxavg*(sdate-1)+tmpx(sdate)) / float(sdate)
        else
          tmp(6) = tmin
          tmpx(6) = tmax
          tmnavg = 0.
          tmxavg = 0.
          do 70 i = 1, 5
            tmp(i) = tmp(i+1)
            tmpx(i) = tmpx(i+1)
            tmnavg = tmnavg + tmp(i)
            tmxavg = tmxavg + tmpx(i)
   70     continue
          tmnavg = tmnavg / 5.
          tmxavg = tmxavg / 5.
        end if
c
      end if
c
      tave = (tmin+tmax) / 2.
c
      if (prcp.le.0.0) then
        do 80 jkl=1,nplane
          norain(jkl) = 0
 80     continue
        return
      end if
c
c
      if (imodel.eq.1) call sumrnf
c
      return
 1000 format ('')
 1100 format ('********************************************************'
     1    /' *** WARNING ***'/
     1    '*Leap year detected with 28 days in February           *'/
     1    '*please modify your climate file by adding 1 day       *'/
     1    '*to February in year ',i4,'.  Leap year annual values    *'/
     1    '*will not be output.                                   *'/
     1    ,' *** WARNING ***',/
     1    '*                                                      *'/
     1    '********************************************************')
 1200 format ('********************************************************'
     1    /'*                       WARNING!!                      *'/
     1    '*Leap year detected with 28 days in February           *'/
     1    '*please modify your climate file by adding 1 day       *'/
     1    '*to February in year ',i4,'.  Simulation will not be     *'/
     1    '*completed.                                            *'/
     1    ,' *** WARNING ***',/
     1    '*                                                      *'/
     1    '********************************************************')
 1300 format ('********************************************************'
     1    /'*                       WARNING!!                      *'/
     1    '* Input precipitation event detected which has a       *'/
     1    '* positive value for precipitation depth and a ZERO    *'/
     1    '* value for precipitation duration.  WEPP has changed  *'/
     1    '* the value for precipitation depth to 0.0.            *'/
     1    '*                                                      *'/
     1    '* PLEASE CHECK YOUR CLIMATE INPUT FILE FOR THE DATE:   *'/
     1    '*        ',i2,'/',i2,'/',i4,'                           ',
     1    '         *'/
     1    ,' *** WARNING ***',/
     1    '*                                                      *'/
     1    '********************************************************')
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
      PROGRAM xbroydn
C     driver for routine broydn
      INTEGER N
      PARAMETER(N=2)
      INTEGER i
      REAL x(N),f(N)
      LOGICAL check
      x(1)=2.
      x(2)=.5
      call broydn(x,N,check)
      call funcv(N,x,f)
      if (check) then
        write(*,*) 'Convergence problems.'
      endif
      write(*,'(1x,a5,t10,a1,t22,a1)') 'Index','x','f'
      do 11 i=1,N
        write(*,'(1x,i2,2x,2f12.6)') i,x(i),f(i)
11    continue
      END
      SUBROUTINE funcv(n,x,f)
      INTEGER n
      REAL x(n),f(n)
      f(1)=x(1)**2+x(2)**2-2.
      f(2)=exp(x(1)-1.)+x(2)**3-2.
      return
      END

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Examples/xbroydn.for


text:
program shape_main

    use :: shape_module

    implicit none

    type(shape_type) :: s1
    real(kind=4)     :: dim1
    real(kind=4)     :: dim2
    real(kind=4)     :: area

    dim1 = 20
    dim2 = 10

    ! Call shape constructor
    s1 = shape_type(dim1, dim2)

    ! Calculate area or volume
    call s1%area(dim1, dim2, area)

    print *, "area = ", area

end program shape_main

max_stars_repo_name: gmengaldo/fortran-gym

max_issues_repo_name: gmengaldo/fortran-gym

max_forks_repo_name: gmengaldo/fortran-gym


text:
!----- LGPL --------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This library is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License as published by the Free Software Foundation version 2.1.
!
!  This library 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
!  Lesser General Public License for more details.
!
!  You should have received a copy of the GNU Lesser General Public
!  License along with this library; 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: ftnunit.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/utils_lgpl/ftnunit/packages/ftnunit/src/ftnunit.f90 $
! ftnunit.f90 --
!     Module that implements part of the "ftnunit" framework:
!     It is inspired by the well-known JUnit framework for
!     integrating unit tests in a Java application.
!
!     The module offers:
!     - a set of common utilities, such as assertion checking
!       routines
!     - a general routine to run the tests if requested
!     - resources that keep track of the status
!
!     Related files:
!     ftnunit_test.f90 -- deprecated
!     runtests.bat
!     runtests.sh
!     runtests.tcl
!
!     TODO:
!     - Test the various assertion failures
!     - HTML output of number of differences
!
!     $Id: ftnunit.f90 5717 2016-01-12 11:35:24Z mourits $
!

! ftnunit_utilities
!     Auxiliary routines (also used in ftnunit_store)
!
module ftnunit_utilities

contains

! ftnunit_get_lun --
!     Auxiliary subroutine to get a free LU-number
! Arguments:
!     lun           The value that can be used
!
subroutine ftnunit_get_lun( lun )
    integer, intent(out) :: lun

    logical       :: opend
    integer, save :: prevlun = 0

    if ( prevlun /= 0 ) then
        inquire( unit = prevlun, opened = opend )
        if ( .not. opend ) then
            lun = prevlun
            return
        endif
    endif

    do prevlun = 10,99
        inquire( unit = prevlun, opened = opend )
        if ( .not. opend ) then
            lun = prevlun
            return
        endif
    enddo

end subroutine ftnunit_get_lun

end module ftnunit_utilities

! ftnunit --
!     Core module for ftnunit
!
module ftnunit
    use ftnunit_utilities
    use ftnunit_hooks

    implicit none

    integer, private, save :: last_test           ! Last test that was started
    integer, private, save :: testno              ! Current test number
    integer, private, save :: nofails             ! Number of assertions that failed
    integer, private, save :: noruns              ! Number of runs so far
    logical, private, save :: call_final = .true. ! Call runtests_final implicitly?
    logical, private, save :: previous   = .false.! Previous test run?
    integer, private, save :: failed_asserts = 0
    logical, private, save :: has_run        = .false.
    character(len=20), private, save :: html_file = 'ftnunit.html'
    character(len=80), private, save :: testname

    interface assert_equal
        module procedure assert_equal_int
        module procedure assert_equal_int1d
        module procedure assert_equal_logical
        module procedure assert_equal_string
    end interface

    interface assert_comparable
        module procedure assert_comparable_real
        module procedure assert_comparable_real1d
        module procedure assert_comparable_real2d
    end interface

contains

! test --
!     Routine to run a unit test
! Arguments:
!     proc          The subroutine implementing the unit test
!     text          Text describing the test
!
subroutine test( proc, text )
    external          :: proc
    character(len=*)  :: text

    integer           :: lun
    integer           :: ierr

    !
    ! Check if the test should run
    !
    testno = testno + 1
    if ( testno <= last_test ) then
        return
    endif

    !
    ! Record the fact that we started the test
    !
    has_run = .true.
    call ftnunit_get_lun( lun )
    open( lun, file = 'ftnunit.lst' )
    write( lun, * ) testno, nofails, noruns, ' ', .true.
    close( lun )

    testname = text
    !
    ! Run the test
    !
    write( *, '(2a)' ) 'Test: ', trim(text)
    call ftnunit_write_html_test_begin( text )
    call ftnunit_hook_test_start( text )

    call proc

    !
    ! No runtime error or premature end of
    ! the program ...
    !
    previous = .true.
    call ftnunit_get_lun( lun )
    open( lun, file = 'ftnunit.lst' )
    write( lun, * ) testno, nofails, noruns, ' ', .true.
    close( lun )

    call ftnunit_hook_test_stop( text )

end subroutine test

! runtests_init --
!     Subroutine to initialise the ftnunit system
! Arguments:
!     None
! Note:
!     Use in conjunction with runtests_final to enable multiple calls
!     to the runtests subroutine. This makes it easier to run tests
!     from different modules, as you have more than one subroutine to
!     do the actual tests.
!
subroutine runtests_init
    call_final = .false.

    if ( .not. ftnunit_file_exists("ftnunit.lst") ) then
        call ftnunit_write_html_header
    endif
end subroutine

! runtests_final --
!     Subroutine to report the overall statistics
! Arguments:
!     None
! Note:
!     Use in conjunction with runtests_init to enable multiple calls
!     to the runtests subroutine. This makes it easier to run tests
!     from different modules, as you have more than one subroutine to
!     do the actual tests.
!
subroutine runtests_final
    if ( ftnunit_file_exists("ftnunit.run") ) then
        write(*,'(a,i5)') 'Number of failed assertions:                ', nofails
        write(*,'(a,i5)') 'Number of runs needed to complete the tests:', noruns
        call ftnunit_remove_file( "ftnunit.lst" )
        call ftnunit_write_html_footer
        call ftnunit_hook_test_completed
        stop
    endif
end subroutine

! runtests_close_report --
!     Subroutine to report the overall statistics
! Arguments:
!     None
! Note:
!     Use in conjunction with runtests_init to enable multiple calls
!     to the runtests subroutine. This makes it easier to run tests
!     from different modules, as you have more than one subroutine to
!     do the actual tests.
!
!     This version does not stop
!
subroutine runtests_close_report
    if ( ftnunit_file_exists("ftnunit.run") ) then
        write(*,'(a,i5)') 'Number of failed assertions:                ', nofails
        write(*,'(a,i5)') 'Number of runs needed to complete the tests:', noruns
        call ftnunit_remove_file( "ftnunit.lst" )
        call ftnunit_write_html_footer
        call ftnunit_hook_test_completed
    endif
end subroutine

! runtests --
!     Subroutine to run the tests if requested
! Arguments:
!     testproc      The test subroutine that actually runs the unit test
!
subroutine runtests( testproc )
    interface
        subroutine testproc
        end subroutine testproc
    end interface

    integer :: lun
    integer :: ierr

    last_test = 0
    nofails   = 0
    noruns    = 0
    testno    = 0
    previous  = .false.
    has_run   = .false.

    if ( ftnunit_file_exists("ftnunit.run") ) then
        if ( ftnunit_file_exists("ftnunit.lst") ) then
            call ftnunit_get_lun( lun )
            open( lun, file = "ftnunit.lst", iostat = ierr )
            if ( ierr == 0 ) then
                read( lun, *, iostat = ierr ) last_test, nofails, noruns, previous
                if ( ierr /= 0 ) then
                    last_test = 0
                    nofails   = 0
                    noruns    = 0
                    previous  = .false.
                endif
                close( lun )
            endif
            if ( previous ) then
                call ftnunit_write_html_previous_failed
            endif
        endif

        noruns = noruns + 1
        if ( noruns == 1 .and. call_final ) then
            call ftnunit_write_html_header
        endif

        call testproc

        if ( call_final ) then
            call runtests_final
        endif

    endif

end subroutine runtests

! assert_true --
!     Subroutine to check if a condition is true
! Arguments:
!     cond          Condition to be checked
!     text          Text describing the assertion
! Side effects:
!     If the assertion fails, this is reported to standard
!     output. Also, nofails is increased by one.
!
subroutine assert_true( cond, text )
    logical, intent(in)          :: cond
    character(len=*), intent(in) :: text

    if ( .not. cond ) then
        nofails = nofails + 1
        write(*,*) '    Condition "',trim(text), '" failed'
        write(*,*) '    It should have been true'
        call ftnunit_write_html_failed_logic( text, .true. )
        call ftnunit_hook_test_assertion_failed( testname, text, "Condition should have been true" )
    endif
end subroutine assert_true

! assert_false --
!     Subroutine to check if a condition is false
! Arguments:
!     cond          Condition to be checked
!     text          Text describing the assertion
! Side effects:
!     If the assertion fails, this is reported to standard
!     output. Also, nofails is increased by one.
!
subroutine assert_false( cond, text )
    logical, intent(in)          :: cond
    character(len=*), intent(in) :: text

    if ( cond ) then
        nofails = nofails + 1
        write(*,*) '    Condition "',trim(text), '" failed'
        write(*,*) '    It should have been false'
        call ftnunit_write_html_failed_logic( text, .false. )
        call ftnunit_hook_test_assertion_failed( testname, text, "Condition should have been false" )
    endif
end subroutine assert_false

! assert_equal_logical --
!     Subroutine to check if two logical values are equivalent
! Arguments:
!     value1        First logical value
!     value2        Second logical value
!     text          Text describing the assertion
! Side effects:
!     If the assertion fails, this is reported to standard
!     output. Also, nofails is increased by one.
!
subroutine assert_equal_logical( value1, value2, text )
    logical, intent(in)          :: value1
    logical, intent(in)          :: value2
    character(len=*), intent(in) :: text

    if ( value1 .neqv. value2 ) then
        nofails = nofails + 1
        write(*,*) '    Assertion "',trim(text), '" failed'
        write(*,*) '    The two logical values are not the same'
        call ftnunit_hook_test_assertion_failed( testname, text, "Logical values are not the same" )
        call ftnunit_write_html_failed_equivalent( text )
    endif
end subroutine assert_equal_logical

! assert_equal_string --
!     Subroutine to check if two string are equal
! Arguments:
!     value1        First string value
!     value2        Second string value
!     text          Text describing the assertion
! Side effects:
!     If the assertion fails, this is reported to standard
!     output. Also, nofails is increased by one.
!
subroutine assert_equal_string( value1, value2, text )
    character(len=*), intent(in) :: value1
    character(len=*), intent(in) :: value2
    character(len=*), intent(in) :: text

    if ( value1 .ne. value2 ) then
        nofails = nofails + 1
        write(*,*) '    Assertion "',trim(text), '" failed'
        write(*,*) '    The two strings are not the same:'
        write(*,*) '    String 1: ' // trim(value1)
        write(*,*) '    String 2: ' // trim(value2)
        call ftnunit_write_html_failed_string( text, value1, value2 )
        call ftnunit_hook_test_assertion_failed( testname, text, "Strings are not the same" )
    endif
end subroutine assert_equal_string

! assert_equal_int --
!     Subroutine to check if two integers are equal
! Arguments:
!     value1        First value
!     value2        Second value
!     text          Text describing the assertion
! Side effects:
!     If the assertion fails, this is reported to standard
!     output. Also, nofails is increased by one.
!
subroutine assert_equal_int( value1, value2, text )
    integer, intent(in)          :: value1
    integer, intent(in)          :: value2
    character(len=*), intent(in) :: text

    if ( value1 /= value2) then
        nofails = nofails + 1
        write(*,*) '    Values not equal: "',trim(text), '" - assertion failed'
        write(*,*) '    Values: ', value1, ' and ', value2
        call ftnunit_write_html_failed_int( text, value1, value2 )
        call ftnunit_hook_test_assertion_failed( testname, text, "Integer values are not the same" )
    endif
end subroutine assert_equal_int

! assert_equal_int1d --
!     Subroutine to check if two integer arrays are equal
! Arguments:
!     array1        First array
!     array2        Second array
!     text          Text describing the assertion
! Side effects:
!     If the assertion fails, this is reported to standard
!     output. Also, nofails is increased by one.
!
subroutine assert_equal_int1d( array1, array2, text )
    integer, dimension(:), intent(in) :: array1
    integer, dimension(:), intent(in) :: array2
    character(len=*), intent(in)      :: text

    integer                           :: i
    integer                           :: count
    logical                           :: addtext

    addtext = .false.

    if ( size(array1) /= size(array2) ) then
        nofails = nofails + 1
        write(*,*) '    Arrays have different sizes: "',trim(text), '" - assertion failed'
        call ftnunit_hook_test_assertion_failed( testname, text, "Arrays have different sizes" )
    else
        if ( any( array1 /= array2 ) ) then
            nofails = nofails + 1
            write(*,*) '    One or more values different: "',trim(text), '" - assertion failed'
            call ftnunit_hook_test_assertion_failed( testname, text, "One or more values are different" )
            count = 0
            do i = 1,size(array1)
                if ( array1(i) /= array2(i) ) then
                    count = count + 1
                    write(*,'(3a10)')    '    Index', '     First', '    Second'
                    if ( count < 50 ) then
                        write(*,'(3i10)')    i, array1(i), array2(i)
                        call ftnunit_write_html_failed_int1d( &
                            text, i, array1(i), array2(i), addtext )
                        addtext = .false.
                    endif
                    write(*,*) 'Number of differences: ', count
                endif
            enddo
        endif
    endif
end subroutine assert_equal_int1d

! assert_comparable_real --
!     Subroutine to check if two reals are approximately equal
! Arguments:
!     value1        First value
!     value2        Second value
!     margin        Allowed margin (relative)
!     text          Text describing the assertion
! Side effects:
!     If the assertion fails, this is reported to standard
!     output. Also, nofails is increased by one.
!
subroutine assert_comparable_real( value1, value2, margin, text )
    real, intent(in)             :: value1
    real, intent(in)             :: value2
    real, intent(in)             :: margin
    character(len=*), intent(in) :: text

    if ( abs(value1-value2) > 0.5 * margin * (abs(value1)+abs(value2)) ) then
        nofails = nofails + 1
        write(*,*) '    Values not comparable: "',trim(text), '" - assertion failed'
        write(*,*) '    Values: ', value1, ' and ', value2
        call ftnunit_hook_test_assertion_failed( testname, text, "One or more values differ more than the margin" )
        call ftnunit_write_html_failed_real( text, value1, value2 )
    endif
end subroutine assert_comparable_real

! assert_comparable_real1d --
!     Subroutine to check if two real arrays are comparable
! Arguments:
!     array1        First array
!     array2        Second array
!     text          Text describing the assertion
! Side effects:
!     If the assertion fails, this is reported to standard
!     output. Also, nofails is increased by one.
!
subroutine assert_comparable_real1d( array1, array2, margin, text )
    real, dimension(:), intent(in)    :: array1
    real, dimension(:), intent(in)    :: array2
    real, intent(in)                  :: margin
    character(len=*), intent(in)      :: text

    integer                           :: i
    integer                           :: count
    logical                           :: addtext

    addtext = .false.

    if ( size(array1) /= size(array2) ) then
        nofails = nofails + 1
        write(*,*) '    Arrays have different sizes: "',trim(text), '" - assertion failed'
        call ftnunit_hook_test_assertion_failed( testname, text, "Arrays have different sizes" )
    else
        if ( any( abs(array1-array2) > 0.5 * margin * (abs(array1)+abs(array2)) ) ) then
            nofails = nofails + 1
            write(*,*) '    One or more values different: "',trim(text), '" - assertion failed'
        call ftnunit_hook_test_assertion_failed( testname, text, "One or more values differ more than the margin" )
            count = 0
            do i = 1,size(array1)
                if ( abs(array1(i)-array2(i)) > &
                         0.5 * margin * (abs(array1(i))+abs(array2(i))) ) then
                    count = count + 1
                    write(*,'(a10,2a15)')    '    Index', '          First', '         Second'
                    if ( count < 50 ) then
                        write(*,'(i10,2e15.5)')    i, array1(i), array2(i)
                        call ftnunit_write_html_failed_real1d( &
                            text, i, array1(i), array2(i), addtext )
                        addtext = .false.
                    endif
                endif
            enddo
            if ( count > 0 ) then
                write(*,*) 'Number of differences: ', count
            endif
        endif
    endif
end subroutine assert_comparable_real1d

! assert_comparable_real2d --
!     Subroutine to check if two two-dimensional real arrays are comparable
! Arguments:
!     array1        First array
!     array2        Second array
!     text          Text describing the assertion
! Side effects:
!     If the assertion fails, this is reported to standard
!     output. Also, nofails is increased by one.
!
subroutine assert_comparable_real2d( array1, array2, margin, text )
    real, dimension(:,:), intent(in)    :: array1
    real, dimension(:,:), intent(in)    :: array2
    real, intent(in)                    :: margin
    character(len=*), intent(in)        :: text

    integer                             :: i
    integer                             :: j
    integer                             :: count
    logical                             :: addtext

    addtext = .false.

    if ( any( shape(array1) /= shape(array2) ) ) then
        nofails = nofails + 1
        write(*,*) '    Arrays have different shapes: "',trim(text), '" - assertion failed'
        call ftnunit_hook_test_assertion_failed( testname, text, "Arrays have different shapes" )
    else
        if ( any( abs(array1-array2) > 0.5 * margin * (abs(array1)+abs(array2)) ) ) then
            nofails = nofails + 1
            write(*,*) '    One or more values different: "',trim(text), '" - assertion failed'
        call ftnunit_hook_test_assertion_failed( testname, text, "One or more values differ more than the margin" )
            count = 0
            do j = 1,size(array1,2)
                do i = 1,size(array1,1)
                    if ( abs(array1(i,j)-array2(i,j)) > &
                             0.5 * margin * (abs(array1(i,j))+abs(array2(i,j))) ) then
                        count = count + 1
                        write(*,'(a10,2a15)')    '    Index', '          First', '         Second'
                        if ( count < 50 ) then
                            write(*,'(2i5,2e15.5)')    i, j, array1(i,j), array2(i,j)
                            call ftnunit_write_html_failed_real2d( &
                                text, i, j, array1(i,j), array2(i,j), addtext )
                            addtext = .false.
                        endif
                    endif
                enddo
            enddo
            if ( count > 0 ) then
                write(*,*) 'Number of differences: ', count
            endif
        endif
    endif
end subroutine assert_comparable_real2d

! assert_files_comparable --
!     Compare two files and establish whether they are equal or not
!
! Arguments:
!     filename1       Name of the first file
!     filename2       Name of the second file
!     text            Text describing the assertion
!     tolerance       Relative tolerance for numbers
!                     (optional; defaults to 1.0e-5)
!
! Note:
!     The files are compared line by line and the items
!     within each line are compared too for numbers that
!     are almost equal.
!     The lines are read into a variable of length line_length
!     Each line may contain up to max_items items (character
!     strings) of length item_length.
!
subroutine assert_files_comparable( filename1, filename2, &
               text, tolerance )

    integer, parameter           :: line_length = 200
    integer, parameter           :: item_length = 20
    integer, parameter           :: max_items   = 100

    character(len=*), intent(in) :: filename1
    character(len=*), intent(in) :: filename2
    character(len=*), intent(in) :: text
    real, intent(in), optional   :: tolerance

    real                         :: tol
    real                         :: value1
    real                         :: value2
    integer                      :: i
    integer                      :: lun1
    integer                      :: lun2
    integer                      :: differences
    integer                      :: ierr1
    integer                      :: ierr2
    logical                      :: header
    character(len=line_length)   :: line1
    character(len=line_length)   :: line2
    character(len=item_length), dimension(max_items) :: item1
    character(len=item_length), dimension(max_items) :: item2

    header      = .true.
    differences = 0
    tol         = 1.0e-5
    if ( present(tolerance) ) tol = tolerance

    call ftnunit_get_lun( lun1 )
    open( lun1, file = filename1, status = 'old', iostat = ierr1 )

    call ftnunit_get_lun( lun2 )
    open( lun2, file = filename2, status = 'old', iostat = ierr2 )

    if ( ierr1 /= 0 .or. ierr2 /= 0 ) then
        nofails = nofails + 1
        call write_header

        if ( ierr1 /= 0 ) then
            write(*, '(a)' ) 'Input file '// trim(filename1) // ' could not be opened'
        endif

        if ( ierr2 /= 0 ) then
            write(*, '(a)' ) 'Input file '// trim(filename2) // ' could not be opened'
        endif

        write(*, '(a)' ) 'Comparison failed'
        call ftnunit_write_html_failed_files( text, filename1, filename2, &
            'One or both files could not be opened' )
        call ftnunit_hook_test_assertion_failed( testname, text, "One or both files could not be opened" )

        return
    endif

    do
        read( lun1, '(a)', iostat = ierr1 ) line1
        read( lun2, '(a)', iostat = ierr2 ) line2

        !
        ! End of file for both?
        !
        if ( ierr1 < 0 .and. ierr2 < 0 ) then
            exit
        endif

        !
        ! End of file or error for only one?
        !
        if ( ierr1 /= 0 .or. ierr2 /= 0 ) then
            call write_header

            if ( ierr1 /= 0 ) then
                write(*, '(a)' ) 'Error/end-of-file reading file '// trim(filename1)
            endif

            if ( ierr2 /= 0 ) then
                write(*, '(a)' ) 'Error/end-of-file reading file '// trim(filename2)
            endif

            write(*, '(a)' ) 'Comparison failed'
            call ftnunit_write_html_failed_files( text, filename1, filename2, &
                'Error/premature end-of-file reading the files' )
            call ftnunit_hook_test_assertion_failed( testname, text, "Error/premature end-of-file reading the files" )
            return
        endif

        !
        ! If the lines are equal, there is nothing to do
        !
        if ( line1 == line2 ) then
            cycle
        endif

        !
        ! The lines are different - but sufficiently different?
        !
        item1 = '?'
        item2 = '?'

        read( line1, *, iostat = ierr1 ) item1
        read( line2, *, iostat = ierr2 ) item2

        !
        ! We are not interested in the error codes, but in
        ! the contents of the items
        !
        do i = 1,max_items
            if ( item1(i) /= item2(i) ) then
                read( item1(i), *, iostat = ierr1 ) value1
                read( item2(i), *, iostat = ierr2 ) value2

                if ( ierr1 == 0 .and. ierr2 == 0 ) then
                    if ( abs(value1-value2) > 0.5 * tol * (abs(value1)+abs(value2)) ) then
                        nofails = nofails + 1
                        call write_header

                        write(*, '(a)' ) 'File 1: ' // trim(line1)
                        write(*, '(a)' ) 'File 2: ' // trim(line2)
                        write(*, '(a,i5,a,e12.4)' ) &
                            'Difference in item ', i, ': ', abs(value1-value2)
                        write(*, '(a,e12.4)' ) &
                            'Mean absolute value       : ', 0.5 * (abs(value1)+abs(value2))

                        differences = differences + 1
                        if ( differences == 1 ) then
                            call ftnunit_write_html_failed_files( text, line1, line2, &
                                'Differences in one or more numerical items' )
                            call ftnunit_hook_test_assertion_failed( testname, text, "Differences in one or more numerical items" )
                        endif
                        exit
                    endif

                elseif ( ierr1 == 0 .neqv. ierr2 == 0 ) then
                    !
                    ! Only one is a number?
                    !
                    nofails = nofails + 1
                    call write_header

                    write(*, '(a)' ) 'File 1: ' // trim(line1)
                    write(*, '(a)' ) 'File 2: ' // trim(line2)
                    write(*, '(a,i5,a,e12.4)' ) &
                        'Item ', i, ' is a number in one file and not in the other'

                    differences = differences + 1
                    if ( differences == 1 ) then
                        call ftnunit_write_html_failed_files( text, line1, line2, &
                             'Corresponding items not both numerical' )
                        call ftnunit_hook_test_assertion_failed( testname, text, "Corresponding items not both numerical" )
                    endif
                    exit
                else
                    !
                    ! Unequal substrings
                    !
                    nofails = nofails + 1
                    call write_header

                    write(*, '(a)' ) 'File 1: ' // trim(line1)
                    write(*, '(a)' ) 'File 2: ' // trim(line2)
                    write(*, '(a,i5,a,e12.4)' ) &
                        'Item ', i, ' differs in the two lines'
                    differences = differences + 1
                    if ( differences == 1 ) then
                        call ftnunit_write_html_failed_files( text, line1, line2, &
                            'Differences in one or more string items' )
                        call ftnunit_hook_test_assertion_failed( testname, text, "Differences in one or more string items" )
                    endif
                    exit
                endif
            endif
        enddo
    enddo

    close( lun1 )
    close( lun2 )
contains
subroutine write_header
    if ( header ) then
        header = .false.
        write(*, '(a)' ) 'Comparing files:'
        write(*, '(a)' ) '    File 1:' // trim(filename1)
        write(*, '(a)' ) '    File 2:' // trim(filename2)
    endif
end subroutine write_header
end subroutine assert_files_comparable

! ftnunit_file_exists --
!     Auxiliary function to see if a file exists
! Arguments:
!     filename      Name of the file to check
! Returns:
!     .true. if the file exists, .false. otherwise
!
logical function ftnunit_file_exists( filename )
    character(len=*), intent(in) :: filename

    inquire( file = filename, exist = ftnunit_file_exists )
end function ftnunit_file_exists

! ftnunit_remove_file --
!     Auxiliary subroutine to remove a file
! Arguments:
!     filename      Name of the file to be removed
!
subroutine ftnunit_remove_file( filename )
    character(len=*), intent(in) :: filename

    integer                      :: lun
    integer                      :: ierr

    call ftnunit_get_lun( lun )
    open( lun, file = filename, iostat = ierr )
    if ( ierr /= 0 ) then
        write(*,'(10a)') '    Could not open file for removal: ', trim(filename)
        nofails = nofails + 1
    else
        close( lun, status = 'delete' )
        if ( ftnunit_file_exists( filename ) ) then
            write(*,'(10a)') '    Removal of file unsuccessful: ', trim(filename)
            nofails = nofails + 1
        endif
    endif

end subroutine ftnunit_remove_file

! ftnunit_make_empty_file --
!     Auxiliary subroutine to make an empty file
! Arguments:
!     filename      Name of the file to be created
!
subroutine ftnunit_make_empty_file( filename )
    character(len=*), intent(in) :: filename

    integer                      :: lun
    integer                      :: ierr

    if ( ftnunit_file_exists( filename ) ) then
        call ftnunit_remove_file( filename )
    endif
    call ftnunit_get_lun( lun )
    open( lun, file = filename, iostat = ierr, status = 'new' )
    if ( ierr /= 0 ) then
        write(*,'(10a)') '    Failed to create empty file: ', trim(filename)
        nofails = nofails + 1
    else
        close( lun )
    endif

end subroutine ftnunit_make_empty_file

! ftnunit_write_html_header --
!     Auxiliary subroutine to write the header of the HTML file
! Arguments:
!     None
!
subroutine ftnunit_write_html_header
    integer :: lun

    call ftnunit_get_lun( lun )
    open( lun, file = html_file )
    write( lun, '(a)' ) &
        '<html>', &
        '<header>', &
        '<title>Results of unit tests</title>', &
        '<style type="text/css">', &
        'span.indent {', &
        '   text-indent: 30;', &
        '}', &
        'span.red {', &
        '   background: red;', &
        '}', &
        'span.green {', &
        '   background: green;', &
        '}', &
        '</style>', &
        '</header>', &
        '<body>', &
        '<h3>Result of unit tests</h3>', &
        '<table>'
    close( lun )

end subroutine ftnunit_write_html_header

! ftnunit_write_html_footer --
!     Auxiliary subroutine to write the footer of the HTML file
! Arguments:
!     None
!
subroutine ftnunit_write_html_footer
    integer :: lun

    if ( .not. has_run ) then
        call ftnunit_write_html_previous_failed
    endif

    call ftnunit_get_lun( lun )
    open( lun, file = html_file, position = 'append' )

    if ( has_run .and. failed_asserts == 0 ) then
        call ftnunit_write_html_close_row( lun )
    endif

    write( lun, '(a)' ) &
        '</tr>',    &
        '</table>'
    write( lun,'(a,i5)') '<p>Number of failed assertions: ', nofails
    write( lun,'(a,i5)') '<br>Number of runs needed to complete the tests:', noruns
    write( lun, '(a)' ) &
        '</body>',  &
        '</html>'
    close( lun )

end subroutine ftnunit_write_html_footer

! ftnunit_write_html_test_begin --
!     Auxiliary subroutine to write the test text to the HTML file
! Arguments:
!     text         Description of the test
!
subroutine ftnunit_write_html_test_begin( text )
    character(len=*)  :: text

    integer           :: lun

    call ftnunit_get_lun( lun )
    open( lun, file = html_file, position = 'append' )

    if ( previous ) then
        if ( failed_asserts == 0 ) then
            write( lun, '(a)' ) &
                '<td><span class="green">OK</span></td></tr>'
        else
            write( lun, '(a)' ) &
                '</tr>'
        endif
    endif

    failed_asserts = 0

    write( lun, '(a)' ) &
        '<tr>', &
        '<td>', trim(text), '</td>'
    close( lun )

    previous = .true.

end subroutine ftnunit_write_html_test_begin

! ftnunit_write_html_previous_failed --
!     Auxiliary subroutine to write the closing of the failed test
!     from the previous run to the HTML file
! Arguments:
!     None
!
! Note:
!     Apparently the test caused a run-time error, so most probably no
!     assertions were checked.
!
subroutine ftnunit_write_html_previous_failed

    integer           :: lun

    call ftnunit_get_lun( lun )
    open( lun, file = html_file, position = 'append' )

    write( lun, '(a)' ) &
        '<td><span class="red">Crashed</span></td></tr>', &
        '<tr><td><span class="indent">Possible run-time failure: check the log file</span></td></tr>'

    close( lun )

    failed_asserts = 1 ! Implicit assertion failed that the test will complete

end subroutine ftnunit_write_html_previous_failed

! ftnunit_write_html_close_row --
!     Auxiliary subroutine to write the closing of a row to the HTML file
! Arguments:
!     lun          LU-number for the HTML file
!
subroutine ftnunit_write_html_close_row( lun )
    integer           :: lun

    if ( failed_asserts > 0 ) then
        if ( failed_asserts == 1 ) then
            write( lun, '(a)' ) &
                '<td><span class="red">Failed</span></td>'
        endif
        write( lun, '(a)' ) &
            '</tr><tr>'
    else
        write( lun, '(a)' ) &
            '<td><span class="green">OK</span></td>'
    endif

end subroutine ftnunit_write_html_close_row

! ftnunit_write_html_failed_logic --
!     Auxiliary subroutine to write a failed logic assertion to the HTML file
! Arguments:
!     text         Description of the test
!     expected     Expected value
!
subroutine ftnunit_write_html_failed_logic( text, expected )
    character(len=*)  :: text
    logical           :: expected

    integer           :: lun

    call ftnunit_get_lun( lun )
    open( lun, file = html_file, position = 'append' )

    failed_asserts = failed_asserts + 1

    call ftnunit_write_html_close_row( lun )

    write( lun, '(a)' ) &
        '<td><span class="indent">', trim(text), '</span></td>', &
        '<td>Value should have been "', trim(merge("true ", "false", expected)), '"</td>'
    close( lun )

end subroutine ftnunit_write_html_failed_logic

! ftnunit_write_html_failed_equivalent --
!     Auxiliary subroutine to write a failed equivalent assertion to the HTML file
! Arguments:
!     text         Description of the test
!     expected     Expected value
!
subroutine ftnunit_write_html_failed_equivalent( text )
    character(len=*)  :: text
    logical           :: expected

    integer           :: lun

    call ftnunit_get_lun( lun )
    open( lun, file = html_file, position = 'append' )

    failed_asserts = failed_asserts + 1

    call ftnunit_write_html_close_row( lun )

    write( lun, '(a)' ) &
        '<td><span class="indent">', trim(text), '</span></td>', &
        '<td>Logical values should have been equal</td>'
    close( lun )

end subroutine ftnunit_write_html_failed_equivalent

! ftnunit_write_html_failed_string --
!     Auxiliary subroutine to write a failed string assertion to the HTML file
! Arguments:
!     text         Description of the test
!     expected     Expected value
!
subroutine ftnunit_write_html_failed_string( text, value1, value2 )
    character(len=*)  :: text
    character(len=*)  :: value1
    character(len=*)  :: value2

    integer           :: lun

    call ftnunit_get_lun( lun )
    open( lun, file = html_file, position = 'append' )

    failed_asserts = failed_asserts + 1

    call ftnunit_write_html_close_row( lun )

    write( lun, '(a)' ) &
        '<td><span class="indent">', trim(text), '</span></td>', &
        '<td>Strings should have been equal:<br>', &
        'String 1: ', value1, '<br>', &
        'String 2: ', value1, '</td>'
    close( lun )

end subroutine ftnunit_write_html_failed_string

! ftnunit_write_html_failed_int --
!     Auxiliary subroutine to write a failed integer assertion to the HTML file
! Arguments:
!     text         Description of the test
!     value1       First value
!     value2       Second value
!
subroutine ftnunit_write_html_failed_int( text, value1, value2 )
    character(len=*)  :: text
    integer           :: value1
    integer           :: value2

    integer           :: lun

    call ftnunit_get_lun( lun )
    open( lun, file = html_file, position = 'append' )

    failed_asserts = failed_asserts + 1

    call ftnunit_write_html_close_row( lun )

    write( lun, '(a)' ) &
        '<td><span class="indent">', trim(text), '</span></td>', &
        '<td>Values: '
    write( lun, '(a,i0,a,i0,a)' ) &
        value1, ' -- ', value2, '</td>'
    close( lun )

end subroutine ftnunit_write_html_failed_int

! ftnunit_write_html_failed_int1d --
!     Auxiliary subroutine to write a failed integer assertion to the HTML file
! Arguments:
!     text         Description of the test
!     idx          Index
!     value1       First value
!     value2       Second value
!     addtext      Add the text or not
!
subroutine ftnunit_write_html_failed_int1d( text, idx, value1, value2, addtext )
    character(len=*)  :: text
    integer           :: idx
    integer           :: value1
    integer           :: value2
    logical           :: addtext

    integer           :: lun

    call ftnunit_get_lun( lun )
    open( lun, file = html_file, position = 'append' )

    failed_asserts = failed_asserts + 1

    call ftnunit_write_html_close_row( lun )

    if ( addtext ) then
        write( lun, '(a)' ) &
            '<td><span class="indent">', trim(text), '</span></td>'
    else
        write( lun, '(a)' ) &
            '<td></td>'
    endif

    write( lun, '(a)' ) &
        '<td>Values at index: '
    write( lun, '(i0,a,i0,a,i0,a)' ) &
        idx, ':', &
        value1, ' -- ', value2, '</td>'
    close( lun )

end subroutine ftnunit_write_html_failed_int1d

! ftnunit_write_html_failed_real --
!     Auxiliary subroutine to write a failed real assertion to the HTML file
! Arguments:
!     text         Description of the test
!     value1       First value
!     value2       Second value
!
subroutine ftnunit_write_html_failed_real( text, value1, value2 )
    character(len=*)  :: text
    real              :: value1
    real              :: value2

    integer           :: lun

    call ftnunit_get_lun( lun )
    open( lun, file = html_file, position = 'append' )

    failed_asserts = failed_asserts + 1

    call ftnunit_write_html_close_row( lun )

    write( lun, '(a)' ) &
        '<td><span class="indent">', trim(text), '</span></td>', &
        '<td>Values: '
    write( lun, '(a,e15.7,a,e15.7,a)' ) &
        value1, ' -- ', value2, '</td>'
    close( lun )

end subroutine ftnunit_write_html_failed_real

! ftnunit_write_html_failed_real1d --
!     Auxiliary subroutine to write a failed real assertion to the HTML file
! Arguments:
!     text         Description of the test
!     idx          Index
!     value1       First value
!     value2       Second value
!     addtext      Add the text or not
!
subroutine ftnunit_write_html_failed_real1d( text, idx, value1, value2, addtext )
    character(len=*)  :: text
    integer           :: idx
    real              :: value1
    real              :: value2
    logical           :: addtext

    integer           :: lun

    call ftnunit_get_lun( lun )
    open( lun, file = html_file, position = 'append' )

    failed_asserts = failed_asserts + 1

    call ftnunit_write_html_close_row( lun )

    if ( addtext ) then
        write( lun, '(a)' ) &
            '<td><span class="indent">', trim(text), '</span></td>'
    else
        write( lun, '(a)' ) &
            '<td></td>'
    endif

    write( lun, '(a)' ) &
        '<td>Values at index: '
    write( lun, '(i0,a,g15.5,a,g15.5,a)' ) &
        idx, ':', &
        value1, ' -- ', value2, '</td>'
    close( lun )
end subroutine ftnunit_write_html_failed_real1d

! ftnunit_write_html_failed_real2d --
!     Auxiliary subroutine to write a failed real assertion to the HTML file
! Arguments:
!     text         Description of the test
!     idx1         Index 1
!     idx2         Index 2
!     value1       First value
!     value2       Second value
!     addtext      Add the text or not
!
subroutine ftnunit_write_html_failed_real2d( text, idx1, idx2, value1, value2, addtext )
    character(len=*)  :: text
    integer           :: idx1
    integer           :: idx2
    real              :: value1
    real              :: value2
    logical           :: addtext

    integer           :: lun

    call ftnunit_get_lun( lun )
    open( lun, file = html_file, position = 'append' )

    failed_asserts = failed_asserts + 1

    call ftnunit_write_html_close_row( lun )

    if ( addtext ) then
        write( lun, '(a)' ) &
            '<td><span class="indent">', trim(text), '</span></td>'
    else
        write( lun, '(a)' ) &
            '<td></td>'
    endif

    write( lun, '(a)' ) &
        '<td>Values at index: '
    write( lun, '(i0,a,i0,a,g15.5,a,g15.5,a)' ) &
        idx1, ',', idx2, ':', &
        value1, ' -- ', value2, '</td>'
    close( lun )
end subroutine ftnunit_write_html_failed_real2d

subroutine ftnunit_write_html_failed_files( text, string1, string2, string3 )
    character(len=*) :: text
    integer          :: type
    character(len=*) :: string1
    character(len=*) :: string2
    character(len=*) :: string3

    integer           :: lun

    call ftnunit_get_lun( lun )
    open( lun, file = html_file, position = 'append' )

    failed_asserts = failed_asserts + 1

    call ftnunit_write_html_close_row( lun )

    write( lun, '(a)' ) &
        '<td><span class="indent">', trim(text), '</span></td>', &
        '<td>Failure: '

    write( lun, '(3a)' ) string3, '<br>'
    write( lun, '(3a)' ) 'File 1: ', trim(string1), '<br>'
    write( lun, '(3a)' ) 'File 2: ', trim(string2), '</td>'

    close( lun )
end subroutine ftnunit_write_html_failed_files

end module ftnunit

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


MODULE dynzad
   !!======================================================================
   !!                       ***  MODULE  dynzad  ***
   !! Ocean dynamics : vertical advection trend
   !!======================================================================
   !! History :  OPA  ! 1991-01  (G. Madec) Original code
   !!   NEMO     0.5  ! 2002-07  (G. Madec) Free form, F90
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   dyn_zad       : vertical advection momentum trend
   !!----------------------------------------------------------------------
   USE oce            ! ocean dynamics and tracers
   USE dom_oce        ! ocean space and time domain
   USE sbc_oce        ! surface boundary condition: ocean
   USE trd_oce        ! trends: ocean variables
   USE trddyn         ! trend manager: dynamics
   !
   USE in_out_manager ! I/O manager
   USE lib_mpp        ! MPP library
   USE prtctl         ! Print control
   USE timing         ! Timing

   IMPLICIT NONE
   PRIVATE

   PUBLIC   dyn_zad       ! routine called by dynadv.F90

   !! * Substitutions
   !!----------------------------------------------------------------------
   !!                   ***  vectopt_loop_substitute  ***
   !!----------------------------------------------------------------------
   !! ** purpose :   substitute the inner loop start/end indices with CPP macro
   !!                allow unrolling of do-loop (useful with vector processors)
   !!----------------------------------------------------------------------
   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id: vectopt_loop_substitute.h90 10068 2018-08-28 14:09:04Z nicolasmartin $
   !! Software governed by the CeCILL license (see ./LICENSE)
   !!----------------------------------------------------------------------
   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id: dynzad.F90 10068 2018-08-28 14:09:04Z nicolasmartin $
   !! Software governed by the CeCILL license (see ./LICENSE)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE dyn_zad ( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE dynzad  ***
      !!
      !! ** Purpose :   Compute the now vertical momentum advection trend and
      !!      add it to the general trend of momentum equation.
      !!
      !! ** Method  :   The now vertical advection of momentum is given by:
      !!         w dz(u) = ua + 1/(e1e2u*e3u) mk+1[ mi(e1e2t*wn) dk(un) ]
      !!         w dz(v) = va + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*wn) dk(vn) ]
      !!      Add this trend to the general trend (ua,va):
      !!         (ua,va) = (ua,va) + w dz(u,v)
      !!
      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
      !!              - Send the trends to trddyn for diagnostics (l_trddyn=T)
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
      !
      INTEGER  ::   ji, jj, jk   ! dummy loop indices
      REAL(wp) ::   zua, zva     ! local scalars
      REAL(wp), DIMENSION(jpi,jpj)     ::   zww
      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwuw, zwvw
      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv
      !!----------------------------------------------------------------------
      !
      IF( ln_timing )   CALL timing_start('dyn_zad')
      !
      IF( kt == nit000 ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'dyn_zad : 2nd order vertical advection scheme'
      ENDIF

      IF( l_trddyn )   THEN         ! Save ua and va trends
         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
         ztrdu(:,:,:) = ua(:,:,:)
         ztrdv(:,:,:) = va(:,:,:)
      ENDIF

      DO jk = 2, jpkm1              ! Vertical momentum advection at level w and u- and v- vertical
         DO jj = 2, jpj                   ! vertical fluxes
            DO ji = 2, jpi             ! vector opt.
               zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
            END DO
         END DO
         DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
            DO ji = 2, jpim1        ! vector opt.
               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )
               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )
            END DO
         END DO
      END DO
      !
      ! Surface and bottom advective fluxes set to zero
      DO jj = 2, jpjm1
         DO ji = 2, jpim1           ! vector opt.
            zwuw(ji,jj, 1 ) = 0._wp
            zwvw(ji,jj, 1 ) = 0._wp
            zwuw(ji,jj,jpk) = 0._wp
            zwvw(ji,jj,jpk) = 0._wp
         END DO
      END DO
      !
      DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points
         DO jj = 2, jpjm1
            DO ji = 2, jpim1       ! vector opt.
               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
               va(ji,jj,jk) = va(ji,jj,jk) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
            END DO
         END DO
      END DO

      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
         DEALLOCATE( ztrdu, ztrdv )
      ENDIF
      !                             ! Control print
      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
      !
      IF( ln_timing )   CALL timing_stop('dyn_zad')
      !
   END SUBROUTINE dyn_zad

   !!======================================================================
END MODULE dynzad


text:
  SUBROUTINE GBSload

    USE basic
    USE fields
    USE constants
    USE futils
    INCLUDE 'mpif.h'
    !
    ! Loads a GBS streamfunction and computes the derivatives which
    ! will be used to find the electric field.
    INTEGER :: k,l,i,tau,tauplus
    DOUBLE PRECISION, DIMENSION(numzpts,numrpts)         :: a
    DOUBLE PRECISION, DIMENSION(numrpts,numzpts)         :: b
    DOUBLE PRECISION, DIMENSION(numrpts)                 :: avg_pot_r
    CHARACTER(LEN=64)			                 :: arraystrmf,tail

  DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: es_potgrid_t_tmp

    INTEGER :: myid, master, numprocs, processor,ierr,stat(MPI_STATUS_SIZE),GBSmaxtimesteps
    CALL mpi_comm_rank(MPI_COMM_WORLD,myid,ierr)
    CALL mpi_comm_size(MPI_COMM_WORLD,numprocs,ierr)
    master = 0
    IF (myid == master) THEN  ! Loading on one processor only, that will then send portion of the perturbation to other processors

!________________________________________________________________________________
  WRITE(*,'(a/)') '=== Load GBS phi(r,z_tilt,t) grid ==='
!________________________________________________________________________________
!
    es_potgrid_t = 0.
    a=0.
    b=0.

    IF (strmfstart == 0) THEN
      write(*,*) 'user must specify a starting record number (strmfstart) from the GBS strmf output -'
      STOP 'leading zeros are built in...stopping'
    ELSE
      CALL openf(GBSfile, fidGBS, 'r')
      CALL getsize(fidGBS, 'data/var2d/time', GBSmaxtimesteps)
    END IF

    IF ((all_bins_same_pot == .FALSE. .AND. GBSmaxtimesteps .LE. strmfstart + numprocs*(numtpts-stalength+stagdist) - stagdist + staglength) .OR. (all_bins_same_pot == .TRUE. .AND. GBSmaxtimesteps .LE. strmfstart + staglength)) THEN
	CALL closef(fidGBS)
	STOP 'GBS simulation is not longer enough for requested simulation...stopping'
    END IF

	    IF (ALLOCATED(es_potgrid_fluct)) THEN
		DEALLOCATE(es_potgrid_fluct)
		ALLOCATE(es_potgrid_fluct(0:numrpts,0:numzpts,0:numtpts))
	    END IF
		ALLOCATE(es_potgrid_t_tmp(0:numrpts,0:numzpts,0:numtpts))

	    tauplus = 0

	    DO tau = 0,numtpts
	      a = 0.0
	      b = 0.0
	      tauplus = strmfstart + tau !+ stagdist*(j1 - 1)
	      write(tail,7001) tauplus
	      7001 FORMAT(i6.6)
	      arraystrmf = 'data/var2d/strmf/'//tail
	      write(*,*) 'arraystrmf ',arraystrmf
	      CALL getarr(fidGBS, arraystrmf, a)
	      ! transpose the array due to the different GBS convention
	      b = transpose(a)
	      ! here, shift the array into the 0:,0: position
	      DO k = 0,numrpts-1
		  DO l = 0,numzpts-1
		    ! multiply by this normalization factor
		    es_potgrid_t_tmp(k,l,tau) = b(k+1,l+1)*phibeamnorm
		  END DO
	      END DO
	    END DO

	    ! adjust the GBS potential to increase the fluctuations
	    ! by a factor of 'GBS_factor' over the background, using
	    ! a file 'avgpotefile' that contains a profile of phi(r) averaged over
	    ! z and time

	    ! adjust the GBS potential to increase the background
	    ! by a factor of 'GBS_bg' while holding fluctuations constant
	    ! using a file 'avgpotefile' that contains phi(r) averaged over
	    ! z and time

	    IF (GBS_factor /= 1 .OR. GBS_bg /= 1) THEN
		open(unit=3,file=avgpotefile, status='old', readonly)
		DO k = 0,numrpts
		  es_potgrid_fluct(k,:,:) = es_potgrid_t_tmp(k,:,:) - phibeamnorm*avg_pot_r(k+1)
		END DO
		es_potgrid_t_tmp = GBS_bg*(es_potgrid_t_tmp - es_potgrid_fluct) + GBS_factor*es_potgrid_fluct
	    END IF

	    ! multiply the e.s. potential by a factor 'GBS_mag'

	    IF (GBS_mag /= 1) THEN
		es_potgrid_t_tmp = GBS_mag*es_potgrid_t_tmp
	    END IF

	    IF (vert_reverse) THEN
		DO l = 0,numzpts-1
		  es_potgrid_t_tmp(:,l,:) = es_potgrid_t_tmp(:,numzpts-1-l,:)
		END DO
	    END IF

	    IF (all_bins_same_pot == .FALSE.) THEN
		strmfstart = strmfstart + numtpts - staglength + stagdist ! Each processor recieve a different portion of the perturbation
	    END IF

    DO processor = master+1,numprocs-1

	    IF ( ALLOCATED(es_potgrid_t) .AND. ALLOCATED(es_potgrid_fluct)) THEN
		DEALLOCATE(es_potgrid_t)
		ALLOCATE(es_potgrid_t(0:numrpts,0:numzpts,0:numtpts))
		DEALLOCATE(es_potgrid_fluct)
		ALLOCATE(es_potgrid_fluct(0:numrpts,0:numzpts,0:numtpts))
	    END IF

	    tauplus = 0

	    DO tau = 0,numtpts
	      a = 0.0
	      b = 0.0
	      tauplus = strmfstart + tau !+ stagdist*(j1 - 1)
	      write(tail,7001) tauplus
	      arraystrmf = 'data/var2d/strmf/'//tail
	      write(*,*) 'arraystrmf ',arraystrmf
	      CALL getarr(fidGBS, arraystrmf, a)
	      ! transpose the array due to the different GBS convention
	      b = transpose(a)
	      ! here, shift the array into the 0:,0: position
	      DO k = 0,numrpts-1
		  DO l = 0,numzpts-1
		    ! multiply by this normalization factor
		    es_potgrid_t(k,l,tau) = b(k+1,l+1)*phibeamnorm
		  END DO
	      END DO
	    END DO

	    ! adjust the GBS potential to increase the fluctuations
	    ! by a factor of 'GBS_factor' over the background, using
	    ! a file 'avgpotefile' that contains a profile of phi(r) averaged over
	    ! z and time

	    ! adjust the GBS potential to increase the background
	    ! by a factor of 'GBS_bg' while holding fluctuations constant
	    ! using a file 'avgpotefile' that contains phi(r) averaged over
	    ! z and time

	    IF (GBS_factor /= 1 .OR. GBS_bg /= 1) THEN
		open(unit=3,file=avgpotefile, status='old', readonly)
		DO k = 0,numrpts
		  es_potgrid_fluct(k,:,:) = es_potgrid_t(k,:,:) - phibeamnorm*avg_pot_r(k+1)
		END DO
		es_potgrid_t = GBS_bg*(es_potgrid_t - es_potgrid_fluct) + GBS_factor*es_potgrid_fluct
	    END IF

	    ! multiply the e.s. potential by a factor 'GBS_mag'

	    IF (GBS_mag /= 1) THEN
		es_potgrid_t = GBS_mag*es_potgrid_t
	    END IF

	    IF (vert_reverse) THEN
		DO l = 0,numzpts-1
		  es_potgrid_t(:,l,:) = es_potgrid_t(:,numzpts-1-l,:)
		END DO
	    END IF

	    IF (all_bins_same_pot == .FALSE.) THEN
		strmfstart = strmfstart + numtpts - staglength + stagdist ! Each processor recieve a different portion of the perturbation
	    END IF

	    ! sending turbulence to other processors

	    IF (processor .NE. master) THEN
		CALL MPI_SEND(es_potgrid_t,(numrpts+1)*(numzpts+1)*(numtpts+1),MPI_DOUBLE_PRECISION,processor,0,MPI_COMM_WORLD,ierr)
	    END IF
    END DO

    CALL closef(fidGBS)
    END IF

    IF (myid .NE. master) THEN
	CALL MPI_RECV(es_potgrid_t,(numrpts+1)*(numzpts+1)*(numtpts+1),MPI_DOUBLE_PRECISION,master,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr)
    END IF

    IF (myid == master) THEN
        es_potgrid_t = es_potgrid_t_tmp
        DEALLOCATE(es_potgrid_t_tmp)
    END IF

    CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)

  END SUBROUTINE GBSload

max_stars_repo_name: alexbovet/fastpart

max_issues_repo_name: alexbovet/fastpart

max_forks_repo_name: alexbovet/fastpart