1       SUBROUTINE SORT03( RC, MU, MV, N, K, U, LDU, V, LDV, WORK, LWORK,
  2      $                   RESULT, INFO )
  3 *
  4 *  -- LAPACK test routine (version 3.1) --
  5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  6 *     November 2006
  7 *
  8 *     .. Scalar Arguments ..
  9       CHARACTER** )    RC
 10       INTEGER            INFO, K, LDU, LDV, LWORK, MU, MV, N
 11       REAL               RESULT
 12 *     ..
 13 *     .. Array Arguments ..
 14       REAL               U( LDU, * ), V( LDV, * ), WORK( * )
 15 *     ..
 16 *
 17 *  Purpose
 18 *  =======
 19 *
 20 *  SORT03 compares two orthogonal matrices U and V to see if their
 21 *  corresponding rows or columns span the same spaces.  The rows are
 22 *  checked if RC = 'R', and the columns are checked if RC = 'C'.
 23 *
 24 *  RESULT is the maximum of
 25 *
 26 *     | V*V' - I | / ( MV ulp ), if RC = 'R', or
 27 *
 28 *     | V'*V - I | / ( MV ulp ), if RC = 'C',
 29 *
 30 *  and the maximum over rows (or columns) 1 to K of
 31 *
 32 *     | U(i) - S*V(i) |/ ( N ulp )
 33 *
 34 *  where S is +-1 (chosen to minimize the expression), U(i) is the i-th
 35 *  row (column) of U, and V(i) is the i-th row (column) of V.
 36 *
 37 *  Arguments
 38 *  ==========
 39 *
 40 *  RC      (input) CHARACTER*1
 41 *          If RC = 'R' the rows of U and V are to be compared.
 42 *          If RC = 'C' the columns of U and V are to be compared.
 43 *
 44 *  MU      (input) INTEGER
 45 *          The number of rows of U if RC = 'R', and the number of
 46 *          columns if RC = 'C'.  If MU = 0 SORT03 does nothing.
 47 *          MU must be at least zero.
 48 *
 49 *  MV      (input) INTEGER
 50 *          The number of rows of V if RC = 'R', and the number of
 51 *          columns if RC = 'C'.  If MV = 0 SORT03 does nothing.
 52 *          MV must be at least zero.
 53 *
 54 *  N       (input) INTEGER
 55 *          If RC = 'R', the number of columns in the matrices U and V,
 56 *          and if RC = 'C', the number of rows in U and V.  If N = 0
 57 *          SORT03 does nothing.  N must be at least zero.
 58 *
 59 *  K       (input) INTEGER
 60 *          The number of rows or columns of U and V to compare.
 61 *          0 <= K <= max(MU,MV).
 62 *
 63 *  U       (input) REAL array, dimension (LDU,N)
 64 *          The first matrix to compare.  If RC = 'R', U is MU by N, and
 65 *          if RC = 'C', U is N by MU.
 66 *
 67 *  LDU     (input) INTEGER
 68 *          The leading dimension of U.  If RC = 'R', LDU >= max(1,MU),
 69 *          and if RC = 'C', LDU >= max(1,N).
 70 *
 71 *  V       (input) REAL array, dimension (LDV,N)
 72 *          The second matrix to compare.  If RC = 'R', V is MV by N, and
 73 *          if RC = 'C', V is N by MV.
 74 *
 75 *  LDV     (input) INTEGER
 76 *          The leading dimension of V.  If RC = 'R', LDV >= max(1,MV),
 77 *          and if RC = 'C', LDV >= max(1,N).
 78 *
 79 *  WORK    (workspace) REAL array, dimension (LWORK)
 80 *
 81 *  LWORK   (input) INTEGER
 82 *          The length of the array WORK.  For best performance, LWORK
 83 *          should be at least N*N if RC = 'C' or M*M if RC = 'R', but
 84 *          the tests will be done even if LWORK is 0.
 85 *
 86 *  RESULT  (output) REAL
 87 *          The value computed by the test described above.  RESULT is
 88 *          limited to 1/ulp to avoid overflow.
 89 *
 90 *  INFO    (output) INTEGER
 91 *          0  indicates a successful exit
 92 *          -k indicates the k-th parameter had an illegal value
 93 *
 94 *  =====================================================================
 95 *
 96 *     .. Parameters ..
 97       REAL               ZERO, ONE
 98       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
 99 *     ..
100 *     .. Local Scalars ..
101       INTEGER            I, IRC, J, LMX
102       REAL               RES1, RES2, S, ULP
103 *     ..
104 *     .. External Functions ..
105       LOGICAL            LSAME
106       INTEGER            ISAMAX
107       REAL               SLAMCH
108       EXTERNAL           LSAME, ISAMAX, SLAMCH
109 *     ..
110 *     .. Intrinsic Functions ..
111       INTRINSIC          ABSMAXMIN, REAL, SIGN
112 *     ..
113 *     .. External Subroutines ..
114       EXTERNAL           SORT01, XERBLA
115 *     ..
116 *     .. Executable Statements ..
117 *
118 *     Check inputs
119 *
120       INFO = 0
121       IF( LSAME( RC, 'R' ) ) THEN
122          IRC = 0
123       ELSE IF( LSAME( RC, 'C' ) ) THEN
124          IRC = 1
125       ELSE
126          IRC = -1
127       END IF
128       IF( IRC.EQ.-1 ) THEN
129          INFO = -1
130       ELSE IF( MU.LT.0 ) THEN
131          INFO = -2
132       ELSE IF( MV.LT.0 ) THEN
133          INFO = -3
134       ELSE IF( N.LT.0 ) THEN
135          INFO = -4
136       ELSE IF( K.LT.0 .OR. K.GT.MAX( MU, MV ) ) THEN
137          INFO = -5
138       ELSE IF( ( IRC.EQ.0 .AND. LDU.LT.MAX1, MU ) ) .OR.
139      $         ( IRC.EQ.1 .AND. LDU.LT.MAX1, N ) ) ) THEN
140          INFO = -7
141       ELSE IF( ( IRC.EQ.0 .AND. LDV.LT.MAX1, MV ) ) .OR.
142      $         ( IRC.EQ.1 .AND. LDV.LT.MAX1, N ) ) ) THEN
143          INFO = -9
144       END IF
145       IF( INFO.NE.0 ) THEN
146          CALL XERBLA( 'SORT03'-INFO )
147          RETURN
148       END IF
149 *
150 *     Initialize result
151 *
152       RESULT = ZERO
153       IF( MU.EQ.0 .OR. MV.EQ.0 .OR. N.EQ.0 )
154      $   RETURN
155 *
156 *     Machine constants
157 *
158       ULP = SLAMCH( 'Precision' )
159 *
160       IF( IRC.EQ.0 ) THEN
161 *
162 *        Compare rows
163 *
164          RES1 = ZERO
165          DO 20 I = 1, K
166             LMX = ISAMAX( N, U( I, 1 ), LDU )
167             S = SIGN( ONE, U( I, LMX ) )*SIGN( ONE, V( I, LMX ) )
168             DO 10 J = 1, N
169                RES1 = MAX( RES1, ABS( U( I, J )-S*V( I, J ) ) )
170    10       CONTINUE
171    20    CONTINUE
172          RES1 = RES1 / ( REAL( N )*ULP )
173 *
174 *        Compute orthogonality of rows of V.
175 *
176          CALL SORT01( 'Rows', MV, N, V, LDV, WORK, LWORK, RES2 )
177 *
178       ELSE
179 *
180 *        Compare columns
181 *
182          RES1 = ZERO
183          DO 40 I = 1, K
184             LMX = ISAMAX( N, U( 1, I ), 1 )
185             S = SIGN( ONE, U( LMX, I ) )*SIGN( ONE, V( LMX, I ) )
186             DO 30 J = 1, N
187                RES1 = MAX( RES1, ABS( U( J, I )-S*V( J, I ) ) )
188    30       CONTINUE
189    40    CONTINUE
190          RES1 = RES1 / ( REAL( N )*ULP )
191 *
192 *        Compute orthogonality of columns of V.
193 *
194          CALL SORT01( 'Columns', N, MV, V, LDV, WORK, LWORK, RES2 )
195       END IF
196 *
197       RESULT = MINMAX( RES1, RES2 ), ONE / ULP )
198       RETURN
199 *
200 *     End of SORT03
201 *
202       END