Program dorcsd_example
! DORCSD Example Program Text
! Copyright 2017, Numerical Algorithms Group Ltd. http://www.nag.com
! .. Use Statements ..
Use blas_interfaces, Only: dgemm
Use lapack_example_aux, Only: nagf_file_print_matrix_real_gen
Use lapack_interfaces, Only: dorcsd
Use lapack_precision, Only: dp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Real (Kind=dp), Parameter :: one = 1.0_dp
Real (Kind=dp), Parameter :: zero = 0.0_dp
Integer, Parameter :: nin = 5, nout = 6, recombine = 1, reprint = 0
! .. Local Scalars ..
Integer :: i, ifail, info, info2, j, ldu, ldu1, ldu2, ldv, ldv1t, ldv2t, &
ldx, ldx11, ldx12, ldx21, ldx22, lwork, m, n11, n12, n21, n22, p, q, r
! .. Local Arrays ..
Real (Kind=dp), Allocatable :: theta(:), u(:, :), u1(:, :), u2(:, :), &
v(:, :), v1t(:, :), v2t(:, :), w(:, :), work(:), x(:, :), x11(:, :), &
x12(:, :), x21(:, :), x22(:, :)
Real (Kind=dp) :: wdum(1)
Integer, Allocatable :: iwork(:)
! .. Intrinsic Procedures ..
Intrinsic :: cos, min, nint, sin
! .. Executable Statements ..
Write (nout, *) 'DORCSD 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 and print orthogonal X from data file
! (as, say, generated by a generalized singular value decomposition).
Read (nin, *)(x(i,1:m), i=1, m)
! Print general matrix using simple matrix printing routine nagf_file_print_matrix_real_gen.
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call nagf_file_print_matrix_real_gen('G', 'N', m, m, x, ldx, &
' Orthogonal matrix X', 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 dorcsd.
lwork = -1
Call dorcsd('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, iwork, info)
lwork = nint(wdum(1))
Allocate (work(lwork))
! 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 dorcsd('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, iwork, info)
If (info/=0) Then
Write (nout, 110) 'Failure in DORCSD. info =', info
Go To 100
End If
! Print Theta, U1, U2, V1T, V2T using matrix printing routine nagf_file_print_matrix_real_gen.
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)
! By changes of sign the first r elements of the first row of U1 can be
! made positive.
Do i = 1, r
If (u(1,i)<0.0_dp) Then
u(1:p, i) = -u(1:p, i)
u(p+1:m, p+i) = -u(p+1:m, p+i)
v(i, 1:q) = -v(i, 1:q)
v(q+i, q+1:m) = -v(q+i, q+1:m)
End If
End Do
ifail = 0
Call nagf_file_print_matrix_real_gen('G', 'N', p, p, u, ldu, ' U1', &
ifail)
Write (nout, *)
Flush (nout)
ifail = 0
Call nagf_file_print_matrix_real_gen('G', 'N', m-p, m-p, u(p+1,p+1), &
ldu, ' U2', ifail)
Write (nout, *)
Flush (nout)
ifail = 0
Call nagf_file_print_matrix_real_gen('G', 'N', q, q, v, ldv, ' V1T', &
ifail)
Write (nout, *)
Flush (nout)
ifail = 0
Call nagf_file_print_matrix_real_gen('G', 'N', m-q, m-q, v(q+1,q+1), &
ldv, ' V2T', ifail)
Write (nout, *)
! And this is how you might pass partitions as separate matrices.
Call dorcsd('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, iwork, &
info2)
If (info2/=0) Then
Write (nout, 110) 'Failure in DORCSD. info =', info
Go To 100
End If
! Print Theta, U1, U2, V1T, V2T using matrix printing routine nagf_file_print_matrix_real_gen.
If (reprint/=0) Then
! By changes of sign the first r elements of the first row of U1 can be
! made positive.
Do i = 1, r
If (u1(1,i)<0.0_dp) Then
u1(1:p, i) = -u1(1:p, i)
u2(1:m-p, i) = -u2(1:m-p, i)
v1t(i, 1:q) = -v1t(i, 1:q)
v2t(i, 1:m-q) = -v2t(i, 1:m-q)
End If
End Do
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)
ifail = 0
Call nagf_file_print_matrix_real_gen('G', 'N', p, p, u1, ldu1, &
' U1', ifail)
Write (nout, *)
Flush (nout)
ifail = 0
Call nagf_file_print_matrix_real_gen('G', 'N', m-p, m-p, u2, ldu2, &
' U2', ifail)
Write (nout, *)
Flush (nout)
ifail = 0
Call nagf_file_print_matrix_real_gen('G', 'N', q, q, v1t, ldv1t, &
' V1T', ifail)
Write (nout, *)
Flush (nout)
ifail = 0
Call nagf_file_print_matrix_real_gen('G', 'N', m-q, m-q, v2t, ldv2t, &
' V2T', 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-1) - r
! Top Half
Do j = 1, n11
x(j, j) = one
End Do
Do j = 1, r
x(j+n11, j+n11) = cos(theta(j))
x(j+n11, j+n11+r+n21+n22) = -sin(theta(j))
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) = sin(theta(j))
x(p+n22+j, j+r+n21+n22) = cos(theta(j))
End Do
Do j = 1, n21
x(p+n22+r+j, n11+r+j) = one
End Do
! multiply U * Sigma_p into w
Call dgemm('n', 'n', m, m, m, one, u, ldu, x, ldx, zero, w, ldx)
! form U * Sigma_p * V^T into u
Call dgemm('n', 'n', m, m, m, one, w, ldx, v, ldv, zero, u, ldu)
Flush (nout)
ifail = 0
Call nagf_file_print_matrix_real_gen('G', 'N', m, m, u, ldu, &
' Recombined matrix X = U * Sigma_p * V^T', ifail)
End If
100 Continue
110 Format (1X, A, I4)
120 Format (1X, A, /)
End Program
