subroutine mat_matfn6()
!
! ident_11="@(#) M_matrix mat_matfn6(3f) evaluate utility functions"
!
integer :: i, j, k
integer :: ia
integer :: ib
integer :: ja
integer :: jb
integer :: location
integer :: la
integer :: lb
integer :: ld
integer :: lj
integer :: ll
integer :: ls
integer :: m
integer :: ma
integer :: mn
integer :: n
integer :: na
integer :: nn
integer,parameter :: unifor(GG_MAX_NAME_LENGTH) = [iachar(['u','n','i','f','o','r','m']),GG_PAD(8:)]
integer,parameter :: normal(GG_MAX_NAME_LENGTH) = [iachar(['n','o','r','m','a','l',' ']),GG_PAD(8:)]
integer,parameter :: seed(GG_MAX_NAME_LENGTH) = [iachar(['s','e','e','d',' ',' ',' ']),GG_PAD(8:)]
integer :: id(GG_MAX_NAME_LENGTH)
doubleprecision :: eps0,eps,s,sr,si,t
character(len=80) :: message
location = G_VAR_DATALOC(G_ARGUMENT_POINTER)
m = G_VAR_ROWS(G_ARGUMENT_POINTER)
n = G_VAR_COLS(G_ARGUMENT_POINTER)
! functions/G_FIN
! magi diag sum prod user eye rand ones chop shape kron tril triu zeros
! 1 2 3 4 5 6 7 8 9 10 11-13 14 15 16
FUN6: select case(G_FIN)
!===================================================================================================================================
case(1) ! COMMAND::MAGIC
N = MAX(int(GM_REALS(location)),0)
IF (N .EQ. 2) N = 0
IF (N .GT. 0) call mat_magic(GM_REALS(location),N,N)
call mat_rset(N*N,0.0D0,GM_IMAGS(location),1)
G_VAR_ROWS(G_ARGUMENT_POINTER) = N
G_VAR_COLS(G_ARGUMENT_POINTER) = N
!===================================================================================================================================
case(11,12,13) ! COMMAND::KRONECKER PRODUCT
if (G_RHS .ne. 2) then
call mat_err(39) ! Incorrect number of arguments
return
endif
G_ARGUMENT_POINTER = G_ARGUMENT_POINTER - 1
location = G_VAR_DATALOC(G_ARGUMENT_POINTER)
MA = G_VAR_ROWS(G_ARGUMENT_POINTER)
NA = G_VAR_COLS(G_ARGUMENT_POINTER)
LA = location + MAX(M*N*MA*NA,M*N+MA*NA)
LB = LA + MA*NA
if(too_much_memory(LB + M*N - G_VAR_DATALOC(G_TOP_OF_SAVED)) )return
! MOVE A AND B ABOVE RESULT
call mat_wcopy(MA*NA+M*N,GM_REALS(location),GM_IMAGS(location),1,GM_REALS(LA),GM_IMAGS(LA),1)
DO JA = 1, NA
DO J = 1, N
LJ = LB + (J-1)*M
DO IA = 1, MA
! GET J-TH COLUMN OF B
call mat_wcopy(M,GM_REALS(LJ),GM_IMAGS(LJ),1,GM_REALS(location),GM_IMAGS(location),1)
! ADDRESS OF A(IA,JA)
LS = LA + IA-1 + (JA-1)*MA
DO I = 1, M
! A(IA,JA) OP B(I,J)
IF (G_FIN .EQ. 11) &
call mat_wmul( GM_REALS(LS), GM_IMAGS(LS), &
GM_REALS(location), GM_IMAGS(location), &
GM_REALS(location), GM_IMAGS(location))
IF (G_FIN .EQ. 12) &
call mat_wdiv( GM_REALS(LS), GM_IMAGS(LS), &
GM_REALS(location), GM_IMAGS(location), &
GM_REALS(location), GM_IMAGS(location))
IF (G_FIN .EQ. 13) &
call mat_wdiv( GM_REALS(location), GM_IMAGS(location), &
GM_REALS(LS), GM_IMAGS(LS), &
GM_REALS(location), GM_IMAGS(location))
IF (G_ERR .GT. 0) return
location = location + 1
enddo
enddo
enddo
enddo
G_VAR_ROWS(G_ARGUMENT_POINTER) = M*MA
G_VAR_COLS(G_ARGUMENT_POINTER) = N*NA
!===================================================================================================================================
case(9) ! COMMAND::CHOP
eps0 = 1.0d0
do ! recalculate epsilon
eps0 = eps0/2.0d0
t = mat_flop(1.0d0 + eps0)
if (t .le. 1.0d0) exit
enddo
eps0 = 2.0d0*eps0
G_FLOP_COUNTER(2) = int(GM_REALS(location))
if (G_SYM .ne. SEMI) then
write(message,'(''CHOP '',I2,'' PLACES.'')') G_FLOP_COUNTER(2)
call journal(message)
endif
eps = 1.0d0
do ! recalculate epsilon
eps = eps/2.0d0
t = mat_flop(1.0d0 + eps)
if (t .le. 1.0d0) exit
enddo
eps = 2.0d0*eps
t = GM_REALS(GM_BIGMEM-4)
if (t.lt.eps .or. t.eq.eps0) GM_REALS(GM_BIGMEM-4) = eps
G_VAR_ROWS(G_ARGUMENT_POINTER) = 0
!===================================================================================================================================
case(3) ! COMMAND::SUM
sr = 0.0d0
si = 0.0d0
mn = m*n
do i = 1, mn
ls = location+i-1
sr = mat_flop(SR+GM_REALS(LS))
si = mat_flop(SI+GM_IMAGS(LS))
enddo
GM_REALS(location) = sr
GM_IMAGS(location) = si
G_VAR_ROWS(G_ARGUMENT_POINTER) = 1
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
!===================================================================================================================================
case(4) ! COMMAND::PROD
SR = 1.0D0
SI = 0.0D0
MN = M*N
DO I = 1, MN
LS = location+I-1
call mat_wmul(GM_REALS(LS),GM_IMAGS(LS),SR,SI,SR,SI)
enddo
GM_REALS(location) = SR
GM_IMAGS(location) = SI
G_VAR_ROWS(G_ARGUMENT_POINTER) = 1
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
!===================================================================================================================================
case(5) ! COMMAND::USER
! The LALA statement "Y = user(X,s,t)" results in a call to the
! subroutine with a copy of the matrix X stored in the argument A,
! its column and row dimensions in M and N, and the scalar parameters
! s and t stored in S and T. If s and t are omitted, they are set
! to 0.0. After the return, A is stored in Y. The dimensions M and
! N may be reset within the subroutine. The statement Y = user(K)"
! results in a call with M = 1, N = 1 and A(1,1) = "float(K)".
! all of the arguments are in a vector that is part of the stack.
! the location points to the last value and M and N are set to the
! the row and column size of the last argument. G_RHS is the number
! of arguments.
s = 0.0d0
t = 0.0d0
if (G_RHS .eq. 2) then
s = GM_REALS(location)
! back up the stack one argument
G_ARGUMENT_POINTER = G_ARGUMENT_POINTER-1
location = G_VAR_DATALOC(G_ARGUMENT_POINTER) ! the end of argument X
m = G_VAR_ROWS(G_ARGUMENT_POINTER) ! the size of X(M,N)
n = G_VAR_COLS(G_ARGUMENT_POINTER)
elseif(G_RHS.gt.2)then
t = GM_REALS(location)
G_ARGUMENT_POINTER = G_ARGUMENT_POINTER-1 ! back up to s
location = G_VAR_DATALOC(G_ARGUMENT_POINTER)
s = GM_REALS(location)
G_ARGUMENT_POINTER = G_ARGUMENT_POINTER-1 ! back up to X
location = G_VAR_DATALOC(G_ARGUMENT_POINTER)
m = G_VAR_ROWS(G_ARGUMENT_POINTER)
n = G_VAR_COLS(G_ARGUMENT_POINTER)
else ! if not 1,2,3 should it be an error???
endif
! ??? if user routine changes size of array and/or should pass vector instead of address ???
! ??? user routine cannot do complex values? Just REAL ???
call usersub(GM_REALS(location:),m,n,s,t)
call mat_rset(m*n,0.0d0,GM_IMAGS(location),1) ! set the imaginary values to zero
G_VAR_COLS(G_ARGUMENT_POINTER) = n ! store the possibly new size
G_VAR_ROWS(G_ARGUMENT_POINTER) = m
!===================================================================================================================================
case(10) ! COMMAND::SHAPE
! store the two output values onto stack
GM_REALS(location) = M
GM_IMAGS(location) = 0.0D0
GM_REALS(location+1) = N
GM_IMAGS(location+1) = 0.0D0
if(G_LHS.eq.1)then
! output is a 1x2 array so store values indicating the shape of the new stack value
G_VAR_ROWS(G_ARGUMENT_POINTER) = 1
G_VAR_COLS(G_ARGUMENT_POINTER) = 2
else
! output is two scalars
G_VAR_ROWS(G_ARGUMENT_POINTER) = 1
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
G_ARGUMENT_POINTER = G_ARGUMENT_POINTER + 1
G_VAR_DATALOC(G_ARGUMENT_POINTER) = location+1
G_VAR_ROWS(G_ARGUMENT_POINTER) = 1
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
endif
!===================================================================================================================================
case(2,14,15) ! COMMAND::DIAG=2
! COMMAND::TRIL=14
! COMMAND::TRIU=15
k = 0
if (G_RHS .eq. 2) then
k = int(GM_REALS(location))
G_ARGUMENT_POINTER = G_ARGUMENT_POINTER-1
location = G_VAR_DATALOC(G_ARGUMENT_POINTER)
m = G_VAR_ROWS(G_ARGUMENT_POINTER)
n = G_VAR_COLS(G_ARGUMENT_POINTER)
endif
if (G_FIN .ge. 14) then ! COMMAND::TRIL, COMMAND::TRIU
do j = 1, n
ld = location + j - k - 1 + (j-1)*m
select case(G_FIN)
case(14)
ll = j - k - 1
ls = ld - ll
case(15)
ll = m - j + k
ls = ld + 1
end select
if (ll .gt. 0) call mat_wset(ll, 0.0d0, 0.0d0, GM_REALS(ls), GM_IMAGS(ls), 1)
enddo
elseif (m .eq. 1 .or. n .eq. 1) then
n = max(m,n)+iabs(k)
if(too_much_memory( location+n*n - G_VAR_DATALOC(G_TOP_OF_SAVED)) )return
G_VAR_ROWS(G_ARGUMENT_POINTER) = n
G_VAR_COLS(G_ARGUMENT_POINTER) = n
do jb = 1, n
do ib = 1, n
j = n+1-jb
i = n+1-ib
sr = 0.0d0
si = 0.0d0
if (k.ge.0) ls = location+i-1
if (k.lt.0) ls = location+j-1
ll = location+i-1+(j-1)*n
if (j-i .eq. k) sr = GM_REALS(ls)
if (j-i .eq. k) si = GM_IMAGS(ls)
GM_REALS(LL) = sr
GM_IMAGS(LL) = si
enddo
enddo
else
if (k.ge.0) mn=min(m,n-k)
if (k.lt.0) mn=min(m+k,n)
G_VAR_ROWS(G_ARGUMENT_POINTER) = max(mn,0)
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
if (mn .le. 0) exit FUN6
do i = 1, mn
if (k.ge.0) ls = location+(i-1)+(i+k-1)*m
if (k.lt.0) ls = location+(i-k-1)+(i-1)*m
ll = location+i-1
GM_REALS(ll) = GM_REALS(ls)
GM_IMAGS(ll) = GM_IMAGS(ls)
enddo
endif
exit FUN6
!-----------------------------------------------------------------------------------------------------------------------------------
case(6,7,8,16) ! COMMAND::EYE,
! COMMAND::RAND,
! COMMAND::ONES,
! COMMAND::ZEROS
if (.not.(m.gt.1 .or. G_RHS.eq.0)) then
if (G_RHS .eq. 2) then
nn = int(GM_REALS(location))
G_ARGUMENT_POINTER = G_ARGUMENT_POINTER-1
location = G_VAR_DATALOC(G_ARGUMENT_POINTER)
n = G_VAR_COLS(G_ARGUMENT_POINTER)
endif
if (G_FIN.eq.7.and.n.lt.GG_MAX_NAME_LENGTH)then ! a call to RAND might be RAND('UNIFORM'|'SEED'|'NORMAL')
id=blank
do i = 1, min(GG_MAX_NAME_LENGTH,n) ! in case it is one of these words store it in the ID array to test if it matches
ls = location+i-1
id(i) = int(GM_REALS(ls))
enddo
if(mat_eqid(id,unifor).or.mat_eqid(id,normal))then ! SWITCH UNIFORM AND NORMAL(if a matrix just happens to match, a bug)
G_CURRENT_RANDOM_TYPE = id(1) - unifor(1) ! set random type to generate by seeing if first letter is a "u"
G_VAR_ROWS(G_ARGUMENT_POINTER) = 0
exit FUN6
elseif (mat_eqid(id,seed)) then ! if a matrix just happens to match "seed" , a bug)
if (G_RHS .eq. 2) G_CURRENT_RANDOM_SEED = nn
GM_REALS(location) = G_CURRENT_RANDOM_SEED
G_VAR_ROWS(G_ARGUMENT_POINTER) = 1
if (G_RHS .eq. 2) G_VAR_ROWS(G_ARGUMENT_POINTER) = 0
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
exit FUN6
endif
endif
if (n .le. 1) then
m = max(int(GM_REALS(location)),0)
if (G_RHS .eq. 2) n = max(nn,0)
if (G_RHS .ne. 2) n = m
if(too_much_memory( location+m*n - G_VAR_DATALOC(G_TOP_OF_SAVED))) return
G_VAR_ROWS(G_ARGUMENT_POINTER) = m
G_VAR_COLS(G_ARGUMENT_POINTER) = n
if (m*n .eq. 0) exit FUN6
endif
endif
do j = 1, n
do i = 1, m
ll = location+i-1+(j-1)*m ! location to place value
GM_IMAGS(ll) = 0.0d0 ! all of these functions set imaginary values to zero
select case(G_FIN)
case( 6 ) !::EYE
if(i.eq.j)then ! on the diagonal
GM_REALS(ll) = 1.0d0
else
GM_REALS(ll) = 0.0d0
endif
case( 7 ) !::RAND
IF(G_CURRENT_RANDOM_TYPE.EQ.0) then
GM_REALS(ll)=mat_flop(mat_urand(G_CURRENT_RANDOM_SEED))
else
do
sr = 2.0d0*mat_urand(G_CURRENT_RANDOM_SEED)-1.0d0
si = 2.0d0*mat_urand(G_CURRENT_RANDOM_SEED)-1.0d0
t = sr*sr + si*si
if (t .le. 1.0d0) exit
enddo
GM_REALS(ll) = mat_flop(sr*dsqrt((-(2.0d0*dlog(t)))/t))
endif
case( 8 ) !::ONES
GM_REALS(ll) = 1.0d0
case( 16) !::ZEROS
GM_REALS(ll) = 0.0d0
case default
call journal('should not get here: internal error')
end select
enddo
enddo
exit FUN6
!===================================================================================================================================
case(17) ! COMMAND::GETENV JSU
GETENV : block
character(len=:),allocatable :: answers(:)
character(len=GG_LINELEN) :: varname
character(len=:),allocatable :: env_value
allocate(character(len=0) :: answers(0) )
! sort out what to do with an array of input later, for now concatenating into one string
if (m.lt.1 .or. G_RHS.eq.0)then
call journal('sc','<ERROR>getenv:needs an argument:rows=',m,' arg_count=',G_RHS)
G_ERR=999
return
endif
if (G_RHS.gt.1)then
call journal('sc','<ERROR>getenv:too many arguments:arg_count=',G_RHS)
G_ERR=999
return
endif
ll=location
do j=1,m
varname=ade2str( int(GM_REALS(ll:ll+n-1)) )
if(.not.mat_is_name(varname))then
call journal('sc',' function name contains unacceptable characters')
return
endif
ll=ll+n
env_value=system_getenv(varname)
! do not leave it undefined or any variable on LHS will not be defined so make sure at least 1
answers=[character(len=max(len(answers),len_trim(env_value),1)) :: answers,env_value]
enddo
m=size(answers,dim=1)
n=len(answers)
if(too_much_memory( location+m*n - G_VAR_DATALOC(G_TOP_OF_SAVED)) )return
G_VAR_ROWS(G_ARGUMENT_POINTER) = m
G_VAR_COLS(G_ARGUMENT_POINTER) = n
if (m*n .eq. 0) exit FUN6
! so starting at GM_REALS(location) convert the characters to numbers and store the M x N number of characters
do j = 1, n
do i = 1, m
ll = location+i-1+(j-1)*m ! location to place value
GM_IMAGS(ll) = 0.0d0 ! all of these functions set imaginary values to zero
nn=iachar(answers(m)(j:j))
if(nn.gt.0)then
GM_REALS(ll) = real(nn)
else
call journal('sc','bad character')
GM_REALS(ll) = 0.0d0
endif
enddo
enddo
endblock GETENV
exit FUN6
!===================================================================================================================================
case(18) ! COMMAND::DAT
DATETIME: block
integer :: time_values(8)
! store the two output values onto stack
call date_and_time(values=time_values)
GM_REALS(location:location+8-1) = dble(time_values)
GM_IMAGS(location:location+8-1) = 0.0D0
! output is a 1x8 array so store values indicating the size of the new stack value
G_VAR_ROWS(G_ARGUMENT_POINTER) = 1
G_VAR_COLS(G_ARGUMENT_POINTER) = 8
endblock DATETIME
!===================================================================================================================================
end select FUN6
end subroutine mat_matfn6