-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathblock_diag.f90
139 lines (130 loc) · 5.6 KB
/
block_diag.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
134
135
136
137
138
139
module block_diag
!
! Simplistic block-diagonalization routines
!
! Note that this set of routines is not complete: only the data types
! We encountered so far are implemented. Sorry.
!
! 2013 Nov 02 - Added block_geev routine
! 2013 Jul 16 - Added sorting of eigenvalues / singular values
! 2014 Jun ?? - Added quad-precision
!
use accuracy
use timer
use lapack
use printing
use sort_tools
implicit none
private
public block_syev, block_gesvd, block_geev
public rcsid_block_diag
!
character(len=clen), save :: rcsid_block_diag = "$Id: block_diag.f90,v 1.8 2021/09/29 13:43:22 ps Exp $"
!
interface block_syev
module procedure block_syev_real
!*qd module procedure block_syev_quad
end interface block_syev
!
interface block_geev
module procedure block_zgeev_real
!*qd module procedure block_zgeev_quad
end interface block_geev
!
interface block_gesvd
module procedure block_gesvd_real
!*qd module procedure block_gesvd_quad
module procedure block_gesvd_complex
!*qd module procedure block_gesvd_quad_complex
end interface block_gesvd
!
integer(ik), parameter :: verbose = 0
contains
subroutine block_syev_real(h,e,eps)
integer, parameter :: ok = rk
include 'block_diag_syev_common.f90'
end subroutine block_syev_real
subroutine block_syev_quad(h,e,eps)
integer, parameter :: ok = xrk
include 'block_diag_syev_common.f90'
end subroutine block_syev_quad
!
! This routine is modeled on lapack_zgeev2
!
subroutine block_zgeev_real(h,e,eps)
integer, parameter :: ok = rk
include "block_diag_zgeev_common.f90"
end subroutine block_zgeev_real
subroutine block_zgeev_quad(h,e,eps)
integer, parameter :: ok = xrk
include "block_diag_zgeev_common.f90"
end subroutine block_zgeev_quad
subroutine block_gesvd_real(a,s,u,vth,eps)
integer, parameter :: ok = rk
! Only the variables where the basic type (not the kind) depend on the
! specific routine are listed below; the remaining variables are defined
! in the include file
real(ok), intent(inout) :: a (:,:) ! In: Matrix to be decomposed
! Out: Matrix destroyed
real(ok), intent(out) :: u (:,:) ! Out: Left singular vectors
real(ok), intent(out) :: vth(:,:) ! Out: Right singular vectors, transposed & conjugated
! The overall result is A = U S VTH
!
real(ok), allocatable :: tmp_a(:,:) ! Temporary used for decomposing the connected block
real(ok), allocatable :: tmp_u(:,:) ! ditto
real(ok), allocatable :: tmp_vth(:,:) ! ditto
!
include 'block_diag_gesvd_common.f90'
end subroutine block_gesvd_real
subroutine block_gesvd_quad(a,s,u,vth,eps)
integer, parameter :: ok = xrk
! Only the variables where the basic type (not the kind) depend on the
! specific routine are listed below; the remaining variables are defined
! in the include file
real(ok), intent(inout) :: a (:,:) ! In: Matrix to be decomposed
! Out: Matrix destroyed
real(ok), intent(out) :: u (:,:) ! Out: Left singular vectors
real(ok), intent(out) :: vth(:,:) ! Out: Right singular vectors, transposed & conjugated
! The overall result is A = U S VTH
!
real(ok), allocatable :: tmp_a(:,:) ! Temporary used for decomposing the connected block
real(ok), allocatable :: tmp_u(:,:) ! ditto
real(ok), allocatable :: tmp_vth(:,:) ! ditto
!
include 'block_diag_gesvd_common.f90'
end subroutine block_gesvd_quad
subroutine block_gesvd_complex(a,s,u,vth,eps)
integer, parameter :: ok = rk
! Only the variables where the basic type (not the kind) depend on the
! specific routine are listed below; the remaining variables are defined
! in the include file
complex(ok), intent(inout) :: a (:,:) ! In: Matrix to be decomposed
! Out: Matrix destroyed
complex(ok), intent(out) :: u (:,:) ! Out: Left singular vectors
complex(ok), intent(out) :: vth(:,:) ! Out: Right singular vectors, transposed & conjugated
! The overall result is A = U S VTH
!
complex(ok), allocatable :: tmp_a(:,:) ! Temporary used for decomposing the connected block
complex(ok), allocatable :: tmp_u(:,:) ! ditto
complex(ok), allocatable :: tmp_vth(:,:) ! ditto
!
include 'block_diag_gesvd_common.f90'
end subroutine block_gesvd_complex
subroutine block_gesvd_quad_complex(a,s,u,vth,eps)
integer, parameter :: ok = xrk
! Only the variables where the basic type (not the kind) depend on the
! specific routine are listed below; the remaining variables are defined
! in the include file
complex(ok), intent(inout) :: a (:,:) ! In: Matrix to be decomposed
! Out: Matrix destroyed
complex(ok), intent(out) :: u (:,:) ! Out: Left singular vectors
complex(ok), intent(out) :: vth(:,:) ! Out: Right singular vectors, transposed & conjugated
! The overall result is A = U S VTH
!
complex(ok), allocatable :: tmp_a(:,:) ! Temporary used for decomposing the connected block
complex(ok), allocatable :: tmp_u(:,:) ! ditto
complex(ok), allocatable :: tmp_vth(:,:) ! ditto
!
include 'block_diag_gesvd_common.f90'
end subroutine block_gesvd_quad_complex
end module block_diag