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