LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zunbdb5.f
Go to the documentation of this file.
1*> \brief \b ZUNBDB5
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download ZUNBDB5 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zunbdb5.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zunbdb5.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zunbdb5.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE ZUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
20* LDQ2, WORK, LWORK, INFO )
21*
22* .. Scalar Arguments ..
23* INTEGER INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2,
24* $ N
25* ..
26* .. Array Arguments ..
27* COMPLEX*16 Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
28* ..
29*
30*
31*> \par Purpose:
32* =============
33*>
34*>\verbatim
35*>
36*> ZUNBDB5 orthogonalizes the column vector
37*> X = [ X1 ]
38*> [ X2 ]
39*> with respect to the columns of
40*> Q = [ Q1 ] .
41*> [ Q2 ]
42*> The columns of Q must be orthonormal.
43*>
44*> If the projection is zero according to Kahan's "twice is enough"
45*> criterion, then some other vector from the orthogonal complement
46*> is returned. This vector is chosen in an arbitrary but deterministic
47*> way.
48*>
49*>\endverbatim
50*
51* Arguments:
52* ==========
53*
54*> \param[in] M1
55*> \verbatim
56*> M1 is INTEGER
57*> The dimension of X1 and the number of rows in Q1. 0 <= M1.
58*> \endverbatim
59*>
60*> \param[in] M2
61*> \verbatim
62*> M2 is INTEGER
63*> The dimension of X2 and the number of rows in Q2. 0 <= M2.
64*> \endverbatim
65*>
66*> \param[in] N
67*> \verbatim
68*> N is INTEGER
69*> The number of columns in Q1 and Q2. 0 <= N.
70*> \endverbatim
71*>
72*> \param[in,out] X1
73*> \verbatim
74*> X1 is COMPLEX*16 array, dimension (M1)
75*> On entry, the top part of the vector to be orthogonalized.
76*> On exit, the top part of the projected vector.
77*> \endverbatim
78*>
79*> \param[in] INCX1
80*> \verbatim
81*> INCX1 is INTEGER
82*> Increment for entries of X1.
83*> \endverbatim
84*>
85*> \param[in,out] X2
86*> \verbatim
87*> X2 is COMPLEX*16 array, dimension (M2)
88*> On entry, the bottom part of the vector to be
89*> orthogonalized. On exit, the bottom part of the projected
90*> vector.
91*> \endverbatim
92*>
93*> \param[in] INCX2
94*> \verbatim
95*> INCX2 is INTEGER
96*> Increment for entries of X2.
97*> \endverbatim
98*>
99*> \param[in] Q1
100*> \verbatim
101*> Q1 is COMPLEX*16 array, dimension (LDQ1, N)
102*> The top part of the orthonormal basis matrix.
103*> \endverbatim
104*>
105*> \param[in] LDQ1
106*> \verbatim
107*> LDQ1 is INTEGER
108*> The leading dimension of Q1. LDQ1 >= M1.
109*> \endverbatim
110*>
111*> \param[in] Q2
112*> \verbatim
113*> Q2 is COMPLEX*16 array, dimension (LDQ2, N)
114*> The bottom part of the orthonormal basis matrix.
115*> \endverbatim
116*>
117*> \param[in] LDQ2
118*> \verbatim
119*> LDQ2 is INTEGER
120*> The leading dimension of Q2. LDQ2 >= M2.
121*> \endverbatim
122*>
123*> \param[out] WORK
124*> \verbatim
125*> WORK is COMPLEX*16 array, dimension (LWORK)
126*> \endverbatim
127*>
128*> \param[in] LWORK
129*> \verbatim
130*> LWORK is INTEGER
131*> The dimension of the array WORK. LWORK >= N.
132*> \endverbatim
133*>
134*> \param[out] INFO
135*> \verbatim
136*> INFO is INTEGER
137*> = 0: successful exit.
138*> < 0: if INFO = -i, the i-th argument had an illegal value.
139*> \endverbatim
140*
141* Authors:
142* ========
143*
144*> \author Univ. of Tennessee
145*> \author Univ. of California Berkeley
146*> \author Univ. of Colorado Denver
147*> \author NAG Ltd.
148*
149*> \ingroup unbdb5
150*
151* =====================================================================
152 SUBROUTINE zunbdb5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1,
153 $ Q2,
154 $ LDQ2, WORK, LWORK, INFO )
155*
156* -- LAPACK computational routine --
157* -- LAPACK is a software package provided by Univ. of Tennessee, --
158* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
159*
160* .. Scalar Arguments ..
161 INTEGER INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2,
162 $ N
163* ..
164* .. Array Arguments ..
165 COMPLEX*16 Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
166* ..
167*
168* =====================================================================
169*
170* .. Parameters ..
171 DOUBLE PRECISION REALZERO
172 PARAMETER ( REALZERO = 0.0d0 )
173 COMPLEX*16 ONE, ZERO
174 parameter( one = (1.0d0,0.0d0), zero = (0.0d0,0.0d0) )
175* ..
176* .. Local Scalars ..
177 INTEGER CHILDINFO, I, J
178 DOUBLE PRECISION EPS, NORM, SCL, SSQ
179* ..
180* .. External Subroutines ..
181 EXTERNAL zlassq, zunbdb6, zscal, xerbla
182* ..
183* .. External Functions ..
184 DOUBLE PRECISION DLAMCH, DZNRM2
185 EXTERNAL DLAMCH, DZNRM2
186* ..
187* .. Intrinsic Function ..
188 INTRINSIC max
189* ..
190* .. Executable Statements ..
191*
192* Test input arguments
193*
194 info = 0
195 IF( m1 .LT. 0 ) THEN
196 info = -1
197 ELSE IF( m2 .LT. 0 ) THEN
198 info = -2
199 ELSE IF( n .LT. 0 ) THEN
200 info = -3
201 ELSE IF( incx1 .LT. 1 ) THEN
202 info = -5
203 ELSE IF( incx2 .LT. 1 ) THEN
204 info = -7
205 ELSE IF( ldq1 .LT. max( 1, m1 ) ) THEN
206 info = -9
207 ELSE IF( ldq2 .LT. max( 1, m2 ) ) THEN
208 info = -11
209 ELSE IF( lwork .LT. n ) THEN
210 info = -13
211 END IF
212*
213 IF( info .NE. 0 ) THEN
214 CALL xerbla( 'ZUNBDB5', -info )
215 RETURN
216 END IF
217*
218 eps = dlamch( 'Precision' )
219*
220* Project X onto the orthogonal complement of Q if X is nonzero
221*
222 scl = realzero
223 ssq = realzero
224 CALL zlassq( m1, x1, incx1, scl, ssq )
225 CALL zlassq( m2, x2, incx2, scl, ssq )
226 norm = scl * sqrt( ssq )
227*
228 IF( norm .GT. n * eps ) THEN
229* Scale vector to unit norm to avoid problems in the caller code.
230* Computing the reciprocal is undesirable but
231* * xLASCL cannot be used because of the vector increments and
232* * the round-off error has a negligible impact on
233* orthogonalization.
234 CALL zscal( m1, one / norm, x1, incx1 )
235 CALL zscal( m2, one / norm, x2, incx2 )
236 CALL zunbdb6( m1, m2, n, x1, incx1, x2, incx2, q1, ldq1, q2,
237 $ ldq2, work, lwork, childinfo )
238*
239* If the projection is nonzero, then return
240*
241 IF( dznrm2(m1,x1,incx1) .NE. realzero
242 $ .OR. dznrm2(m2,x2,incx2) .NE. realzero ) THEN
243 RETURN
244 END IF
245 END IF
246*
247* Project each standard basis vector e_1,...,e_M1 in turn, stopping
248* when a nonzero projection is found
249*
250 DO i = 1, m1
251 DO j = 1, m1
252 x1(j) = zero
253 END DO
254 x1(i) = one
255 DO j = 1, m2
256 x2(j) = zero
257 END DO
258 CALL zunbdb6( m1, m2, n, x1, incx1, x2, incx2, q1, ldq1, q2,
259 $ ldq2, work, lwork, childinfo )
260 IF( dznrm2(m1,x1,incx1) .NE. realzero
261 $ .OR. dznrm2(m2,x2,incx2) .NE. realzero ) THEN
262 RETURN
263 END IF
264 END DO
265*
266* Project each standard basis vector e_(M1+1),...,e_(M1+M2) in turn,
267* stopping when a nonzero projection is found
268*
269 DO i = 1, m2
270 DO j = 1, m1
271 x1(j) = zero
272 END DO
273 DO j = 1, m2
274 x2(j) = zero
275 END DO
276 x2(i) = one
277 CALL zunbdb6( m1, m2, n, x1, incx1, x2, incx2, q1, ldq1, q2,
278 $ ldq2, work, lwork, childinfo )
279 IF( dznrm2(m1,x1,incx1) .NE. realzero
280 $ .OR. dznrm2(m2,x2,incx2) .NE. realzero ) THEN
281 RETURN
282 END IF
283 END DO
284*
285 RETURN
286*
287* End of ZUNBDB5
288*
289 END
290
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zlassq
subroutine zlassq(n, x, incx, scale, sumsq)
ZLASSQ updates a sum of squares represented in scaled form.
Definition zlassq.f90:122
zscal
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
zunbdb5
subroutine zunbdb5(m1, m2, n, x1, incx1, x2, incx2, q1, ldq1, q2, ldq2, work, lwork, info)
ZUNBDB5
Definition zunbdb5.f:155
zunbdb6
subroutine zunbdb6(m1, m2, n, x1, incx1, x2, incx2, q1, ldq1, q2, ldq2, work, lwork, info)
ZUNBDB6
Definition zunbdb6.f:158