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.
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 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 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 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 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 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 complexOTHERcomputational
150*
151*  =====================================================================
152      SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
153     $                    LDQ2, WORK, LWORK, INFO )
154*
155*  -- LAPACK computational routine --
156*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
157*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
158*
159*     .. Scalar Arguments ..
160      INTEGER            INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2,
161     $                   N
162*     ..
163*     .. Array Arguments ..
164      COMPLEX            Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
165*     ..
166*
167*  =====================================================================
168*
169*     .. Parameters ..
170      REAL               ALPHASQ, REALONE, REALZERO
171      PARAMETER          ( ALPHASQ = 0.01E0, REALONE = 1.0E0,
172     $                     REALZERO = 0.0E0 )
173      COMPLEX            NEGONE, ONE, ZERO
174      PARAMETER          ( NEGONE = (-1.0E0,0.0E0), ONE = (1.0E0,0.0E0),
175     $                     ZERO = (0.0E0,0.0E0) )
176*     ..
177*     .. Local Scalars ..
178      INTEGER            I
179      REAL               NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2
180*     ..
181*     .. External Subroutines ..
182      EXTERNAL           CGEMV, CLASSQ, XERBLA
183*     ..
184*     .. Intrinsic Function ..
185      INTRINSIC          MAX
186*     ..
187*     .. Executable Statements ..
188*
189*     Test input arguments
190*
191      INFO = 0
192      IF( M1 .LT. 0 ) THEN
193         INFO = -1
194      ELSE IF( M2 .LT. 0 ) THEN
195         INFO = -2
196      ELSE IF( N .LT. 0 ) THEN
197         INFO = -3
198      ELSE IF( INCX1 .LT. 1 ) THEN
199         INFO = -5
200      ELSE IF( INCX2 .LT. 1 ) THEN
201         INFO = -7
202      ELSE IF( LDQ1 .LT. MAX( 1, M1 ) ) THEN
203         INFO = -9
204      ELSE IF( LDQ2 .LT. MAX( 1, M2 ) ) THEN
205         INFO = -11
206      ELSE IF( LWORK .LT. N ) THEN
207         INFO = -13
208      END IF
209*
210      IF( INFO .NE. 0 ) THEN
211         CALL XERBLA( 'CUNBDB6', -INFO )
212         RETURN
213      END IF
214*
215*     First, project X onto the orthogonal complement of Q's column
216*     space
217*
218      SCL1 = REALZERO
219      SSQ1 = REALONE
220      CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
221      SCL2 = REALZERO
222      SSQ2 = REALONE
223      CALL CLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
224      NORMSQ1 = SCL1**2*SSQ1 + SCL2**2*SSQ2
225*
226      IF( M1 .EQ. 0 ) THEN
227         DO I = 1, N
228            WORK(I) = ZERO
229         END DO
230      ELSE
231         CALL CGEMV( 'C', M1, N, ONE, Q1, LDQ1, X1, INCX1, ZERO, WORK,
232     $               1 )
233      END IF
234*
235      CALL CGEMV( 'C', M2, N, ONE, Q2, LDQ2, X2, INCX2, ONE, WORK, 1 )
236*
237      CALL CGEMV( 'N', M1, N, NEGONE, Q1, LDQ1, WORK, 1, ONE, X1,
238     $            INCX1 )
239      CALL CGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
240     $            INCX2 )
241*
242      SCL1 = REALZERO
243      SSQ1 = REALONE
244      CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
245      SCL2 = REALZERO
246      SSQ2 = REALONE
247      CALL CLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
248      NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
249*
250*     If projection is sufficiently large in norm, then stop.
251*     If projection is zero, then stop.
252*     Otherwise, project again.
253*
254      IF( NORMSQ2 .GE. ALPHASQ*NORMSQ1 ) THEN
255         RETURN
256      END IF
257*
258      IF( NORMSQ2 .EQ. ZERO ) THEN
259         RETURN
260      END IF
261*
262      NORMSQ1 = NORMSQ2
263*
264      DO I = 1, N
265         WORK(I) = ZERO
266      END DO
267*
268      IF( M1 .EQ. 0 ) THEN
269         DO I = 1, N
270            WORK(I) = ZERO
271         END DO
272      ELSE
273         CALL CGEMV( 'C', M1, N, ONE, Q1, LDQ1, X1, INCX1, ZERO, WORK,
274     $               1 )
275      END IF
276*
277      CALL CGEMV( 'C', M2, N, ONE, Q2, LDQ2, X2, INCX2, ONE, WORK, 1 )
278*
279      CALL CGEMV( 'N', M1, N, NEGONE, Q1, LDQ1, WORK, 1, ONE, X1,
280     $            INCX1 )
281      CALL CGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
282     $            INCX2 )
283*
284      SCL1 = REALZERO
285      SSQ1 = REALONE
286      CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
287      SCL2 = REALZERO
288      SSQ2 = REALONE
289      CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
290      NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
291*
292*     If second projection is sufficiently large in norm, then do
293*     nothing more. Alternatively, if it shrunk significantly, then
294*     truncate it to zero.
295*
296      IF( NORMSQ2 .LT. ALPHASQ*NORMSQ1 ) THEN
297         DO I = 1, M1
298            X1(I) = ZERO
299         END DO
300         DO I = 1, M2
301            X2(I) = ZERO
302         END DO
303      END IF
304*
305      RETURN
306*
307*     End of CUNBDB6
308*
309      END
310
311