! 19. Array Types.					File: arraytype19.f90
!------------------------------------------------------------------------------!

MODULE small
	IMPLICIT NONE
	SAVE

	! make the small number epsilon available for use in testing whether 
	! real numbers are effectively zero.

	REAL, PARAMETER :: epsilon=1.0e-5
	REAL :: test
	INTEGER :: i

END MODULE small


! module for real matrix tests
MODULE matrixtestreal
	USE small
	IMPLICIT NONE
CONTAINS

	SUBROUTINE inputreal(nnn,aa,aainv)
		IMPLICIT NONE

		! declarations
		INTEGER, INTENT(IN) :: nnn
		! declare aa, aainv as explicit-shape arrays
		REAL, DIMENSION(nnn,nnn), INTENT(OUT) :: aa, aainv

		! input
		PRINT *
		PRINT *,'Type in the elements of the matrix, column by column'
		PRINT *
		READ *, aa
		PRINT *
		PRINT *,'Now type in the elements of the inverse of the '
		PRINT *,'matrix, column by column'
		PRINT *
		READ *, aainv
		PRINT *

		! check that aa IS inverse of a
		CALL checkreal(aa,aainv)

	END SUBROUTINE inputreal


	SUBROUTINE checkreal(aaa,aaainv)
		IMPLICIT NONE

		! evaluate aaa*aaainv and make sure equals identity matrix

		! declarations

		! declare aaa,aaainv as assumed-shape rank-2 arrays
		REAL, DIMENSION(:,:), INTENT(IN) :: aaa, aaainv

		! declare identity matrix and a temporary 'work' array 
		! as automatic arrays
		REAL, DIMENSION(SIZE(aaa,1),SIZE(aaa,1)) :: lambda, temparr

		! form the identity matrix, lambda
		lambda = 0.0
		! make diagonal elements = 1.0
		DO i=1,SIZE(aaa,1)
			lambda(i,i) = 1.0
		END DO

		! evaluate aaa*aaainv
		temparr = MATMUL(aaa,aaainv)

		! test
		test = MAXVAL(ABS(temparr-lambda))

		IF (test > epsilon) THEN
			PRINT *
			PRINT *,'The matrices you have input are not inverses'
			PRINT *,'of one another. Start again.'
			PRINT *

			STOP
		END IF

	END SUBROUTINE checkreal


	SUBROUTINE symmetric(a)
		IMPLICIT NONE

		! declarations

		! declare the assumed-shape rank-2 array - a
		REAL, DIMENSION(:,:), INTENT(IN) :: a

		! declare the automatic rank-2 array - atrans
		REAL, DIMENSION(SIZE(a,1),SIZE(a,1)) :: atrans

		atrans = TRANSPOSE(a)

		! test
		test = MAXVAL(ABS(atrans - a))

		IF (test < epsilon) THEN
			PRINT *
			PRINT *,'This matrix is symmetric'
			PRINT *
		ELSE
			PRINT *
			PRINT *,'This matrix is NOT symmetric'
			PRINT *
		END IF

	END SUBROUTINE symmetric


	SUBROUTINE orthogonal(a,ainv)
		IMPLICIT NONE

		! declarations

		! declare the assumed shape rank-2 arrays - a,ainv
		REAL, DIMENSION(:,:), INTENT(IN) :: a,ainv

		! declare the automatic rank-2 array - atrans
		REAL, DIMENSION(SIZE(a,1),SIZE(a,1)) :: atrans

		atrans = TRANSPOSE(a)

		! test
		test = MAXVAL(ABS(atrans - ainv))

		IF (test < epsilon) THEN
			PRINT *
			PRINT *,'This matrix is orthogonal'
			PRINT *
		ELSE
			PRINT *
			PRINT *,'This matrix is NOT orthogonal'
			PRINT *
		END IF

	END SUBROUTINE orthogonal

END MODULE matrixtestreal


! Module for complex-valued tests
MODULE matrixtestcomplex
	USE small
	IMPLICIT NONE
CONTAINS

	SUBROUTINE inputcomplex(n,aa,aainv)
		IMPLICIT NONE

		! declarations
		INTEGER, INTENT(IN) :: n
		COMPLEX, DIMENSION(n,n), INTENT(OUT) :: aa, aainv

		! input
		PRINT *
		PRINT *,'Type in the elements of the matrix, column by column'
		PRINT *,'Remember: type in a complex number as a pair of reals'
		PRINT *,'in the form (x,y)'
		PRINT *
		READ *, aa
		PRINT *
		PRINT *,'Now type in the elements of the inverse of the '
		PRINT *,'matrix, column by column'
		PRINT *
		READ *, aainv
		PRINT *

		! check that aa IS inverse of a
		CALL checkcomplex(aa,aainv)

	END SUBROUTINE inputcomplex


	SUBROUTINE checkcomplex(aaa,aaainv)
		IMPLICIT NONE

		! evaluate aaa*aaainv and make sure equals identity matrix

		! declarations
		! declare aaa,aaainv as assumed shape rank-2 arrays
		COMPLEX, DIMENSION(:,:), INTENT(IN) :: aaa, aaainv
		! declare identity matrix and a temporary 'work' array 
		! as automatic arrays
		COMPLEX, DIMENSION(SIZE(aaa,1),SIZE(aaa,1)) :: lambda
		COMPLEX, DIMENSION(SIZE(aaa,1),SIZE(aaa,1)) :: temparr

		! form the identity matrix, lambda
		lambda = (0.0,0.0)
		! make diagonal elements = 1.0
		DO i=1,SIZE(aaa,1)
			lambda(i,i) = (1.0,0.0)
		END DO

		! evaluate aaa*aaainv
		temparr = MATMUL(aaa,aaainv)

		! test
		test = MAXVAL(ABS(temparr-lambda))

		IF (test > epsilon) THEN
			PRINT *
			PRINT *,'The matrices you have input are not inverses'
			PRINT *,'of one another. Start again.'
			PRINT *

			STOP
		END IF

	END SUBROUTINE checkcomplex


	SUBROUTINE hermitian(a)
		IMPLICIT NONE

		! declarations
		! declare the assumed shape rank-2 array - a
		COMPLEX, DIMENSION(:,:), INTENT(IN) :: a
		! declare the automatic rank-2 array - adagger
		COMPLEX, DIMENSION(SIZE(a,1),SIZE(a,1)) :: adagger

		adagger = TRANSPOSE(CONJG(a))

		! test
		test = MAXVAL(ABS(adagger - a))

		IF (test < epsilon) THEN
			PRINT *
			PRINT *,'This matrix is Hermitian'
			PRINT *
		ELSE
			PRINT *
			PRINT *,'This matrix is NOT Hermitian'
			PRINT *
		END IF

	END SUBROUTINE hermitian


	SUBROUTINE unitary(a,ainv)
		IMPLICIT NONE

		! declarations
		! declare the assumed shape rank-2 arrays - a,ainv
		COMPLEX, DIMENSION(:,:), INTENT(IN) :: a,ainv
		! declare the automatic rank-2 array - adagger
		COMPLEX, DIMENSION(SIZE(a,1),SIZE(a,1)) :: adagger, temparr

		adagger = TRANSPOSE(CONJG(a))

		! test
		test = MAXVAL(ABS(adagger - ainv))

		IF (test < epsilon) THEN
			PRINT *
			PRINT *,'This matrix is unitary'
			PRINT *
		ELSE
			PRINT *
			PRINT *,'This matrix is NOT unitary'
			PRINT *
		END IF

	END SUBROUTINE unitary

END MODULE matrixtestcomplex


SUBROUTINE matrixreal(nn)
	USE matrixtestreal
	IMPLICIT NONE

	! declarations
	INTEGER, INTENT(IN) :: nn
	! declare the explicit-shape arrays, a and ainverse
	REAL, DIMENSION(nn,nn) :: a,ainverse

	CALL inputreal(nn,a,ainverse)
	CALL symmetric(a)
	CALL orthogonal(a,ainverse)

END SUBROUTINE matrixreal


SUBROUTINE matrixcomplex(nn)
	USE matrixtestcomplex
	IMPLICIT NONE

	! declarations
	INTEGER, INTENT(IN) :: nn
	! declare the explicit-shape arrays, a and ainverse
	COMPLEX, DIMENSION(nn,nn) :: a,ainverse

	CALL inputcomplex(nn,a,ainverse)
	CALL hermitian(a)
	CALL unitary(a,ainverse)

END SUBROUTINE matrixcomplex


PROGRAM matrix
	IMPLICIT NONE 

	! declarations
	INTEGER :: choice,n
	
	PRINT *
	PRINT *,'This program tests whether a square matrix A is '
	PRINT *,'symmetric, orthogonal, Hermitian, unitary'
	PRINT *
	PRINT *,'If A real then symmetric if A = Atranspose'
	PRINT *,'If A real then orthogonal if Ainverse = Atranspose'
	PRINT *,'If A complex then Hermitian if A = Adagger'
	PRINT *,'If A complex then unitary if Ainverse = Adagger'
	PRINT *,'[Adagger = transpose of complex conjugate of A]'
	PRINT *
	PRINT *,'It requires you to specify '
	PRINT *,'	(i) whether the matrix is real or complex '
	PRINT *,'	(ii) the extent of any dimension of the matrix '
	PRINT *,'	(iii) the matrix A '
	PRINT *,'	(iv) the inverse of the matrix A'
	PRINT *
	PRINT *
	PRINT *,'Type 1 if the matrix is real'
	PRINT *,'Type 2 if the matrix is complex'
	PRINT *,'Type 3 if you want to quit'
	PRINT *
	READ *, choice
	PRINT *


	SELECT CASE(choice)
	CASE(1)
		PRINT *
		PRINT *,'Specify the extent of any dimension of the '
		PRINT *,'real square matrix'
		PRINT *
		READ *, n

		CALL matrixreal(n)
	CASE(2)
		PRINT *
		PRINT *,'Specify the extent of any dimension of the '
		PRINT *,'real square matrix'
		PRINT *
		READ *, n

		CALL matrixcomplex(n)
	CASE(3)
		STOP
	CASE DEFAULT
		PRINT *,'Invalid entry. Try again'
		STOP
	END SELECT

END PROGRAM matrix


!! 1.	Hang on to your hats. The above program involves three modules, ten 
!!	subroutines and the main program unit. It is primarily designed to 
!!	illustrate the different types of arrays which Fortran90 allows. 
!!	Situations arise where procedures have arrays as dummy arguments or, 
!!	within a procedure, a local 'work' array is needed yet, until they are 
!!	called/referenced, an array shape cannot be set. How is a procedure to 
!!	decide on the shape of an array?
!!
!!	Three distinct types of arrays are used in the above program:
!!
!!			(a) explicit-shape arrays
!!			(b) assumed-shape arrays
!!			(c) automatic arrays
!!
!!	EXPLICIT-SHAPE arrays are those where the bounds of each dimension are 
!!	declared explicitly, OR may be calculated from information available 
!!	when the arrays are declared. For example, in the subroutine 
!!	'matrixreal', 'a' and 'ainverse' are explicit-shape since 'n' is known: 
!!
!!		REAL, DIMENSION(nn,nn) :: a,ainverse
!!
!!	ASSUMED-SHAPE arrays MUST be dummy arguments in procedures. Then the 
!!	array assumes the shape of the actual array when the subroutine is 
!!	called or function referenced. Note, in this regard, that it is only 
!!	the shape that must agree and NOT the dimension bounds. The procedure 
!!	MUST have an explicit interface and should be placed in a module. The 
!!	declaration for, say, a real assumed-shape rank-2 array would be
!!
!!		REAL, DIMENSION(lower:,:) :: array
!!
!!	Note that the first dimension specifies a lower index bound which 
!!	can be any integer you desire. When the procedure is called, the 
!!	upper bound will be determined. If 'lower' is absent then it is 
!!	equal to 1 by default, as in the second dimension specification.
!!
!!	AUTOMATIC ARRAYS are explicit-shape arrays placed within procedures 
!!	but which are NOT dummy arguments of the procedure. They are arrays 
!!	local to that procedure and sometimes referred to as 'work' arrays. 
!!	One index bound must not be constant. The array is created dynamically 
!!	on entry to the procedure and 'deallocated' on exit.
!!
!!	DEFERRED-SHAPE (or ALLOCATABLE) arrays may be used instead of automatic 
!!	arrays if more flexibility is required over the creation of temporary 
!!	arrays. See text book for details.

!! 2.	If index bounds are needed within a procedure, they can be passed as 
!!	arguments of the procedure, or be made available by USE of a module.

!! 3.	Note the use of the array intrinsic function MAXVAL to identify the 
!!	maximum array element.
!!
!!	Note the use of SIZE(array,m) to establish the extent of the m-th 
!!	dimension of the array 'array'.

!! 4.	Exercises.
!!
!!	(a) Write a subroutine which swaps two matrices. Use assumed-shape 
!!	arrays as the arguments of the subroutine and an automatic array as 
!!	a local work array to hold data while the swap is being carried out. 
!!	Place the subroutine in a module. Why is this necessary?
!!
!!	(b) Write a function subprogram which evaluates the determinant of a 
!!	real 2x2 matrix. Then write a second function subprogram which obtains 
!!	the inverse of a 2x2 matrix. Include a check to ensure that the 
!!	determinant is not zero. In this the inverse is not defined.

!! 5.	Read chapter 7 and chapter 13.1-13.7 of Ellis.

!!	END OF FILE: arraytype19.f90
