LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
sort03.f
Go to the documentation of this file.
1 *> \brief \b SORT03
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 SORT03( 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 * REAL RESULT
18 * ..
19 * .. Array Arguments ..
20 * REAL U( LDU, * ), V( LDV, * ), WORK( * )
21 * ..
22 *
23 *
24 *> \par Purpose:
25 * =============
26 *>
27 *> \verbatim
28 *>
29 *> SORT03 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 SORT03 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 SORT03 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 *> SORT03 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 REAL 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 REAL 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 REAL 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 REAL
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 single_eig
154 *
155 * =====================================================================
156  SUBROUTINE sort03( 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  REAL result
168 * ..
169 * .. Array Arguments ..
170  REAL u( ldu, * ), v( ldv, * ), work( * )
171 * ..
172 *
173 * =====================================================================
174 *
175 * .. Parameters ..
176  REAL zero, one
177  parameter( zero = 0.0e0, one = 1.0e0 )
178 * ..
179 * .. Local Scalars ..
180  INTEGER i, irc, j, lmx
181  REAL res1, res2, s, ulp
182 * ..
183 * .. External Functions ..
184  LOGICAL lsame
185  INTEGER isamax
186  REAL slamch
187  EXTERNAL lsame, isamax, slamch
188 * ..
189 * .. Intrinsic Functions ..
190  INTRINSIC abs, max, min, REAL, sign
191 * ..
192 * .. External Subroutines ..
193  EXTERNAL sort01, 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( 'SORT03', -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 = slamch( 'Precision' )
238 *
239  IF( irc.EQ.0 ) THEN
240 *
241 * Compare rows
242 *
243  res1 = zero
244  DO 20 i = 1, k
245  lmx = isamax( 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 / ( REAL( n )*ulp )
252 *
253 * Compute orthogonality of rows of V.
254 *
255  CALL sort01( '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 = isamax( 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 / ( REAL( n )*ulp )
270 *
271 * Compute orthogonality of columns of V.
272 *
273  CALL sort01( '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 SORT03
280 *
281  END