LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cunbdb6.f
Go to the documentation of this file.
1*> \brief \b CUNBDB6
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CUNBDB6 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cunbdb6.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cunbdb6.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cunbdb6.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
22* LDQ2, WORK, LWORK, INFO )
23*
24* .. Scalar Arguments ..
25* INTEGER INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2,
26* $ N
27* ..
28* .. Array Arguments ..
29* COMPLEX Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*>\verbatim
37*>
38*> CUNBDB6 orthogonalizes the column vector
39*> X = [ X1 ]
40*> [ X2 ]
41*> with respect to the columns of
42*> Q = [ Q1 ] .
43*> [ Q2 ]
44*> The columns of Q must be orthonormal. The orthogonalized vector will
45*> be zero if and only if it lies entirely in the range of Q.
46*>
47*> The projection is computed with at most two iterations of the
48*> classical Gram-Schmidt algorithm, see
49*> * L. Giraud, J. Langou, M. Rozložník. "On the round-off error
50*> analysis of the Gram-Schmidt algorithm with reorthogonalization."
51*> 2002. CERFACS Technical Report No. TR/PA/02/33. URL:
52*> https://www.cerfacs.fr/algor/reports/2002/TR_PA_02_33.pdf
53*>
54*>\endverbatim
55*
56* Arguments:
57* ==========
58*
59*> \param[in] M1
60*> \verbatim
61*> M1 is INTEGER
62*> The dimension of X1 and the number of rows in Q1. 0 <= M1.
63*> \endverbatim
64*>
65*> \param[in] M2
66*> \verbatim
67*> M2 is INTEGER
68*> The dimension of X2 and the number of rows in Q2. 0 <= M2.
69*> \endverbatim
70*>
71*> \param[in] N
72*> \verbatim
73*> N is INTEGER
74*> The number of columns in Q1 and Q2. 0 <= N.
75*> \endverbatim
76*>
77*> \param[in,out] X1
78*> \verbatim
79*> X1 is COMPLEX array, dimension (M1)
80*> On entry, the top part of the vector to be orthogonalized.
81*> On exit, the top part of the projected vector.
82*> \endverbatim
83*>
84*> \param[in] INCX1
85*> \verbatim
86*> INCX1 is INTEGER
87*> Increment for entries of X1.
88*> \endverbatim
89*>
90*> \param[in,out] X2
91*> \verbatim
92*> X2 is COMPLEX array, dimension (M2)
93*> On entry, the bottom part of the vector to be
94*> orthogonalized. On exit, the bottom part of the projected
95*> vector.
96*> \endverbatim
97*>
98*> \param[in] INCX2
99*> \verbatim
100*> INCX2 is INTEGER
101*> Increment for entries of X2.
102*> \endverbatim
103*>
104*> \param[in] Q1
105*> \verbatim
106*> Q1 is COMPLEX array, dimension (LDQ1, N)
107*> The top part of the orthonormal basis matrix.
108*> \endverbatim
109*>
110*> \param[in] LDQ1
111*> \verbatim
112*> LDQ1 is INTEGER
113*> The leading dimension of Q1. LDQ1 >= M1.
114*> \endverbatim
115*>
116*> \param[in] Q2
117*> \verbatim
118*> Q2 is COMPLEX array, dimension (LDQ2, N)
119*> The bottom part of the orthonormal basis matrix.
120*> \endverbatim
121*>
122*> \param[in] LDQ2
123*> \verbatim
124*> LDQ2 is INTEGER
125*> The leading dimension of Q2. LDQ2 >= M2.
126*> \endverbatim
127*>
128*> \param[out] WORK
129*> \verbatim
130*> WORK is COMPLEX array, dimension (LWORK)
131*> \endverbatim
132*>
133*> \param[in] LWORK
134*> \verbatim
135*> LWORK is INTEGER
136*> The dimension of the array WORK. LWORK >= N.
137*> \endverbatim
138*>
139*> \param[out] INFO
140*> \verbatim
141*> INFO is INTEGER
142*> = 0: successful exit.
143*> < 0: if INFO = -i, the i-th argument had an illegal value.
144*> \endverbatim
145*
146* Authors:
147* ========
148*
149*> \author Univ. of Tennessee
150*> \author Univ. of California Berkeley
151*> \author Univ. of Colorado Denver
152*> \author NAG Ltd.
153*
154*> \ingroup unbdb6
155*
156* =====================================================================
157 SUBROUTINE cunbdb6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
158 $ LDQ2, WORK, LWORK, INFO )
159*
160* -- LAPACK computational routine --
161* -- LAPACK is a software package provided by Univ. of Tennessee, --
162* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
163*
164* .. Scalar Arguments ..
165 INTEGER INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2,
166 $ n
167* ..
168* .. Array Arguments ..
169 COMPLEX Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
170* ..
171*
172* =====================================================================
173*
174* .. Parameters ..
175 REAL ALPHA, REALONE, REALZERO
176 parameter( alpha = 0.83e0, realone = 1.0e0,
177 $ realzero = 0.0e0 )
178 COMPLEX NEGONE, ONE, ZERO
179 parameter( negone = (-1.0e0,0.0e0), one = (1.0e0,0.0e0),
180 $ zero = (0.0e0,0.0e0) )
181* ..
182* .. Local Scalars ..
183 INTEGER I, IX
184 REAL EPS, NORM, NORM_NEW, SCL, SSQ
185* ..
186* .. External Functions ..
187 REAL SLAMCH
188* ..
189* .. External Subroutines ..
190 EXTERNAL cgemv, classq, xerbla
191* ..
192* .. Intrinsic Function ..
193 INTRINSIC max
194* ..
195* .. Executable Statements ..
196*
197* Test input arguments
198*
199 info = 0
200 IF( m1 .LT. 0 ) THEN
201 info = -1
202 ELSE IF( m2 .LT. 0 ) THEN
203 info = -2
204 ELSE IF( n .LT. 0 ) THEN
205 info = -3
206 ELSE IF( incx1 .LT. 1 ) THEN
207 info = -5
208 ELSE IF( incx2 .LT. 1 ) THEN
209 info = -7
210 ELSE IF( ldq1 .LT. max( 1, m1 ) ) THEN
211 info = -9
212 ELSE IF( ldq2 .LT. max( 1, m2 ) ) THEN
213 info = -11
214 ELSE IF( lwork .LT. n ) THEN
215 info = -13
216 END IF
217*
218 IF( info .NE. 0 ) THEN
219 CALL xerbla( 'CUNBDB6', -info )
220 RETURN
221 END IF
222*
223 eps = slamch( 'Precision' )
224*
225* Compute the Euclidean norm of X
226*
227 scl = realzero
228 ssq = realzero
229 CALL classq( m1, x1, incx1, scl, ssq )
230 CALL classq( m2, x2, incx2, scl, ssq )
231 norm = scl * sqrt( ssq )
232*
233* First, project X onto the orthogonal complement of Q's column
234* space
235*
236 IF( m1 .EQ. 0 ) THEN
237 DO i = 1, n
238 work(i) = zero
239 END DO
240 ELSE
241 CALL cgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
242 $ 1 )
243 END IF
244*
245 CALL cgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
246*
247 CALL cgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
248 $ incx1 )
249 CALL cgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
250 $ incx2 )
251*
252 scl = realzero
253 ssq = realzero
254 CALL classq( m1, x1, incx1, scl, ssq )
255 CALL classq( m2, x2, incx2, scl, ssq )
256 norm_new = scl * sqrt(ssq)
257*
258* If projection is sufficiently large in norm, then stop.
259* If projection is zero, then stop.
260* Otherwise, project again.
261*
262 IF( norm_new .GE. alpha * norm ) THEN
263 RETURN
264 END IF
265*
266 IF( norm_new .LE. n * eps * norm ) THEN
267 DO ix = 1, 1 + (m1-1)*incx1, incx1
268 x1( ix ) = zero
269 END DO
270 DO ix = 1, 1 + (m2-1)*incx2, incx2
271 x2( ix ) = zero
272 END DO
273 RETURN
274 END IF
275*
276 norm = norm_new
277*
278 DO i = 1, n
279 work(i) = zero
280 END DO
281*
282 IF( m1 .EQ. 0 ) THEN
283 DO i = 1, n
284 work(i) = zero
285 END DO
286 ELSE
287 CALL cgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
288 $ 1 )
289 END IF
290*
291 CALL cgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
292*
293 CALL cgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
294 $ incx1 )
295 CALL cgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
296 $ incx2 )
297*
298 scl = realzero
299 ssq = realzero
300 CALL classq( m1, x1, incx1, scl, ssq )
301 CALL classq( m2, x2, incx2, scl, ssq )
302 norm_new = scl * sqrt(ssq)
303*
304* If second projection is sufficiently large in norm, then do
305* nothing more. Alternatively, if it shrunk significantly, then
306* truncate it to zero.
307*
308 IF( norm_new .LT. alpha * norm ) THEN
309 DO ix = 1, 1 + (m1-1)*incx1, incx1
310 x1(ix) = zero
311 END DO
312 DO ix = 1, 1 + (m2-1)*incx2, incx2
313 x2(ix) = zero
314 END DO
315 END IF
316*
317 RETURN
318*
319* End of CUNBDB6
320*
321 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
subroutine classq(n, x, incx, scale, sumsq)
CLASSQ updates a sum of squares represented in scaled form.
Definition classq.f90:124
subroutine cunbdb6(m1, m2, n, x1, incx1, x2, incx2, q1, ldq1, q2, ldq2, work, lwork, info)
CUNBDB6
Definition cunbdb6.f:159