1: SUBROUTINE DGBBRD( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
2: $ LDQ, PT, LDPT, C, LDC, WORK, INFO )
3: *
4: * -- LAPACK routine (version 3.3.1) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * -- April 2011 --
8: *
9: * .. Scalar Arguments ..
10: CHARACTER VECT
11: INTEGER INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION AB( LDAB, * ), C( LDC, * ), D( * ), E( * ),
15: $ PT( LDPT, * ), Q( LDQ, * ), WORK( * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DGBBRD reduces a real general m-by-n band matrix A to upper
22: * bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
23: *
24: * The routine computes B, and optionally forms Q or P**T, or computes
25: * Q**T*C for a given matrix C.
26: *
27: * Arguments
28: * =========
29: *
30: * VECT (input) CHARACTER*1
31: * Specifies whether or not the matrices Q and P**T are to be
32: * formed.
33: * = 'N': do not form Q or P**T;
34: * = 'Q': form Q only;
35: * = 'P': form P**T only;
36: * = 'B': form both.
37: *
38: * M (input) INTEGER
39: * The number of rows of the matrix A. M >= 0.
40: *
41: * N (input) INTEGER
42: * The number of columns of the matrix A. N >= 0.
43: *
44: * NCC (input) INTEGER
45: * The number of columns of the matrix C. NCC >= 0.
46: *
47: * KL (input) INTEGER
48: * The number of subdiagonals of the matrix A. KL >= 0.
49: *
50: * KU (input) INTEGER
51: * The number of superdiagonals of the matrix A. KU >= 0.
52: *
53: * AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
54: * On entry, the m-by-n band matrix A, stored in rows 1 to
55: * KL+KU+1. The j-th column of A is stored in the j-th column of
56: * the array AB as follows:
57: * AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
58: * On exit, A is overwritten by values generated during the
59: * reduction.
60: *
61: * LDAB (input) INTEGER
62: * The leading dimension of the array A. LDAB >= KL+KU+1.
63: *
64: * D (output) DOUBLE PRECISION array, dimension (min(M,N))
65: * The diagonal elements of the bidiagonal matrix B.
66: *
67: * E (output) DOUBLE PRECISION array, dimension (min(M,N)-1)
68: * The superdiagonal elements of the bidiagonal matrix B.
69: *
70: * Q (output) DOUBLE PRECISION array, dimension (LDQ,M)
71: * If VECT = 'Q' or 'B', the m-by-m orthogonal matrix Q.
72: * If VECT = 'N' or 'P', the array Q is not referenced.
73: *
74: * LDQ (input) INTEGER
75: * The leading dimension of the array Q.
76: * LDQ >= max(1,M) if VECT = 'Q' or 'B'; LDQ >= 1 otherwise.
77: *
78: * PT (output) DOUBLE PRECISION array, dimension (LDPT,N)
79: * If VECT = 'P' or 'B', the n-by-n orthogonal matrix P'.
80: * If VECT = 'N' or 'Q', the array PT is not referenced.
81: *
82: * LDPT (input) INTEGER
83: * The leading dimension of the array PT.
84: * LDPT >= max(1,N) if VECT = 'P' or 'B'; LDPT >= 1 otherwise.
85: *
86: * C (input/output) DOUBLE PRECISION array, dimension (LDC,NCC)
87: * On entry, an m-by-ncc matrix C.
88: * On exit, C is overwritten by Q**T*C.
89: * C is not referenced if NCC = 0.
90: *
91: * LDC (input) INTEGER
92: * The leading dimension of the array C.
93: * LDC >= max(1,M) if NCC > 0; LDC >= 1 if NCC = 0.
94: *
95: * WORK (workspace) DOUBLE PRECISION array, dimension (2*max(M,N))
96: *
97: * INFO (output) INTEGER
98: * = 0: successful exit.
99: * < 0: if INFO = -i, the i-th argument had an illegal value.
100: *
101: * =====================================================================
102: *
103: * .. Parameters ..
104: DOUBLE PRECISION ZERO, ONE
105: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
106: * ..
107: * .. Local Scalars ..
108: LOGICAL WANTB, WANTC, WANTPT, WANTQ
109: INTEGER I, INCA, J, J1, J2, KB, KB1, KK, KLM, KLU1,
110: $ KUN, L, MINMN, ML, ML0, MN, MU, MU0, NR, NRT
111: DOUBLE PRECISION RA, RB, RC, RS
112: * ..
113: * .. External Subroutines ..
114: EXTERNAL DLARGV, DLARTG, DLARTV, DLASET, DROT, XERBLA
115: * ..
116: * .. Intrinsic Functions ..
117: INTRINSIC MAX, MIN
118: * ..
119: * .. External Functions ..
120: LOGICAL LSAME
121: EXTERNAL LSAME
122: * ..
123: * .. Executable Statements ..
124: *
125: * Test the input parameters
126: *
127: WANTB = LSAME( VECT, 'B' )
128: WANTQ = LSAME( VECT, 'Q' ) .OR. WANTB
129: WANTPT = LSAME( VECT, 'P' ) .OR. WANTB
130: WANTC = NCC.GT.0
131: KLU1 = KL + KU + 1
132: INFO = 0
133: IF( .NOT.WANTQ .AND. .NOT.WANTPT .AND. .NOT.LSAME( VECT, 'N' ) )
134: $ THEN
135: INFO = -1
136: ELSE IF( M.LT.0 ) THEN
137: INFO = -2
138: ELSE IF( N.LT.0 ) THEN
139: INFO = -3
140: ELSE IF( NCC.LT.0 ) THEN
141: INFO = -4
142: ELSE IF( KL.LT.0 ) THEN
143: INFO = -5
144: ELSE IF( KU.LT.0 ) THEN
145: INFO = -6
146: ELSE IF( LDAB.LT.KLU1 ) THEN
147: INFO = -8
148: ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. LDQ.LT.MAX( 1, M ) ) THEN
149: INFO = -12
150: ELSE IF( LDPT.LT.1 .OR. WANTPT .AND. LDPT.LT.MAX( 1, N ) ) THEN
151: INFO = -14
152: ELSE IF( LDC.LT.1 .OR. WANTC .AND. LDC.LT.MAX( 1, M ) ) THEN
153: INFO = -16
154: END IF
155: IF( INFO.NE.0 ) THEN
156: CALL XERBLA( 'DGBBRD', -INFO )
157: RETURN
158: END IF
159: *
160: * Initialize Q and P**T to the unit matrix, if needed
161: *
162: IF( WANTQ )
163: $ CALL DLASET( 'Full', M, M, ZERO, ONE, Q, LDQ )
164: IF( WANTPT )
165: $ CALL DLASET( 'Full', N, N, ZERO, ONE, PT, LDPT )
166: *
167: * Quick return if possible.
168: *
169: IF( M.EQ.0 .OR. N.EQ.0 )
170: $ RETURN
171: *
172: MINMN = MIN( M, N )
173: *
174: IF( KL+KU.GT.1 ) THEN
175: *
176: * Reduce to upper bidiagonal form if KU > 0; if KU = 0, reduce
177: * first to lower bidiagonal form and then transform to upper
178: * bidiagonal
179: *
180: IF( KU.GT.0 ) THEN
181: ML0 = 1
182: MU0 = 2
183: ELSE
184: ML0 = 2
185: MU0 = 1
186: END IF
187: *
188: * Wherever possible, plane rotations are generated and applied in
189: * vector operations of length NR over the index set J1:J2:KLU1.
190: *
191: * The sines of the plane rotations are stored in WORK(1:max(m,n))
192: * and the cosines in WORK(max(m,n)+1:2*max(m,n)).
193: *
194: MN = MAX( M, N )
195: KLM = MIN( M-1, KL )
196: KUN = MIN( N-1, KU )
197: KB = KLM + KUN
198: KB1 = KB + 1
199: INCA = KB1*LDAB
200: NR = 0
201: J1 = KLM + 2
202: J2 = 1 - KUN
203: *
204: DO 90 I = 1, MINMN
205: *
206: * Reduce i-th column and i-th row of matrix to bidiagonal form
207: *
208: ML = KLM + 1
209: MU = KUN + 1
210: DO 80 KK = 1, KB
211: J1 = J1 + KB
212: J2 = J2 + KB
213: *
214: * generate plane rotations to annihilate nonzero elements
215: * which have been created below the band
216: *
217: IF( NR.GT.0 )
218: $ CALL DLARGV( NR, AB( KLU1, J1-KLM-1 ), INCA,
219: $ WORK( J1 ), KB1, WORK( MN+J1 ), KB1 )
220: *
221: * apply plane rotations from the left
222: *
223: DO 10 L = 1, KB
224: IF( J2-KLM+L-1.GT.N ) THEN
225: NRT = NR - 1
226: ELSE
227: NRT = NR
228: END IF
229: IF( NRT.GT.0 )
230: $ CALL DLARTV( NRT, AB( KLU1-L, J1-KLM+L-1 ), INCA,
231: $ AB( KLU1-L+1, J1-KLM+L-1 ), INCA,
232: $ WORK( MN+J1 ), WORK( J1 ), KB1 )
233: 10 CONTINUE
234: *
235: IF( ML.GT.ML0 ) THEN
236: IF( ML.LE.M-I+1 ) THEN
237: *
238: * generate plane rotation to annihilate a(i+ml-1,i)
239: * within the band, and apply rotation from the left
240: *
241: CALL DLARTG( AB( KU+ML-1, I ), AB( KU+ML, I ),
242: $ WORK( MN+I+ML-1 ), WORK( I+ML-1 ),
243: $ RA )
244: AB( KU+ML-1, I ) = RA
245: IF( I.LT.N )
246: $ CALL DROT( MIN( KU+ML-2, N-I ),
247: $ AB( KU+ML-2, I+1 ), LDAB-1,
248: $ AB( KU+ML-1, I+1 ), LDAB-1,
249: $ WORK( MN+I+ML-1 ), WORK( I+ML-1 ) )
250: END IF
251: NR = NR + 1
252: J1 = J1 - KB1
253: END IF
254: *
255: IF( WANTQ ) THEN
256: *
257: * accumulate product of plane rotations in Q
258: *
259: DO 20 J = J1, J2, KB1
260: CALL DROT( M, Q( 1, J-1 ), 1, Q( 1, J ), 1,
261: $ WORK( MN+J ), WORK( J ) )
262: 20 CONTINUE
263: END IF
264: *
265: IF( WANTC ) THEN
266: *
267: * apply plane rotations to C
268: *
269: DO 30 J = J1, J2, KB1
270: CALL DROT( NCC, C( J-1, 1 ), LDC, C( J, 1 ), LDC,
271: $ WORK( MN+J ), WORK( J ) )
272: 30 CONTINUE
273: END IF
274: *
275: IF( J2+KUN.GT.N ) THEN
276: *
277: * adjust J2 to keep within the bounds of the matrix
278: *
279: NR = NR - 1
280: J2 = J2 - KB1
281: END IF
282: *
283: DO 40 J = J1, J2, KB1
284: *
285: * create nonzero element a(j-1,j+ku) above the band
286: * and store it in WORK(n+1:2*n)
287: *
288: WORK( J+KUN ) = WORK( J )*AB( 1, J+KUN )
289: AB( 1, J+KUN ) = WORK( MN+J )*AB( 1, J+KUN )
290: 40 CONTINUE
291: *
292: * generate plane rotations to annihilate nonzero elements
293: * which have been generated above the band
294: *
295: IF( NR.GT.0 )
296: $ CALL DLARGV( NR, AB( 1, J1+KUN-1 ), INCA,
297: $ WORK( J1+KUN ), KB1, WORK( MN+J1+KUN ),
298: $ KB1 )
299: *
300: * apply plane rotations from the right
301: *
302: DO 50 L = 1, KB
303: IF( J2+L-1.GT.M ) THEN
304: NRT = NR - 1
305: ELSE
306: NRT = NR
307: END IF
308: IF( NRT.GT.0 )
309: $ CALL DLARTV( NRT, AB( L+1, J1+KUN-1 ), INCA,
310: $ AB( L, J1+KUN ), INCA,
311: $ WORK( MN+J1+KUN ), WORK( J1+KUN ),
312: $ KB1 )
313: 50 CONTINUE
314: *
315: IF( ML.EQ.ML0 .AND. MU.GT.MU0 ) THEN
316: IF( MU.LE.N-I+1 ) THEN
317: *
318: * generate plane rotation to annihilate a(i,i+mu-1)
319: * within the band, and apply rotation from the right
320: *
321: CALL DLARTG( AB( KU-MU+3, I+MU-2 ),
322: $ AB( KU-MU+2, I+MU-1 ),
323: $ WORK( MN+I+MU-1 ), WORK( I+MU-1 ),
324: $ RA )
325: AB( KU-MU+3, I+MU-2 ) = RA
326: CALL DROT( MIN( KL+MU-2, M-I ),
327: $ AB( KU-MU+4, I+MU-2 ), 1,
328: $ AB( KU-MU+3, I+MU-1 ), 1,
329: $ WORK( MN+I+MU-1 ), WORK( I+MU-1 ) )
330: END IF
331: NR = NR + 1
332: J1 = J1 - KB1
333: END IF
334: *
335: IF( WANTPT ) THEN
336: *
337: * accumulate product of plane rotations in P**T
338: *
339: DO 60 J = J1, J2, KB1
340: CALL DROT( N, PT( J+KUN-1, 1 ), LDPT,
341: $ PT( J+KUN, 1 ), LDPT, WORK( MN+J+KUN ),
342: $ WORK( J+KUN ) )
343: 60 CONTINUE
344: END IF
345: *
346: IF( J2+KB.GT.M ) THEN
347: *
348: * adjust J2 to keep within the bounds of the matrix
349: *
350: NR = NR - 1
351: J2 = J2 - KB1
352: END IF
353: *
354: DO 70 J = J1, J2, KB1
355: *
356: * create nonzero element a(j+kl+ku,j+ku-1) below the
357: * band and store it in WORK(1:n)
358: *
359: WORK( J+KB ) = WORK( J+KUN )*AB( KLU1, J+KUN )
360: AB( KLU1, J+KUN ) = WORK( MN+J+KUN )*AB( KLU1, J+KUN )
361: 70 CONTINUE
362: *
363: IF( ML.GT.ML0 ) THEN
364: ML = ML - 1
365: ELSE
366: MU = MU - 1
367: END IF
368: 80 CONTINUE
369: 90 CONTINUE
370: END IF
371: *
372: IF( KU.EQ.0 .AND. KL.GT.0 ) THEN
373: *
374: * A has been reduced to lower bidiagonal form
375: *
376: * Transform lower bidiagonal form to upper bidiagonal by applying
377: * plane rotations from the left, storing diagonal elements in D
378: * and off-diagonal elements in E
379: *
380: DO 100 I = 1, MIN( M-1, N )
381: CALL DLARTG( AB( 1, I ), AB( 2, I ), RC, RS, RA )
382: D( I ) = RA
383: IF( I.LT.N ) THEN
384: E( I ) = RS*AB( 1, I+1 )
385: AB( 1, I+1 ) = RC*AB( 1, I+1 )
386: END IF
387: IF( WANTQ )
388: $ CALL DROT( M, Q( 1, I ), 1, Q( 1, I+1 ), 1, RC, RS )
389: IF( WANTC )
390: $ CALL DROT( NCC, C( I, 1 ), LDC, C( I+1, 1 ), LDC, RC,
391: $ RS )
392: 100 CONTINUE
393: IF( M.LE.N )
394: $ D( M ) = AB( 1, M )
395: ELSE IF( KU.GT.0 ) THEN
396: *
397: * A has been reduced to upper bidiagonal form
398: *
399: IF( M.LT.N ) THEN
400: *
401: * Annihilate a(m,m+1) by applying plane rotations from the
402: * right, storing diagonal elements in D and off-diagonal
403: * elements in E
404: *
405: RB = AB( KU, M+1 )
406: DO 110 I = M, 1, -1
407: CALL DLARTG( AB( KU+1, I ), RB, RC, RS, RA )
408: D( I ) = RA
409: IF( I.GT.1 ) THEN
410: RB = -RS*AB( KU, I )
411: E( I-1 ) = RC*AB( KU, I )
412: END IF
413: IF( WANTPT )
414: $ CALL DROT( N, PT( I, 1 ), LDPT, PT( M+1, 1 ), LDPT,
415: $ RC, RS )
416: 110 CONTINUE
417: ELSE
418: *
419: * Copy off-diagonal elements to E and diagonal elements to D
420: *
421: DO 120 I = 1, MINMN - 1
422: E( I ) = AB( KU, I+1 )
423: 120 CONTINUE
424: DO 130 I = 1, MINMN
425: D( I ) = AB( KU+1, I )
426: 130 CONTINUE
427: END IF
428: ELSE
429: *
430: * A is diagonal. Set elements of E to zero and copy diagonal
431: * elements to D.
432: *
433: DO 140 I = 1, MINMN - 1
434: E( I ) = ZERO
435: 140 CONTINUE
436: DO 150 I = 1, MINMN
437: D( I ) = AB( 1, I )
438: 150 CONTINUE
439: END IF
440: RETURN
441: *
442: * End of DGBBRD
443: *
444: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgbbrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack