Program zuncsd_example
! ZUNCSD Example Program Text
! Copyright 2017, Numerical Algorithms Group Ltd. http://www.nag.com
! .. Use Statements ..
Use blas_interfaces, Only: zgemm
Use lapack_example_aux, Only: nagf_file_print_matrix_complex_gen_comp, &
nagf_file_print_matrix_real_gen
Use lapack_interfaces, Only: zuncsd
Use lapack_precision, Only: dp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Complex (Kind=dp), Parameter :: one = (1.0_dp, 0.0_dp)
Complex (Kind=dp), Parameter :: zero = (0.0_dp, 0.0_dp)
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=dp), Allocatable :: u(:, :), u1(:, :), u2(:, :), v(:, :), &
v1t(:, :), v2t(:, :), w(:, :), work(:), x(:, :), x11(:, :), x12(:, :), &
x21(:, :), x22(:, :)
Complex (Kind=dp) :: wdum(1)
Real (Kind=dp) :: rwdum(1)
Real (Kind=dp), 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, *) 'ZUNCSD 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 nagf_file_print_matrix_complex_gen_comp.
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call nagf_file_print_matrix_complex_gen_comp('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.
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, 110) 'Failure in ZUNCSD. info =', info
Go To 100
End If
! Print Theta using real matrix printing routine nagf_file_print_matrix_real_gen
! 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, 120) 'Theta Component of CS factorization of X:'
Flush (nout)
ifail = 0
Call nagf_file_print_matrix_real_gen('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, 110) 'Failure in ZUNCSD. info =', info
Go To 100
End If
! Reprint Theta using matrix printing routine nagf_file_print_matrix_real_gen.
If (reprint/=0) Then
Write (nout, 120) 'Components of CS factorization of X:'
Flush (nout)
ifail = 0
Call nagf_file_print_matrix_real_gen('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_dp, kind=dp)
x(j+n11, j+n11+r+n21+n22) = cmplx(-sin(theta(j)), 0.0_dp, kind=dp)
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_dp, kind=dp)
x(p+n22+j, j+r+n21+n22) = cmplx(cos(theta(j)), 0.0_dp, kind=dp)
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 nagf_file_print_matrix_complex_gen_comp.
Write (nout, *)
Flush (nout)
ifail = 0
Call nagf_file_print_matrix_complex_gen_comp('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
110 Format (1X, A, I4)
120 Format (/, 1X, A, /)
End Program
