Program f08rnfe
! F08RNF Example Program Text
! Mark 28.3 Release. NAG Copyright 2022.
! .. Use Statements ..
Use nag_library, Only: nag_wp, x04caf, x04dbf, zgemm, zuncsd
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Complex (Kind=nag_wp), Parameter :: one = (1.0_nag_wp,0.0_nag_wp)
Complex (Kind=nag_wp), Parameter :: zero = (0.0_nag_wp,0.0_nag_wp)
Integer, Parameter :: nin = 5, nout = 6, recombine = 1, &
reprint = 1
! .. Local Scalars ..
Integer :: i, ifail, info, info2, j, ldu, ldu1, &
ldu2, ldv, ldv1t, ldv2t, ldx, ldx11, &
ldx12, ldx21, ldx22, lrwork, lwork, &
m, n11, n12, n21, n22, p, q, r
! .. Local Arrays ..
Complex (Kind=nag_wp), Allocatable :: u(:,:), u1(:,:), u2(:,:), v(:,:), &
v1t(:,:), v2t(:,:), w(:,:), work(:), &
x(:,:), x11(:,:), x12(:,:), &
x21(:,:), x22(:,:)
Complex (Kind=nag_wp) :: wdum(1)
Real (Kind=nag_wp) :: rwdum(1)
Real (Kind=nag_wp), Allocatable :: rwork(:), theta(:)
Integer, Allocatable :: iwork(:)
Character (1) :: clabs(1), rlabs(1)
! .. Intrinsic Procedures ..
Intrinsic :: cmplx, cos, min, nint, real, sin
! .. Executable Statements ..
Write (nout,*) 'F08RNF Example Program Results'
Write (nout,*)
Flush (nout)
! Skip heading in data file
Read (nin,*)
Read (nin,*) m, p, q
r = min(min(p,q),min(m-p,m-q))
ldx = m
ldx11 = p
ldx12 = p
ldx21 = m - p
ldx22 = m - p
ldu = m
ldu1 = p
ldu2 = m - p
ldv = m
ldv1t = q
ldv2t = m - q
Allocate (x(ldx,m),u(ldu,m),v(ldv,m),theta(r),iwork(m),w(ldx,m))
Allocate (x11(ldx11,q),x12(ldx12,m-q),x21(ldx21,q),x22(ldx22,m-q))
Allocate (u1(ldu1,p),u2(ldu2,m-p),v1t(ldv1t,q),v2t(ldv2t,m-q))
! Read (by column) and print unitary X from data file
! (as, say, generated by a generalized singular value decomposition).
Do i = 1, m
Read (nin,*) x(1:m,i)
End Do
! Print general complex matrix using matrix printing routine x04dbf.
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call x04dbf('General','N',m,m,x,ldx,'Bracketed','F7.4', &
' Unitary matrix X','Integer',rlabs,'Integer',clabs,80,0,ifail)
Write (nout,*)
! Compute the complete CS factorization of X:
! X11 is stored in X(1:p, 1:q), X12 is stored in X(1:p, q+1:m)
! X21 is stored in X(p+1:m, 1:q), X22 is stored in X(p+1:m, q+1:m)
! U1 is stored in U(1:p, 1:p), U2 is stored in U(p+1:m, p+1:m)
! V1 is stored in V(1:q, 1:q), V2 is stored in V(q+1:m, q+1:m)
x11(1:p,1:q) = x(1:p,1:q)
x12(1:p,1:m-q) = x(1:p,q+1:m)
x21(1:m-p,1:q) = x(p+1:m,1:q)
x22(1:m-p,1:m-q) = x(p+1:m,q+1:m)
! Workspace query first to get length of work array for complete CS
! factorization routine zuncsd/f08rnf.
lwork = -1
lrwork = -1
Call zuncsd('Yes U1','Yes U2','Yes V1T','Yes V2T','Column','Default',m, &
p,q,x,ldx,x(1,q+1),ldx,x(p+1,1),ldx,x(p+1,q+1),ldx,theta,u,ldu, &
u(p+1,p+1),ldu,v,ldv,v(q+1,q+1),ldv,wdum,lwork,rwdum,lrwork,iwork, &
info)
lwork = nint(real(wdum(1)))
lrwork = nint(rwdum(1))
Allocate (work(lwork),rwork(lrwork))
! Initialize all of u, v to zero.
u(1:m,1:m) = zero
v(1:m,1:m) = zero
! This is how you might pass partitions as sub-matrices
Call zuncsd('Yes U1','Yes U2','Yes V1T','Yes V2T','Column','Default',m, &
p,q,x,ldx,x(1,q+1),ldx,x(p+1,1),ldx,x(p+1,q+1),ldx,theta,u,ldu, &
u(p+1,p+1),ldu,v,ldv,v(q+1,q+1),ldv,work,lwork,rwork,lrwork,iwork, &
info)
If (info/=0) Then
Write (nout,99999) 'Failure in ZUNCSD/F08RNF. info =', info
Go To 100
End If
! Print Theta using real matrix printing routine x04caf
! Note: U1, U2, V1T, V2T not printed since these may differ by a sign
! change in columns of U1, U2 and corresponding rows of V1T, V2T.
Write (nout,99998) 'Theta Component of CS factorization of X:'
Flush (nout)
ifail = 0
Call x04caf('G','N',r,1,theta,r,' Theta',ifail)
Write (nout,*)
Flush (nout)
! And this is how you might pass partitions as separate matrices.
Call zuncsd('Yes U1','Yes U2','Yes V1T','Yes V2T','Column','Default',m, &
p,q,x11,ldx11,x12,ldx12,x21,ldx21,x22,ldx22,theta,u1,ldu1,u2,ldu2,v1t, &
ldv1t,v2t,ldv2t,work,lwork,rwork,lrwork,iwork,info2)
If (info/=0) Then
Write (nout,99999) 'Failure in ZUNCSD/F08RNF. info =', info
Go To 100
End If
! Reprint Theta using matrix printing routine x04caf.
If (reprint/=0) Then
Write (nout,99998) 'Components of CS factorization of X:'
Flush (nout)
ifail = 0
Call x04caf('G','N',r,1,theta,r,' Theta',ifail)
Write (nout,*)
Flush (nout)
End If
If (recombine/=0) Then
! Recombining should return the original matrix
! Assemble Sigma_p into X
x(1:m,1:m) = zero
n11 = min(p,q) - r
n12 = min(p,m-q) - r
n21 = min(m-p,q) - r
n22 = min(m-p,m-q) - r
! Top Half
Do j = 1, n11
x(j,j) = one
End Do
Do j = 1, r
x(j+n11,j+n11) = cmplx(cos(theta(j)),0.0_nag_wp,kind=nag_wp)
x(j+n11,j+n11+r+n21+n22) = cmplx(-sin(theta(j)),0.0_nag_wp, &
kind=nag_wp)
End Do
Do j = 1, n12
x(j+n11+r,j+n11+r+n21+n22+r) = -one
End Do
! Bottom half
Do j = 1, n22
x(p+j,q+j) = one
End Do
Do j = 1, r
x(p+n22+j,j+n11) = cmplx(sin(theta(j)),0.0_nag_wp,kind=nag_wp)
x(p+n22+j,j+r+n21+n22) = cmplx(cos(theta(j)),0.0_nag_wp,kind=nag_wp)
End Do
Do j = 1, n21
x(p+n22+r+j,n11+r+j) = one
End Do
! multiply U * Sigma_p into w
Call zgemm('n','n',m,m,m,one,u,ldu,x,ldx,zero,w,ldx)
! form U * Sigma_p * V^T into u
Call zgemm('n','n',m,m,m,one,w,ldx,v,ldv,zero,u,ldu)
! Print recombined matrix using complex matrix printing routine x04dbf.
Write (nout,*)
Flush (nout)
ifail = 0
Call x04dbf('General','N',m,m,u,ldu,'Bracketed','F7.4', &
' Recombined matrix X = U * Sigma_p * V^H','Integer',rlabs, &
'Integer',clabs,80,0,ifail)
End If
100 Continue
99999 Format (1X,A,I4)
99998 Format (/,1X,A,/)
End Program f08rnfe