41 lines
1.1 KiB
Forth
41 lines
1.1 KiB
Forth
program qrtask
|
|
implicit none
|
|
integer, parameter :: n = 4
|
|
real(8) :: durer(n, n) = reshape([16d0, 5d0, 9d0, 4d0, &
|
|
3d0, 10d0, 6d0, 15d0, &
|
|
2d0, 11d0, 7d0, 14d0, &
|
|
13d0, 8d0, 12d0, 1d0], [n, n])
|
|
real(8) :: q(n, n), r(n, n), qr(n, n), id(n, n), tau(n)
|
|
integer, parameter :: lwork = 1024
|
|
real(8) :: work(lwork)
|
|
integer :: info, i, j
|
|
|
|
q = durer
|
|
call dgeqrf(n, n, q, n, tau, work, lwork, info)
|
|
|
|
r = 0d0
|
|
forall (i = 1:n, j = 1:n, j >= i) r(i, j) = q(i, j)
|
|
|
|
call dorgqr(n, n, n, q, n, tau, work, lwork, info)
|
|
|
|
qr = matmul(q, r)
|
|
id = matmul(q, transpose(q))
|
|
|
|
call show(4, durer, "A")
|
|
call show(4, q, "Q")
|
|
call show(4, r, "R")
|
|
call show(4, qr, "Q*R")
|
|
call show(4, id, "Q*Q'")
|
|
contains
|
|
subroutine show(n, a, s)
|
|
character(*) :: s
|
|
integer :: n, i
|
|
real(8) :: a(n, n)
|
|
|
|
print *, s
|
|
do i = 1, n
|
|
print 1, a(i, :)
|
|
1 format (*(F12.6,:,' '))
|
|
end do
|
|
end subroutine
|
|
end program
|