LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sort01.f
Go to the documentation of this file.
00001 *> \brief \b SORT01
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 SORT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RESID )
00012 * 
00013 *       .. Scalar Arguments ..
00014 *       CHARACTER          ROWCOL
00015 *       INTEGER            LDU, LWORK, M, N
00016 *       REAL               RESID
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       REAL               U( LDU, * ), WORK( * )
00020 *       ..
00021 *  
00022 *
00023 *> \par Purpose:
00024 *  =============
00025 *>
00026 *> \verbatim
00027 *>
00028 *> SORT01 checks that the matrix U is orthogonal by computing the ratio
00029 *>
00030 *>    RESID = norm( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
00031 *> or
00032 *>    RESID = norm( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
00033 *>
00034 *> Alternatively, if there isn't sufficient workspace to form
00035 *> I - U*U' or I - U'*U, the ratio is computed as
00036 *>
00037 *>    RESID = abs( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
00038 *> or
00039 *>    RESID = abs( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
00040 *>
00041 *> where EPS is the machine precision.  ROWCOL is used only if m = n;
00042 *> if m > n, ROWCOL is assumed to be 'C', and if m < n, ROWCOL is
00043 *> assumed to be 'R'.
00044 *> \endverbatim
00045 *
00046 *  Arguments:
00047 *  ==========
00048 *
00049 *> \param[in] ROWCOL
00050 *> \verbatim
00051 *>          ROWCOL is CHARACTER
00052 *>          Specifies whether the rows or columns of U should be checked
00053 *>          for orthogonality.  Used only if M = N.
00054 *>          = 'R':  Check for orthogonal rows of U
00055 *>          = 'C':  Check for orthogonal columns of U
00056 *> \endverbatim
00057 *>
00058 *> \param[in] M
00059 *> \verbatim
00060 *>          M is INTEGER
00061 *>          The number of rows of the matrix U.
00062 *> \endverbatim
00063 *>
00064 *> \param[in] N
00065 *> \verbatim
00066 *>          N is INTEGER
00067 *>          The number of columns of the matrix U.
00068 *> \endverbatim
00069 *>
00070 *> \param[in] U
00071 *> \verbatim
00072 *>          U is REAL array, dimension (LDU,N)
00073 *>          The orthogonal matrix U.  U is checked for orthogonal columns
00074 *>          if m > n or if m = n and ROWCOL = 'C'.  U is checked for
00075 *>          orthogonal rows if m < n or if m = n and ROWCOL = 'R'.
00076 *> \endverbatim
00077 *>
00078 *> \param[in] LDU
00079 *> \verbatim
00080 *>          LDU is INTEGER
00081 *>          The leading dimension of the array U.  LDU >= max(1,M).
00082 *> \endverbatim
00083 *>
00084 *> \param[out] WORK
00085 *> \verbatim
00086 *>          WORK is REAL array, dimension (LWORK)
00087 *> \endverbatim
00088 *>
00089 *> \param[in] LWORK
00090 *> \verbatim
00091 *>          LWORK is INTEGER
00092 *>          The length of the array WORK.  For best performance, LWORK
00093 *>          should be at least N*(N+1) if ROWCOL = 'C' or M*(M+1) if
00094 *>          ROWCOL = 'R', but the test will be done even if LWORK is 0.
00095 *> \endverbatim
00096 *>
00097 *> \param[out] RESID
00098 *> \verbatim
00099 *>          RESID is REAL
00100 *>          RESID = norm( I - U * U' ) / ( n * EPS ), if ROWCOL = 'R', or
00101 *>          RESID = norm( I - U' * U ) / ( m * EPS ), if ROWCOL = 'C'.
00102 *> \endverbatim
00103 *
00104 *  Authors:
00105 *  ========
00106 *
00107 *> \author Univ. of Tennessee 
00108 *> \author Univ. of California Berkeley 
00109 *> \author Univ. of Colorado Denver 
00110 *> \author NAG Ltd. 
00111 *
00112 *> \date November 2011
00113 *
00114 *> \ingroup single_eig
00115 *
00116 *  =====================================================================
00117       SUBROUTINE SORT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RESID )
00118 *
00119 *  -- LAPACK test routine (version 3.4.0) --
00120 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00121 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00122 *     November 2011
00123 *
00124 *     .. Scalar Arguments ..
00125       CHARACTER          ROWCOL
00126       INTEGER            LDU, LWORK, M, N
00127       REAL               RESID
00128 *     ..
00129 *     .. Array Arguments ..
00130       REAL               U( LDU, * ), WORK( * )
00131 *     ..
00132 *
00133 *  =====================================================================
00134 *
00135 *     .. Parameters ..
00136       REAL               ZERO, ONE
00137       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00138 *     ..
00139 *     .. Local Scalars ..
00140       CHARACTER          TRANSU
00141       INTEGER            I, J, K, LDWORK, MNMIN
00142       REAL               EPS, TMP
00143 *     ..
00144 *     .. External Functions ..
00145       LOGICAL            LSAME
00146       REAL               SDOT, SLAMCH, SLANSY
00147       EXTERNAL           LSAME, SDOT, SLAMCH, SLANSY
00148 *     ..
00149 *     .. External Subroutines ..
00150       EXTERNAL           SLASET, SSYRK
00151 *     ..
00152 *     .. Intrinsic Functions ..
00153       INTRINSIC          MAX, MIN, REAL
00154 *     ..
00155 *     .. Executable Statements ..
00156 *
00157       RESID = ZERO
00158 *
00159 *     Quick return if possible
00160 *
00161       IF( M.LE.0 .OR. N.LE.0 )
00162      $   RETURN
00163 *
00164       EPS = SLAMCH( 'Precision' )
00165       IF( M.LT.N .OR. ( M.EQ.N .AND. LSAME( ROWCOL, 'R' ) ) ) THEN
00166          TRANSU = 'N'
00167          K = N
00168       ELSE
00169          TRANSU = 'T'
00170          K = M
00171       END IF
00172       MNMIN = MIN( M, N )
00173 *
00174       IF( ( MNMIN+1 )*MNMIN.LE.LWORK ) THEN
00175          LDWORK = MNMIN
00176       ELSE
00177          LDWORK = 0
00178       END IF
00179       IF( LDWORK.GT.0 ) THEN
00180 *
00181 *        Compute I - U*U' or I - U'*U.
00182 *
00183          CALL SLASET( 'Upper', MNMIN, MNMIN, ZERO, ONE, WORK, LDWORK )
00184          CALL SSYRK( 'Upper', TRANSU, MNMIN, K, -ONE, U, LDU, ONE, WORK,
00185      $               LDWORK )
00186 *
00187 *        Compute norm( I - U*U' ) / ( K * EPS ) .
00188 *
00189          RESID = SLANSY( '1', 'Upper', MNMIN, WORK, LDWORK,
00190      $           WORK( LDWORK*MNMIN+1 ) )
00191          RESID = ( RESID / REAL( K ) ) / EPS
00192       ELSE IF( TRANSU.EQ.'T' ) THEN
00193 *
00194 *        Find the maximum element in abs( I - U'*U ) / ( m * EPS )
00195 *
00196          DO 20 J = 1, N
00197             DO 10 I = 1, J
00198                IF( I.NE.J ) THEN
00199                   TMP = ZERO
00200                ELSE
00201                   TMP = ONE
00202                END IF
00203                TMP = TMP - SDOT( M, U( 1, I ), 1, U( 1, J ), 1 )
00204                RESID = MAX( RESID, ABS( TMP ) )
00205    10       CONTINUE
00206    20    CONTINUE
00207          RESID = ( RESID / REAL( M ) ) / EPS
00208       ELSE
00209 *
00210 *        Find the maximum element in abs( I - U*U' ) / ( n * EPS )
00211 *
00212          DO 40 J = 1, M
00213             DO 30 I = 1, J
00214                IF( I.NE.J ) THEN
00215                   TMP = ZERO
00216                ELSE
00217                   TMP = ONE
00218                END IF
00219                TMP = TMP - SDOT( N, U( J, 1 ), LDU, U( I, 1 ), LDU )
00220                RESID = MAX( RESID, ABS( TMP ) )
00221    30       CONTINUE
00222    40    CONTINUE
00223          RESID = ( RESID / REAL( N ) ) / EPS
00224       END IF
00225       RETURN
00226 *
00227 *     End of SORT01
00228 *
00229       END
 All Files Functions