LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
zunbdb6.f
Go to the documentation of this file.
1 *> \brief \b ZUNBDB6
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZUNBDB6 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zunbdb6.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zunbdb6.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zunbdb6.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZUNBDB6( 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*16 Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *>\verbatim
37 *>
38 *> ZUNBDB6 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.
45 *>
46 *> If the projection is zero according to Kahan's "twice is enough"
47 *> criterion, then the zero vector is returned.
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 *> \date July 2012
150 *
151 *> \ingroup complex16OTHERcomputational
152 *
153 * =====================================================================
154  SUBROUTINE zunbdb6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
155  $ LDQ2, WORK, LWORK, INFO )
156 *
157 * -- LAPACK computational routine (version 3.7.1) --
158 * -- LAPACK is a software package provided by Univ. of Tennessee, --
159 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
160 * July 2012
161 *
162 * .. Scalar Arguments ..
163  INTEGER INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2,
164  $ n
165 * ..
166 * .. Array Arguments ..
167  COMPLEX*16 Q1(ldq1,*), Q2(ldq2,*), WORK(*), X1(*), X2(*)
168 * ..
169 *
170 * =====================================================================
171 *
172 * .. Parameters ..
173  DOUBLE PRECISION ALPHASQ, REALONE, REALZERO
174  parameter( alphasq = 0.01d0, realone = 1.0d0,
175  $ realzero = 0.0d0 )
176  COMPLEX*16 NEGONE, ONE, ZERO
177  parameter( negone = (-1.0d0,0.0d0), one = (1.0d0,0.0d0),
178  $ zero = (0.0d0,0.0d0) )
179 * ..
180 * .. Local Scalars ..
181  INTEGER I
182  DOUBLE PRECISION NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2
183 * ..
184 * .. External Subroutines ..
185  EXTERNAL zgemv, zlassq, xerbla
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( 'ZUNBDB6', -info )
215  RETURN
216  END IF
217 *
218 * First, project X onto the orthogonal complement of Q's column
219 * space
220 *
221  scl1 = realzero
222  ssq1 = realone
223  CALL zlassq( m1, x1, incx1, scl1, ssq1 )
224  scl2 = realzero
225  ssq2 = realone
226  CALL zlassq( m2, x2, incx2, scl2, ssq2 )
227  normsq1 = scl1**2*ssq1 + scl2**2*ssq2
228 *
229  IF( m1 .EQ. 0 ) THEN
230  DO i = 1, n
231  work(i) = zero
232  END DO
233  ELSE
234  CALL zgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
235  $ 1 )
236  END IF
237 *
238  CALL zgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
239 *
240  CALL zgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
241  $ incx1 )
242  CALL zgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
243  $ incx2 )
244 *
245  scl1 = realzero
246  ssq1 = realone
247  CALL zlassq( m1, x1, incx1, scl1, ssq1 )
248  scl2 = realzero
249  ssq2 = realone
250  CALL zlassq( m2, x2, incx2, scl2, ssq2 )
251  normsq2 = scl1**2*ssq1 + scl2**2*ssq2
252 *
253 * If projection is sufficiently large in norm, then stop.
254 * If projection is zero, then stop.
255 * Otherwise, project again.
256 *
257  IF( normsq2 .GE. alphasq*normsq1 ) THEN
258  RETURN
259  END IF
260 *
261  IF( normsq2 .EQ. zero ) THEN
262  RETURN
263  END IF
264 *
265  normsq1 = normsq2
266 *
267  DO i = 1, n
268  work(i) = zero
269  END DO
270 *
271  IF( m1 .EQ. 0 ) THEN
272  DO i = 1, n
273  work(i) = zero
274  END DO
275  ELSE
276  CALL zgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
277  $ 1 )
278  END IF
279 *
280  CALL zgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
281 *
282  CALL zgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
283  $ incx1 )
284  CALL zgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
285  $ incx2 )
286 *
287  scl1 = realzero
288  ssq1 = realone
289  CALL zlassq( m1, x1, incx1, scl1, ssq1 )
290  scl2 = realzero
291  ssq2 = realone
292  CALL zlassq( m1, x1, incx1, scl1, ssq1 )
293  normsq2 = scl1**2*ssq1 + scl2**2*ssq2
294 *
295 * If second projection is sufficiently large in norm, then do
296 * nothing more. Alternatively, if it shrunk significantly, then
297 * truncate it to zero.
298 *
299  IF( normsq2 .LT. alphasq*normsq1 ) THEN
300  DO i = 1, m1
301  x1(i) = zero
302  END DO
303  DO i = 1, m2
304  x2(i) = zero
305  END DO
306  END IF
307 *
308  RETURN
309 *
310 * End of ZUNBDB6
311 *
312  END
313 
subroutine zlassq(N, X, INCX, SCALE, SUMSQ)
ZLASSQ updates a sum of squares represented in scaled form.
Definition: zlassq.f:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zunbdb6(M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, LDQ2, WORK, LWORK, INFO)
ZUNBDB6
Definition: zunbdb6.f:156
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160