subroutine mat_matfn1()
! ident_25="@(#) M_matrix mat_matfn1(3fp) evaluate functions involving gaussian elimination"
doubleprecision :: dtr(2)
doubleprecision :: dti(2)
doubleprecision :: sr(1)
doubleprecision :: si(1)
doubleprecision :: rcond
doubleprecision :: t
doubleprecision :: t0
doubleprecision :: t1
doubleprecision :: eps
character(len=80) :: mline
integer :: i
integer :: info
integer :: j
integer :: k
integer :: ka
integer :: kb
integer :: location
integer :: l2
integer :: l3
integer :: li
integer :: lj
integer :: lk
integer :: ll
integer :: ls
integer :: lu
integer :: m
integer :: m2
integer :: n
integer :: n2
integer :: nn
!
location = G_VAR_DATALOC(G_ARGUMENT_POINTER)
M = G_VAR_ROWS(G_ARGUMENT_POINTER)
N = G_VAR_COLS(G_ARGUMENT_POINTER)
!===================================================================================================================================
select case(G_FIN)
!===================================================================================================================================
case(-1) ! MATRIX RIGHT DIVISION, A/A2
l2 = G_VAR_DATALOC(G_ARGUMENT_POINTER+1)
m2 = G_VAR_ROWS(G_ARGUMENT_POINTER+1)
n2 = G_VAR_COLS(G_ARGUMENT_POINTER+1)
if (m2 .ne. n2) then
call mat_err(20)
return
endif
if (m*n .ne. 1) then
if (n .ne. n2) then
call mat_err(11)
return
endif
l3 = l2 + m2*n2
if(too_much_memory( l3+n2 - G_VAR_DATALOC(G_TOP_OF_SAVED) ) )return
call ml_wgeco(GM_REALS(l2),GM_IMAGS(l2),m2,n2,G_BUF,rcond,GM_REALS(l3),GM_IMAGS(l3))
if (rcond .eq. 0.0d0) then
call mat_err(19)
return
endif
t = mat_flop(1.0d0 + rcond)
if (t.eq.1.0d0 .and. G_FUN.ne.21)then
call journal('WARNING:')
call journal('MATRIX IS CLOSE TO SINGULAR OR BADLY SCALED.')
WRITE(mline,'(''RESULTS MAY BE INACCURATE. RCOND='',1PD13.4)') RCOND
call journal(mline)
endif
if (t.eq.1.0d0 .and. G_FUN.eq.21)then
call journal('WARNING')
call journal('EIGENVECTORS ARE BADLY CONDITIONED.')
WRITE(mline,'(''RESULTS MAY BE INACCURATE. RCOND='',1PD13.4)') RCOND
call journal(mline)
endif
do i = 1, m
do j = 1, n
ls = location+i-1+(j-1)*m
ll = l3+j-1
GM_REALS(ll) = GM_REALS(ls)
GM_IMAGS(ll) = -GM_IMAGS(ls)
enddo
call ml_wgesl(GM_REALS(l2),GM_IMAGS(l2),m2,n2,G_BUF,GM_REALS(l3),GM_IMAGS(l3),1)
do j = 1, n
ll = location+i-1+(j-1)*m
ls = l3+j-1
GM_REALS(ll) = GM_REALS(ls)
GM_IMAGS(ll) = -GM_IMAGS(ls)
enddo
enddo
if (G_FUN .ne. 21) goto 99
!
! CHECK FOR IMAGINARY ROUNDOFF IN MATRIX FUNCTIONS
sr(1) = mat_wasum(n*n,GM_REALS(location),GM_REALS(location),1)
si(1) = mat_wasum(n*n,GM_IMAGS(location),GM_IMAGS(location),1)
eps = GM_REALS(GM_BIGMEM-4)
t = eps*sr(1)
if (si(1) .le. eps*sr(1)) call mat_rset(n*n,0.0d0,GM_IMAGS(location),1)
goto 99
!
endif
sr(1) = GM_REALS(location)
si(1) = GM_IMAGS(location)
n = n2
m = n
G_VAR_ROWS(G_ARGUMENT_POINTER) = n
G_VAR_COLS(G_ARGUMENT_POINTER) = n
call mat_wcopy(n*n,GM_REALS(l2),GM_IMAGS(l2),1,GM_REALS(location),GM_IMAGS(location),1)
!===================================================================================================================================
case(-2) ! MATRIX LEFT DIVISION A BACKSLASH A2
l2 = G_VAR_DATALOC(G_ARGUMENT_POINTER+1)
m2 = G_VAR_ROWS(G_ARGUMENT_POINTER+1)
n2 = G_VAR_COLS(G_ARGUMENT_POINTER+1)
if (m .ne. n) then
call mat_err(20)
return
endif
if (m2*n2 .ne. 1) then
l3 = l2 + m2*n2
if(too_much_memory( l3+n - G_VAR_DATALOC(G_TOP_OF_SAVED) ) )return
call ml_wgeco(GM_REALS(location),GM_IMAGS(location),m,n,G_BUF,rcond,GM_REALS(l3),GM_IMAGS(l3))
if (rcond .eq. 0.0d0) then
call mat_err(19)
return
endif
t = mat_flop(1.0d0 + rcond)
if (t .eq. 1.0d0) then
call journal('WARNING:')
call journal('MATRIX IS CLOSE TO SINGULAR OR BADLY SCALED.')
WRITE(mline,'(''RESULTS MAY BE INACCURATE. RCOND='',1PD13.4)') RCOND
call journal(mline)
endif
if (m2 .ne. n) then
call mat_err(12)
return
endif
do j = 1, n2
lj = l2+(j-1)*m2
call ml_wgesl(GM_REALS(location),GM_IMAGS(location),m,n,G_BUF,GM_REALS(lj),GM_IMAGS(lj),0)
enddo
G_VAR_COLS(G_ARGUMENT_POINTER) = n2
call mat_wcopy(m2*n2,GM_REALS(l2),GM_IMAGS(l2),1,GM_REALS(location),GM_IMAGS(location),1)
goto 99
endif
sr(1) = GM_REALS(l2)
si(1) = GM_IMAGS(l2)
!===================================================================================================================================
end select
!===================================================================================================================================
select case(G_FIN)
!===================================================================================================================================
case(1) ! COMMAND::INV
if (m .ne. n) then
call mat_err(20)
return
endif
do j = 1, n
do i = 1, n
ls = location+i-1+(j-1)*n
t0 = GM_REALS(ls)
t1 = mat_flop(1.0d0/(dble(i+j-1)))
if (t0 .ne. t1) goto 32
enddo
enddo
call mat_inverse_hilbert(GM_REALS(location),n,n)
call mat_rset(n*n,0.0d0,GM_IMAGS(location),1)
if (G_FIN .lt. 0) call mat_wscal(n*n,sr(1),si(1),GM_REALS(location),GM_IMAGS(location),1)
goto 99
32 continue
l3 = location + n*n
if(too_much_memory( l3+n - G_VAR_DATALOC(G_TOP_OF_SAVED) ) )return
call ml_wgeco(GM_REALS(location),GM_IMAGS(location),m,n,G_BUF,rcond,GM_REALS(l3),GM_IMAGS(l3))
if (rcond .eq. 0.0d0) then
call mat_err(19)
return
endif
t = mat_flop(1.0d0 + rcond)
if (t .eq. 1.0d0) then
call journal('warning:')
call journal('matrix is close to singular or badly scaled.')
write(mline,'(''results may be inaccurate. rcond='',1pd13.4)') rcond
call journal(mline)
endif
call ml_wgedi(GM_REALS(location),GM_IMAGS(location),m,n,G_BUF,dtr,dti,GM_REALS(l3),GM_IMAGS(l3),1)
if (G_FIN .lt. 0) call mat_wscal(n*n,sr(1),si(1),GM_REALS(location),GM_IMAGS(location),1)
!===================================================================================================================================
case (2) ! COMMAND::DET
if (m .ne. n) then
call mat_err(20)
return
endif
call ml_wgefa(GM_REALS(location),GM_IMAGS(location),m,n,G_BUF,info)
!SUBROUTINE ML_WGEDI(ar,ai,LDA,N,ipvt,detr,deti,workr,worki,JOB)
call ml_wgedi(GM_REALS(location),GM_IMAGS(location),m,n,G_BUF,dtr,dti,sr(1),si(1),10)
k = int(dtr(2))
ka = iabs(k)+2
t = 1.0d0
do i = 1, ka
t = t/10.0d0
if (t .ne. 0.0d0) goto 42
enddo
GM_REALS(location) = dtr(1)*10.d0**k
GM_IMAGS(location) = dti(1)*10.d0**k
G_VAR_ROWS(G_ARGUMENT_POINTER) = 1
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
goto 99
42 continue
if (dti(1) .eq. 0.0d0)then
write(mline,43) dtr(1),k
call journal(mline)
else
write(mline,44) dtr(1),dti(1),k
call journal(mline)
endif
GM_REALS(location) = dtr(1)
GM_IMAGS(location) = dti(1)
GM_REALS(location+1) = dtr(2)
GM_IMAGS(location+1) = 0.0d0
G_VAR_ROWS(G_ARGUMENT_POINTER) = 1
G_VAR_COLS(G_ARGUMENT_POINTER) = 2
43 format(' det = ',f7.4,' * 10**',i4)
44 format(' det = ',f7.4,' + ',f7.4,' i ',' * 10**',i4)
!===================================================================================================================================
case(3) ! COMMAND::RCOND
if (m .ne. n) then
call mat_err(20)
return
endif
l3 = location + n*n
if(too_much_memory( l3+n - G_VAR_DATALOC(G_TOP_OF_SAVED) ) )return
call ml_wgeco(GM_REALS(location),GM_IMAGS(location),m,n,G_BUF,rcond,GM_REALS(l3),GM_IMAGS(l3))
GM_REALS(location) = rcond
GM_IMAGS(location) = 0.0d0
G_VAR_ROWS(G_ARGUMENT_POINTER) = 1
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
if (G_lhs .ne. 1)then
location = location + 1
call mat_wcopy(n,GM_REALS(l3),GM_IMAGS(l3),1,GM_REALS(location),GM_IMAGS(location),1)
G_ARGUMENT_POINTER = G_ARGUMENT_POINTER + 1
G_VAR_DATALOC(G_ARGUMENT_POINTER) = location
G_VAR_ROWS(G_ARGUMENT_POINTER) = n
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
endif
!===================================================================================================================================
case(4) ! COMMAND::LU
if (m .ne. n) then
call mat_err(20)
return
endif
call ml_wgefa(GM_REALS(location),GM_IMAGS(location),m,n,G_BUF,info)
if (G_lhs .ne. 2) goto 99
nn = n*n
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) = location + nn
G_VAR_ROWS(G_ARGUMENT_POINTER) = n
G_VAR_COLS(G_ARGUMENT_POINTER) = n
if(too_much_memory( location+nn+nn - G_VAR_DATALOC(G_TOP_OF_SAVED) ) )return
do kb = 1, n
k = n+1-kb
do i = 1, n
ll = location+i-1+(k-1)*n
lu = ll + nn
if (i .le. k) GM_REALS(lu) = GM_REALS(ll)
if (i .le. k) GM_IMAGS(lu) = GM_IMAGS(ll)
if (i .gt. k) GM_REALS(lu) = 0.0d0
if (i .gt. k) GM_IMAGS(lu) = 0.0d0
if (i .lt. k) GM_REALS(ll) = 0.0d0
if (i .lt. k) GM_IMAGS(ll) = 0.0d0
if (i .eq. k) GM_REALS(ll) = 1.0d0
if (i .eq. k) GM_IMAGS(ll) = 0.0d0
if (i .gt. k) GM_REALS(ll) = -GM_REALS(ll)
if (i .gt. k) GM_IMAGS(ll) = -GM_IMAGS(ll)
enddo
i = G_BUF(k)
if (i .eq. k) cycle
li = location+i-1+(k-1)*n
lk = location+k-1+(k-1)*n
call mat_wswap(n-k+1,GM_REALS(li),GM_IMAGS(li),n,GM_REALS(lk),GM_IMAGS(lk),n)
enddo
!===================================================================================================================================
case(5) ! COMMAND::inverse_hilbert
n = int(GM_REALS(location))
G_VAR_ROWS(G_ARGUMENT_POINTER) = n
G_VAR_COLS(G_ARGUMENT_POINTER) = n
call mat_inverse_hilbert(GM_REALS(location),n,n)
call mat_rset(n*n,0.0d0,GM_IMAGS(location),1)
if (G_FIN .lt. 0) call mat_wscal(n*n,sr(1),si(1),GM_REALS(location),GM_IMAGS(location),1)
!===================================================================================================================================
case(6) ! COMMAND::CHOLESKY
if (m .ne. n) then
call mat_err(20)
return
endif
call mat_wpofa(GM_REALS(location),GM_IMAGS(location),m,n,G_err)
if (G_err .ne. 0) then
call mat_err(29)
return
endif
do j = 1, n
ll = location+j+(j-1)*m
call mat_wset(m-j,0.0d0,0.0d0,GM_REALS(ll),GM_IMAGS(ll),1)
enddo
!===================================================================================================================================
case(7) ! COMMAND::RREF
if (G_RHS .ge. 2)then
G_ARGUMENT_POINTER = G_ARGUMENT_POINTER-1
location = G_VAR_DATALOC(G_ARGUMENT_POINTER)
if (G_VAR_ROWS(G_ARGUMENT_POINTER) .ne. m) then
call mat_err(5)
return
endif
n = n + G_VAR_COLS(G_ARGUMENT_POINTER)
endif
call mat_rref(GM_REALS(location),GM_IMAGS(location),m,m,n,GM_REALS(GM_BIGMEM-4))
G_VAR_COLS(G_ARGUMENT_POINTER) = n
!===================================================================================================================================
end select
!
99 continue
end subroutine mat_matfn1