LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dort01.f
Go to the documentation of this file.
1 *> \brief \b DORT01
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 DORT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RESID )
12 *
13 * .. Scalar Arguments ..
14 * CHARACTER ROWCOL
15 * INTEGER LDU, LWORK, M, N
16 * DOUBLE PRECISION RESID
17 * ..
18 * .. Array Arguments ..
19 * DOUBLE PRECISION U( LDU, * ), WORK( * )
20 * ..
21 *
22 *
23 *> \par Purpose:
24 * =============
25 *>
26 *> \verbatim
27 *>
28 *> DORT01 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 double_eig
115 *
116 * =====================================================================
117  SUBROUTINE dort01( 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  DOUBLE PRECISION resid
128 * ..
129 * .. Array Arguments ..
130  DOUBLE PRECISION u( ldu, * ), work( * )
131 * ..
132 *
133 * =====================================================================
134 *
135 * .. Parameters ..
136  DOUBLE PRECISION zero, one
137  parameter( zero = 0.0d+0, one = 1.0d+0 )
138 * ..
139 * .. Local Scalars ..
140  CHARACTER transu
141  INTEGER i, j, k, ldwork, mnmin
142  DOUBLE PRECISION eps, tmp
143 * ..
144 * .. External Functions ..
145  LOGICAL lsame
146  DOUBLE PRECISION ddot, dlamch, dlansy
147  EXTERNAL lsame, ddot, dlamch, dlansy
148 * ..
149 * .. External Subroutines ..
150  EXTERNAL dlaset, dsyrk
151 * ..
152 * .. Intrinsic Functions ..
153  INTRINSIC abs, dble, max, min
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 = dlamch( '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 dlaset( 'Upper', mnmin, mnmin, zero, one, work, ldwork )
184  CALL dsyrk( 'Upper', transu, mnmin, k, -one, u, ldu, one, work,
185  $ ldwork )
186 *
187 * Compute norm( I - U*U' ) / ( K * EPS ) .
188 *
189  resid = dlansy( '1', 'Upper', mnmin, work, ldwork,
190  $ work( ldwork*mnmin+1 ) )
191  resid = ( resid / dble( 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 - ddot( m, u( 1, i ), 1, u( 1, j ), 1 )
204  resid = max( resid, abs( tmp ) )
205  10 continue
206  20 continue
207  resid = ( resid / dble( 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 - ddot( n, u( j, 1 ), ldu, u( i, 1 ), ldu )
220  resid = max( resid, abs( tmp ) )
221  30 continue
222  40 continue
223  resid = ( resid / dble( n ) ) / eps
224  END IF
225  return
226 *
227 * End of DORT01
228 *
229  END