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