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