LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sget54.f
Go to the documentation of this file.
00001 *> \brief \b SGET54
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 SGET54( 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 *       REAL               RESULT
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       REAL               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 *> SGET54 checks a generalized decomposition of the form
00031 *>
00032 *>          A = U*S*V'  and B = U*T* V'
00033 *>
00034 *> where ' means transpose and U and V are orthogonal.
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, SGET54 does nothing.
00048 *>          It must be at least zero.
00049 *> \endverbatim
00050 *>
00051 *> \param[in] A
00052 *> \verbatim
00053 *>          A is REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL array, dimension (3*N**2)
00134 *> \endverbatim
00135 *>
00136 *> \param[out] RESULT
00137 *> \verbatim
00138 *>          RESULT is REAL
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 single_eig
00154 *
00155 *  =====================================================================
00156       SUBROUTINE SGET54( 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       REAL               RESULT
00167 *     ..
00168 *     .. Array Arguments ..
00169       REAL               A( LDA, * ), B( LDB, * ), S( LDS, * ),
00170      $                   T( LDT, * ), U( LDU, * ), V( LDV, * ),
00171      $                   WORK( * )
00172 *     ..
00173 *
00174 *  =====================================================================
00175 *
00176 *     .. Parameters ..
00177       REAL               ZERO, ONE
00178       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00179 *     ..
00180 *     .. Local Scalars ..
00181       REAL               ABNORM, ULP, UNFL, WNORM
00182 *     ..
00183 *     .. Local Arrays ..
00184       REAL               DUM( 1 )
00185 *     ..
00186 *     .. External Functions ..
00187       REAL               SLAMCH, SLANGE
00188       EXTERNAL           SLAMCH, SLANGE
00189 *     ..
00190 *     .. External Subroutines ..
00191       EXTERNAL           SGEMM, SLACPY
00192 *     ..
00193 *     .. Intrinsic Functions ..
00194       INTRINSIC          MAX, MIN, REAL
00195 *     ..
00196 *     .. Executable Statements ..
00197 *
00198       RESULT = ZERO
00199       IF( N.LE.0 )
00200      $   RETURN
00201 *
00202 *     Constants
00203 *
00204       UNFL = SLAMCH( 'Safe minimum' )
00205       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
00206 *
00207 *     compute the norm of (A,B)
00208 *
00209       CALL SLACPY( 'Full', N, N, A, LDA, WORK, N )
00210       CALL SLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
00211       ABNORM = MAX( SLANGE( '1', N, 2*N, WORK, N, DUM ), UNFL )
00212 *
00213 *     Compute W1 = A - U*S*V', and put in the array WORK(1:N*N)
00214 *
00215       CALL SLACPY( ' ', N, N, A, LDA, WORK, N )
00216       CALL SGEMM( 'N', 'N', N, N, N, ONE, U, LDU, S, LDS, ZERO,
00217      $            WORK( N*N+1 ), N )
00218 *
00219       CALL SGEMM( 'N', 'C', N, N, N, -ONE, WORK( N*N+1 ), N, V, LDV,
00220      $            ONE, WORK, N )
00221 *
00222 *     Compute W2 = B - U*T*V', and put in the workarray W(N*N+1:2*N*N)
00223 *
00224       CALL SLACPY( ' ', N, N, B, LDB, WORK( N*N+1 ), N )
00225       CALL SGEMM( 'N', 'N', N, N, N, ONE, U, LDU, T, LDT, ZERO,
00226      $            WORK( 2*N*N+1 ), N )
00227 *
00228       CALL SGEMM( 'N', 'C', N, N, N, -ONE, WORK( 2*N*N+1 ), N, V, LDV,
00229      $            ONE, WORK( N*N+1 ), N )
00230 *
00231 *     Compute norm(W)/ ( ulp*norm((A,B)) )
00232 *
00233       WNORM = SLANGE( '1', N, 2*N, WORK, N, DUM )
00234 *
00235       IF( ABNORM.GT.WNORM ) THEN
00236          RESULT = ( WNORM / ABNORM ) / ( 2*N*ULP )
00237       ELSE
00238          IF( ABNORM.LT.ONE ) THEN
00239             RESULT = ( MIN( WNORM, 2*N*ABNORM ) / ABNORM ) / ( 2*N*ULP )
00240          ELSE
00241             RESULT = MIN( WNORM / ABNORM, REAL( 2*N ) ) / ( 2*N*ULP )
00242          END IF
00243       END IF
00244 *
00245       RETURN
00246 *
00247 *     End of SGET54
00248 *
00249       END
 All Files Functions