OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
zgbbrd.f
Go to the documentation of this file.
1*> \brief \b ZGBBRD
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZGBBRD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgbbrd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgbbrd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgbbrd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZGBBRD( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
22* LDQ, PT, LDPT, C, LDC, WORK, RWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER VECT
26* INTEGER INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
27* ..
28* .. Array Arguments ..
29* DOUBLE PRECISION D( * ), E( * ), RWORK( * )
30* COMPLEX*16 AB( LDAB, * ), C( LDC, * ), PT( LDPT, * ),
31* $ Q( LDQ, * ), WORK( * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> ZGBBRD reduces a complex general m-by-n band matrix A to real upper
41*> bidiagonal form B by a unitary transformation: Q**H * A * P = B.
42*>
43*> The routine computes B, and optionally forms Q or P**H, or computes
44*> Q**H*C for a given matrix C.
45*> \endverbatim
46*
47* Arguments:
48* ==========
49*
50*> \param[in] VECT
51*> \verbatim
52*> VECT is CHARACTER*1
53*> Specifies whether or not the matrices Q and P**H are to be
54*> formed.
55*> = 'N': do not form Q or P**H;
56*> = 'Q': form Q only;
57*> = 'P': form P**H only;
58*> = 'B': form both.
59*> \endverbatim
60*>
61*> \param[in] M
62*> \verbatim
63*> M is INTEGER
64*> The number of rows of the matrix A. M >= 0.
65*> \endverbatim
66*>
67*> \param[in] N
68*> \verbatim
69*> N is INTEGER
70*> The number of columns of the matrix A. N >= 0.
71*> \endverbatim
72*>
73*> \param[in] NCC
74*> \verbatim
75*> NCC is INTEGER
76*> The number of columns of the matrix C. NCC >= 0.
77*> \endverbatim
78*>
79*> \param[in] KL
80*> \verbatim
81*> KL is INTEGER
82*> The number of subdiagonals of the matrix A. KL >= 0.
83*> \endverbatim
84*>
85*> \param[in] KU
86*> \verbatim
87*> KU is INTEGER
88*> The number of superdiagonals of the matrix A. KU >= 0.
89*> \endverbatim
90*>
91*> \param[in,out] AB
92*> \verbatim
93*> AB is COMPLEX*16 array, dimension (LDAB,N)
94*> On entry, the m-by-n band matrix A, stored in rows 1 to
95*> KL+KU+1. The j-th column of A is stored in the j-th column of
96*> the array AB as follows:
97*> AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
98*> On exit, A is overwritten by values generated during the
99*> reduction.
100*> \endverbatim
101*>
102*> \param[in] LDAB
103*> \verbatim
104*> LDAB is INTEGER
105*> The leading dimension of the array A. LDAB >= KL+KU+1.
106*> \endverbatim
107*>
108*> \param[out] D
109*> \verbatim
110*> D is DOUBLE PRECISION array, dimension (min(M,N))
111*> The diagonal elements of the bidiagonal matrix B.
112*> \endverbatim
113*>
114*> \param[out] E
115*> \verbatim
116*> E is DOUBLE PRECISION array, dimension (min(M,N)-1)
117*> The superdiagonal elements of the bidiagonal matrix B.
118*> \endverbatim
119*>
120*> \param[out] Q
121*> \verbatim
122*> Q is COMPLEX*16 array, dimension (LDQ,M)
123*> If VECT = 'Q' or 'B', the m-by-m unitary matrix Q.
124*> If VECT = 'N' or 'P', the array Q is not referenced.
125*> \endverbatim
126*>
127*> \param[in] LDQ
128*> \verbatim
129*> LDQ is INTEGER
130*> The leading dimension of the array Q.
131*> LDQ >= max(1,M) if VECT = 'Q' or 'B'; LDQ >= 1 otherwise.
132*> \endverbatim
133*>
134*> \param[out] PT
135*> \verbatim
136*> PT is COMPLEX*16 array, dimension (LDPT,N)
137*> If VECT = 'P' or 'B', the n-by-n unitary matrix P'.
138*> If VECT = 'N' or 'Q', the array PT is not referenced.
139*> \endverbatim
140*>
141*> \param[in] LDPT
142*> \verbatim
143*> LDPT is INTEGER
144*> The leading dimension of the array PT.
145*> LDPT >= max(1,N) if VECT = 'P' or 'B'; LDPT >= 1 otherwise.
146*> \endverbatim
147*>
148*> \param[in,out] C
149*> \verbatim
150*> C is COMPLEX*16 array, dimension (LDC,NCC)
151*> On entry, an m-by-ncc matrix C.
152*> On exit, C is overwritten by Q**H*C.
153*> C is not referenced if NCC = 0.
154*> \endverbatim
155*>
156*> \param[in] LDC
157*> \verbatim
158*> LDC is INTEGER
159*> The leading dimension of the array C.
160*> LDC >= max(1,M) if NCC > 0; LDC >= 1 if NCC = 0.
161*> \endverbatim
162*>
163*> \param[out] WORK
164*> \verbatim
165*> WORK is COMPLEX*16 array, dimension (max(M,N))
166*> \endverbatim
167*>
168*> \param[out] RWORK
169*> \verbatim
170*> RWORK is DOUBLE PRECISION array, dimension (max(M,N))
171*> \endverbatim
172*>
173*> \param[out] INFO
174*> \verbatim
175*> INFO is INTEGER
176*> = 0: successful exit.
177*> < 0: if INFO = -i, the i-th argument had an illegal value.
178*> \endverbatim
179*
180* Authors:
181* ========
182*
183*> \author Univ. of Tennessee
184*> \author Univ. of California Berkeley
185*> \author Univ. of Colorado Denver
186*> \author NAG Ltd.
187*
188*> \ingroup complex16GBcomputational
189*
190* =====================================================================
191 SUBROUTINE zgbbrd( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
192 $ LDQ, PT, LDPT, C, LDC, WORK, RWORK, INFO )
193*
194* -- LAPACK computational routine --
195* -- LAPACK is a software package provided by Univ. of Tennessee, --
196* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
197*
198* .. Scalar Arguments ..
199 CHARACTER VECT
200 INTEGER INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
201* ..
202* .. Array Arguments ..
203 DOUBLE PRECISION D( * ), E( * ), RWORK( * )
204 COMPLEX*16 AB( LDAB, * ), C( LDC, * ), PT( LDPT, * ),
205 $ q( ldq, * ), work( * )
206* ..
207*
208* =====================================================================
209*
210* .. Parameters ..
211 DOUBLE PRECISION ZERO
212 parameter( zero = 0.0d+0 )
213 COMPLEX*16 CZERO, CONE
214 parameter( czero = ( 0.0d+0, 0.0d+0 ),
215 $ cone = ( 1.0d+0, 0.0d+0 ) )
216* ..
217* .. Local Scalars ..
218 LOGICAL WANTB, WANTC, WANTPT, WANTQ
219 INTEGER I, INCA, J, J1, J2, KB, KB1, KK, KLM, KLU1,
220 $ kun, l, minmn, ml, ml0, mu, mu0, nr, nrt
221 DOUBLE PRECISION ABST, RC
222 COMPLEX*16 RA, RB, RS, T
223* ..
224* .. External Subroutines ..
225 EXTERNAL xerbla, zlargv, zlartg, zlartv, zlaset, zrot,
226 $ zscal
227* ..
228* .. Intrinsic Functions ..
229 INTRINSIC abs, dconjg, max, min
230* ..
231* .. External Functions ..
232 LOGICAL LSAME
233 EXTERNAL lsame
234* ..
235* .. Executable Statements ..
236*
237* Test the input parameters
238*
239 wantb = lsame( vect, 'B' )
240 wantq = lsame( vect, 'Q' ) .OR. wantb
241 wantpt = lsame( vect, 'P' ) .OR. wantb
242 wantc = ncc.GT.0
243 klu1 = kl + ku + 1
244 info = 0
245 IF( .NOT.wantq .AND. .NOT.wantpt .AND. .NOT.lsame( vect, 'N' ) )
246 $ THEN
247 info = -1
248 ELSE IF( m.LT.0 ) THEN
249 info = -2
250 ELSE IF( n.LT.0 ) THEN
251 info = -3
252 ELSE IF( ncc.LT.0 ) THEN
253 info = -4
254 ELSE IF( kl.LT.0 ) THEN
255 info = -5
256 ELSE IF( ku.LT.0 ) THEN
257 info = -6
258 ELSE IF( ldab.LT.klu1 ) THEN
259 info = -8
260 ELSE IF( ldq.LT.1 .OR. wantq .AND. ldq.LT.max( 1, m ) ) THEN
261 info = -12
262 ELSE IF( ldpt.LT.1 .OR. wantpt .AND. ldpt.LT.max( 1, n ) ) THEN
263 info = -14
264 ELSE IF( ldc.LT.1 .OR. wantc .AND. ldc.LT.max( 1, m ) ) THEN
265 info = -16
266 END IF
267 IF( info.NE.0 ) THEN
268 CALL xerbla( 'ZGBBRD', -info )
269 RETURN
270 END IF
271*
272* Initialize Q and P**H to the unit matrix, if needed
273*
274 IF( wantq )
275 $ CALL zlaset( 'Full', m, m, czero, cone, q, ldq )
276 IF( wantpt )
277 $ CALL zlaset( 'Full', n, n, czero, cone, pt, ldpt )
278*
279* Quick return if possible.
280*
281 IF( m.EQ.0 .OR. n.EQ.0 )
282 $ RETURN
283*
284 minmn = min( m, n )
285*
286 IF( kl+ku.GT.1 ) THEN
287*
288* Reduce to upper bidiagonal form if KU > 0; if KU = 0, reduce
289* first to lower bidiagonal form and then transform to upper
290* bidiagonal
291*
292 IF( ku.GT.0 ) THEN
293 ml0 = 1
294 mu0 = 2
295 ELSE
296 ml0 = 2
297 mu0 = 1
298 END IF
299*
300* Wherever possible, plane rotations are generated and applied in
301* vector operations of length NR over the index set J1:J2:KLU1.
302*
303* The complex sines of the plane rotations are stored in WORK,
304* and the real cosines in RWORK.
305*
306 klm = min( m-1, kl )
307 kun = min( n-1, ku )
308 kb = klm + kun
309 kb1 = kb + 1
310 inca = kb1*ldab
311 nr = 0
312 j1 = klm + 2
313 j2 = 1 - kun
314*
315 DO 90 i = 1, minmn
316*
317* Reduce i-th column and i-th row of matrix to bidiagonal form
318*
319 ml = klm + 1
320 mu = kun + 1
321 DO 80 kk = 1, kb
322 j1 = j1 + kb
323 j2 = j2 + kb
324*
325* generate plane rotations to annihilate nonzero elements
326* which have been created below the band
327*
328 IF( nr.GT.0 )
329 $ CALL zlargv( nr, ab( klu1, j1-klm-1 ), inca,
330 $ work( j1 ), kb1, rwork( j1 ), kb1 )
331*
332* apply plane rotations from the left
333*
334 DO 10 l = 1, kb
335 IF( j2-klm+l-1.GT.n ) THEN
336 nrt = nr - 1
337 ELSE
338 nrt = nr
339 END IF
340 IF( nrt.GT.0 )
341 $ CALL zlartv( nrt, ab( klu1-l, j1-klm+l-1 ), inca,
342 $ ab( klu1-l+1, j1-klm+l-1 ), inca,
343 $ rwork( j1 ), work( j1 ), kb1 )
344 10 CONTINUE
345*
346 IF( ml.GT.ml0 ) THEN
347 IF( ml.LE.m-i+1 ) THEN
348*
349* generate plane rotation to annihilate a(i+ml-1,i)
350* within the band, and apply rotation from the left
351*
352 CALL zlartg( ab( ku+ml-1, i ), ab( ku+ml, i ),
353 $ rwork( i+ml-1 ), work( i+ml-1 ), ra )
354 ab( ku+ml-1, i ) = ra
355 IF( i.LT.n )
356 $ CALL zrot( min( ku+ml-2, n-i ),
357 $ ab( ku+ml-2, i+1 ), ldab-1,
358 $ ab( ku+ml-1, i+1 ), ldab-1,
359 $ rwork( i+ml-1 ), work( i+ml-1 ) )
360 END IF
361 nr = nr + 1
362 j1 = j1 - kb1
363 END IF
364*
365 IF( wantq ) THEN
366*
367* accumulate product of plane rotations in Q
368*
369 DO 20 j = j1, j2, kb1
370 CALL zrot( m, q( 1, j-1 ), 1, q( 1, j ), 1,
371 $ rwork( j ), dconjg( work( j ) ) )
372 20 CONTINUE
373 END IF
374*
375 IF( wantc ) THEN
376*
377* apply plane rotations to C
378*
379 DO 30 j = j1, j2, kb1
380 CALL zrot( ncc, c( j-1, 1 ), ldc, c( j, 1 ), ldc,
381 $ rwork( j ), work( j ) )
382 30 CONTINUE
383 END IF
384*
385 IF( j2+kun.GT.n ) THEN
386*
387* adjust J2 to keep within the bounds of the matrix
388*
389 nr = nr - 1
390 j2 = j2 - kb1
391 END IF
392*
393 DO 40 j = j1, j2, kb1
394*
395* create nonzero element a(j-1,j+ku) above the band
396* and store it in WORK(n+1:2*n)
397*
398 work( j+kun ) = work( j )*ab( 1, j+kun )
399 ab( 1, j+kun ) = rwork( j )*ab( 1, j+kun )
400 40 CONTINUE
401*
402* generate plane rotations to annihilate nonzero elements
403* which have been generated above the band
404*
405 IF( nr.GT.0 )
406 $ CALL zlargv( nr, ab( 1, j1+kun-1 ), inca,
407 $ work( j1+kun ), kb1, rwork( j1+kun ),
408 $ kb1 )
409*
410* apply plane rotations from the right
411*
412 DO 50 l = 1, kb
413 IF( j2+l-1.GT.m ) THEN
414 nrt = nr - 1
415 ELSE
416 nrt = nr
417 END IF
418 IF( nrt.GT.0 )
419 $ CALL zlartv( nrt, ab( l+1, j1+kun-1 ), inca,
420 $ ab( l, j1+kun ), inca,
421 $ rwork( j1+kun ), work( j1+kun ), kb1 )
422 50 CONTINUE
423*
424 IF( ml.EQ.ml0 .AND. mu.GT.mu0 ) THEN
425 IF( mu.LE.n-i+1 ) THEN
426*
427* generate plane rotation to annihilate a(i,i+mu-1)
428* within the band, and apply rotation from the right
429*
430 CALL zlartg( ab( ku-mu+3, i+mu-2 ),
431 $ ab( ku-mu+2, i+mu-1 ),
432 $ rwork( i+mu-1 ), work( i+mu-1 ), ra )
433 ab( ku-mu+3, i+mu-2 ) = ra
434 CALL zrot( min( kl+mu-2, m-i ),
435 $ ab( ku-mu+4, i+mu-2 ), 1,
436 $ ab( ku-mu+3, i+mu-1 ), 1,
437 $ rwork( i+mu-1 ), work( i+mu-1 ) )
438 END IF
439 nr = nr + 1
440 j1 = j1 - kb1
441 END IF
442*
443 IF( wantpt ) THEN
444*
445* accumulate product of plane rotations in P**H
446*
447 DO 60 j = j1, j2, kb1
448 CALL zrot( n, pt( j+kun-1, 1 ), ldpt,
449 $ pt( j+kun, 1 ), ldpt, rwork( j+kun ),
450 $ dconjg( work( j+kun ) ) )
451 60 CONTINUE
452 END IF
453*
454 IF( j2+kb.GT.m ) THEN
455*
456* adjust J2 to keep within the bounds of the matrix
457*
458 nr = nr - 1
459 j2 = j2 - kb1
460 END IF
461*
462 DO 70 j = j1, j2, kb1
463*
464* create nonzero element a(j+kl+ku,j+ku-1) below the
465* band and store it in WORK(1:n)
466*
467 work( j+kb ) = work( j+kun )*ab( klu1, j+kun )
468 ab( klu1, j+kun ) = rwork( j+kun )*ab( klu1, j+kun )
469 70 CONTINUE
470*
471 IF( ml.GT.ml0 ) THEN
472 ml = ml - 1
473 ELSE
474 mu = mu - 1
475 END IF
476 80 CONTINUE
477 90 CONTINUE
478 END IF
479*
480 IF( ku.EQ.0 .AND. kl.GT.0 ) THEN
481*
482* A has been reduced to complex lower bidiagonal form
483*
484* Transform lower bidiagonal form to upper bidiagonal by applying
485* plane rotations from the left, overwriting superdiagonal
486* elements on subdiagonal elements
487*
488 DO 100 i = 1, min( m-1, n )
489 CALL zlartg( ab( 1, i ), ab( 2, i ), rc, rs, ra )
490 ab( 1, i ) = ra
491 IF( i.LT.n ) THEN
492 ab( 2, i ) = rs*ab( 1, i+1 )
493 ab( 1, i+1 ) = rc*ab( 1, i+1 )
494 END IF
495 IF( wantq )
496 $ CALL zrot( m, q( 1, i ), 1, q( 1, i+1 ), 1, rc,
497 $ dconjg( rs ) )
498 IF( wantc )
499 $ CALL zrot( ncc, c( i, 1 ), ldc, c( i+1, 1 ), ldc, rc,
500 $ rs )
501 100 CONTINUE
502 ELSE
503*
504* A has been reduced to complex upper bidiagonal form or is
505* diagonal
506*
507 IF( ku.GT.0 .AND. m.LT.n ) THEN
508*
509* Annihilate a(m,m+1) by applying plane rotations from the
510* right
511*
512 rb = ab( ku, m+1 )
513 DO 110 i = m, 1, -1
514 CALL zlartg( ab( ku+1, i ), rb, rc, rs, ra )
515 ab( ku+1, i ) = ra
516 IF( i.GT.1 ) THEN
517 rb = -dconjg( rs )*ab( ku, i )
518 ab( ku, i ) = rc*ab( ku, i )
519 END IF
520 IF( wantpt )
521 $ CALL zrot( n, pt( i, 1 ), ldpt, pt( m+1, 1 ), ldpt,
522 $ rc, dconjg( rs ) )
523 110 CONTINUE
524 END IF
525 END IF
526*
527* Make diagonal and superdiagonal elements real, storing them in D
528* and E
529*
530 t = ab( ku+1, 1 )
531 DO 120 i = 1, minmn
532 abst = abs( t )
533 d( i ) = abst
534 IF( abst.NE.zero ) THEN
535 t = t / abst
536 ELSE
537 t = cone
538 END IF
539 IF( wantq )
540 $ CALL zscal( m, t, q( 1, i ), 1 )
541 IF( wantc )
542 $ CALL zscal( ncc, dconjg( t ), c( i, 1 ), ldc )
543 IF( i.LT.minmn ) THEN
544 IF( ku.EQ.0 .AND. kl.EQ.0 ) THEN
545 e( i ) = zero
546 t = ab( 1, i+1 )
547 ELSE
548 IF( ku.EQ.0 ) THEN
549 t = ab( 2, i )*dconjg( t )
550 ELSE
551 t = ab( ku, i+1 )*dconjg( t )
552 END IF
553 abst = abs( t )
554 e( i ) = abst
555 IF( abst.NE.zero ) THEN
556 t = t / abst
557 ELSE
558 t = cone
559 END IF
560 IF( wantpt )
561 $ CALL zscal( n, t, pt( i+1, 1 ), ldpt )
562 t = ab( ku+1, i+1 )*dconjg( t )
563 END IF
564 END IF
565 120 CONTINUE
566 RETURN
567*
568* End of ZGBBRD
569*
570 END
zlartg
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition zlartg.f90:118
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
zgbbrd
subroutine zgbbrd(vect, m, n, ncc, kl, ku, ab, ldab, d, e, q, ldq, pt, ldpt, c, ldc, work, rwork, info)
ZGBBRD
Definition zgbbrd.f:193
zlargv
subroutine zlargv(n, x, incx, y, incy, c, incc)
ZLARGV generates a vector of plane rotations with real cosines and complex sines.
Definition zlargv.f:122
zrot
subroutine zrot(n, cx, incx, cy, incy, c, s)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition zrot.f:103
zlartv
subroutine zlartv(n, x, incx, y, incy, c, s, incc)
ZLARTV applies a vector of plane rotations with real cosines and complex sines to the elements of a p...
Definition zlartv.f:107
zlaset
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:106
zscal
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21