LAPACK: TESTING/EIG/sort01.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SORT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RESID )
00002 *
00003 *  -- LAPACK test routine (version 3.1) --
00004 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00005 *     November 2006
00006 *
00007 *     .. Scalar Arguments ..
00008       CHARACTER          ROWCOL
00009       INTEGER            LDU, LWORK, M, N
00010       REAL               RESID
00011 *     ..
00012 *     .. Array Arguments ..
00013       REAL               U( LDU, * ), WORK( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  SORT01 checks that the matrix U is orthogonal by computing the ratio
00020 *
00021 *     RESID = norm( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
00022 *  or
00023 *     RESID = norm( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
00024 *
00025 *  Alternatively, if there isn't sufficient workspace to form
00026 *  I - U*U' or I - U'*U, the ratio is computed as
00027 *
00028 *     RESID = abs( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
00029 *  or
00030 *     RESID = abs( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
00031 *
00032 *  where EPS is the machine precision.  ROWCOL is used only if m = n;
00033 *  if m > n, ROWCOL is assumed to be 'C', and if m < n, ROWCOL is
00034 *  assumed to be 'R'.
00035 *
00036 *  Arguments
00037 *  =========
00038 *
00039 *  ROWCOL  (input) CHARACTER
00040 *          Specifies whether the rows or columns of U should be checked
00041 *          for orthogonality.  Used only if M = N.
00042 *          = 'R':  Check for orthogonal rows of U
00043 *          = 'C':  Check for orthogonal columns of U
00044 *
00045 *  M       (input) INTEGER
00046 *          The number of rows of the matrix U.
00047 *
00048 *  N       (input) INTEGER
00049 *          The number of columns of the matrix U.
00050 *
00051 *  U       (input) REAL array, dimension (LDU,N)
00052 *          The orthogonal matrix U.  U is checked for orthogonal columns
00053 *          if m > n or if m = n and ROWCOL = 'C'.  U is checked for
00054 *          orthogonal rows if m < n or if m = n and ROWCOL = 'R'.
00055 *
00056 *  LDU     (input) INTEGER
00057 *          The leading dimension of the array U.  LDU >= max(1,M).
00058 *
00059 *  WORK    (workspace) REAL array, dimension (LWORK)
00060 *
00061 *  LWORK   (input) INTEGER
00062 *          The length of the array WORK.  For best performance, LWORK
00063 *          should be at least N*(N+1) if ROWCOL = 'C' or M*(M+1) if
00064 *          ROWCOL = 'R', but the test will be done even if LWORK is 0.
00065 *
00066 *  RESID   (output) REAL
00067 *          RESID = norm( I - U * U' ) / ( n * EPS ), if ROWCOL = 'R', or
00068 *          RESID = norm( I - U' * U ) / ( m * EPS ), if ROWCOL = 'C'.
00069 *
00070 *  =====================================================================
00071 *
00072 *     .. Parameters ..
00073       REAL               ZERO, ONE
00074       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00075 *     ..
00076 *     .. Local Scalars ..
00077       CHARACTER          TRANSU
00078       INTEGER            I, J, K, LDWORK, MNMIN
00079       REAL               EPS, TMP
00080 *     ..
00081 *     .. External Functions ..
00082       LOGICAL            LSAME
00083       REAL               SDOT, SLAMCH, SLANSY
00084       EXTERNAL           LSAME, SDOT, SLAMCH, SLANSY
00085 *     ..
00086 *     .. External Subroutines ..
00087       EXTERNAL           SLASET, SSYRK
00088 *     ..
00089 *     .. Intrinsic Functions ..
00090       INTRINSIC          MAX, MIN, REAL
00091 *     ..
00092 *     .. Executable Statements ..
00093 *
00094       RESID = ZERO
00095 *
00096 *     Quick return if possible
00097 *
00098       IF( M.LE.0 .OR. N.LE.0 )
00099      $   RETURN
00100 *
00101       EPS = SLAMCH( 'Precision' )
00102       IF( M.LT.N .OR. ( M.EQ.N .AND. LSAME( ROWCOL, 'R' ) ) ) THEN
00103          TRANSU = 'N'
00104          K = N
00105       ELSE
00106          TRANSU = 'T'
00107          K = M
00108       END IF
00109       MNMIN = MIN( M, N )
00110 *
00111       IF( ( MNMIN+1 )*MNMIN.LE.LWORK ) THEN
00112          LDWORK = MNMIN
00113       ELSE
00114          LDWORK = 0
00115       END IF
00116       IF( LDWORK.GT.0 ) THEN
00117 *
00118 *        Compute I - U*U' or I - U'*U.
00119 *
00120          CALL SLASET( 'Upper', MNMIN, MNMIN, ZERO, ONE, WORK, LDWORK )
00121          CALL SSYRK( 'Upper', TRANSU, MNMIN, K, -ONE, U, LDU, ONE, WORK,
00122      $               LDWORK )
00123 *
00124 *        Compute norm( I - U*U' ) / ( K * EPS ) .
00125 *
00126          RESID = SLANSY( '1', 'Upper', MNMIN, WORK, LDWORK,
00127      $           WORK( LDWORK*MNMIN+1 ) )
00128          RESID = ( RESID / REAL( K ) ) / EPS
00129       ELSE IF( TRANSU.EQ.'T' ) THEN
00130 *
00131 *        Find the maximum element in abs( I - U'*U ) / ( m * EPS )
00132 *
00133          DO 20 J = 1, N
00134             DO 10 I = 1, J
00135                IF( I.NE.J ) THEN
00136                   TMP = ZERO
00137                ELSE
00138                   TMP = ONE
00139                END IF
00140                TMP = TMP - SDOT( M, U( 1, I ), 1, U( 1, J ), 1 )
00141                RESID = MAX( RESID, ABS( TMP ) )
00142    10       CONTINUE
00143    20    CONTINUE
00144          RESID = ( RESID / REAL( M ) ) / EPS
00145       ELSE
00146 *
00147 *        Find the maximum element in abs( I - U*U' ) / ( n * EPS )
00148 *
00149          DO 40 J = 1, M
00150             DO 30 I = 1, J
00151                IF( I.NE.J ) THEN
00152                   TMP = ZERO
00153                ELSE
00154                   TMP = ONE
00155                END IF
00156                TMP = TMP - SDOT( N, U( J, 1 ), LDU, U( I, 1 ), LDU )
00157                RESID = MAX( RESID, ABS( TMP ) )
00158    30       CONTINUE
00159    40    CONTINUE
00160          RESID = ( RESID / REAL( N ) ) / EPS
00161       END IF
00162       RETURN
00163 *
00164 *     End of SORT01
00165 *
00166       END