LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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

Download CUNBDB6 + dependencies [TGZ] [ZIP] [TXT]

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. The orthogonalized vector will
 be zero if and only if it lies entirely in the range of Q.

 The projection is computed with at most two iterations of the
 classical Gram-Schmidt algorithm, see
 * L. Giraud, J. Langou, M. Rozložník. "On the round-off error
   analysis of the Gram-Schmidt algorithm with reorthogonalization."
   2002. CERFACS Technical Report No. TR/PA/02/33. URL:
   https://www.cerfacs.fr/algor/reports/2002/TR_PA_02_33.pdf
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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 157 of file cunbdb6.f.

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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
classq
subroutine classq(n, x, incx, scale, sumsq)
CLASSQ updates a sum of squares represented in scaled form.
Definition classq.f90:124
Here is the call graph for this function:
Here is the caller graph for this function: