 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ cunbdb6()

 subroutine cunbdb6 ( integer M1, integer M2, integer N, complex, dimension(*) X1, integer INCX1, complex, dimension(*) X2, integer INCX2, complex, dimension(ldq1,*) Q1, integer LDQ1, complex, dimension(ldq2,*) Q2, integer LDQ2, complex, dimension(*) WORK, integer LWORK, integer INFO )

CUNBDB6

Purpose:
CUNBDB6 orthogonalizes the column vector
X = [ X1 ]
[ X2 ]
with respect to the columns of
Q = [ Q1 ] .
[ Q2 ]
The columns of Q must be orthonormal.

If the projection is zero according to Kahan's "twice is enough"
criterion, then the zero vector is returned.
Parameters
 [in] M1 M1 is INTEGER The dimension of X1 and the number of rows in Q1. 0 <= M1. [in] M2 M2 is INTEGER The dimension of X2 and the number of rows in Q2. 0 <= M2. [in] N N is INTEGER The number of columns in Q1 and Q2. 0 <= N. [in,out] X1 X1 is COMPLEX array, dimension (M1) On entry, the top part of the vector to be orthogonalized. On exit, the top part of the projected vector. [in] INCX1 INCX1 is INTEGER Increment for entries of X1. [in,out] X2 X2 is COMPLEX array, dimension (M2) On entry, the bottom part of the vector to be orthogonalized. On exit, the bottom part of the projected vector. [in] INCX2 INCX2 is INTEGER Increment for entries of X2. [in] Q1 Q1 is COMPLEX array, dimension (LDQ1, N) The top part of the orthonormal basis matrix. [in] LDQ1 LDQ1 is INTEGER The leading dimension of Q1. LDQ1 >= M1. [in] Q2 Q2 is COMPLEX array, dimension (LDQ2, N) The bottom part of the orthonormal basis matrix. [in] LDQ2 LDQ2 is INTEGER The leading dimension of Q2. LDQ2 >= M2. [out] WORK WORK is COMPLEX array, dimension (LWORK) [in] LWORK LWORK is INTEGER The dimension of the array WORK. LWORK >= N. [out] INFO INFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value.

Definition at line 152 of file cunbdb6.f.

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 *
subroutine classq(n, x, incx, scl, sumsq)
CLASSQ updates a sum of squares represented in scaled form.
Definition: classq.f90:126
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
Here is the call graph for this function:
Here is the caller graph for this function: