cunt03.f

Go to the documentation of this file.

*> \brief \b CUNT03
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE CUNT03( RC, MU, MV, N, K, U, LDU, V, LDV, WORK, LWORK,
*                          RWORK, RESULT, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER*( * )    RC
*       INTEGER            INFO, K, LDU, LDV, LWORK, MU, MV, N
*       REAL               RESULT
*       ..
*       .. Array Arguments ..
*       REAL               RWORK( * )
*       COMPLEX            U( LDU, * ), V( LDV, * ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> CUNT03 compares two unitary matrices U and V to see if their
*> corresponding rows or columns span the same spaces.  The rows are
*> checked if RC = 'R', and the columns are checked if RC = 'C'.
*>
*> RESULT is the maximum of
*>
*>    | V*V' - I | / ( MV ulp ), if RC = 'R', or
*>
*>    | V'*V - I | / ( MV ulp ), if RC = 'C',
*>
*> and the maximum over rows (or columns) 1 to K of
*>
*>    | U(i) - S*V(i) |/ ( N ulp )
*>
*> where abs(S) = 1 (chosen to minimize the expression), U(i) is the
*> i-th row (column) of U, and V(i) is the i-th row (column) of V.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] RC
*> \verbatim
*>          RC is CHARACTER*1
*>          If RC = 'R' the rows of U and V are to be compared.
*>          If RC = 'C' the columns of U and V are to be compared.
*> \endverbatim
*>
*> \param[in] MU
*> \verbatim
*>          MU is INTEGER
*>          The number of rows of U if RC = 'R', and the number of
*>          columns if RC = 'C'.  If MU = 0 CUNT03 does nothing.
*>          MU must be at least zero.
*> \endverbatim
*>
*> \param[in] MV
*> \verbatim
*>          MV is INTEGER
*>          The number of rows of V if RC = 'R', and the number of
*>          columns if RC = 'C'.  If MV = 0 CUNT03 does nothing.
*>          MV must be at least zero.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          If RC = 'R', the number of columns in the matrices U and V,
*>          and if RC = 'C', the number of rows in U and V.  If N = 0
*>          CUNT03 does nothing.  N must be at least zero.
*> \endverbatim
*>
*> \param[in] K
*> \verbatim
*>          K is INTEGER
*>          The number of rows or columns of U and V to compare.
*>          0 <= K <= max(MU,MV).
*> \endverbatim
*>
*> \param[in] U
*> \verbatim
*>          U is COMPLEX array, dimension (LDU,N)
*>          The first matrix to compare.  If RC = 'R', U is MU by N, and
*>          if RC = 'C', U is N by MU.
*> \endverbatim
*>
*> \param[in] LDU
*> \verbatim
*>          LDU is INTEGER
*>          The leading dimension of U.  If RC = 'R', LDU >= max(1,MU),
*>          and if RC = 'C', LDU >= max(1,N).
*> \endverbatim
*>
*> \param[in] V
*> \verbatim
*>          V is COMPLEX array, dimension (LDV,N)
*>          The second matrix to compare.  If RC = 'R', V is MV by N, and
*>          if RC = 'C', V is N by MV.
*> \endverbatim
*>
*> \param[in] LDV
*> \verbatim
*>          LDV is INTEGER
*>          The leading dimension of V.  If RC = 'R', LDV >= max(1,MV),
*>          and if RC = 'C', LDV >= max(1,N).
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX array, dimension (LWORK)
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The length of the array WORK.  For best performance, LWORK
*>          should be at least N*N if RC = 'C' or M*M if RC = 'R', but
*>          the tests will be done even if LWORK is 0.
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*>          RWORK is REAL array, dimension (max(MV,N))
*> \endverbatim
*>
*> \param[out] RESULT
*> \verbatim
*>          RESULT is REAL
*>          The value computed by the test described above.  RESULT is
*>          limited to 1/ulp to avoid overflow.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          0  indicates a successful exit
*>          -k indicates the k-th parameter had an illegal value
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup complex_eig
*
*  =====================================================================
      SUBROUTINE cunt03( RC, MU, MV, N, K, U, LDU, V, LDV, WORK, LWORK,
     $                   RWORK, RESULT, INFO )
*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER*( * )    RC
      INTEGER            INFO, K, LDU, LDV, LWORK, MU, MV, N
      REAL               RESULT
*     ..
*     .. Array Arguments ..
      REAL               RWORK( * )
      COMPLEX            U( LDU, * ), V( LDV, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter( zero = 0.0e0, one = 1.0e0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, IRC, J, LMX
      REAL               RES1, RES2, ULP
      COMPLEX            S, SU, SV
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ICAMAX
      REAL               SLAMCH
      EXTERNAL           lsame, icamax, slamch
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, cmplx, max, min, real
*     ..
*     .. External Subroutines ..
      EXTERNAL           cunt01, xerbla
*     ..
*     .. Executable Statements ..
*
*     Check inputs
*
      info = 0
      IF( lsame( rc, 'R' ) ) THEN
         irc = 0
      ELSE IF( lsame( rc, 'C' ) ) THEN
         irc = 1
      ELSE
         irc = -1
      END IF
      IF( irc.EQ.-1 ) THEN
         info = -1
      ELSE IF( mu.LT.0 ) THEN
         info = -2
      ELSE IF( mv.LT.0 ) THEN
         info = -3
      ELSE IF( n.LT.0 ) THEN
         info = -4
      ELSE IF( k.LT.0 .OR. k.GT.max( mu, mv ) ) THEN
         info = -5
      ELSE IF( ( irc.EQ.0 .AND. ldu.LT.max( 1, mu ) ) .OR.
     $         ( irc.EQ.1 .AND. ldu.LT.max( 1, n ) ) ) THEN
         info = -7
      ELSE IF( ( irc.EQ.0 .AND. ldv.LT.max( 1, mv ) ) .OR.
     $         ( irc.EQ.1 .AND. ldv.LT.max( 1, n ) ) ) THEN
         info = -9
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'CUNT03', -info )
         RETURN
      END IF
*
*     Initialize result
*
      result = zero
      IF( mu.EQ.0 .OR. mv.EQ.0 .OR. n.EQ.0 )
     $   RETURN
*
*     Machine constants
*
      ulp = slamch( 'Precision' )
*
      IF( irc.EQ.0 ) THEN
*
*        Compare rows
*
         res1 = zero
         DO 20 i = 1, k
            lmx = icamax( n, u( i, 1 ), ldu )
            IF( v( i, lmx ).EQ.cmplx( zero ) ) THEN
               sv = one
            ELSE
               sv = abs( v( i, lmx ) ) / v( i, lmx )
            END IF
            IF( u( i, lmx ).EQ.cmplx( zero ) ) THEN
               su = one
            ELSE
               su = abs( u( i, lmx ) ) / u( i, lmx )
            END IF
            s = sv / su
            DO 10 j = 1, n
               res1 = max( res1, abs( u( i, j )-s*v( i, j ) ) )
   10       CONTINUE
   20    CONTINUE
         res1 = res1 / ( real( n )*ulp )
*
*        Compute orthogonality of rows of V.
*
         CALL cunt01( 'Rows', mv, n, v, ldv, work, lwork, rwork, res2 )
*
      ELSE
*
*        Compare columns
*
         res1 = zero
         DO 40 i = 1, k
            lmx = icamax( n, u( 1, i ), 1 )
            IF( v( lmx, i ).EQ.cmplx( zero ) ) THEN
               sv = one
            ELSE
               sv = abs( v( lmx, i ) ) / v( lmx, i )
            END IF
            IF( u( lmx, i ).EQ.cmplx( zero ) ) THEN
               su = one
            ELSE
               su = abs( u( lmx, i ) ) / u( lmx, i )
            END IF
            s = sv / su
            DO 30 j = 1, n
               res1 = max( res1, abs( u( j, i )-s*v( j, i ) ) )
   30       CONTINUE
   40    CONTINUE
         res1 = res1 / ( real( n )*ulp )
*
*        Compute orthogonality of columns of V.
*
         CALL cunt01( 'Columns', n, mv, v, ldv, work, lwork, rwork,
     $                res2 )
      END IF
*
      result = min( max( res1, res2 ), one / ulp )
      RETURN
*
*     End of CUNT03
*
      END