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)
261 bsmall=floordiv(n,nblocks)
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