LAPACK: TESTING/EIG/zunt01.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZUNT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK,
00002      $                   RESID )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          ROWCOL
00010       INTEGER            LDU, LWORK, M, N
00011       DOUBLE PRECISION   RESID
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   RWORK( * )
00015       COMPLEX*16         U( LDU, * ), WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  ZUNT01 checks that the matrix U is unitary by computing the ratio
00022 *
00023 *     RESID = norm( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
00024 *  or
00025 *     RESID = norm( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
00026 *
00027 *  Alternatively, if there isn't sufficient workspace to form
00028 *  I - U*U' or I - U'*U, the ratio is computed as
00029 *
00030 *     RESID = abs( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
00031 *  or
00032 *     RESID = abs( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
00033 *
00034 *  where EPS is the machine precision.  ROWCOL is used only if m = n;
00035 *  if m > n, ROWCOL is assumed to be 'C', and if m < n, ROWCOL is
00036 *  assumed to be 'R'.
00037 *
00038 *  Arguments
00039 *  =========
00040 *
00041 *  ROWCOL  (input) CHARACTER
00042 *          Specifies whether the rows or columns of U should be checked
00043 *          for orthogonality.  Used only if M = N.
00044 *          = 'R':  Check for orthogonal rows of U
00045 *          = 'C':  Check for orthogonal columns of U
00046 *
00047 *  M       (input) INTEGER
00048 *          The number of rows of the matrix U.
00049 *
00050 *  N       (input) INTEGER
00051 *          The number of columns of the matrix U.
00052 *
00053 *  U       (input) COMPLEX*16 array, dimension (LDU,N)
00054 *          The unitary matrix U.  U is checked for orthogonal columns
00055 *          if m > n or if m = n and ROWCOL = 'C'.  U is checked for
00056 *          orthogonal rows if m < n or if m = n and ROWCOL = 'R'.
00057 *
00058 *  LDU     (input) INTEGER
00059 *          The leading dimension of the array U.  LDU >= max(1,M).
00060 *
00061 *  WORK    (workspace) COMPLEX*16 array, dimension (LWORK)
00062 *
00063 *  LWORK   (input) INTEGER
00064 *          The length of the array WORK.  For best performance, LWORK
00065 *          should be at least N*N if ROWCOL = 'C' or M*M if
00066 *          ROWCOL = 'R', but the test will be done even if LWORK is 0.
00067 *
00068 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (min(M,N))
00069 *          Used only if LWORK is large enough to use the Level 3 BLAS
00070 *          code.
00071 *
00072 *  RESID   (output) DOUBLE PRECISION
00073 *          RESID = norm( I - U * U' ) / ( n * EPS ), if ROWCOL = 'R', or
00074 *          RESID = norm( I - U' * U ) / ( m * EPS ), if ROWCOL = 'C'.
00075 *
00076 *  =====================================================================
00077 *
00078 *     .. Parameters ..
00079       DOUBLE PRECISION   ZERO, ONE
00080       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00081 *     ..
00082 *     .. Local Scalars ..
00083       CHARACTER          TRANSU
00084       INTEGER            I, J, K, LDWORK, MNMIN
00085       DOUBLE PRECISION   EPS
00086       COMPLEX*16         TMP, ZDUM
00087 *     ..
00088 *     .. External Functions ..
00089       LOGICAL            LSAME
00090       DOUBLE PRECISION   DLAMCH, ZLANSY
00091       COMPLEX*16         ZDOTC
00092       EXTERNAL           LSAME, DLAMCH, ZLANSY, ZDOTC
00093 *     ..
00094 *     .. External Subroutines ..
00095       EXTERNAL           ZHERK, ZLASET
00096 *     ..
00097 *     .. Intrinsic Functions ..
00098       INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, MAX, MIN
00099 *     ..
00100 *     .. Statement Functions ..
00101       DOUBLE PRECISION   CABS1
00102 *     ..
00103 *     .. Statement Function definitions ..
00104       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
00105 *     ..
00106 *     .. Executable Statements ..
00107 *
00108       RESID = ZERO
00109 *
00110 *     Quick return if possible
00111 *
00112       IF( M.LE.0 .OR. N.LE.0 )
00113      $   RETURN
00114 *
00115       EPS = DLAMCH( 'Precision' )
00116       IF( M.LT.N .OR. ( M.EQ.N .AND. LSAME( ROWCOL, 'R' ) ) ) THEN
00117          TRANSU = 'N'
00118          K = N
00119       ELSE
00120          TRANSU = 'C'
00121          K = M
00122       END IF
00123       MNMIN = MIN( M, N )
00124 *
00125       IF( ( MNMIN+1 )*MNMIN.LE.LWORK ) THEN
00126          LDWORK = MNMIN
00127       ELSE
00128          LDWORK = 0
00129       END IF
00130       IF( LDWORK.GT.0 ) THEN
00131 *
00132 *        Compute I - U*U' or I - U'*U.
00133 *
00134          CALL ZLASET( 'Upper', MNMIN, MNMIN, DCMPLX( ZERO ),
00135      $                DCMPLX( ONE ), WORK, LDWORK )
00136          CALL ZHERK( 'Upper', TRANSU, MNMIN, K, -ONE, U, LDU, ONE, WORK,
00137      $               LDWORK )
00138 *
00139 *        Compute norm( I - U*U' ) / ( K * EPS ) .
00140 *
00141          RESID = ZLANSY( '1', 'Upper', MNMIN, WORK, LDWORK, RWORK )
00142          RESID = ( RESID / DBLE( K ) ) / EPS
00143       ELSE IF( TRANSU.EQ.'C' ) THEN
00144 *
00145 *        Find the maximum element in abs( I - U'*U ) / ( m * EPS )
00146 *
00147          DO 20 J = 1, N
00148             DO 10 I = 1, J
00149                IF( I.NE.J ) THEN
00150                   TMP = ZERO
00151                ELSE
00152                   TMP = ONE
00153                END IF
00154                TMP = TMP - ZDOTC( M, U( 1, I ), 1, U( 1, J ), 1 )
00155                RESID = MAX( RESID, CABS1( TMP ) )
00156    10       CONTINUE
00157    20    CONTINUE
00158          RESID = ( RESID / DBLE( M ) ) / EPS
00159       ELSE
00160 *
00161 *        Find the maximum element in abs( I - U*U' ) / ( n * EPS )
00162 *
00163          DO 40 J = 1, M
00164             DO 30 I = 1, J
00165                IF( I.NE.J ) THEN
00166                   TMP = ZERO
00167                ELSE
00168                   TMP = ONE
00169                END IF
00170                TMP = TMP - ZDOTC( N, U( J, 1 ), LDU, U( I, 1 ), LDU )
00171                RESID = MAX( RESID, CABS1( TMP ) )
00172    30       CONTINUE
00173    40    CONTINUE
00174          RESID = ( RESID / DBLE( N ) ) / EPS
00175       END IF
00176       RETURN
00177 *
00178 *     End of ZUNT01
00179 *
00180       END