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

◆ zunbdb6()

subroutine zunbdb6 ( integer  m1,
integer  m2,
integer  n,
complex*16, dimension(*)  x1,
integer  incx1,
complex*16, dimension(*)  x2,
integer  incx2,
complex*16, dimension(ldq1,*)  q1,
integer  ldq1,
complex*16, dimension(ldq2,*)  q2,
integer  ldq2,
complex*16, dimension(*)  work,
integer  lwork,
integer  info 
)

ZUNBDB6

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

Purpose:
 ZUNBDB6 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*16 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*16 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*16 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*16 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*16 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 zunbdb6.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*16 Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
170* ..
171*
172* =====================================================================
173*
174* .. Parameters ..
175 DOUBLE PRECISION ALPHA, REALONE, REALZERO
176 parameter( alpha = 0.83d0, realone = 1.0d0,
177 $ realzero = 0.0d0 )
178 COMPLEX*16 NEGONE, ONE, ZERO
179 parameter( negone = (-1.0d0,0.0d0), one = (1.0d0,0.0d0),
180 $ zero = (0.0d0,0.0d0) )
181* ..
182* .. Local Scalars ..
183 INTEGER I, IX
184 DOUBLE PRECISION EPS, NORM, NORM_NEW, SCL, SSQ
185* ..
186* .. External Functions ..
187 DOUBLE PRECISION DLAMCH
188* ..
189* .. External Subroutines ..
190 EXTERNAL zgemv, zlassq, 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( 'ZUNBDB6', -info )
220 RETURN
221 END IF
222*
223 eps = dlamch( 'Precision' )
224*
225* Compute the Euclidean norm of X
226*
227 scl = realzero
228 ssq = realzero
229 CALL zlassq( m1, x1, incx1, scl, ssq )
230 CALL zlassq( 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 zgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
242 $ 1 )
243 END IF
244*
245 CALL zgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
246*
247 CALL zgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
248 $ incx1 )
249 CALL zgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
250 $ incx2 )
251*
252 scl = realzero
253 ssq = realzero
254 CALL zlassq( m1, x1, incx1, scl, ssq )
255 CALL zlassq( 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 zgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
288 $ 1 )
289 END IF
290*
291 CALL zgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
292*
293 CALL zgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
294 $ incx1 )
295 CALL zgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
296 $ incx2 )
297*
298 scl = realzero
299 ssq = realzero
300 CALL zlassq( m1, x1, incx1, scl, ssq )
301 CALL zlassq( 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 ZUNBDB6
320*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
zlassq
subroutine zlassq(n, x, incx, scale, sumsq)
ZLASSQ updates a sum of squares represented in scaled form.
Definition zlassq.f90:124
Here is the call graph for this function:
Here is the caller graph for this function: