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