subroutine mat_matfn2()
integer :: i
integer :: inc
integer :: j
integer :: job
integer :: k
integer :: location
integer :: l1
integer :: l2
integer :: ld
integer :: le
integer :: lj
integer :: ll
integer :: ls
integer :: lw
integer :: m
integer :: n
integer :: nn
!
! evaluate elementary functions and functions involving eigenvalues and eigenvectors
!
doubleprecision tr(1),ti(1),sr,si,powr,powi
logical herm,schur,vect,hess
!
! functions/G_FIN
! ** SIN COS ATAN EXP SQRT LOG
! 0 1 2 3 4 5 6
! EIG SCHU HESS POLY ROOT
! 11 12 13 14 15
! ABS ROUN REAL IMAG CONJ
! 21 22 23 24 25
if (G_FIN .ne. 0) goto 05
location = G_VAR_DATALOC(G_ARGUMENT_POINTER+1)
powr = GM_REALS(location)
powi = GM_IMAGS(location)
05 continue
location = G_VAR_DATALOC(G_ARGUMENT_POINTER)
m = G_VAR_ROWS(G_ARGUMENT_POINTER)
n = G_VAR_COLS(G_ARGUMENT_POINTER)
if (G_FIN .ge. 11 .and. G_FIN .le. 13) goto 10
if (G_FIN .eq. 14 .and. (m.eq.1 .or. n.eq.1))then
goto 50
endif
if (G_FIN .eq. 14) goto 10
if (G_FIN .eq. 15) goto 60
if (G_FIN .gt. 20) goto 40
if (m .eq. 1 .or. n .eq. 1) goto 40
! what about fall-though?
!===================================================================================================================================
! EIGENVALUES AND VECTORS
10 continue
IF (M .NE. N) then
call mat_err(20)
return
endif
SCHUR = G_FIN .EQ. 12
HESS = G_FIN .EQ. 13
VECT = G_LHS.EQ.2 .OR. G_FIN.LT.10
NN = N*N
L2 = location + NN
LD = L2 + NN
LE = LD + N
LW = LE + N
if(too_much_memory( LW+N - G_VAR_DATALOC(G_TOP_OF_SAVED) ) )return
call mat_wcopy(NN,GM_REALS(location),GM_IMAGS(location),1,GM_REALS(L2),GM_IMAGS(L2),1)
!
! CHECK IF HERMITIAN
HERM=.FALSE.
DO J = 1, N
DO I = 1, J
LS = location+I-1+(J-1)*N
LL = location+(I-1)*N+J-1
HERM = GM_REALS(LL).EQ.GM_REALS(LS) .AND. GM_IMAGS(LL).EQ.-GM_IMAGS(LS)
IF (.NOT. HERM) goto 30
enddo
enddo
!
! HERMITIAN EIGENVALUE PROBLEM
call mat_wset(NN,0.0D0,0.0D0,GM_REALS(location),GM_IMAGS(location),1)
call mat_wset(N,1.0D0,0.0D0,GM_REALS(location),GM_IMAGS(location),N+1)
call mat_wset(N,0.0D0,0.0D0,GM_IMAGS(LD),GM_IMAGS(LE),1)
job = 0
IF (VECT) JOB = 1
call ML_HTRIDI(N,N, &
GM_REALS(L2),GM_IMAGS(L2), &
GM_REALS(LD),GM_REALS(LE), &
GM_REALS(LE),GM_REALS(LW))
IF(.NOT.HESS)call ML_IMTQL2(N,N,GM_REALS(LD),GM_REALS(LE),GM_REALS(location),G_ERR,JOB)
IF (G_ERR .GT. 0) then
call mat_err(24)
return
endif
IF (JOB .NE. 0) call ML_HTRIBK(N,N,GM_REALS(L2),GM_IMAGS(L2), &
GM_REALS(LW),N,GM_REALS(location), &
GM_IMAGS(location))
goto 31
!
! NON-HERMITIAN EIGENVALUE PROBLEM
30 continue
call ML_CORTH(N,N,1,N,GM_REALS(L2),GM_IMAGS(L2), &
GM_REALS(LW),GM_IMAGS(LW))
IF (.NOT.VECT .AND. HESS) goto 31
JOB = 0
IF (VECT) JOB = 2
IF (VECT .AND. SCHUR) JOB = 1
IF (HESS) JOB = 3
call ML_COMQR3(N,N,1,N,GM_REALS(LW),GM_IMAGS(LW), &
GM_REALS(L2),GM_IMAGS(L2), &
GM_REALS(LD),GM_IMAGS(LD), &
GM_REALS(location),GM_IMAGS(location), &
G_ERR,JOB)
IF (G_ERR .GT. 0) then
call mat_err(24)
return
endif
!
! VECTORS
31 continue
IF (.NOT.VECT) goto 34
IF (G_ARGUMENT_POINTER+1 .GE. G_TOP_OF_SAVED) then
call mat_err(18)
return
endif
G_ARGUMENT_POINTER = G_ARGUMENT_POINTER+1
G_VAR_DATALOC(G_ARGUMENT_POINTER) = L2
G_VAR_ROWS(G_ARGUMENT_POINTER) = N
G_VAR_COLS(G_ARGUMENT_POINTER) = N
!
! DIAGONAL OF VALUES OR CANONICAL FORMS
34 continue
IF (.NOT.VECT .AND. .NOT.SCHUR .AND. .NOT.HESS) goto 37
DO J = 1, N
LJ = L2+(J-1)*N
IF (SCHUR .AND. (.NOT.HERM)) LJ = LJ+J
IF (HESS .AND. (.NOT.HERM)) LJ = LJ+J+1
LL = L2+J*N-LJ
call mat_wset(LL,0.0D0,0.0D0,GM_REALS(LJ),GM_IMAGS(LJ),1)
enddo
IF (.NOT.HESS .OR. HERM) call mat_wcopy(N,GM_REALS(LD),GM_IMAGS(LD),1,GM_REALS(L2),GM_IMAGS(L2),N+1)
LL = L2+1
IF (HESS .AND. HERM)call mat_wcopy(N-1,GM_REALS(LE+1),GM_IMAGS(LE+1),1,GM_REALS(LL),GM_IMAGS(LL),N+1)
LL = L2+N
IF (HESS .AND. HERM)call mat_wcopy(N-1,GM_REALS(LE+1),GM_IMAGS(LE+1),1,GM_REALS(LL),GM_IMAGS(LL),N+1)
IF (G_FIN .LT. 10) goto 42
IF (VECT .OR. .NOT.(SCHUR.OR.HESS)) goto 99
call mat_wcopy(NN,GM_REALS(L2),GM_IMAGS(L2),1,GM_REALS(location),GM_IMAGS(location),1)
goto 99
!
! VECTOR OF EIGENVALUES
37 continue
IF (G_FIN .EQ. 14) goto 52
call mat_wcopy(N,GM_REALS(LD),GM_IMAGS(LD),1,GM_REALS(location),GM_IMAGS(location),1)
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
goto 99
!===================================================================================================================================
! elementary functions
! for matrices.. x,d = eig(a), fun(a) = x*fun(d)/x
40 continue
inc = 1
n = m*n
l2 = location
goto 44
42 continue
INC = N+1
44 continue
do j = 1, n
ls = l2+(j-1)*inc
sr = GM_REALS(ls)
si = GM_IMAGS(ls)
ti = 0.0d0
if (G_FIN .eq. 0) then
call mat_wlog(sr,si,sr,si)
call mat_wmul(sr,si,powr,powi,sr,si)
tr(1) = dexp(sr)*dcos(si)
ti(1) = dexp(sr)*dsin(si)
endif
select case(G_FIN)
case( 1) ! sin
tr(1) = dsin(sr)*dcosh(si)
ti(1) = dcos(sr)*dsinh(si)
case( 2) ! cos
tr(1) = dcos(sr)*dcosh(si)
ti(1) = (-dsin(sr))*dsinh(si)
case( 3) ! atan
call mat_watan(sr,si,tr(1),ti(1))
case( 4) ! exp
tr(1) = dexp(sr)*dcos(si)
ti(1) = dexp(sr)*dsin(si)
case( 5) ! sqrt
call mat_wsqrt(sr,si,tr(1),ti(1))
case( 6) ! log
call mat_wlog(sr,si,tr(1),ti(1))
case( 21)
tr(1) = mat_pythag(sr,si)
case( 22)
tr(1) = mat_round(sr)
case( 23)
tr(1) = sr
case( 24)
tr(1) = si
case( 25)
tr(1) = sr
ti(1) = -si
end select
if (G_ERR .gt. 0) return
GM_REALS(ls) = mat_flop(tr(1))
GM_IMAGS(ls) = 0.0d0
if (ti(1) .ne. 0.0d0) GM_IMAGS(ls) = mat_flop(ti(1))
enddo
if (inc .eq. 1) goto 99
do j = 1, n
ls = l2+(j-1)*inc
sr = GM_REALS(ls)
si = GM_IMAGS(ls)
ls = location+(j-1)*n
ll = l2+(j-1)*n
call mat_wcopy(n,GM_REALS(ls),GM_IMAGS(ls),1,GM_REALS(ll),GM_IMAGS(ll),1)
call mat_wscal(n,sr,si,GM_REALS(ls),GM_IMAGS(ls),1)
enddo
! signal matfn1 to divide by eigenvectors
G_FUN = 21
G_FIN = -1
G_ARGUMENT_POINTER = G_ARGUMENT_POINTER-1
goto 99
!===================================================================================================================================
! POLY
! form polynomial with given vector as roots
50 continue
N = MAX(M,N)
LD = location+N+1
call mat_wcopy(N,GM_REALS(location),GM_IMAGS(location),1,GM_REALS(LD),GM_IMAGS(LD),1)
goto 52
!===================================================================================================================================
! FORM CHARACTERISTIC POLYNOMIAL
52 continue
call mat_wset(N+1,0.0D0,0.0D0,GM_REALS(location),GM_IMAGS(location),1)
GM_REALS(location) = 1.0D0
DO J = 1, N
call matX_waxpy(J,-GM_REALS(LD),-GM_IMAGS(LD), &
GM_REALS(location),GM_IMAGS(location), &
-1, &
GM_REALS(location+1),GM_IMAGS(location+1), &
-1)
LD = LD+1
enddo
G_VAR_ROWS(G_ARGUMENT_POINTER) = N+1
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
goto 99
!===================================================================================================================================
! ROOTS
60 continue
LL = location+M*N
GM_REALS(LL) = -1.0D0
GM_IMAGS(LL) = 0.0D0
K = -1
61 continue
K = K+1
L1 = location+K
IF (DABS(GM_REALS(L1))+DABS(GM_IMAGS(L1)) .EQ. 0.0D0) goto 61
N = MAX(M*N - K-1, 0)
IF (N .LE. 0) goto 65
L2 = L1+N+1
LW = L2+N*N
if(too_much_memory( LW+N - G_VAR_DATALOC(G_TOP_OF_SAVED) ) )return
call mat_wset(N*N+N,0.0D0,0.0D0,GM_REALS(L2),GM_IMAGS(L2),1)
DO J = 1, N
LL = L2+J+(J-1)*N
GM_REALS(LL) = 1.0D0
LS = L1+J
LL = L2+(J-1)*N
call mat_wdiv(-GM_REALS(LS),-GM_IMAGS(LS), &
GM_REALS(L1),GM_IMAGS(L1), &
GM_REALS(LL),GM_IMAGS(LL))
IF (G_ERR .GT. 0) return
enddo
call ML_COMQR3(N,N,1,N,GM_REALS(LW),GM_IMAGS(LW), &
GM_REALS(L2),GM_IMAGS(L2), &
GM_REALS(location),GM_IMAGS(location), &
TR,TI,G_ERR,0)
IF (G_ERR .GT. 0) then
call mat_err(24)
return
endif
65 continue
G_VAR_ROWS(G_ARGUMENT_POINTER) = N
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
goto 99
!===================================================================================================================================
99 continue
end subroutine mat_matfn2