LAPACK 3.12.1
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. 
!>   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 155 of file zunbdb6.f.

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