LAPACK 3.3.0

dget51.f

Go to the documentation of this file.
00001       SUBROUTINE DGET51( ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK,
00002      $                   RESULT )
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       INTEGER            ITYPE, LDA, LDB, LDU, LDV, N
00010       DOUBLE PRECISION   RESULT
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), U( LDU, * ),
00014      $                   V( LDV, * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *       DGET51  generally checks a decomposition of the form
00021 *
00022 *               A = U B V'
00023 *
00024 *       where ' means transpose and U and V are orthogonal.
00025 *
00026 *       Specifically, if ITYPE=1
00027 *
00028 *               RESULT = | A - U B V' | / ( |A| n ulp )
00029 *
00030 *       If ITYPE=2, then:
00031 *
00032 *               RESULT = | A - B | / ( |A| n ulp )
00033 *
00034 *       If ITYPE=3, then:
00035 *
00036 *               RESULT = | I - UU' | / ( n ulp )
00037 *
00038 *  Arguments
00039 *  =========
00040 *
00041 *  ITYPE   (input) INTEGER
00042 *          Specifies the type of tests to be performed.
00043 *          =1: RESULT = | A - U B V' | / ( |A| n ulp )
00044 *          =2: RESULT = | A - B | / ( |A| n ulp )
00045 *          =3: RESULT = | I - UU' | / ( n ulp )
00046 *
00047 *  N       (input) INTEGER
00048 *          The size of the matrix.  If it is zero, DGET51 does nothing.
00049 *          It must be at least zero.
00050 *
00051 *  A       (input) DOUBLE PRECISION array, dimension (LDA, N)
00052 *          The original (unfactored) matrix.
00053 *
00054 *  LDA     (input) INTEGER
00055 *          The leading dimension of A.  It must be at least 1
00056 *          and at least N.
00057 *
00058 *  B       (input) DOUBLE PRECISION array, dimension (LDB, N)
00059 *          The factored matrix.
00060 *
00061 *  LDB     (input) INTEGER
00062 *          The leading dimension of B.  It must be at least 1
00063 *          and at least N.
00064 *
00065 *  U       (input) DOUBLE PRECISION array, dimension (LDU, N)
00066 *          The orthogonal matrix on the left-hand side in the
00067 *          decomposition.
00068 *          Not referenced if ITYPE=2
00069 *
00070 *  LDU     (input) INTEGER
00071 *          The leading dimension of U.  LDU must be at least N and
00072 *          at least 1.
00073 *
00074 *  V       (input) DOUBLE PRECISION array, dimension (LDV, N)
00075 *          The orthogonal matrix on the left-hand side in the
00076 *          decomposition.
00077 *          Not referenced if ITYPE=2
00078 *
00079 *  LDV     (input) INTEGER
00080 *          The leading dimension of V.  LDV must be at least N and
00081 *          at least 1.
00082 *
00083 *  WORK    (workspace) DOUBLE PRECISION array, dimension (2*N**2)
00084 *
00085 *  RESULT  (output) DOUBLE PRECISION
00086 *          The values computed by the test specified by ITYPE.  The
00087 *          value is currently limited to 1/ulp, to avoid overflow.
00088 *          Errors are flagged by RESULT=10/ulp.
00089 *
00090 *  =====================================================================
00091 *
00092 *     .. Parameters ..
00093       DOUBLE PRECISION   ZERO, ONE, TEN
00094       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0 )
00095 *     ..
00096 *     .. Local Scalars ..
00097       INTEGER            JCOL, JDIAG, JROW
00098       DOUBLE PRECISION   ANORM, ULP, UNFL, WNORM
00099 *     ..
00100 *     .. External Functions ..
00101       DOUBLE PRECISION   DLAMCH, DLANGE
00102       EXTERNAL           DLAMCH, DLANGE
00103 *     ..
00104 *     .. External Subroutines ..
00105       EXTERNAL           DGEMM, DLACPY
00106 *     ..
00107 *     .. Intrinsic Functions ..
00108       INTRINSIC          DBLE, MAX, MIN
00109 *     ..
00110 *     .. Executable Statements ..
00111 *
00112       RESULT = ZERO
00113       IF( N.LE.0 )
00114      $   RETURN
00115 *
00116 *     Constants
00117 *
00118       UNFL = DLAMCH( 'Safe minimum' )
00119       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
00120 *
00121 *     Some Error Checks
00122 *
00123       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
00124          RESULT = TEN / ULP
00125          RETURN
00126       END IF
00127 *
00128       IF( ITYPE.LE.2 ) THEN
00129 *
00130 *        Tests scaled by the norm(A)
00131 *
00132          ANORM = MAX( DLANGE( '1', N, N, A, LDA, WORK ), UNFL )
00133 *
00134          IF( ITYPE.EQ.1 ) THEN
00135 *
00136 *           ITYPE=1: Compute W = A - UBV'
00137 *
00138             CALL DLACPY( ' ', N, N, A, LDA, WORK, N )
00139             CALL DGEMM( 'N', 'N', N, N, N, ONE, U, LDU, B, LDB, ZERO,
00140      $                  WORK( N**2+1 ), N )
00141 *
00142             CALL DGEMM( 'N', 'C', N, N, N, -ONE, WORK( N**2+1 ), N, V,
00143      $                  LDV, ONE, WORK, N )
00144 *
00145          ELSE
00146 *
00147 *           ITYPE=2: Compute W = A - B
00148 *
00149             CALL DLACPY( ' ', N, N, B, LDB, WORK, N )
00150 *
00151             DO 20 JCOL = 1, N
00152                DO 10 JROW = 1, N
00153                   WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) )
00154      $                - A( JROW, JCOL )
00155    10          CONTINUE
00156    20       CONTINUE
00157          END IF
00158 *
00159 *        Compute norm(W)/ ( ulp*norm(A) )
00160 *
00161          WNORM = DLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) )
00162 *
00163          IF( ANORM.GT.WNORM ) THEN
00164             RESULT = ( WNORM / ANORM ) / ( N*ULP )
00165          ELSE
00166             IF( ANORM.LT.ONE ) THEN
00167                RESULT = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
00168             ELSE
00169                RESULT = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP )
00170             END IF
00171          END IF
00172 *
00173       ELSE
00174 *
00175 *        Tests not scaled by norm(A)
00176 *
00177 *        ITYPE=3: Compute  UU' - I
00178 *
00179          CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK,
00180      $               N )
00181 *
00182          DO 30 JDIAG = 1, N
00183             WORK( ( N+1 )*( JDIAG-1 )+1 ) = WORK( ( N+1 )*( JDIAG-1 )+
00184      $         1 ) - ONE
00185    30    CONTINUE
00186 *
00187          RESULT = MIN( DLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) ),
00188      $            DBLE( N ) ) / ( N*ULP )
00189       END IF
00190 *
00191       RETURN
00192 *
00193 *     End of DGET51
00194 *
00195       END
 All Files Functions