subroutine mat_matfn3()
! ident_26="@(#) M_matrix mat_matfn3(3fp) evaluate functions involving singular value decomposition"
integer :: i
integer :: j
integer :: jb
integer :: job
integer :: k
integer :: location
integer :: l1
integer :: l2
integer :: ld
integer :: li
integer :: lj
integer :: ll
integer :: ls
integer :: lu
integer :: lv
integer :: m
integer :: mn
integer :: n
logical :: fro,inf
doubleprecision :: p,s,t(1,1),tol,eps
!
if (G_FIN.eq.1 .and. G_RHS.eq.2) 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)
mn = m*n
! SVD PINV COND NORM RANK
! 1 2 3 4 5
FUN3: select case(G_FIN)
!===================================================================================================================================
case(3) ! COMMAND::COND
ld = location + m*n
l1 = ld + min(m+1,n)
l2 = l1 + n
if(too_much_memory( l2+min(m,n) - G_VAR_DATALOC(G_TOP_OF_SAVED) ) )return
call ml_wsvdc(GM_REALS(location),GM_IMAGS(location), &
& m,m,n, &
& GM_REALS(ld),GM_IMAGS(ld), &
& GM_REALS(l1),GM_IMAGS(l1), &
& t,t,1,t,t,1, &
& GM_REALS(l2),GM_IMAGS(l2), &
& 0,G_err)
if (G_err .ne. 0) then
call mat_err(24)
return
endif
s = GM_REALS(ld)
ld = ld + min(m,n) - 1
t(1,1) = GM_REALS(ld)
if (t(1,1) .ne. 0.0d0) then
GM_REALS(location) = mat_flop(s/t(1,1))
GM_IMAGS(location) = 0.0d0
G_VAR_ROWS(G_ARGUMENT_POINTER) = 1
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
else
call journal(' CONDITION IS INFINITE')
G_VAR_ROWS(G_ARGUMENT_POINTER) = 0
endif
!===================================================================================================================================
case(4) ! command::norm
p = 2.0d0
inf = .false.
if (G_RHS .eq. 2)then
fro = int(GM_REALS(location)).eq.iachar('f') .and. mn.gt.1
inf = int(GM_REALS(location)).eq.iachar('i') .and. mn.gt.1
if (.not. fro) then
p = GM_REALS(location)
endif
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)
mn = m*n
if (fro) then
m = mn
n = 1
endif
endif
if (m .gt. 1 .and. n .gt. 1) then
! matrix norm
if (inf)then
s = 0.0d0
do i = 1, m
li = location+i-1
t(1,1) = mat_wasum(n,GM_REALS(LI),GM_IMAGS(li),m)
s = dmax1(s,t(1,1))
enddo
elseif (p .eq. 1.0d0) then
s = 0.0d0
do j = 1, n
lj = location+(j-1)*m
t(1,1) = mat_wasum(m,GM_REALS(LJ),GM_IMAGS(lj),1)
s = dmax1(s,t(1,1))
enddo
elseif (p .ne. 2.0d0) then
call mat_err(23) ! Only 1, 2 or INF norm of matrix
return
else
ld = location + m*n
l1 = ld + min(m+1,n)
l2 = l1 + n
if(too_much_memory( l2+min(m,n) - G_VAR_DATALOC(G_TOP_OF_SAVED) ) )then
return
endif
call ml_wsvdc(GM_REALS(location),GM_IMAGS(location), &
& m,m,n, &
& GM_REALS(ld),GM_IMAGS(ld), &
& GM_REALS(l1),GM_IMAGS(l1), &
& t,t,1,t,t,1, &
& GM_REALS(l2),GM_IMAGS(l2), &
& 0,G_err)
if (G_ERR .ne. 0)then
call mat_err(24)
return
endif
s = GM_REALS(LD)
endif
elseif (p .eq. 1.0d0)then
s = mat_wasum(MN,GM_REALS(location),GM_IMAGS(location),1)
elseif (p .eq. 2.0d0) then
s = mat_wnrm2(MN,GM_REALS(location),GM_IMAGS(location),1)
else
i = mat_iwamax(mn,GM_REALS(location),GM_IMAGS(location),1) + location - 1
s = dabs(GM_REALS(i)) + dabs(GM_IMAGS(i))
if (.not.(inf .or. s .eq. 0.0d0))then
t(1,1) = 0.0d0
do i = 1, mn
ls = location+i-1
t(1,1) = mat_flop(t(1,1) + (mat_pythag(GM_REALS(ls),GM_IMAGS(ls))/s)**p)
enddo
if (p .ne. 0.0d0) then
p = 1.0d0/p
endif
s = mat_flop(s*t(1,1)**p)
endif
endif
GM_REALS(location) = s
GM_IMAGS(location) = 0.0d0
G_VAR_ROWS(G_ARGUMENT_POINTER) = 1
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
!===================================================================================================================================
case(1) ! COMMAND::SVD
IF (G_LHS .EQ. 3)then
K = M
IF (G_RHS .EQ. 2) K = MIN(M,N)
LU = location + M*N
LD = LU + M*K
LV = LD + K*N
L1 = LV + N*N
L2 = L1 + N
if(too_much_memory( L2+MIN(M,N) - G_VAR_DATALOC(G_TOP_OF_SAVED) ) )return
JOB = 11
IF (G_RHS .EQ. 2) JOB = 21
call ml_wsvdc(GM_REALS(location),GM_IMAGS(location), &
& m,m,n, &
& GM_REALS(ld),GM_IMAGS(ld), &
& GM_REALS(l1),GM_IMAGS(l1), &
& GM_REALS(lu),GM_IMAGS(lu), &
& m, &
& GM_REALS(lv),GM_IMAGS(lv), &
& n, &
& GM_REALS(l2),GM_IMAGS(l2), &
& job,G_err)
DO JB = 1, N
DO I = 1, K
J = N+1-JB
LL = LD+I-1+(J-1)*K
IF (I.NE.J) GM_REALS(LL) = 0.0D0
GM_IMAGS(LL) = 0.0D0
LS = LD+I-1
IF (I.EQ.J) GM_REALS(LL) = GM_REALS(LS)
LS = L1+I-1
IF (G_ERR.NE.0 .AND. I.EQ.J-1) GM_REALS(LL) = GM_REALS(LS)
enddo
enddo
IF (G_ERR .NE. 0) call mat_err(24)
G_ERR = 0
call mat_wcopy(M*K+K*N+N*N, &
& GM_REALS(LU),GM_IMAGS(LU), &
& 1, &
& GM_REALS(location),GM_IMAGS(location), &
& 1)
G_VAR_ROWS(G_ARGUMENT_POINTER) = M
G_VAR_COLS(G_ARGUMENT_POINTER) = K
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 + M*K
G_VAR_ROWS(G_ARGUMENT_POINTER) = K
G_VAR_COLS(G_ARGUMENT_POINTER) = 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 + M*K + K*N
G_VAR_ROWS(G_ARGUMENT_POINTER) = N
G_VAR_COLS(G_ARGUMENT_POINTER) = N
else
LD = location + M*N
L1 = LD + MIN(M+1,N)
L2 = L1 + N
if(too_much_memory( L2+MIN(M,N) - G_VAR_DATALOC(G_TOP_OF_SAVED) ) )return
call ml_wsvdc(GM_REALS(location),GM_IMAGS(location),m,m,n, &
& GM_REALS(ld),GM_IMAGS(ld),GM_REALS(l1),GM_IMAGS(l1), &
& t,t,1,t,t,1,GM_REALS(l2),GM_IMAGS(l2),0,G_err)
IF (G_ERR .NE. 0) then
call mat_err(24)
return
endif
K = MIN(M,N)
call mat_wcopy(K,GM_REALS(LD),GM_IMAGS(LD),1,GM_REALS(location),GM_IMAGS(location),1)
G_VAR_ROWS(G_ARGUMENT_POINTER) = K
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
endif
!===================================================================================================================================
case(2,5) ! COMMAND::PINV AND RANK
TOL = -1.0D0
IF (G_RHS .EQ. 2) then
TOL = 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
LU = location + M*N
LD = LU + M*M
IF (G_FIN .EQ. 5) LD = location + M*N
LV = LD + M*N
L1 = LV + N*N
IF (G_FIN .EQ. 5) L1 = LD + N
L2 = L1 + N
if(too_much_memory( L2+MIN(M,N) - G_VAR_DATALOC(G_TOP_OF_SAVED) ) )return
IF (G_FIN .EQ. 2) JOB = 11
IF (G_FIN .EQ. 5) JOB = 0
call ML_WSVDC(GM_REALS(location),GM_IMAGS(location),M,M,N, &
& GM_REALS(LD),GM_IMAGS(LD), &
& GM_REALS(L1),GM_IMAGS(L1), &
& GM_REALS(LU),GM_IMAGS(LU), &
& M, &
& GM_REALS(LV),GM_IMAGS(LV), &
& N, &
& GM_REALS(L2),GM_IMAGS(L2), &
& JOB,G_ERR)
IF (G_ERR .NE. 0) then
call mat_err(24)
return
endif
EPS = GM_REALS(GM_BIGMEM-4)
IF (TOL .LT. 0.0D0) TOL = mat_flop(dble(MAX(M,N))*EPS*GM_REALS(LD))
MN = MIN(M,N)
K = 0
DO J = 1, MN
LS = LD+J-1
S = GM_REALS(LS)
IF (S .LE. TOL) exit
K = J
LL = LV+(J-1)*N
IF (G_FIN .EQ. 2) call mat_wrscal(N,1.0D0/S,GM_REALS(LL),GM_IMAGS(LL),1)
enddo
if (G_FIN .ne. 5) then
do j = 1, m
do i = 1, n
ll = location+i-1+(j-1)*n
l1 = lv+i-1
l2 = lu+j-1
GM_REALS(ll) = mat_wdotcr(k,GM_REALS(l2),GM_IMAGS(l2),m,GM_REALS(l1),GM_IMAGS(l1),n)
GM_IMAGS(ll) = mat_wdotci(k,GM_REALS(l2),GM_IMAGS(l2),m,GM_REALS(l1),GM_IMAGS(l1),n)
enddo
enddo
G_VAR_ROWS(G_ARGUMENT_POINTER) = n
G_VAR_COLS(G_ARGUMENT_POINTER) = m
else
GM_REALS(location) = dble(k)
GM_IMAGS(location) = 0.0d0
G_VAR_ROWS(G_ARGUMENT_POINTER) = 1
G_VAR_COLS(G_ARGUMENT_POINTER) = 1
endif
!===================================================================================================================================
end select FUN3
!
end subroutine mat_matfn3