NAG Library Manual, Mark 30
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
    Program f08rnfe

!     F08RNF Example Program Text

!     Mark 30.0 Release. NAG Copyright 2024.

!     .. 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