LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zget54.f
Go to the documentation of this file.
00001 *> \brief \b ZGET54
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 ZGET54( N, A, LDA, B, LDB, S, LDS, T, LDT, U, LDU, V,
00012 *                          LDV, WORK, RESULT )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       INTEGER            LDA, LDB, LDS, LDT, LDU, LDV, N
00016 *       DOUBLE PRECISION   RESULT
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       COMPLEX*16         A( LDA, * ), B( LDB, * ), S( LDS, * ),
00020 *      $                   T( LDT, * ), U( LDU, * ), V( LDV, * ),
00021 *      $                   WORK( * )
00022 *       ..
00023 *  
00024 *
00025 *> \par Purpose:
00026 *  =============
00027 *>
00028 *> \verbatim
00029 *>
00030 *> ZGET54 checks a generalized decomposition of the form
00031 *>
00032 *>          A = U*S*V'  and B = U*T* V'
00033 *>
00034 *> where ' means conjugate transpose and U and V are unitary.
00035 *>
00036 *> Specifically,
00037 *>
00038 *>   RESULT = ||( A - U*S*V', B - U*T*V' )|| / (||( A, B )||*n*ulp )
00039 *> \endverbatim
00040 *
00041 *  Arguments:
00042 *  ==========
00043 *
00044 *> \param[in] N
00045 *> \verbatim
00046 *>          N is INTEGER
00047 *>          The size of the matrix.  If it is zero, DGET54 does nothing.
00048 *>          It must be at least zero.
00049 *> \endverbatim
00050 *>
00051 *> \param[in] A
00052 *> \verbatim
00053 *>          A is COMPLEX*16 array, dimension (LDA, N)
00054 *>          The original (unfactored) matrix A.
00055 *> \endverbatim
00056 *>
00057 *> \param[in] LDA
00058 *> \verbatim
00059 *>          LDA is INTEGER
00060 *>          The leading dimension of A.  It must be at least 1
00061 *>          and at least N.
00062 *> \endverbatim
00063 *>
00064 *> \param[in] B
00065 *> \verbatim
00066 *>          B is COMPLEX*16 array, dimension (LDB, N)
00067 *>          The original (unfactored) matrix B.
00068 *> \endverbatim
00069 *>
00070 *> \param[in] LDB
00071 *> \verbatim
00072 *>          LDB is INTEGER
00073 *>          The leading dimension of B.  It must be at least 1
00074 *>          and at least N.
00075 *> \endverbatim
00076 *>
00077 *> \param[in] S
00078 *> \verbatim
00079 *>          S is COMPLEX*16 array, dimension (LDS, N)
00080 *>          The factored matrix S.
00081 *> \endverbatim
00082 *>
00083 *> \param[in] LDS
00084 *> \verbatim
00085 *>          LDS is INTEGER
00086 *>          The leading dimension of S.  It must be at least 1
00087 *>          and at least N.
00088 *> \endverbatim
00089 *>
00090 *> \param[in] T
00091 *> \verbatim
00092 *>          T is COMPLEX*16 array, dimension (LDT, N)
00093 *>          The factored matrix T.
00094 *> \endverbatim
00095 *>
00096 *> \param[in] LDT
00097 *> \verbatim
00098 *>          LDT is INTEGER
00099 *>          The leading dimension of T.  It must be at least 1
00100 *>          and at least N.
00101 *> \endverbatim
00102 *>
00103 *> \param[in] U
00104 *> \verbatim
00105 *>          U is COMPLEX*16 array, dimension (LDU, N)
00106 *>          The orthogonal matrix on the left-hand side in the
00107 *>          decomposition.
00108 *> \endverbatim
00109 *>
00110 *> \param[in] LDU
00111 *> \verbatim
00112 *>          LDU is INTEGER
00113 *>          The leading dimension of U.  LDU must be at least N and
00114 *>          at least 1.
00115 *> \endverbatim
00116 *>
00117 *> \param[in] V
00118 *> \verbatim
00119 *>          V is COMPLEX*16 array, dimension (LDV, N)
00120 *>          The orthogonal matrix on the left-hand side in the
00121 *>          decomposition.
00122 *> \endverbatim
00123 *>
00124 *> \param[in] LDV
00125 *> \verbatim
00126 *>          LDV is INTEGER
00127 *>          The leading dimension of V.  LDV must be at least N and
00128 *>          at least 1.
00129 *> \endverbatim
00130 *>
00131 *> \param[out] WORK
00132 *> \verbatim
00133 *>          WORK is COMPLEX*16 array, dimension (3*N**2)
00134 *> \endverbatim
00135 *>
00136 *> \param[out] RESULT
00137 *> \verbatim
00138 *>          RESULT is DOUBLE PRECISION
00139 *>          The value RESULT, It is currently limited to 1/ulp, to
00140 *>          avoid overflow. Errors are flagged by RESULT=10/ulp.
00141 *> \endverbatim
00142 *
00143 *  Authors:
00144 *  ========
00145 *
00146 *> \author Univ. of Tennessee 
00147 *> \author Univ. of California Berkeley 
00148 *> \author Univ. of Colorado Denver 
00149 *> \author NAG Ltd. 
00150 *
00151 *> \date November 2011
00152 *
00153 *> \ingroup complex16_eig
00154 *
00155 *  =====================================================================
00156       SUBROUTINE ZGET54( N, A, LDA, B, LDB, S, LDS, T, LDT, U, LDU, V,
00157      $                   LDV, WORK, RESULT )
00158 *
00159 *  -- LAPACK test routine (version 3.4.0) --
00160 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00161 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00162 *     November 2011
00163 *
00164 *     .. Scalar Arguments ..
00165       INTEGER            LDA, LDB, LDS, LDT, LDU, LDV, N
00166       DOUBLE PRECISION   RESULT
00167 *     ..
00168 *     .. Array Arguments ..
00169       COMPLEX*16         A( LDA, * ), B( LDB, * ), S( LDS, * ),
00170      $                   T( LDT, * ), U( LDU, * ), V( LDV, * ),
00171      $                   WORK( * )
00172 *     ..
00173 *
00174 *  =====================================================================
00175 *
00176 *     .. Parameters ..
00177       DOUBLE PRECISION   ZERO, ONE
00178       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00179       COMPLEX*16         CZERO, CONE
00180       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00181      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00182 *     ..
00183 *     .. Local Scalars ..
00184       DOUBLE PRECISION   ABNORM, ULP, UNFL, WNORM
00185 *     ..
00186 *     .. Local Arrays ..
00187       DOUBLE PRECISION   DUM( 1 )
00188 *     ..
00189 *     .. External Functions ..
00190       DOUBLE PRECISION   DLAMCH, ZLANGE
00191       EXTERNAL           DLAMCH, ZLANGE
00192 *     ..
00193 *     .. External Subroutines ..
00194       EXTERNAL           ZGEMM, ZLACPY
00195 *     ..
00196 *     .. Intrinsic Functions ..
00197       INTRINSIC          DBLE, MAX, MIN
00198 *     ..
00199 *     .. Executable Statements ..
00200 *
00201       RESULT = ZERO
00202       IF( N.LE.0 )
00203      $   RETURN
00204 *
00205 *     Constants
00206 *
00207       UNFL = DLAMCH( 'Safe minimum' )
00208       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
00209 *
00210 *     compute the norm of (A,B)
00211 *
00212       CALL ZLACPY( 'Full', N, N, A, LDA, WORK, N )
00213       CALL ZLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
00214       ABNORM = MAX( ZLANGE( '1', N, 2*N, WORK, N, DUM ), UNFL )
00215 *
00216 *     Compute W1 = A - U*S*V', and put in the array WORK(1:N*N)
00217 *
00218       CALL ZLACPY( ' ', N, N, A, LDA, WORK, N )
00219       CALL ZGEMM( 'N', 'N', N, N, N, CONE, U, LDU, S, LDS, CZERO,
00220      $            WORK( N*N+1 ), N )
00221 *
00222       CALL ZGEMM( 'N', 'C', N, N, N, -CONE, WORK( N*N+1 ), N, V, LDV,
00223      $            CONE, WORK, N )
00224 *
00225 *     Compute W2 = B - U*T*V', and put in the workarray W(N*N+1:2*N*N)
00226 *
00227       CALL ZLACPY( ' ', N, N, B, LDB, WORK( N*N+1 ), N )
00228       CALL ZGEMM( 'N', 'N', N, N, N, CONE, U, LDU, T, LDT, CZERO,
00229      $            WORK( 2*N*N+1 ), N )
00230 *
00231       CALL ZGEMM( 'N', 'C', N, N, N, -CONE, WORK( 2*N*N+1 ), N, V, LDV,
00232      $            CONE, WORK( N*N+1 ), N )
00233 *
00234 *     Compute norm(W)/ ( ulp*norm((A,B)) )
00235 *
00236       WNORM = ZLANGE( '1', N, 2*N, WORK, N, DUM )
00237 *
00238       IF( ABNORM.GT.WNORM ) THEN
00239          RESULT = ( WNORM / ABNORM ) / ( 2*N*ULP )
00240       ELSE
00241          IF( ABNORM.LT.ONE ) THEN
00242             RESULT = ( MIN( WNORM, 2*N*ABNORM ) / ABNORM ) / ( 2*N*ULP )
00243          ELSE
00244             RESULT = MIN( WNORM / ABNORM, DBLE( 2*N ) ) / ( 2*N*ULP )
00245          END IF
00246       END IF
00247 *
00248       RETURN
00249 *
00250 *     End of ZGET54
00251 *
00252       END
 All Files Functions