LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zunt01.f
Go to the documentation of this file.
00001 *> \brief \b ZUNT01
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE ZUNT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK,
00012 *                          RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          ROWCOL
00016 *       INTEGER            LDU, LWORK, M, N
00017 *       DOUBLE PRECISION   RESID
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       DOUBLE PRECISION   RWORK( * )
00021 *       COMPLEX*16         U( LDU, * ), WORK( * )
00022 *       ..
00023 *  
00024 *
00025 *> \par Purpose:
00026 *  =============
00027 *>
00028 *> \verbatim
00029 *>
00030 *> ZUNT01 checks that the matrix U is unitary by computing the ratio
00031 *>
00032 *>    RESID = norm( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
00033 *> or
00034 *>    RESID = norm( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
00035 *>
00036 *> Alternatively, if there isn't sufficient workspace to form
00037 *> I - U*U' or I - U'*U, the ratio is computed as
00038 *>
00039 *>    RESID = abs( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
00040 *> or
00041 *>    RESID = abs( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
00042 *>
00043 *> where EPS is the machine precision.  ROWCOL is used only if m = n;
00044 *> if m > n, ROWCOL is assumed to be 'C', and if m < n, ROWCOL is
00045 *> assumed to be 'R'.
00046 *> \endverbatim
00047 *
00048 *  Arguments:
00049 *  ==========
00050 *
00051 *> \param[in] ROWCOL
00052 *> \verbatim
00053 *>          ROWCOL is CHARACTER
00054 *>          Specifies whether the rows or columns of U should be checked
00055 *>          for orthogonality.  Used only if M = N.
00056 *>          = 'R':  Check for orthogonal rows of U
00057 *>          = 'C':  Check for orthogonal columns of U
00058 *> \endverbatim
00059 *>
00060 *> \param[in] M
00061 *> \verbatim
00062 *>          M is INTEGER
00063 *>          The number of rows of the matrix U.
00064 *> \endverbatim
00065 *>
00066 *> \param[in] N
00067 *> \verbatim
00068 *>          N is INTEGER
00069 *>          The number of columns of the matrix U.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] U
00073 *> \verbatim
00074 *>          U is COMPLEX*16 array, dimension (LDU,N)
00075 *>          The unitary matrix U.  U is checked for orthogonal columns
00076 *>          if m > n or if m = n and ROWCOL = 'C'.  U is checked for
00077 *>          orthogonal rows if m < n or if m = n and ROWCOL = 'R'.
00078 *> \endverbatim
00079 *>
00080 *> \param[in] LDU
00081 *> \verbatim
00082 *>          LDU is INTEGER
00083 *>          The leading dimension of the array U.  LDU >= max(1,M).
00084 *> \endverbatim
00085 *>
00086 *> \param[out] WORK
00087 *> \verbatim
00088 *>          WORK is COMPLEX*16 array, dimension (LWORK)
00089 *> \endverbatim
00090 *>
00091 *> \param[in] LWORK
00092 *> \verbatim
00093 *>          LWORK is INTEGER
00094 *>          The length of the array WORK.  For best performance, LWORK
00095 *>          should be at least N*N if ROWCOL = 'C' or M*M if
00096 *>          ROWCOL = 'R', but the test will be done even if LWORK is 0.
00097 *> \endverbatim
00098 *>
00099 *> \param[out] RWORK
00100 *> \verbatim
00101 *>          RWORK is DOUBLE PRECISION array, dimension (min(M,N))
00102 *>          Used only if LWORK is large enough to use the Level 3 BLAS
00103 *>          code.
00104 *> \endverbatim
00105 *>
00106 *> \param[out] RESID
00107 *> \verbatim
00108 *>          RESID is DOUBLE PRECISION
00109 *>          RESID = norm( I - U * U' ) / ( n * EPS ), if ROWCOL = 'R', or
00110 *>          RESID = norm( I - U' * U ) / ( m * EPS ), if ROWCOL = 'C'.
00111 *> \endverbatim
00112 *
00113 *  Authors:
00114 *  ========
00115 *
00116 *> \author Univ. of Tennessee 
00117 *> \author Univ. of California Berkeley 
00118 *> \author Univ. of Colorado Denver 
00119 *> \author NAG Ltd. 
00120 *
00121 *> \date November 2011
00122 *
00123 *> \ingroup complex16_eig
00124 *
00125 *  =====================================================================
00126       SUBROUTINE ZUNT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK,
00127      $                   RESID )
00128 *
00129 *  -- LAPACK test routine (version 3.4.0) --
00130 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00131 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00132 *     November 2011
00133 *
00134 *     .. Scalar Arguments ..
00135       CHARACTER          ROWCOL
00136       INTEGER            LDU, LWORK, M, N
00137       DOUBLE PRECISION   RESID
00138 *     ..
00139 *     .. Array Arguments ..
00140       DOUBLE PRECISION   RWORK( * )
00141       COMPLEX*16         U( LDU, * ), WORK( * )
00142 *     ..
00143 *
00144 *  =====================================================================
00145 *
00146 *     .. Parameters ..
00147       DOUBLE PRECISION   ZERO, ONE
00148       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00149 *     ..
00150 *     .. Local Scalars ..
00151       CHARACTER          TRANSU
00152       INTEGER            I, J, K, LDWORK, MNMIN
00153       DOUBLE PRECISION   EPS
00154       COMPLEX*16         TMP, ZDUM
00155 *     ..
00156 *     .. External Functions ..
00157       LOGICAL            LSAME
00158       DOUBLE PRECISION   DLAMCH, ZLANSY
00159       COMPLEX*16         ZDOTC
00160       EXTERNAL           LSAME, DLAMCH, ZLANSY, ZDOTC
00161 *     ..
00162 *     .. External Subroutines ..
00163       EXTERNAL           ZHERK, ZLASET
00164 *     ..
00165 *     .. Intrinsic Functions ..
00166       INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, MAX, MIN
00167 *     ..
00168 *     .. Statement Functions ..
00169       DOUBLE PRECISION   CABS1
00170 *     ..
00171 *     .. Statement Function definitions ..
00172       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
00173 *     ..
00174 *     .. Executable Statements ..
00175 *
00176       RESID = ZERO
00177 *
00178 *     Quick return if possible
00179 *
00180       IF( M.LE.0 .OR. N.LE.0 )
00181      $   RETURN
00182 *
00183       EPS = DLAMCH( 'Precision' )
00184       IF( M.LT.N .OR. ( M.EQ.N .AND. LSAME( ROWCOL, 'R' ) ) ) THEN
00185          TRANSU = 'N'
00186          K = N
00187       ELSE
00188          TRANSU = 'C'
00189          K = M
00190       END IF
00191       MNMIN = MIN( M, N )
00192 *
00193       IF( ( MNMIN+1 )*MNMIN.LE.LWORK ) THEN
00194          LDWORK = MNMIN
00195       ELSE
00196          LDWORK = 0
00197       END IF
00198       IF( LDWORK.GT.0 ) THEN
00199 *
00200 *        Compute I - U*U' or I - U'*U.
00201 *
00202          CALL ZLASET( 'Upper', MNMIN, MNMIN, DCMPLX( ZERO ),
00203      $                DCMPLX( ONE ), WORK, LDWORK )
00204          CALL ZHERK( 'Upper', TRANSU, MNMIN, K, -ONE, U, LDU, ONE, WORK,
00205      $               LDWORK )
00206 *
00207 *        Compute norm( I - U*U' ) / ( K * EPS ) .
00208 *
00209          RESID = ZLANSY( '1', 'Upper', MNMIN, WORK, LDWORK, RWORK )
00210          RESID = ( RESID / DBLE( K ) ) / EPS
00211       ELSE IF( TRANSU.EQ.'C' ) THEN
00212 *
00213 *        Find the maximum element in abs( I - U'*U ) / ( m * EPS )
00214 *
00215          DO 20 J = 1, N
00216             DO 10 I = 1, J
00217                IF( I.NE.J ) THEN
00218                   TMP = ZERO
00219                ELSE
00220                   TMP = ONE
00221                END IF
00222                TMP = TMP - ZDOTC( M, U( 1, I ), 1, U( 1, J ), 1 )
00223                RESID = MAX( RESID, CABS1( TMP ) )
00224    10       CONTINUE
00225    20    CONTINUE
00226          RESID = ( RESID / DBLE( M ) ) / EPS
00227       ELSE
00228 *
00229 *        Find the maximum element in abs( I - U*U' ) / ( n * EPS )
00230 *
00231          DO 40 J = 1, M
00232             DO 30 I = 1, J
00233                IF( I.NE.J ) THEN
00234                   TMP = ZERO
00235                ELSE
00236                   TMP = ONE
00237                END IF
00238                TMP = TMP - ZDOTC( N, U( J, 1 ), LDU, U( I, 1 ), LDU )
00239                RESID = MAX( RESID, CABS1( TMP ) )
00240    30       CONTINUE
00241    40    CONTINUE
00242          RESID = ( RESID / DBLE( N ) ) / EPS
00243       END IF
00244       RETURN
00245 *
00246 *     End of ZUNT01
00247 *
00248       END
 All Files Functions