cunt01.f

Go to the documentation of this file.

*> \brief \b CUNT01
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*  Definition:
*  ===========
*
*       SUBROUTINE CUNT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK,
*                          RESID )
* 
*       .. Scalar Arguments ..
*       CHARACTER          ROWCOL
*       INTEGER            LDU, LWORK, M, N
*       REAL               RESID
*       ..
*       .. Array Arguments ..
*       REAL               RWORK( * )
*       COMPLEX            U( LDU, * ), WORK( * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> CUNT01 checks that the matrix U is unitary by computing the ratio
*>
*>    RESID = norm( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
*> or
*>    RESID = norm( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
*>
*> Alternatively, if there isn't sufficient workspace to form
*> I - U*U' or I - U'*U, the ratio is computed as
*>
*>    RESID = abs( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
*> or
*>    RESID = abs( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
*>
*> where EPS is the machine precision.  ROWCOL is used only if m = n;
*> if m > n, ROWCOL is assumed to be 'C', and if m < n, ROWCOL is
*> assumed to be 'R'.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] ROWCOL
*> \verbatim
*>          ROWCOL is CHARACTER
*>          Specifies whether the rows or columns of U should be checked
*>          for orthogonality.  Used only if M = N.
*>          = 'R':  Check for orthogonal rows of U
*>          = 'C':  Check for orthogonal columns of U
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The number of rows of the matrix U.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of columns of the matrix U.
*> \endverbatim
*>
*> \param[in] U
*> \verbatim
*>          U is COMPLEX array, dimension (LDU,N)
*>          The unitary matrix U.  U is checked for orthogonal columns
*>          if m > n or if m = n and ROWCOL = 'C'.  U is checked for
*>          orthogonal rows if m < n or if m = n and ROWCOL = 'R'.
*> \endverbatim
*>
*> \param[in] LDU
*> \verbatim
*>          LDU is INTEGER
*>          The leading dimension of the array U.  LDU >= max(1,M).
*> \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 ROWCOL = 'C' or M*M if
*>          ROWCOL = 'R', but the test will be done even if LWORK is 0.
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*>          RWORK is REAL array, dimension (min(M,N))
*>          Used only if LWORK is large enough to use the Level 3 BLAS
*>          code.
*> \endverbatim
*>
*> \param[out] RESID
*> \verbatim
*>          RESID is REAL
*>          RESID = norm( I - U * U' ) / ( n * EPS ), if ROWCOL = 'R', or
*>          RESID = norm( I - U' * U ) / ( m * EPS ), if ROWCOL = 'C'.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date November 2011
*
*> \ingroup complex_eig
*
*  =====================================================================
      SUBROUTINE cunt01( ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK,
     $                   resid )
*
*  -- LAPACK test routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      CHARACTER          ROWCOL
      INTEGER            LDU, LWORK, M, N
      REAL               RESID
*     ..
*     .. Array Arguments ..
      REAL               RWORK( * )
      COMPLEX            U( ldu, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter                ( zero = 0.0e+0, one = 1.0e+0 )
*     ..
*     .. Local Scalars ..
      CHARACTER          TRANSU
      INTEGER            I, J, K, LDWORK, MNMIN
      REAL               EPS
      COMPLEX            TMP, ZDUM
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      REAL               CLANSY, SLAMCH
      COMPLEX            CDOTC
      EXTERNAL           lsame, clansy, slamch, cdotc
*     ..
*     .. External Subroutines ..
      EXTERNAL           cherk, claset
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, aimag, cmplx, max, min, real
*     ..
*     .. Statement Functions ..
      REAL               CABS1
*     ..
*     .. Statement Function definitions ..
      cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( AIMAG( zdum ) )
*     ..
*     .. Executable Statements ..
*
      resid = zero
*
*     Quick return if possible
*
      IF( m.LE.0 .OR. n.LE.0 )
     $   RETURN
*
      eps = slamch( 'Precision' )
      IF( m.LT.n .OR. ( m.EQ.n .AND. lsame( rowcol, 'R' ) ) ) THEN
         transu = 'N'
         k = n
      ELSE
         transu = 'C'
         k = m
      END IF
      mnmin = min( m, n )
*
      IF( ( mnmin+1 )*mnmin.LE.lwork ) THEN
         ldwork = mnmin
      ELSE
         ldwork = 0
      END IF
      IF( ldwork.GT.0 ) THEN
*
*        Compute I - U*U' or I - U'*U.
*
         CALL claset( 'Upper', mnmin, mnmin, cmplx( zero ),
     $                cmplx( one ), work, ldwork )
         CALL cherk( 'Upper', transu, mnmin, k, -one, u, ldu, one, work,
     $               ldwork )
*
*        Compute norm( I - U*U' ) / ( K * EPS ) .
*
         resid = clansy( '1', 'Upper', mnmin, work, ldwork, rwork )
         resid = ( resid / REAL( K ) ) / eps
      ELSE IF( transu.EQ.'C' ) THEN
*
*        Find the maximum element in abs( I - U'*U ) / ( m * EPS )
*
         DO 20 j = 1, n
            DO 10 i = 1, j
               IF( i.NE.j ) THEN
                  tmp = zero
               ELSE
                  tmp = one
               END IF
               tmp = tmp - cdotc( m, u( 1, i ), 1, u( 1, j ), 1 )
               resid = max( resid, cabs1( tmp ) )
   10       CONTINUE
   20    CONTINUE
         resid = ( resid / REAL( M ) ) / eps
      ELSE
*
*        Find the maximum element in abs( I - U*U' ) / ( n * EPS )
*
         DO 40 j = 1, m
            DO 30 i = 1, j
               IF( i.NE.j ) THEN
                  tmp = zero
               ELSE
                  tmp = one
               END IF
               tmp = tmp - cdotc( n, u( j, 1 ), ldu, u( i, 1 ), ldu )
               resid = max( resid, cabs1( tmp ) )
   30       CONTINUE
   40    CONTINUE
         resid = ( resid / REAL( N ) ) / eps
      END IF
      RETURN
*
*     End of CUNT01
*
      END