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