 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
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

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.```

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*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine zlassq(n, x, incx, scl, sumsq)
ZLASSQ updates a sum of squares represented in scaled form.
Definition: zlassq.f90:137
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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: