TESTING/EIG/sort01.f

   1 *> \brief \b SORT01
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE SORT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RESID )
  12 *
  13 *       .. Scalar Arguments ..
  14 *       CHARACTER          ROWCOL
  15 *       INTEGER            LDU, LWORK, M, N
  16 *       REAL               RESID
  17 *       ..
  18 *       .. Array Arguments ..
  19 *       REAL               U( LDU, * ), WORK( * )
  20 *       ..
  21 *
  22 *
  23 *> \par Purpose:
  24 *  =============
  25 *>
  26 *> \verbatim
  27 *>
  28 *> SORT01 checks that the matrix U is orthogonal by computing the ratio
  29 *>
  30 *>    RESID = norm( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
  31 *> or
  32 *>    RESID = norm( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
  33 *>
  34 *> Alternatively, if there isn't sufficient workspace to form
  35 *> I - U*U' or I - U'*U, the ratio is computed as
  36 *>
  37 *>    RESID = abs( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
  38 *> or
  39 *>    RESID = abs( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
  40 *>
  41 *> where EPS is the machine precision.  ROWCOL is used only if m = n;
  42 *> if m > n, ROWCOL is assumed to be 'C', and if m < n, ROWCOL is
  43 *> assumed to be 'R'.
  44 *> \endverbatim
  45 *
  46 *  Arguments:
  47 *  ==========
  48 *
  49 *> \param[in] ROWCOL
  50 *> \verbatim
  51 *>          ROWCOL is CHARACTER
  52 *>          Specifies whether the rows or columns of U should be checked
  53 *>          for orthogonality.  Used only if M = N.
  54 *>          = 'R':  Check for orthogonal rows of U
  55 *>          = 'C':  Check for orthogonal columns of U
  56 *> \endverbatim
  57 *>
  58 *> \param[in] M
  59 *> \verbatim
  60 *>          M is INTEGER
  61 *>          The number of rows of the matrix U.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] N
  65 *> \verbatim
  66 *>          N is INTEGER
  67 *>          The number of columns of the matrix U.
  68 *> \endverbatim
  69 *>
  70 *> \param[in] U
  71 *> \verbatim
  72 *>          U is REAL array, dimension (LDU,N)
  73 *>          The orthogonal matrix U.  U is checked for orthogonal columns
  74 *>          if m > n or if m = n and ROWCOL = 'C'.  U is checked for
  75 *>          orthogonal rows if m < n or if m = n and ROWCOL = 'R'.
  76 *> \endverbatim
  77 *>
  78 *> \param[in] LDU
  79 *> \verbatim
  80 *>          LDU is INTEGER
  81 *>          The leading dimension of the array U.  LDU >= max(1,M).
  82 *> \endverbatim
  83 *>
  84 *> \param[out] WORK
  85 *> \verbatim
  86 *>          WORK is REAL array, dimension (LWORK)
  87 *> \endverbatim
  88 *>
  89 *> \param[in] LWORK
  90 *> \verbatim
  91 *>          LWORK is INTEGER
  92 *>          The length of the array WORK.  For best performance, LWORK
  93 *>          should be at least N*(N+1) if ROWCOL = 'C' or M*(M+1) if
  94 *>          ROWCOL = 'R', but the test will be done even if LWORK is 0.
  95 *> \endverbatim
  96 *>
  97 *> \param[out] RESID
  98 *> \verbatim
  99 *>          RESID is REAL
 100 *>          RESID = norm( I - U * U' ) / ( n * EPS ), if ROWCOL = 'R', or
 101 *>          RESID = norm( I - U' * U ) / ( m * EPS ), if ROWCOL = 'C'.
 102 *> \endverbatim
 103 *
 104 *  Authors:
 105 *  ========
 106 *
 107 *> \author Univ. of Tennessee
 108 *> \author Univ. of California Berkeley
 109 *> \author Univ. of Colorado Denver
 110 *> \author NAG Ltd.
 111 *
 112 *> \date November 2011
 113 *
 114 *> \ingroup single_eig
 115 *
 116 *  =====================================================================
 117       SUBROUTINE SORT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RESID )
 118 *
 119 *  -- LAPACK test routine (version 3.4.0) --
 120 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 121 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 122 *     November 2011
 123 *
 124 *     .. Scalar Arguments ..
 125       CHARACTER          ROWCOL
 126       INTEGER            LDU, LWORK, M, N
 127       REAL               RESID
 128 *     ..
 129 *     .. Array Arguments ..
 130       REAL               U( LDU, * ), WORK( * )
 131 *     ..
 132 *
 133 *  =====================================================================
 134 *
 135 *     .. Parameters ..
 136       REAL               ZERO, ONE
 137       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 138 *     ..
 139 *     .. Local Scalars ..
 140       CHARACTER          TRANSU
 141       INTEGER            I, J, K, LDWORK, MNMIN
 142       REAL               EPS, TMP
 143 *     ..
 144 *     .. External Functions ..
 145       LOGICAL            LSAME
 146       REAL               SDOT, SLAMCH, SLANSY
 147       EXTERNAL           LSAME, SDOT, SLAMCH, SLANSY
 148 *     ..
 149 *     .. External Subroutines ..
 150       EXTERNAL           SLASET, SSYRK
 151 *     ..
 152 *     .. Intrinsic Functions ..
 153       INTRINSIC          MAX, MIN, REAL
 154 *     ..
 155 *     .. Executable Statements ..
 156 *
 157       RESID = ZERO
 158 *
 159 *     Quick return if possible
 160 *
 161       IF( M.LE.0 .OR. N.LE.0 )
 162      $   RETURN
 163 *
 164       EPS = SLAMCH( 'Precision' )
 165       IF( M.LT.N .OR. ( M.EQ.N .AND. LSAME( ROWCOL, 'R' ) ) ) THEN
 166          TRANSU = 'N'
 167          K = N
 168       ELSE
 169          TRANSU = 'T'
 170          K = M
 171       END IF
 172       MNMIN = MIN( M, N )
 173 *
 174       IF( ( MNMIN+1 )*MNMIN.LE.LWORK ) THEN
 175          LDWORK = MNMIN
 176       ELSE
 177          LDWORK = 0
 178       END IF
 179       IF( LDWORK.GT.0 ) THEN
 180 *
 181 *        Compute I - U*U' or I - U'*U.
 182 *
 183          CALL SLASET( 'Upper', MNMIN, MNMIN, ZERO, ONE, WORK, LDWORK )
 184          CALL SSYRK( 'Upper', TRANSU, MNMIN, K, -ONE, U, LDU, ONE, WORK,
 185      $               LDWORK )
 186 *
 187 *        Compute norm( I - U*U' ) / ( K * EPS ) .
 188 *
 189          RESID = SLANSY( '1', 'Upper', MNMIN, WORK, LDWORK,
 190      $           WORK( LDWORK*MNMIN+1 ) )
 191          RESID = ( RESID / REAL( K ) ) / EPS
 192       ELSE IF( TRANSU.EQ.'T' ) THEN
 193 *
 194 *        Find the maximum element in abs( I - U'*U ) / ( m * EPS )
 195 *
 196          DO 20 J = 1, N
 197             DO 10 I = 1, J
 198                IF( I.NE.J ) THEN
 199                   TMP = ZERO
 200                ELSE
 201                   TMP = ONE
 202                END IF
 203                TMP = TMP - SDOT( M, U( 1, I ), 1, U( 1, J ), 1 )
 204                RESID = MAX( RESID, ABS( TMP ) )
 205    10       CONTINUE
 206    20    CONTINUE
 207          RESID = ( RESID / REAL( M ) ) / EPS
 208       ELSE
 209 *
 210 *        Find the maximum element in abs( I - U*U' ) / ( n * EPS )
 211 *
 212          DO 40 J = 1, M
 213             DO 30 I = 1, J
 214                IF( I.NE.J ) THEN
 215                   TMP = ZERO
 216                ELSE
 217                   TMP = ONE
 218                END IF
 219                TMP = TMP - SDOT( N, U( J, 1 ), LDU, U( I, 1 ), LDU )
 220                RESID = MAX( RESID, ABS( TMP ) )
 221    30       CONTINUE
 222    40    CONTINUE
 223          RESID = ( RESID / REAL( N ) ) / EPS
 224       END IF
 225       RETURN
 226 *
 227 *     End of SORT01
 228 *
 229       END