27 REAL(KIND=8),
DIMENSION(:),
ALLOCATABLE ::
work 32 REAL(KIND=8),
PARAMETER ::
real_eps = 2.2204460492503131e-16
45 IF (
ALLOCATED(
work))
THEN 47 IF (allocstat /= 0) stop
"*** Deallocation problem ! ***" 50 IF (allocstat /= 0) stop
"*** Not enough memory ! ***" 62 SUBROUTINE sqrtm(A,sqA,info,info_triu,bs)
63 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: A
64 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: sqA
65 INTEGER,
INTENT(IN),
OPTIONAL :: bs
66 INTEGER,
INTENT(OUT) :: info,info_triu
67 REAL(KIND=8),
DIMENSION(SIZE(A,1),SIZE(A,1)) :: T,Z,R
68 COMPLEX(KIND=16),
DIMENSION(SIZE(A,1),SIZE(A,1)) :: Tz,Zz,Rz
69 REAL(KIND=8),
DIMENSION(SIZE(A,1)) :: wr,wi
70 LOGICAL,
DIMENSION(SIZE(A,1)) :: bwork
77 CALL dgees(
'v',
'n',selectev,n,t,n,sdim,wr,wi,z,n,
work,
lwork,bwork,info)
104 rz=matmul(zz,matmul(rz,conjg(transpose(zz))))
111 IF (
PRESENT(bs))
THEN 116 sqa=matmul(z,matmul(r,transpose(z)))
131 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: A
132 INTEGER,
INTENT(IN),
OPTIONAL :: bs
133 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: sqA
134 INTEGER,
INTENT(OUT) :: info
135 REAL(KIND=8),
DIMENSION(SIZE(A,1)) :: A_diag
136 REAL(KIND=8),
DIMENSION(SIZE(A,1),SIZE(A,1)) :: R,Sm,Rii,Rjj
137 INTEGER,
DIMENSION(2*SIZE(A,1),2) :: start_stop_pairs
138 REAL(KIND=8) :: s,denom,scale
139 INTEGER :: i,j,k,start,n,sstop,m
140 INTEGER :: istart,istop,jstart,jstop
141 INTEGER :: nblocks,blocksize
142 INTEGER :: bsmall,blarge,nlarge,nsmall
145 IF (
PRESENT(bs)) blocksize=bs
152 r(i,i)=sqrt(a_diag(i))
156 nblocks=max(
floordiv(n,blocksize),1)
158 nlarge=mod(n,nblocks)
160 nsmall=nblocks-nlarge
161 IF (nsmall*bsmall + nlarge*blarge /= n) stop
'Sqrtm: Internal inconsistency' 167 start_stop_pairs(i,1)=start
168 start_stop_pairs(i,2)=start+bsmall-1
171 DO i=nsmall+1,nsmall+nlarge
172 start_stop_pairs(i,1)=start
173 start_stop_pairs(i,2)=start+blarge-1
183 start=start_stop_pairs(k,1)
184 sstop=start_stop_pairs(k,2)
188 IF (j-i>1) s= dot_product(r(i,i+1:j-1),r(i+1:j-1,j))
190 IF (denom==0.d0) stop
'Sqrtm: Failed to find the matrix square root' 191 r(i,j)=(a(i,j)-s)/denom
200 jstart=start_stop_pairs(j,1)
201 jstop=start_stop_pairs(j,2)
203 istart=start_stop_pairs(i,1)
204 istop=start_stop_pairs(i,2)
206 sm(istart:istop,jstart:jstop)=a(istart:istop,jstart:jstop)
207 IF (j-i>1) sm(istart:istop,jstart:jstop) = sm(istart:istop&
208 &,jstart:jstop) - matmul(r(istart:istop,istop:jstart)&
209 &,r(istop:jstart,jstart:jstop))
211 rii = r(istart:istop, istart:istop)
213 rjj = r(jstart:jstop, jstart:jstop)
225 CALL dtrsyl(
'N',
'N',k,m,n,rii(istart:istop, istart:istop),m&
226 &,rjj(jstart:jstop,jstart:jstop),n,sm(istart:istop&
227 &,jstart:jstop),m,scale,info)
228 r(istart:istop,jstart:jstop)=sm(istart:istop,jstart:jstop)*scale
235 COMPLEX(KIND=16),
DIMENSION(:,:),
INTENT(IN) :: A
236 INTEGER,
INTENT(IN),
OPTIONAL :: bs
237 COMPLEX(KIND=16),
DIMENSION(:,:),
INTENT(OUT) :: sqA
238 INTEGER,
INTENT(OUT) :: info
239 COMPLEX(KIND=16),
DIMENSION(SIZE(A,1)) :: A_diag
240 COMPLEX(KIND=16),
DIMENSION(SIZE(A,1),SIZE(A,1)) :: R,Sm,Rii,Rjj
241 INTEGER,
DIMENSION(2*SIZE(A,1),2) :: start_stop_pairs
242 COMPLEX(KIND=16) :: s,denom,scale
243 INTEGER :: i,j,k,start,n,sstop,m
244 INTEGER :: istart,istop,jstart,jstop
245 INTEGER :: nblocks,blocksize
246 INTEGER :: bsmall,blarge,nlarge,nsmall
249 IF (
PRESENT(bs)) blocksize=bs
256 r(i,i)=sqrt(a_diag(i))
260 nblocks=max(
floordiv(n,blocksize),1)
262 nlarge=mod(n,nblocks)
264 nsmall=nblocks-nlarge
265 IF (nsmall*bsmall + nlarge*blarge /= n) stop
'Sqrtm: Internal inconsistency' 271 start_stop_pairs(i,1)=start
272 start_stop_pairs(i,2)=start+bsmall-1
275 DO i=nsmall+1,nsmall+nlarge
276 start_stop_pairs(i,1)=start
277 start_stop_pairs(i,2)=start+blarge-1
287 start=start_stop_pairs(k,1)
288 sstop=start_stop_pairs(k,2)
292 IF (j-i>1) s= dot_product(r(i,i+1:j-1),r(i+1:j-1,j))
294 IF (denom==0.d0) stop
'Sqrtm: Failed to find the matrix square root' 295 r(i,j)=(a(i,j)-s)/denom
304 jstart=start_stop_pairs(j,1)
305 jstop=start_stop_pairs(j,2)
307 istart=start_stop_pairs(i,1)
308 istop=start_stop_pairs(i,2)
310 sm(istart:istop,jstart:jstop)=a(istart:istop,jstart:jstop)
311 IF (j-i>1) sm(istart:istop,jstart:jstop) = sm(istart:istop&
312 &,jstart:jstop) - matmul(r(istart:istop,istop:jstart)&
313 &,r(istop:jstart,jstart:jstop))
315 rii = r(istart:istop, istart:istop)
317 rjj = r(jstart:jstop, jstart:jstop)
329 CALL ztrsyl(
'N',
'N',k,m,n,rii(istart:istop, istart:istop),m&
330 &,rjj(jstart:jstop,jstart:jstop),n,sm(istart:istop&
331 &,jstart:jstop),m,scale,info)
332 r(istart:istop,jstart:jstop)=sm(istart:istop,jstart:jstop)*scale
339 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: T,Z
340 COMPLEX(KIND=16),
DIMENSION(:,:),
INTENT(OUT) :: Tz,Zz
341 INTEGER,
PARAMETER :: nb=2
342 COMPLEX(KIND=16),
DIMENSION(nb) :: w
344 COMPLEX(KIND=16) :: r,c,s,mu
345 COMPLEX(KIND=16),
DIMENSION(nb,nb) :: G,Gc
355 IF (abs(tz(m,m-1)) >
real_eps*(abs(tz(m-1,m-1)) + abs(tz(m,m))))
THEN 362 c=g(1,1)*g(2,2)-g(1,2)*g(2,1)
363 w(1)=s/2+sqrt(s**2/4-c)
365 r=sqrt(mu*conjg(mu)+tz(m,m-1)*conjg(tz(m,m-1)))
372 gc=conjg(transpose(g))
373 tz(m-1:m,m-1:n)=matmul(g,tz(m-1:m,m-1:n))
374 tz(1:m,m-1:m)=matmul(tz(1:m,m-1:m),gc)
375 zz(:,m-1:m)=matmul(zz(:,m-1:m),gc)
377 tz(m,m-1)=cmplx(0.d0,kind=16)
385 SUBROUTINE chol(A,sqA,info)
386 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: A
387 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: sqA
388 INTEGER,
INTENT(OUT) :: info
391 CALL dpotrf(
'L',
SIZE(sqa,1),sqa,
SIZE(sqa,1),info)
400 SUBROUTINE sqrtm_svd(A,sqA,info,info_triu,bs)
401 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: A
402 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: sqA
403 INTEGER,
INTENT(IN),
OPTIONAL :: bs
404 INTEGER,
INTENT(OUT) :: info,info_triu
405 REAL(KIND=8),
DIMENSION(SIZE(A,1)) :: S
406 REAL(KIND=8),
DIMENSION(SIZE(A,1),SIZE(A,1)) :: Sq,U,VT
411 CALL dgesvd(
'A',
'A',n,n,sqa,n,s,u,n,vt,n,
work,
lwork,info)
416 sqa=matmul(u,matmul(sq,vt))
integer ndim
Number of variables (dimension of the model)
subroutine, public chol(A, sqA, info)
Routine to perform a Cholesky decomposition.
integer function, public floordiv(i, j)
subroutine, public cdiag(A, d)
subroutine, public sqrtm(A, sqA, info, info_triu, bs)
Routine to compute a real square-root of a matrix.
real(kind=8), dimension(:), allocatable work
subroutine, public diag(A, d)
subroutine, public triu(A, T)
subroutine sqrtm_triu(A, sqA, info, bs)
subroutine rsf2csf(T, Z, Tz, Zz)
subroutine, public init_sqrt
logical function selectev(a, b)
The model parameters module.
real(kind=8), parameter real_eps
Utility module with various routine to compute matrix square root.
subroutine csqrtm_triu(A, sqA, info, bs)
subroutine, public sqrtm_svd(A, sqA, info, info_triu, bs)
Routine to compute a real square-root of a matrix via a SVD decomposition.