TESTING/EIG/dort03.f

   1 *> \brief \b DORT03
   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 DORT03( RC, MU, MV, N, K, U, LDU, V, LDV, WORK, LWORK,
  12 *                          RESULT, INFO )
  13 *
  14 *       .. Scalar Arguments ..
  15 *       CHARACTER*( * )    RC
  16 *       INTEGER            INFO, K, LDU, LDV, LWORK, MU, MV, N
  17 *       DOUBLE PRECISION   RESULT
  18 *       ..
  19 *       .. Array Arguments ..
  20 *       DOUBLE PRECISION   U( LDU, * ), V( LDV, * ), WORK( * )
  21 *       ..
  22 *
  23 *
  24 *> \par Purpose:
  25 *  =============
  26 *>
  27 *> \verbatim
  28 *>
  29 *> DORT03 compares two orthogonal matrices U and V to see if their
  30 *> corresponding rows or columns span the same spaces.  The rows are
  31 *> checked if RC = 'R', and the columns are checked if RC = 'C'.
  32 *>
  33 *> RESULT is the maximum of
  34 *>
  35 *>    | V*V' - I | / ( MV ulp ), if RC = 'R', or
  36 *>
  37 *>    | V'*V - I | / ( MV ulp ), if RC = 'C',
  38 *>
  39 *> and the maximum over rows (or columns) 1 to K of
  40 *>
  41 *>    | U(i) - S*V(i) |/ ( N ulp )
  42 *>
  43 *> where S is +-1 (chosen to minimize the expression), U(i) is the i-th
  44 *> row (column) of U, and V(i) is the i-th row (column) of V.
  45 *> \endverbatim
  46 *
  47 *  Arguments:
  48 *  ==========
  49 *
  50 *> \param[in] RC
  51 *> \verbatim
  52 *>          RC is CHARACTER*1
  53 *>          If RC = 'R' the rows of U and V are to be compared.
  54 *>          If RC = 'C' the columns of U and V are to be compared.
  55 *> \endverbatim
  56 *>
  57 *> \param[in] MU
  58 *> \verbatim
  59 *>          MU is INTEGER
  60 *>          The number of rows of U if RC = 'R', and the number of
  61 *>          columns if RC = 'C'.  If MU = 0 DORT03 does nothing.
  62 *>          MU must be at least zero.
  63 *> \endverbatim
  64 *>
  65 *> \param[in] MV
  66 *> \verbatim
  67 *>          MV is INTEGER
  68 *>          The number of rows of V if RC = 'R', and the number of
  69 *>          columns if RC = 'C'.  If MV = 0 DORT03 does nothing.
  70 *>          MV must be at least zero.
  71 *> \endverbatim
  72 *>
  73 *> \param[in] N
  74 *> \verbatim
  75 *>          N is INTEGER
  76 *>          If RC = 'R', the number of columns in the matrices U and V,
  77 *>          and if RC = 'C', the number of rows in U and V.  If N = 0
  78 *>          DORT03 does nothing.  N must be at least zero.
  79 *> \endverbatim
  80 *>
  81 *> \param[in] K
  82 *> \verbatim
  83 *>          K is INTEGER
  84 *>          The number of rows or columns of U and V to compare.
  85 *>          0 <= K <= max(MU,MV).
  86 *> \endverbatim
  87 *>
  88 *> \param[in] U
  89 *> \verbatim
  90 *>          U is DOUBLE PRECISION array, dimension (LDU,N)
  91 *>          The first matrix to compare.  If RC = 'R', U is MU by N, and
  92 *>          if RC = 'C', U is N by MU.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] LDU
  96 *> \verbatim
  97 *>          LDU is INTEGER
  98 *>          The leading dimension of U.  If RC = 'R', LDU >= max(1,MU),
  99 *>          and if RC = 'C', LDU >= max(1,N).
 100 *> \endverbatim
 101 *>
 102 *> \param[in] V
 103 *> \verbatim
 104 *>          V is DOUBLE PRECISION array, dimension (LDV,N)
 105 *>          The second matrix to compare.  If RC = 'R', V is MV by N, and
 106 *>          if RC = 'C', V is N by MV.
 107 *> \endverbatim
 108 *>
 109 *> \param[in] LDV
 110 *> \verbatim
 111 *>          LDV is INTEGER
 112 *>          The leading dimension of V.  If RC = 'R', LDV >= max(1,MV),
 113 *>          and if RC = 'C', LDV >= max(1,N).
 114 *> \endverbatim
 115 *>
 116 *> \param[out] WORK
 117 *> \verbatim
 118 *>          WORK is DOUBLE PRECISION array, dimension (LWORK)
 119 *> \endverbatim
 120 *>
 121 *> \param[in] LWORK
 122 *> \verbatim
 123 *>          LWORK is INTEGER
 124 *>          The length of the array WORK.  For best performance, LWORK
 125 *>          should be at least N*N if RC = 'C' or M*M if RC = 'R', but
 126 *>          the tests will be done even if LWORK is 0.
 127 *> \endverbatim
 128 *>
 129 *> \param[out] RESULT
 130 *> \verbatim
 131 *>          RESULT is DOUBLE PRECISION
 132 *>          The value computed by the test described above.  RESULT is
 133 *>          limited to 1/ulp to avoid overflow.
 134 *> \endverbatim
 135 *>
 136 *> \param[out] INFO
 137 *> \verbatim
 138 *>          INFO is INTEGER
 139 *>          0  indicates a successful exit
 140 *>          -k indicates the k-th parameter had an illegal value
 141 *> \endverbatim
 142 *
 143 *  Authors:
 144 *  ========
 145 *
 146 *> \author Univ. of Tennessee
 147 *> \author Univ. of California Berkeley
 148 *> \author Univ. of Colorado Denver
 149 *> \author NAG Ltd.
 150 *
 151 *> \date November 2011
 152 *
 153 *> \ingroup double_eig
 154 *
 155 *  =====================================================================
 156       SUBROUTINE DORT03( RC, MU, MV, N, K, U, LDU, V, LDV, WORK, LWORK,
 157      $                   RESULT, INFO )
 158 *
 159 *  -- LAPACK test routine (version 3.4.0) --
 160 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 161 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 162 *     November 2011
 163 *
 164 *     .. Scalar Arguments ..
 165       CHARACTER*( * )    RC
 166       INTEGER            INFO, K, LDU, LDV, LWORK, MU, MV, N
 167       DOUBLE PRECISION   RESULT
 168 *     ..
 169 *     .. Array Arguments ..
 170       DOUBLE PRECISION   U( LDU, * ), V( LDV, * ), WORK( * )
 171 *     ..
 172 *
 173 *  =====================================================================
 174 *
 175 *     .. Parameters ..
 176       DOUBLE PRECISION   ZERO, ONE
 177       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 178 *     ..
 179 *     .. Local Scalars ..
 180       INTEGER            I, IRC, J, LMX
 181       DOUBLE PRECISION   RES1, RES2, S, ULP
 182 *     ..
 183 *     .. External Functions ..
 184       LOGICAL            LSAME
 185       INTEGER            IDAMAX
 186       DOUBLE PRECISION   DLAMCH
 187       EXTERNAL           LSAME, IDAMAX, DLAMCH
 188 *     ..
 189 *     .. Intrinsic Functions ..
 190       INTRINSIC          ABS, DBLE, MAX, MIN, SIGN
 191 *     ..
 192 *     .. External Subroutines ..
 193       EXTERNAL           DORT01, XERBLA
 194 *     ..
 195 *     .. Executable Statements ..
 196 *
 197 *     Check inputs
 198 *
 199       INFO = 0
 200       IF( LSAME( RC, 'R' ) ) THEN
 201          IRC = 0
 202       ELSE IF( LSAME( RC, 'C' ) ) THEN
 203          IRC = 1
 204       ELSE
 205          IRC = -1
 206       END IF
 207       IF( IRC.EQ.-1 ) THEN
 208          INFO = -1
 209       ELSE IF( MU.LT.0 ) THEN
 210          INFO = -2
 211       ELSE IF( MV.LT.0 ) THEN
 212          INFO = -3
 213       ELSE IF( N.LT.0 ) THEN
 214          INFO = -4
 215       ELSE IF( K.LT.0 .OR. K.GT.MAX( MU, MV ) ) THEN
 216          INFO = -5
 217       ELSE IF( ( IRC.EQ.0 .AND. LDU.LT.MAX( 1, MU ) ) .OR.
 218      $         ( IRC.EQ.1 .AND. LDU.LT.MAX( 1, N ) ) ) THEN
 219          INFO = -7
 220       ELSE IF( ( IRC.EQ.0 .AND. LDV.LT.MAX( 1, MV ) ) .OR.
 221      $         ( IRC.EQ.1 .AND. LDV.LT.MAX( 1, N ) ) ) THEN
 222          INFO = -9
 223       END IF
 224       IF( INFO.NE.0 ) THEN
 225          CALL XERBLA( 'DORT03', -INFO )
 226          RETURN
 227       END IF
 228 *
 229 *     Initialize result
 230 *
 231       RESULT = ZERO
 232       IF( MU.EQ.0 .OR. MV.EQ.0 .OR. N.EQ.0 )
 233      $   RETURN
 234 *
 235 *     Machine constants
 236 *
 237       ULP = DLAMCH( 'Precision' )
 238 *
 239       IF( IRC.EQ.0 ) THEN
 240 *
 241 *        Compare rows
 242 *
 243          RES1 = ZERO
 244          DO 20 I = 1, K
 245             LMX = IDAMAX( N, U( I, 1 ), LDU )
 246             S = SIGN( ONE, U( I, LMX ) )*SIGN( ONE, V( I, LMX ) )
 247             DO 10 J = 1, N
 248                RES1 = MAX( RES1, ABS( U( I, J )-S*V( I, J ) ) )
 249    10       CONTINUE
 250    20    CONTINUE
 251          RES1 = RES1 / ( DBLE( N )*ULP )
 252 *
 253 *        Compute orthogonality of rows of V.
 254 *
 255          CALL DORT01( 'Rows', MV, N, V, LDV, WORK, LWORK, RES2 )
 256 *
 257       ELSE
 258 *
 259 *        Compare columns
 260 *
 261          RES1 = ZERO
 262          DO 40 I = 1, K
 263             LMX = IDAMAX( N, U( 1, I ), 1 )
 264             S = SIGN( ONE, U( LMX, I ) )*SIGN( ONE, V( LMX, I ) )
 265             DO 30 J = 1, N
 266                RES1 = MAX( RES1, ABS( U( J, I )-S*V( J, I ) ) )
 267    30       CONTINUE
 268    40    CONTINUE
 269          RES1 = RES1 / ( DBLE( N )*ULP )
 270 *
 271 *        Compute orthogonality of columns of V.
 272 *
 273          CALL DORT01( 'Columns', N, MV, V, LDV, WORK, LWORK, RES2 )
 274       END IF
 275 *
 276       RESULT = MIN( MAX( RES1, RES2 ), ONE / ULP )
 277       RETURN
 278 *
 279 *     End of DORT03
 280 *
 281       END