This repository has been archived by the owner on Oct 3, 2018. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
/
linalg.f90
133 lines (99 loc) · 3.89 KB
/
linalg.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
module linalg
use types
use matrix_routines
use utils
implicit none
private
public :: qdr, rdq
public :: stable_mm_qdr, stable_mm_rdq
contains
subroutine qdr(a, q, d, r, info) !{{{
real(dp), dimension(:,:), intent(in) :: a
real(dp), dimension(:,:), allocatable, intent(out) :: q,d,r
integer, intent(out) :: info
real(dp), dimension(:,:), allocatable :: tr
real(dp), dimension(:), allocatable :: vd,vdi,tau, work, work_q
integer, dimension(:), allocatable :: sh
integer :: lwork, m, ws
integer :: i,s
q = a
sh = shape(a)
m = sh(1)
if (m /= sh(2)) then
call stop_error("qdr: non-square matrix detected.")
end if
allocate(tau(m))
allocate(work_q(1))
! Compute the optimal work array: query
call dgeqrf(m, m, q, m, tau, work_q, -1, info)
ws = int(abs(work_q(1)))
allocate(work(ws))
! Compute the QR factorization
call dgeqrf(m, m, q, m, tau, work, ws, info)
! R = D R'
tr = triangular_upper(q)
vd = diagonal(tr)
s = size(vd)
allocate(vdi(s))
do i=1,s
vdi(i) = 1/vd(i)
end do
r = tr
call dgemm('n','n',m,m,m,1.0_dp,diagonal(vdi),m,tr,m,0.0_dp,r,m)
d = diagonal(vd)
! Compute the orthogonal matrix from the reflectors resulting from dgeqrf
call dorgqr(m, m, m, q, m, tau, work_q, -1, info)
! Reallocate the work array if the optimal size is different from the QR factorization above.
if (work_q(1) .ne. ws) then
ws = work_q(1)
deallocate(work)
allocate(work(ws))
end if
call dorgqr(m, m, m, q, m, tau, work, ws, info)
end subroutine qdr !}}}
subroutine rdq(a, r, d, q, info) !{{{
! RDQ decomposition of the matrix A. Assumes that A is quadratic!
real(dp), dimension(:,:), intent(in) :: a
real(dp), dimension(:,:), allocatable, intent(out) :: r, d, q
integer, intent(out) :: info
real(dp), dimension(:,:), allocatable :: tr
real(dp), dimension(:), allocatable :: vd,vdi,tau, work, work_q
integer, dimension(:), allocatable :: sh
integer :: lwork, m, n, mn, ws
integer :: i,s
q = a
sh = shape(a)
m = sh(1)
n = sh(2)
if (m /= n) then
call stop_error("rdq: non-square matrix detected.")
end if
allocate(tau(m))
allocate(work_q(1))
! Compute the optimal work array: query
call dgerqf(m, m, q, m, tau, work_q, -1, info)
ws = int(abs(work_q(1)))
allocate(work(ws))
! Compute the RQ factorization
call dgerqf(m, m, q, m, tau, work, ws, info)
tr = triangular_upper(q)
vd = diagonal(tr)
s = size(vd)
allocate(vdi(s))
do i=1,s
vdi(i) = 1/vd(i)
end do
r = tr
call dgemm('n','n',m,m,m,1.0_dp,tr,m,diagonal(vdi),m,0.0_dp,r,m)
d = diagonal(vd)
! Compute the orthogonal matrix from the reflectors resulting from dgerqf
call dorgrq(m, m, m, q, m, tau, work_q, -1, info)
! Reallocate the work array if the optimal size is different from the RQ factorization above.
if (work_q(1) .ne. ws) then
ws = work_q(1)
deallocate(work)
allocate(work(ws))
end if
call dorgrq(m, m, m, q, m, tau, work, ws, info)
end subroutine rdq !}}}
end module linalg