LAPACK 3.11.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 Euclidean norm of X must be one and 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 158 of file zunbdb6.f.

160*
161* -- LAPACK computational routine --
162* -- LAPACK is a software package provided by Univ. of Tennessee, --
163* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
164*
165* .. Scalar Arguments ..
166 INTEGER INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2,
167 $ N
168* ..
169* .. Array Arguments ..
170 COMPLEX*16 Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
171* ..
172*
173* =====================================================================
174*
175* .. Parameters ..
176 DOUBLE PRECISION ALPHA, REALONE, REALZERO
177 parameter( alpha = 0.01d0, realone = 1.0d0,
178 $ realzero = 0.0d0 )
179 COMPLEX*16 NEGONE, ONE, ZERO
180 parameter( negone = (-1.0d0,0.0d0), one = (1.0d0,0.0d0),
181 $ zero = (0.0d0,0.0d0) )
182* ..
183* .. Local Scalars ..
184 INTEGER I, IX
185 DOUBLE PRECISION EPS, NORM, NORM_NEW, SCL, SSQ
186* ..
187* .. External Functions ..
188 DOUBLE PRECISION DLAMCH
189* ..
190* .. External Subroutines ..
191 EXTERNAL zgemv, zlassq, xerbla
192* ..
193* .. Intrinsic Function ..
194 INTRINSIC max
195* ..
196* .. Executable Statements ..
197*
198* Test input arguments
199*
200 info = 0
201 IF( m1 .LT. 0 ) THEN
202 info = -1
203 ELSE IF( m2 .LT. 0 ) THEN
204 info = -2
205 ELSE IF( n .LT. 0 ) THEN
206 info = -3
207 ELSE IF( incx1 .LT. 1 ) THEN
208 info = -5
209 ELSE IF( incx2 .LT. 1 ) THEN
210 info = -7
211 ELSE IF( ldq1 .LT. max( 1, m1 ) ) THEN
212 info = -9
213 ELSE IF( ldq2 .LT. max( 1, m2 ) ) THEN
214 info = -11
215 ELSE IF( lwork .LT. n ) THEN
216 info = -13
217 END IF
218*
219 IF( info .NE. 0 ) THEN
220 CALL xerbla( 'ZUNBDB6', -info )
221 RETURN
222 END IF
223*
224 eps = dlamch( 'Precision' )
225*
226* First, project X onto the orthogonal complement of Q's column
227* space
228*
229* Christoph Conrads: In debugging mode the norm should be computed
230* and an assertion added comparing the norm with one. Alas, Fortran
231* never made it into 1989 when assert() was introduced into the C
232* programming language.
233 norm = realone
234*
235 IF( m1 .EQ. 0 ) THEN
236 DO i = 1, n
237 work(i) = zero
238 END DO
239 ELSE
240 CALL zgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
241 $ 1 )
242 END IF
243*
244 CALL zgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
245*
246 CALL zgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
247 $ incx1 )
248 CALL zgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
249 $ incx2 )
250*
251 scl = realzero
252 ssq = realzero
253 CALL zlassq( m1, x1, incx1, scl, ssq )
254 CALL zlassq( m2, x2, incx2, scl, ssq )
255 norm_new = scl * sqrt(ssq)
256*
257* If projection is sufficiently large in norm, then stop.
258* If projection is zero, then stop.
259* Otherwise, project again.
260*
261 IF( norm_new .GE. alpha * norm ) THEN
262 RETURN
263 END IF
264*
265 IF( norm_new .LE. n * eps * norm ) THEN
266 DO ix = 1, 1 + (m1-1)*incx1, incx1
267 x1( ix ) = zero
268 END DO
269 DO ix = 1, 1 + (m2-1)*incx2, incx2
270 x2( ix ) = zero
271 END DO
272 RETURN
273 END IF
274*
275 norm = norm_new
276*
277 DO i = 1, n
278 work(i) = zero
279 END DO
280*
281 IF( m1 .EQ. 0 ) THEN
282 DO i = 1, n
283 work(i) = zero
284 END DO
285 ELSE
286 CALL zgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
287 $ 1 )
288 END IF
289*
290 CALL zgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
291*
292 CALL zgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
293 $ incx1 )
294 CALL zgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
295 $ incx2 )
296*
297 scl = realzero
298 ssq = realzero
299 CALL zlassq( m1, x1, incx1, scl, ssq )
300 CALL zlassq( m2, x2, incx2, scl, ssq )
301 norm_new = scl * sqrt(ssq)
302*
303* If second projection is sufficiently large in norm, then do
304* nothing more. Alternatively, if it shrunk significantly, then
305* truncate it to zero.
306*
307 IF( norm_new .LT. alpha * norm ) THEN
308 DO ix = 1, 1 + (m1-1)*incx1, incx1
309 x1(ix) = zero
310 END DO
311 DO ix = 1, 1 + (m2-1)*incx2, incx2
312 x2(ix) = zero
313 END DO
314 END IF
315*
316 RETURN
317*
318* End of ZUNBDB6
319*
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
zlassq
subroutine zlassq(n, x, incx, scl, sumsq)
ZLASSQ updates a sum of squares represented in scaled form.
Definition: zlassq.f90:137
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
zgemv
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
Here is the call graph for this function:
Here is the caller graph for this function: