001:       SUBROUTINE CGBBRD( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
002:      $                   LDQ, PT, LDPT, C, LDC, WORK, RWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          VECT
010:       INTEGER            INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
011: *     ..
012: *     .. Array Arguments ..
013:       REAL               D( * ), E( * ), RWORK( * )
014:       COMPLEX            AB( LDAB, * ), C( LDC, * ), PT( LDPT, * ),
015:      $                   Q( LDQ, * ), WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  CGBBRD reduces a complex general m-by-n band matrix A to real upper
022: *  bidiagonal form B by a unitary transformation: Q' * A * P = B.
023: *
024: *  The routine computes B, and optionally forms Q or P', or computes
025: *  Q'*C for a given matrix C.
026: *
027: *  Arguments
028: *  =========
029: *
030: *  VECT    (input) CHARACTER*1
031: *          Specifies whether or not the matrices Q and P' are to be
032: *          formed.
033: *          = 'N': do not form Q or P';
034: *          = 'Q': form Q only;
035: *          = 'P': form P' only;
036: *          = 'B': form both.
037: *
038: *  M       (input) INTEGER
039: *          The number of rows of the matrix A.  M >= 0.
040: *
041: *  N       (input) INTEGER
042: *          The number of columns of the matrix A.  N >= 0.
043: *
044: *  NCC     (input) INTEGER
045: *          The number of columns of the matrix C.  NCC >= 0.
046: *
047: *  KL      (input) INTEGER
048: *          The number of subdiagonals of the matrix A. KL >= 0.
049: *
050: *  KU      (input) INTEGER
051: *          The number of superdiagonals of the matrix A. KU >= 0.
052: *
053: *  AB      (input/output) COMPLEX array, dimension (LDAB,N)
054: *          On entry, the m-by-n band matrix A, stored in rows 1 to
055: *          KL+KU+1. The j-th column of A is stored in the j-th column of
056: *          the array AB as follows:
057: *          AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
058: *          On exit, A is overwritten by values generated during the
059: *          reduction.
060: *
061: *  LDAB    (input) INTEGER
062: *          The leading dimension of the array A. LDAB >= KL+KU+1.
063: *
064: *  D       (output) REAL array, dimension (min(M,N))
065: *          The diagonal elements of the bidiagonal matrix B.
066: *
067: *  E       (output) REAL array, dimension (min(M,N)-1)
068: *          The superdiagonal elements of the bidiagonal matrix B.
069: *
070: *  Q       (output) COMPLEX array, dimension (LDQ,M)
071: *          If VECT = 'Q' or 'B', the m-by-m unitary matrix Q.
072: *          If VECT = 'N' or 'P', the array Q is not referenced.
073: *
074: *  LDQ     (input) INTEGER
075: *          The leading dimension of the array Q.
076: *          LDQ >= max(1,M) if VECT = 'Q' or 'B'; LDQ >= 1 otherwise.
077: *
078: *  PT      (output) COMPLEX array, dimension (LDPT,N)
079: *          If VECT = 'P' or 'B', the n-by-n unitary matrix P'.
080: *          If VECT = 'N' or 'Q', the array PT is not referenced.
081: *
082: *  LDPT    (input) INTEGER
083: *          The leading dimension of the array PT.
084: *          LDPT >= max(1,N) if VECT = 'P' or 'B'; LDPT >= 1 otherwise.
085: *
086: *  C       (input/output) COMPLEX array, dimension (LDC,NCC)
087: *          On entry, an m-by-ncc matrix C.
088: *          On exit, C is overwritten by Q'*C.
089: *          C is not referenced if NCC = 0.
090: *
091: *  LDC     (input) INTEGER
092: *          The leading dimension of the array C.
093: *          LDC >= max(1,M) if NCC > 0; LDC >= 1 if NCC = 0.
094: *
095: *  WORK    (workspace) COMPLEX array, dimension (max(M,N))
096: *
097: *  RWORK   (workspace) REAL array, dimension (max(M,N))
098: *
099: *  INFO    (output) INTEGER
100: *          = 0:  successful exit.
101: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
102: *
103: *  =====================================================================
104: *
105: *     .. Parameters ..
106:       REAL               ZERO
107:       PARAMETER          ( ZERO = 0.0E+0 )
108:       COMPLEX            CZERO, CONE
109:       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
110:      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
111: *     ..
112: *     .. Local Scalars ..
113:       LOGICAL            WANTB, WANTC, WANTPT, WANTQ
114:       INTEGER            I, INCA, J, J1, J2, KB, KB1, KK, KLM, KLU1,
115:      $                   KUN, L, MINMN, ML, ML0, MU, MU0, NR, NRT
116:       REAL               ABST, RC
117:       COMPLEX            RA, RB, RS, T
118: *     ..
119: *     .. External Subroutines ..
120:       EXTERNAL           CLARGV, CLARTG, CLARTV, CLASET, CROT, CSCAL,
121:      $                   XERBLA
122: *     ..
123: *     .. Intrinsic Functions ..
124:       INTRINSIC          ABS, CONJG, MAX, MIN
125: *     ..
126: *     .. External Functions ..
127:       LOGICAL            LSAME
128:       EXTERNAL           LSAME
129: *     ..
130: *     .. Executable Statements ..
131: *
132: *     Test the input parameters
133: *
134:       WANTB = LSAME( VECT, 'B' )
135:       WANTQ = LSAME( VECT, 'Q' ) .OR. WANTB
136:       WANTPT = LSAME( VECT, 'P' ) .OR. WANTB
137:       WANTC = NCC.GT.0
138:       KLU1 = KL + KU + 1
139:       INFO = 0
140:       IF( .NOT.WANTQ .AND. .NOT.WANTPT .AND. .NOT.LSAME( VECT, 'N' ) )
141:      $     THEN
142:          INFO = -1
143:       ELSE IF( M.LT.0 ) THEN
144:          INFO = -2
145:       ELSE IF( N.LT.0 ) THEN
146:          INFO = -3
147:       ELSE IF( NCC.LT.0 ) THEN
148:          INFO = -4
149:       ELSE IF( KL.LT.0 ) THEN
150:          INFO = -5
151:       ELSE IF( KU.LT.0 ) THEN
152:          INFO = -6
153:       ELSE IF( LDAB.LT.KLU1 ) THEN
154:          INFO = -8
155:       ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. LDQ.LT.MAX( 1, M ) ) THEN
156:          INFO = -12
157:       ELSE IF( LDPT.LT.1 .OR. WANTPT .AND. LDPT.LT.MAX( 1, N ) ) THEN
158:          INFO = -14
159:       ELSE IF( LDC.LT.1 .OR. WANTC .AND. LDC.LT.MAX( 1, M ) ) THEN
160:          INFO = -16
161:       END IF
162:       IF( INFO.NE.0 ) THEN
163:          CALL XERBLA( 'CGBBRD', -INFO )
164:          RETURN
165:       END IF
166: *
167: *     Initialize Q and P' to the unit matrix, if needed
168: *
169:       IF( WANTQ )
170:      $   CALL CLASET( 'Full', M, M, CZERO, CONE, Q, LDQ )
171:       IF( WANTPT )
172:      $   CALL CLASET( 'Full', N, N, CZERO, CONE, PT, LDPT )
173: *
174: *     Quick return if possible.
175: *
176:       IF( M.EQ.0 .OR. N.EQ.0 )
177:      $   RETURN
178: *
179:       MINMN = MIN( M, N )
180: *
181:       IF( KL+KU.GT.1 ) THEN
182: *
183: *        Reduce to upper bidiagonal form if KU > 0; if KU = 0, reduce
184: *        first to lower bidiagonal form and then transform to upper
185: *        bidiagonal
186: *
187:          IF( KU.GT.0 ) THEN
188:             ML0 = 1
189:             MU0 = 2
190:          ELSE
191:             ML0 = 2
192:             MU0 = 1
193:          END IF
194: *
195: *        Wherever possible, plane rotations are generated and applied in
196: *        vector operations of length NR over the index set J1:J2:KLU1.
197: *
198: *        The complex sines of the plane rotations are stored in WORK,
199: *        and the real cosines in RWORK.
200: *
201:          KLM = MIN( M-1, KL )
202:          KUN = MIN( N-1, KU )
203:          KB = KLM + KUN
204:          KB1 = KB + 1
205:          INCA = KB1*LDAB
206:          NR = 0
207:          J1 = KLM + 2
208:          J2 = 1 - KUN
209: *
210:          DO 90 I = 1, MINMN
211: *
212: *           Reduce i-th column and i-th row of matrix to bidiagonal form
213: *
214:             ML = KLM + 1
215:             MU = KUN + 1
216:             DO 80 KK = 1, KB
217:                J1 = J1 + KB
218:                J2 = J2 + KB
219: *
220: *              generate plane rotations to annihilate nonzero elements
221: *              which have been created below the band
222: *
223:                IF( NR.GT.0 )
224:      $            CALL CLARGV( NR, AB( KLU1, J1-KLM-1 ), INCA,
225:      $                         WORK( J1 ), KB1, RWORK( J1 ), KB1 )
226: *
227: *              apply plane rotations from the left
228: *
229:                DO 10 L = 1, KB
230:                   IF( J2-KLM+L-1.GT.N ) THEN
231:                      NRT = NR - 1
232:                   ELSE
233:                      NRT = NR
234:                   END IF
235:                   IF( NRT.GT.0 )
236:      $               CALL CLARTV( NRT, AB( KLU1-L, J1-KLM+L-1 ), INCA,
237:      $                            AB( KLU1-L+1, J1-KLM+L-1 ), INCA,
238:      $                            RWORK( J1 ), WORK( J1 ), KB1 )
239:    10          CONTINUE
240: *
241:                IF( ML.GT.ML0 ) THEN
242:                   IF( ML.LE.M-I+1 ) THEN
243: *
244: *                    generate plane rotation to annihilate a(i+ml-1,i)
245: *                    within the band, and apply rotation from the left
246: *
247:                      CALL CLARTG( AB( KU+ML-1, I ), AB( KU+ML, I ),
248:      $                            RWORK( I+ML-1 ), WORK( I+ML-1 ), RA )
249:                      AB( KU+ML-1, I ) = RA
250:                      IF( I.LT.N )
251:      $                  CALL CROT( MIN( KU+ML-2, N-I ),
252:      $                             AB( KU+ML-2, I+1 ), LDAB-1,
253:      $                             AB( KU+ML-1, I+1 ), LDAB-1,
254:      $                             RWORK( I+ML-1 ), WORK( I+ML-1 ) )
255:                   END IF
256:                   NR = NR + 1
257:                   J1 = J1 - KB1
258:                END IF
259: *
260:                IF( WANTQ ) THEN
261: *
262: *                 accumulate product of plane rotations in Q
263: *
264:                   DO 20 J = J1, J2, KB1
265:                      CALL CROT( M, Q( 1, J-1 ), 1, Q( 1, J ), 1,
266:      $                          RWORK( J ), CONJG( WORK( J ) ) )
267:    20             CONTINUE
268:                END IF
269: *
270:                IF( WANTC ) THEN
271: *
272: *                 apply plane rotations to C
273: *
274:                   DO 30 J = J1, J2, KB1
275:                      CALL CROT( NCC, C( J-1, 1 ), LDC, C( J, 1 ), LDC,
276:      $                          RWORK( J ), WORK( J ) )
277:    30             CONTINUE
278:                END IF
279: *
280:                IF( J2+KUN.GT.N ) THEN
281: *
282: *                 adjust J2 to keep within the bounds of the matrix
283: *
284:                   NR = NR - 1
285:                   J2 = J2 - KB1
286:                END IF
287: *
288:                DO 40 J = J1, J2, KB1
289: *
290: *                 create nonzero element a(j-1,j+ku) above the band
291: *                 and store it in WORK(n+1:2*n)
292: *
293:                   WORK( J+KUN ) = WORK( J )*AB( 1, J+KUN )
294:                   AB( 1, J+KUN ) = RWORK( J )*AB( 1, J+KUN )
295:    40          CONTINUE
296: *
297: *              generate plane rotations to annihilate nonzero elements
298: *              which have been generated above the band
299: *
300:                IF( NR.GT.0 )
301:      $            CALL CLARGV( NR, AB( 1, J1+KUN-1 ), INCA,
302:      $                         WORK( J1+KUN ), KB1, RWORK( J1+KUN ),
303:      $                         KB1 )
304: *
305: *              apply plane rotations from the right
306: *
307:                DO 50 L = 1, KB
308:                   IF( J2+L-1.GT.M ) THEN
309:                      NRT = NR - 1
310:                   ELSE
311:                      NRT = NR
312:                   END IF
313:                   IF( NRT.GT.0 )
314:      $               CALL CLARTV( NRT, AB( L+1, J1+KUN-1 ), INCA,
315:      $                            AB( L, J1+KUN ), INCA,
316:      $                            RWORK( J1+KUN ), WORK( J1+KUN ), KB1 )
317:    50          CONTINUE
318: *
319:                IF( ML.EQ.ML0 .AND. MU.GT.MU0 ) THEN
320:                   IF( MU.LE.N-I+1 ) THEN
321: *
322: *                    generate plane rotation to annihilate a(i,i+mu-1)
323: *                    within the band, and apply rotation from the right
324: *
325:                      CALL CLARTG( AB( KU-MU+3, I+MU-2 ),
326:      $                            AB( KU-MU+2, I+MU-1 ),
327:      $                            RWORK( I+MU-1 ), WORK( I+MU-1 ), RA )
328:                      AB( KU-MU+3, I+MU-2 ) = RA
329:                      CALL CROT( MIN( KL+MU-2, M-I ),
330:      $                          AB( KU-MU+4, I+MU-2 ), 1,
331:      $                          AB( KU-MU+3, I+MU-1 ), 1,
332:      $                          RWORK( I+MU-1 ), WORK( I+MU-1 ) )
333:                   END IF
334:                   NR = NR + 1
335:                   J1 = J1 - KB1
336:                END IF
337: *
338:                IF( WANTPT ) THEN
339: *
340: *                 accumulate product of plane rotations in P'
341: *
342:                   DO 60 J = J1, J2, KB1
343:                      CALL CROT( N, PT( J+KUN-1, 1 ), LDPT,
344:      $                          PT( J+KUN, 1 ), LDPT, RWORK( J+KUN ),
345:      $                          CONJG( WORK( J+KUN ) ) )
346:    60             CONTINUE
347:                END IF
348: *
349:                IF( J2+KB.GT.M ) THEN
350: *
351: *                 adjust J2 to keep within the bounds of the matrix
352: *
353:                   NR = NR - 1
354:                   J2 = J2 - KB1
355:                END IF
356: *
357:                DO 70 J = J1, J2, KB1
358: *
359: *                 create nonzero element a(j+kl+ku,j+ku-1) below the
360: *                 band and store it in WORK(1:n)
361: *
362:                   WORK( J+KB ) = WORK( J+KUN )*AB( KLU1, J+KUN )
363:                   AB( KLU1, J+KUN ) = RWORK( J+KUN )*AB( KLU1, J+KUN )
364:    70          CONTINUE
365: *
366:                IF( ML.GT.ML0 ) THEN
367:                   ML = ML - 1
368:                ELSE
369:                   MU = MU - 1
370:                END IF
371:    80       CONTINUE
372:    90    CONTINUE
373:       END IF
374: *
375:       IF( KU.EQ.0 .AND. KL.GT.0 ) THEN
376: *
377: *        A has been reduced to complex lower bidiagonal form
378: *
379: *        Transform lower bidiagonal form to upper bidiagonal by applying
380: *        plane rotations from the left, overwriting superdiagonal
381: *        elements on subdiagonal elements
382: *
383:          DO 100 I = 1, MIN( M-1, N )
384:             CALL CLARTG( AB( 1, I ), AB( 2, I ), RC, RS, RA )
385:             AB( 1, I ) = RA
386:             IF( I.LT.N ) THEN
387:                AB( 2, I ) = RS*AB( 1, I+1 )
388:                AB( 1, I+1 ) = RC*AB( 1, I+1 )
389:             END IF
390:             IF( WANTQ )
391:      $         CALL CROT( M, Q( 1, I ), 1, Q( 1, I+1 ), 1, RC,
392:      $                    CONJG( RS ) )
393:             IF( WANTC )
394:      $         CALL CROT( NCC, C( I, 1 ), LDC, C( I+1, 1 ), LDC, RC,
395:      $                    RS )
396:   100    CONTINUE
397:       ELSE
398: *
399: *        A has been reduced to complex upper bidiagonal form or is
400: *        diagonal
401: *
402:          IF( KU.GT.0 .AND. M.LT.N ) THEN
403: *
404: *           Annihilate a(m,m+1) by applying plane rotations from the
405: *           right
406: *
407:             RB = AB( KU, M+1 )
408:             DO 110 I = M, 1, -1
409:                CALL CLARTG( AB( KU+1, I ), RB, RC, RS, RA )
410:                AB( KU+1, I ) = RA
411:                IF( I.GT.1 ) THEN
412:                   RB = -CONJG( RS )*AB( KU, I )
413:                   AB( KU, I ) = RC*AB( KU, I )
414:                END IF
415:                IF( WANTPT )
416:      $            CALL CROT( N, PT( I, 1 ), LDPT, PT( M+1, 1 ), LDPT,
417:      $                       RC, CONJG( RS ) )
418:   110       CONTINUE
419:          END IF
420:       END IF
421: *
422: *     Make diagonal and superdiagonal elements real, storing them in D
423: *     and E
424: *
425:       T = AB( KU+1, 1 )
426:       DO 120 I = 1, MINMN
427:          ABST = ABS( T )
428:          D( I ) = ABST
429:          IF( ABST.NE.ZERO ) THEN
430:             T = T / ABST
431:          ELSE
432:             T = CONE
433:          END IF
434:          IF( WANTQ )
435:      $      CALL CSCAL( M, T, Q( 1, I ), 1 )
436:          IF( WANTC )
437:      $      CALL CSCAL( NCC, CONJG( T ), C( I, 1 ), LDC )
438:          IF( I.LT.MINMN ) THEN
439:             IF( KU.EQ.0 .AND. KL.EQ.0 ) THEN
440:                E( I ) = ZERO
441:                T = AB( 1, I+1 )
442:             ELSE
443:                IF( KU.EQ.0 ) THEN
444:                   T = AB( 2, I )*CONJG( T )
445:                ELSE
446:                   T = AB( KU, I+1 )*CONJG( T )
447:                END IF
448:                ABST = ABS( T )
449:                E( I ) = ABST
450:                IF( ABST.NE.ZERO ) THEN
451:                   T = T / ABST
452:                ELSE
453:                   T = CONE
454:                END IF
455:                IF( WANTPT )
456:      $            CALL CSCAL( N, T, PT( I+1, 1 ), LDPT )
457:                T = AB( KU+1, I+1 )*CONJG( T )
458:             END IF
459:          END IF
460:   120 CONTINUE
461:       RETURN
462: *
463: *     End of CGBBRD
464: *
465:       END
466: