subroutine mat_stack2(op)
! ident_30="@(#) M_matrix ml_stackp(3fp) binary and ternary operations"
integer :: op
doubleprecision :: sr,si,e1,st,e2
integer :: i
integer :: j
integer :: k
integer :: k1
integer :: k2
integer :: kexp
integer :: location
integer :: l1
integer :: l2
integer :: l3
integer :: ll
integer :: ls
integer :: m
integer :: m2
integer :: mn
integer :: n
integer :: n2
integer :: nexp
integer :: op_select
l2 = G_VAR_DATALOC(G_ARGUMENT_POINTER)
m2 = G_VAR_ROWS(G_ARGUMENT_POINTER)
n2 = G_VAR_COLS(G_ARGUMENT_POINTER)
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)
G_FUN = 0
if(op.eq.DSTAR)then
op_select=-op
else
op_select=op
endif
DO_OP: select case(op_select)
!-----------------------------------------------------------------------------------------------------------------------------------
case (PLUS) ! ADDITION
if (m .lt. 0) then
if (m2 .ne. n2) then
call mat_err(8)
exit DO_OP
endif
m = m2
n = n2
G_VAR_ROWS(G_ARGUMENT_POINTER) = m
G_VAR_COLS(G_ARGUMENT_POINTER) = n
sr = GM_REALS(location)
si = GM_IMAGS(location)
call mat_wcopy(m*n,GM_REALS(location+1),GM_IMAGS(location+1),1,GM_REALS(location),GM_IMAGS(location),1)
call finish()
exit DO_OP
endif
if (m2 .lt. 0)then
if (m .ne. n) then
call mat_err(8)
exit DO_OP
endif
sr = GM_REALS(l2)
si = GM_IMAGS(l2)
call finish()
exit DO_OP
endif
if (m .ne. m2) then
call mat_err(8)
exit DO_OP
endif
if (n .ne. n2) then
call mat_err(8)
exit DO_OP
endif
call matX_waxpy(m*n,1.0d0,0.0d0,GM_REALS(l2),GM_IMAGS(l2),1,GM_REALS(location),GM_IMAGS(location),1)
!-----------------------------------------------------------------------------------------------------------------------------------
case (MINUS) ! SUBTRACTION
if (m .lt. 0) then
if (m2 .ne. n2)then
call mat_err(9)
exit do_op
endif
m = m2
n = n2
G_VAR_ROWS(G_ARGUMENT_POINTER) = m
G_VAR_COLS(G_ARGUMENT_POINTER) = n
sr = GM_REALS(location)
si = GM_IMAGS(location)
call mat_wcopy(m*n,GM_REALS(location+1),GM_IMAGS(location+1),1,GM_REALS(location),GM_IMAGS(location),1)
call mat_wrscal(m*n,-1.0d0,GM_REALS(location),GM_IMAGS(location),1)
call finish()
exit DO_OP
endif
if (m2 .lt. 0) then
! add or subtract scalar
if (m .ne. n) then
call mat_err(9)
exit DO_OP
endif
sr = -GM_REALS(l2)
si = -GM_IMAGS(l2)
call finish()
exit DO_OP
endif
if (m .ne. m2)then
call mat_err(9)
exit DO_OP
endif
if (n .ne. n2) then
call mat_err(9)
exit DO_OP
endif
call matX_waxpy(M*N,-1.0D0,0.0D0,GM_REALS(L2),GM_IMAGS(L2),1,GM_REALS(location),GM_IMAGS(location),1)
!-----------------------------------------------------------------------------------------------------------------------------------
case (STAR) ! MULTIPLICATION
if (m2*m2*n2 .eq. 1) goto 10
if (m*n .eq. 1) goto 11
if (m2*n2 .eq. 1) goto 10
if (n .ne. m2) then
call mat_err(10)
exit do_op
endif
mn = m*n2
ll = location + mn
if(too_much_memory( ll+m*n+m2*n2 - G_VAR_DATALOC(G_TOP_OF_SAVED)) ) exit do_op
call mat_wcopy(m*n+m2*n2,GM_REALS(location),GM_IMAGS(location),-1,GM_REALS(ll),GM_IMAGS(ll),-1)
do j = 1, n2
do i = 1, m
k1 = location + mn + (i-1)
k2 = l2 + mn + (j-1)*m2
k = location + (i-1) + (j-1)*m
GM_REALS(k) = mat_wdotur(N,GM_REALS(k1),GM_IMAGS(k1),m,GM_REALS(k2),GM_IMAGS(k2),1)
GM_IMAGS(k) = mat_wdotui(N,GM_REALS(k1),GM_IMAGS(k1),m,GM_REALS(k2),GM_IMAGS(k2),1)
enddo
enddo
G_VAR_COLS(G_ARGUMENT_POINTER) = n2
exit do_op
!-----------------------------------------------------------------------------------------------------------------------------------
! multiplication by scalar
10 continue
sr = GM_REALS(l2)
si = GM_IMAGS(l2)
l1 = location
goto 13
11 continue
sr = GM_REALS(location)
si = GM_IMAGS(location)
l1 = location+1
G_VAR_ROWS(G_ARGUMENT_POINTER) = m2
G_VAR_COLS(G_ARGUMENT_POINTER) = n2
13 continue
mn = G_VAR_ROWS(G_ARGUMENT_POINTER)*G_VAR_COLS(G_ARGUMENT_POINTER)
call mat_wscal(mn,sr,si,GM_REALS(l1),GM_IMAGS(l1),1)
if (l1.ne.location) call mat_wcopy(mn,GM_REALS(l1),GM_IMAGS(l1),1,GM_REALS(location),GM_IMAGS(location),1)
!-----------------------------------------------------------------------------------------------------------------------------------
case (-DSTAR) ! POWER
IF (M2*N2 .NE. 1) then
call mat_err(30)
exit do_op
endif
IF (M .NE. N) then
call mat_err(20)
exit do_op
endif
NEXP = int(GM_REALS(L2))
IF ( (GM_REALS(L2) .NE. dble(NEXP)) .or. (GM_IMAGS(L2) .NE. 0.0D0) .or. (NEXP .LT. 2) )then
! NONINTEGER OR NONPOSITIVE POWER, USE EIGENVECTORS
G_FUN = 2
G_FIN = 0
exit DO_OP
endif
MN = M*N
if(too_much_memory( L2+MN+N - G_VAR_DATALOC(G_TOP_OF_SAVED)) ) exit do_op
call mat_wcopy(MN,GM_REALS(location),GM_IMAGS(location),1,GM_REALS(L2),GM_IMAGS(L2),1)
L3 = L2+MN
DO KEXP = 2, NEXP
DO J = 1, N
LS = location+(J-1)*N
call mat_wcopy(N,GM_REALS(LS),GM_IMAGS(LS),1,GM_REALS(L3),GM_IMAGS(L3),1)
DO I = 1, N
LS = L2+I-1
LL = location+I-1+(J-1)*N
GM_REALS(LL)=mat_wdotur(N,GM_REALS(LS),GM_IMAGS(LS),N,GM_REALS(L3),GM_IMAGS(L3),1)
GM_IMAGS(LL)=mat_wdotui(N,GM_REALS(LS),GM_IMAGS(LS),N,GM_REALS(L3),GM_IMAGS(L3),1)
enddo
enddo
enddo
!-----------------------------------------------------------------------------------------------------------------------------------
case (SLASH) ! right division
if (m2*n2 .ne. 1) then
if (m2 .eq. n2) G_FUN = 1
if (m2 .ne. n2) G_FUN = 4
G_FIN = -1
G_RHS = 2
exit DO_OP
endif
sr = GM_REALS(l2)
si = GM_IMAGS(l2)
mn = m*n
do i = 1, mn
ll = location+i-1
call mat_wdiv(GM_REALS(ll),GM_IMAGS(ll),sr,si,GM_REALS(ll),GM_IMAGS(ll))
if (G_ERR .gt. 0) exit
enddo
!-----------------------------------------------------------------------------------------------------------------------------------
case (BSLASH) ! LEFT DIVISION
if (m*n .ne. 1) then
if (m .eq. n) G_FUN = 1
if (m .ne. n) G_FUN = 4
G_FIN = -2
G_RHS = 2
exit DO_OP
endif
SR = GM_REALS(location)
SI = GM_IMAGS(location)
G_VAR_ROWS(G_ARGUMENT_POINTER) = M2
G_VAR_COLS(G_ARGUMENT_POINTER) = N2
MN = M2*N2
DO I = 1, MN
LL = location+I-1
call mat_wdiv(GM_REALS(LL+1),GM_IMAGS(LL+1),SR,SI,GM_REALS(LL),GM_IMAGS(LL))
IF (G_ERR .GT. 0) exit
enddo
!-----------------------------------------------------------------------------------------------------------------------------------
case (COLON) ! COLON
E2 = GM_REALS(L2)
ST = 1.0D0
N = 0
IF (G_RHS .GE. 3) then
ST = GM_REALS(location)
G_ARGUMENT_POINTER = G_ARGUMENT_POINTER-1
location = G_VAR_DATALOC(G_ARGUMENT_POINTER)
IF (ST .EQ. 0.0D0) goto 63
endif
E1 = GM_REALS(location)
! CHECK FOR CLAUSE
IF (G_RSTK(G_PT) .EQ. 3) then
! FOR CLAUSE
GM_REALS(location) = E1
GM_REALS(location+1) = ST
GM_REALS(location+2) = E2
G_VAR_ROWS(G_ARGUMENT_POINTER) = -3
G_VAR_COLS(G_ARGUMENT_POINTER) = -1
exit DO_OP
endif
if(too_much_memory( location + MAX(3,int((E2-E1)/ST)) - G_VAR_DATALOC(G_TOP_OF_SAVED) ) ) exit do_op
do
IF (ST .GT. 0.0D0 .AND. GM_REALS(location) .GT. E2) exit
IF (ST .LT. 0.0D0 .AND. GM_REALS(location) .LT. E2) exit
N = N+1
location = location+1
GM_REALS(location) = E1 + dble(N)*ST
GM_IMAGS(location) = 0.0D0
enddo
63 continue
G_VAR_COLS(G_ARGUMENT_POINTER) = N
G_VAR_ROWS(G_ARGUMENT_POINTER) = 1
IF (N .EQ. 0) G_VAR_ROWS(G_ARGUMENT_POINTER) = 0
!-----------------------------------------------------------------------------------------------------------------------------------
case (1000:2000-1) ! element-wise operations
op = op -1000
if (m.ne.m2 .or. n.ne.n2) then
call mat_err(10)
exit do_op
endif
mn = m*n
do i = 1, mn
j = location+i-1
k = l2+i-1
select case(op)
case(STAR)
call mat_wmul(GM_REALS(J),GM_IMAGS(J), &
GM_REALS(K),GM_IMAGS(K), &
GM_REALS(J),GM_IMAGS(J))
case(SLASH)
call mat_wdiv(GM_REALS(J),GM_IMAGS(J), &
GM_REALS(K),GM_IMAGS(K), &
GM_REALS(J),GM_IMAGS(J))
case(BSLASH)
call mat_wdiv(GM_REALS(K),GM_IMAGS(K), &
GM_REALS(J),GM_IMAGS(J), &
GM_REALS(J),GM_IMAGS(J))
end select
IF (G_ERR .GT. 0) exit
enddo
!-----------------------------------------------------------------------------------------------------------------------------------
case (2000:) ! kronecker
G_FIN = op - 2000 - star + 11
G_FUN = 6
G_ARGUMENT_POINTER = G_ARGUMENT_POINTER + 1
G_RHS = 2
!-----------------------------------------------------------------------------------------------------------------------------------
case default
write(*,*)'<ERROR> unknown operator ',op
stop
end select DO_OP
!-----------------------------------------------------------------------------------------------------------------------------------
contains
subroutine finish()
do i = 1, n
ll = location + (i-1)*(n+1)
GM_REALS(ll) = mat_flop(GM_REALS(LL)+sr)
GM_IMAGS(ll) = mat_flop(GM_IMAGS(LL)+si)
enddo
end subroutine finish
end subroutine mat_stack2