SUBROUTINE mat_matfn4()
! ident_27="@(#) M_matrix mat_matfn4(3fp) evaluate functions involving qr decomposition (least squares)"
integer :: info
integer :: j
integer :: jb
integer :: job
integer :: k
integer :: location
integer :: l2
integer :: l3
integer :: l4
integer :: le
integer :: ll
integer :: ls
integer :: m
integer :: m2
integer :: mm
integer :: mn
integer :: n
integer :: n2
character(len=81) :: message
DOUBLEPRECISION :: T(1),TOL,EPS
!
location = G_VAR_DATALOC(G_ARGUMENT_POINTER)
M = G_VAR_ROWS(G_ARGUMENT_POINTER)
N = G_VAR_COLS(G_ARGUMENT_POINTER)
IF (G_FIN .EQ. -1) then
goto 10
elseIF (G_FIN .EQ. -2) then
goto 20
else
goto 40
endif
!
! RECTANGULAR MATRIX RIGHT DIVISION, A/A2
10 continue
L2 = G_VAR_DATALOC(G_ARGUMENT_POINTER+1)
M2 = G_VAR_ROWS(G_ARGUMENT_POINTER+1)
N2 = G_VAR_COLS(G_ARGUMENT_POINTER+1)
G_ARGUMENT_POINTER = G_ARGUMENT_POINTER + 1
IF (N.GT.1 .AND. N.NE.N2) then
call mat_err(11)
return
endif
call mat_stack1(QUOTE)
IF (G_ERR .GT. 0) return
LL = L2+M2*N2
call mat_wcopy(M*N,GM_REALS(location),GM_IMAGS(location),1,GM_REALS(LL),GM_IMAGS(LL),1)
call mat_wcopy(M*N+M2*N2,GM_REALS(L2),GM_IMAGS(L2),1,GM_REALS(location),GM_IMAGS(location),1)
G_VAR_DATALOC(G_ARGUMENT_POINTER) = location+M2*N2
G_VAR_ROWS(G_ARGUMENT_POINTER) = M
G_VAR_COLS(G_ARGUMENT_POINTER) = N
call mat_stack1(QUOTE)
IF (G_ERR .GT. 0) return
G_ARGUMENT_POINTER = G_ARGUMENT_POINTER - 1
M = N2
N = M2
goto 20
!
! RECTANGULAR MATRIX LEFT DIVISION A BACKSLASH A2
!
20 continue
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*N2 .GT. 1) goto 21
M2 = M
N2 = M
if(too_much_memory( L2+M*M - G_VAR_DATALOC(G_TOP_OF_SAVED) ) )return
call mat_wset(M*M-1,0.0D0,0.0D0,GM_REALS(L2+1),GM_IMAGS(L2+1),1)
call mat_wcopy(M,GM_REALS(L2),GM_IMAGS(L2),0,GM_REALS(L2),GM_IMAGS(L2),M+1)
21 continue
IF (M2 .NE. M) then
call mat_err(12)
return
endif
L3 = L2 + MAX(M,N)*N2
L4 = L3 + N
if(too_much_memory( L4 + N - G_VAR_DATALOC(G_TOP_OF_SAVED) ) )return
IF (M .GT. N) goto 23
DO JB = 1, N2
J = N+1-JB
LS = L2 + (J-1)*M
LL = L2 + (J-1)*N
call mat_wcopy(M,GM_REALS(LS),GM_IMAGS(LS),-1,GM_REALS(LL),GM_IMAGS(LL),-1)
enddo
23 continue
DO J = 1, N
G_BUF(J) = 0
enddo
call ML_WQRDC(GM_REALS(location),GM_IMAGS(location), &
& M,M,N, &
& GM_REALS(L4),GM_IMAGS(L4), &
& G_BUF, &
& GM_REALS(L3),GM_IMAGS(L3), &
& 1)
K = 0
EPS = GM_REALS(GM_BIGMEM-4)
T(1) = DABS(GM_REALS(location))+DABS(GM_IMAGS(location))
TOL = mat_flop(dble(MAX(M,N))*EPS*T(1))
MN = MIN(M,N)
DO J = 1, MN
LS = location+J-1+(J-1)*M
T(1) = DABS(GM_REALS(LS)) + DABS(GM_IMAGS(LS))
IF (T(1) .GT. TOL) K = J
enddo
IF (K .LT. MN) then
WRITE(message,'(" RANK DEFICIENT, RANK =",I4,", TOL =",1PD13.4)') K,TOL
call journal(message)
endif
MN = MAX(M,N)
DO J = 1, N2
LS = L2+(J-1)*MN
call ML_WQRSL(GM_REALS(location),GM_IMAGS(location), &
& M,M,K, &
& GM_REALS(L4),GM_IMAGS(L4), &
& GM_REALS(LS),GM_IMAGS(LS), &
& T,T, &
& GM_REALS(LS),GM_IMAGS(LS), &
& GM_REALS(LS),GM_IMAGS(LS), &
& T,T,T,T,100,INFO)
LL = LS+K
call mat_wset(N-K,0.0D0,0.0D0,GM_REALS(LL),GM_IMAGS(LL),1)
enddo
DO J = 1, N
G_BUF(J) = -G_BUF(J)
enddo
DO J = 1, N
IF (G_BUF(J) .GT. 0) cycle
K = -G_BUF(J)
G_BUF(J) = K
33 CONTINUE
IF (K .EQ. J) cycle
LS = L2+J-1
LL = L2+K-1
call mat_wswap(N2,GM_REALS(LS),GM_IMAGS(LS),MN,GM_REALS(LL),GM_IMAGS(LL),MN)
G_BUF(K) = -G_BUF(K)
K = G_BUF(K)
goto 33
enddo
DO J = 1, N2
LS = L2+(J-1)*MN
LL = location+(J-1)*N
call mat_wcopy(N,GM_REALS(LS),GM_IMAGS(LS),1,GM_REALS(LL),GM_IMAGS(LL),1)
enddo
G_VAR_ROWS(G_ARGUMENT_POINTER) = N
G_VAR_COLS(G_ARGUMENT_POINTER) = N2
IF (G_FIN .EQ. -1) call mat_stack1(QUOTE)
IF (G_ERR .GT. 0) return
goto 99
!===================================================================================================================================
! QR
!
40 continue
mm = max(m,n)
ls = location + mm*mm
if (G_LHS.eq.1 .and. G_FIN.eq.1) ls = location
le = ls + m*n
l4 = le + mm
if(too_much_memory( l4+mm - G_VAR_DATALOC(G_TOP_OF_SAVED) ) )return
if (ls.ne.location) then
call mat_wcopy(m*n,GM_REALS(location),GM_IMAGS(location),1,GM_REALS(ls),GM_IMAGS(ls),1)
endif
job = 1
if (G_LHS.lt.3) job = 0
do j = 1, n
G_BUF(j) = 0
enddo
call ml_wqrdc(GM_REALS(ls),GM_IMAGS(ls), &
& m,m,n, &
& GM_REALS(l4),GM_IMAGS(l4), &
& G_BUF, &
& GM_REALS(le),GM_IMAGS(le), &
& job)
if (G_LHS.eq.1 .and. G_FIN.eq.1) goto 99
call mat_wset(m*m,0.0d0,0.0d0,GM_REALS(location),GM_IMAGS(location),1)
call mat_wset(m,1.0d0,0.0d0,GM_REALS(location),GM_IMAGS(location),m+1)
do j = 1, m
ll = location+(j-1)*m
call ml_wqrsl(GM_REALS(ls),GM_IMAGS(ls),m,m,n,GM_REALS(l4),GM_IMAGS(l4), &
& GM_REALS(ll),GM_IMAGS(ll),GM_REALS(ll),GM_IMAGS(ll),t,t, &
& t,t,t,t,t,t,10000,info)
enddo
if (G_FIN .eq. 2) goto 99
G_VAR_COLS(G_ARGUMENT_POINTER) = M
do j = 1, n
ll = ls+j+(j-1)*m
call mat_wset(m-j,0.0d0,0.0d0,GM_REALS(ll),GM_IMAGS(ll),1)
enddo
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) = ls
G_VAR_ROWS(G_ARGUMENT_POINTER) = m
G_VAR_COLS(G_ARGUMENT_POINTER) = n
if (G_LHS .eq. 2) goto 99
call mat_wset(N*N,0.0D0,0.0D0,GM_REALS(le),GM_IMAGS(le),1)
do j = 1, n
ll = le+G_BUF(j)-1+(j-1)*n
GM_REALS(ll) = 1.0d0
enddo
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) = le
G_VAR_ROWS(G_ARGUMENT_POINTER) = n
G_VAR_COLS(G_ARGUMENT_POINTER) = n
goto 99
!===================================================================================================================================
!
99 continue
END SUBROUTINE mat_matfn4