Lots of trailing whitespaces in the files of Syd. Cleaning this. No big deal.
[platform/upstream/lapack.git] / SRC / sgbbrd.f
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 *> \date November 2011
183 *
184 *> \ingroup realGBcomputational
185 *
186 *  =====================================================================
187       SUBROUTINE SGBBRD( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
188      $                   LDQ, PT, LDPT, C, LDC, WORK, INFO )
189 *
190 *  -- LAPACK computational routine (version 3.4.0) --
191 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
192 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
193 *     November 2011
194 *
195 *     .. Scalar Arguments ..
196       CHARACTER          VECT
197       INTEGER            INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
198 *     ..
199 *     .. Array Arguments ..
200       REAL               AB( LDAB, * ), C( LDC, * ), D( * ), E( * ),
201      $                   PT( LDPT, * ), Q( LDQ, * ), WORK( * )
202 *     ..
203 *
204 *  =====================================================================
205 *
206 *     .. Parameters ..
207       REAL               ZERO, ONE
208       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
209 *     ..
210 *     .. Local Scalars ..
211       LOGICAL            WANTB, WANTC, WANTPT, WANTQ
212       INTEGER            I, INCA, J, J1, J2, KB, KB1, KK, KLM, KLU1,
213      $                   KUN, L, MINMN, ML, ML0, MN, MU, MU0, NR, NRT
214       REAL               RA, RB, RC, RS
215 *     ..
216 *     .. External Subroutines ..
217       EXTERNAL           SLARGV, SLARTG, SLARTV, SLASET, SROT, XERBLA
218 *     ..
219 *     .. Intrinsic Functions ..
220       INTRINSIC          MAX, MIN
221 *     ..
222 *     .. External Functions ..
223       LOGICAL            LSAME
224       EXTERNAL           LSAME
225 *     ..
226 *     .. Executable Statements ..
227 *
228 *     Test the input parameters
229 *
230       WANTB = LSAME( VECT, 'B' )
231       WANTQ = LSAME( VECT, 'Q' ) .OR. WANTB
232       WANTPT = LSAME( VECT, 'P' ) .OR. WANTB
233       WANTC = NCC.GT.0
234       KLU1 = KL + KU + 1
235       INFO = 0
236       IF( .NOT.WANTQ .AND. .NOT.WANTPT .AND. .NOT.LSAME( VECT, 'N' ) )
237      $     THEN
238          INFO = -1
239       ELSE IF( M.LT.0 ) THEN
240          INFO = -2
241       ELSE IF( N.LT.0 ) THEN
242          INFO = -3
243       ELSE IF( NCC.LT.0 ) THEN
244          INFO = -4
245       ELSE IF( KL.LT.0 ) THEN
246          INFO = -5
247       ELSE IF( KU.LT.0 ) THEN
248          INFO = -6
249       ELSE IF( LDAB.LT.KLU1 ) THEN
250          INFO = -8
251       ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. LDQ.LT.MAX( 1, M ) ) THEN
252          INFO = -12
253       ELSE IF( LDPT.LT.1 .OR. WANTPT .AND. LDPT.LT.MAX( 1, N ) ) THEN
254          INFO = -14
255       ELSE IF( LDC.LT.1 .OR. WANTC .AND. LDC.LT.MAX( 1, M ) ) THEN
256          INFO = -16
257       END IF
258       IF( INFO.NE.0 ) THEN
259          CALL XERBLA( 'SGBBRD', -INFO )
260          RETURN
261       END IF
262 *
263 *     Initialize Q and P**T to the unit matrix, if needed
264 *
265       IF( WANTQ )
266      $   CALL SLASET( 'Full', M, M, ZERO, ONE, Q, LDQ )
267       IF( WANTPT )
268      $   CALL SLASET( 'Full', N, N, ZERO, ONE, PT, LDPT )
269 *
270 *     Quick return if possible.
271 *
272       IF( M.EQ.0 .OR. N.EQ.0 )
273      $   RETURN
274 *
275       MINMN = MIN( M, N )
276 *
277       IF( KL+KU.GT.1 ) THEN
278 *
279 *        Reduce to upper bidiagonal form if KU > 0; if KU = 0, reduce
280 *        first to lower bidiagonal form and then transform to upper
281 *        bidiagonal
282 *
283          IF( KU.GT.0 ) THEN
284             ML0 = 1
285             MU0 = 2
286          ELSE
287             ML0 = 2
288             MU0 = 1
289          END IF
290 *
291 *        Wherever possible, plane rotations are generated and applied in
292 *        vector operations of length NR over the index set J1:J2:KLU1.
293 *
294 *        The sines of the plane rotations are stored in WORK(1:max(m,n))
295 *        and the cosines in WORK(max(m,n)+1:2*max(m,n)).
296 *
297          MN = MAX( M, N )
298          KLM = MIN( M-1, KL )
299          KUN = MIN( N-1, KU )
300          KB = KLM + KUN
301          KB1 = KB + 1
302          INCA = KB1*LDAB
303          NR = 0
304          J1 = KLM + 2
305          J2 = 1 - KUN
306 *
307          DO 90 I = 1, MINMN
308 *
309 *           Reduce i-th column and i-th row of matrix to bidiagonal form
310 *
311             ML = KLM + 1
312             MU = KUN + 1
313             DO 80 KK = 1, KB
314                J1 = J1 + KB
315                J2 = J2 + KB
316 *
317 *              generate plane rotations to annihilate nonzero elements
318 *              which have been created below the band
319 *
320                IF( NR.GT.0 )
321      $            CALL SLARGV( NR, AB( KLU1, J1-KLM-1 ), INCA,
322      $                         WORK( J1 ), KB1, WORK( MN+J1 ), KB1 )
323 *
324 *              apply plane rotations from the left
325 *
326                DO 10 L = 1, KB
327                   IF( J2-KLM+L-1.GT.N ) THEN
328                      NRT = NR - 1
329                   ELSE
330                      NRT = NR
331                   END IF
332                   IF( NRT.GT.0 )
333      $               CALL SLARTV( NRT, AB( KLU1-L, J1-KLM+L-1 ), INCA,
334      $                            AB( KLU1-L+1, J1-KLM+L-1 ), INCA,
335      $                            WORK( MN+J1 ), WORK( J1 ), KB1 )
336    10          CONTINUE
337 *
338                IF( ML.GT.ML0 ) THEN
339                   IF( ML.LE.M-I+1 ) THEN
340 *
341 *                    generate plane rotation to annihilate a(i+ml-1,i)
342 *                    within the band, and apply rotation from the left
343 *
344                      CALL SLARTG( AB( KU+ML-1, I ), AB( KU+ML, I ),
345      $                            WORK( MN+I+ML-1 ), WORK( I+ML-1 ),
346      $                            RA )
347                      AB( KU+ML-1, I ) = RA
348                      IF( I.LT.N )
349      $                  CALL SROT( MIN( KU+ML-2, N-I ),
350      $                             AB( KU+ML-2, I+1 ), LDAB-1,
351      $                             AB( KU+ML-1, I+1 ), LDAB-1,
352      $                             WORK( MN+I+ML-1 ), WORK( I+ML-1 ) )
353                   END IF
354                   NR = NR + 1
355                   J1 = J1 - KB1
356                END IF
357 *
358                IF( WANTQ ) THEN
359 *
360 *                 accumulate product of plane rotations in Q
361 *
362                   DO 20 J = J1, J2, KB1
363                      CALL SROT( M, Q( 1, J-1 ), 1, Q( 1, J ), 1,
364      $                          WORK( MN+J ), WORK( J ) )
365    20             CONTINUE
366                END IF
367 *
368                IF( WANTC ) THEN
369 *
370 *                 apply plane rotations to C
371 *
372                   DO 30 J = J1, J2, KB1
373                      CALL SROT( NCC, C( J-1, 1 ), LDC, C( J, 1 ), LDC,
374      $                          WORK( MN+J ), WORK( J ) )
375    30             CONTINUE
376                END IF
377 *
378                IF( J2+KUN.GT.N ) THEN
379 *
380 *                 adjust J2 to keep within the bounds of the matrix
381 *
382                   NR = NR - 1
383                   J2 = J2 - KB1
384                END IF
385 *
386                DO 40 J = J1, J2, KB1
387 *
388 *                 create nonzero element a(j-1,j+ku) above the band
389 *                 and store it in WORK(n+1:2*n)
390 *
391                   WORK( J+KUN ) = WORK( J )*AB( 1, J+KUN )
392                   AB( 1, J+KUN ) = WORK( MN+J )*AB( 1, J+KUN )
393    40          CONTINUE
394 *
395 *              generate plane rotations to annihilate nonzero elements
396 *              which have been generated above the band
397 *
398                IF( NR.GT.0 )
399      $            CALL SLARGV( NR, AB( 1, J1+KUN-1 ), INCA,
400      $                         WORK( J1+KUN ), KB1, WORK( MN+J1+KUN ),
401      $                         KB1 )
402 *
403 *              apply plane rotations from the right
404 *
405                DO 50 L = 1, KB
406                   IF( J2+L-1.GT.M ) THEN
407                      NRT = NR - 1
408                   ELSE
409                      NRT = NR
410                   END IF
411                   IF( NRT.GT.0 )
412      $               CALL SLARTV( NRT, AB( L+1, J1+KUN-1 ), INCA,
413      $                            AB( L, J1+KUN ), INCA,
414      $                            WORK( MN+J1+KUN ), WORK( J1+KUN ),
415      $                            KB1 )
416    50          CONTINUE
417 *
418                IF( ML.EQ.ML0 .AND. MU.GT.MU0 ) THEN
419                   IF( MU.LE.N-I+1 ) THEN
420 *
421 *                    generate plane rotation to annihilate a(i,i+mu-1)
422 *                    within the band, and apply rotation from the right
423 *
424                      CALL SLARTG( AB( KU-MU+3, I+MU-2 ),
425      $                            AB( KU-MU+2, I+MU-1 ),
426      $                            WORK( MN+I+MU-1 ), WORK( I+MU-1 ),
427      $                            RA )
428                      AB( KU-MU+3, I+MU-2 ) = RA
429                      CALL SROT( MIN( KL+MU-2, M-I ),
430      $                          AB( KU-MU+4, I+MU-2 ), 1,
431      $                          AB( KU-MU+3, I+MU-1 ), 1,
432      $                          WORK( MN+I+MU-1 ), WORK( I+MU-1 ) )
433                   END IF
434                   NR = NR + 1
435                   J1 = J1 - KB1
436                END IF
437 *
438                IF( WANTPT ) THEN
439 *
440 *                 accumulate product of plane rotations in P**T
441 *
442                   DO 60 J = J1, J2, KB1
443                      CALL SROT( N, PT( J+KUN-1, 1 ), LDPT,
444      $                          PT( J+KUN, 1 ), LDPT, WORK( MN+J+KUN ),
445      $                          WORK( J+KUN ) )
446    60             CONTINUE
447                END IF
448 *
449                IF( J2+KB.GT.M ) THEN
450 *
451 *                 adjust J2 to keep within the bounds of the matrix
452 *
453                   NR = NR - 1
454                   J2 = J2 - KB1
455                END IF
456 *
457                DO 70 J = J1, J2, KB1
458 *
459 *                 create nonzero element a(j+kl+ku,j+ku-1) below the
460 *                 band and store it in WORK(1:n)
461 *
462                   WORK( J+KB ) = WORK( J+KUN )*AB( KLU1, J+KUN )
463                   AB( KLU1, J+KUN ) = WORK( MN+J+KUN )*AB( KLU1, J+KUN )
464    70          CONTINUE
465 *
466                IF( ML.GT.ML0 ) THEN
467                   ML = ML - 1
468                ELSE
469                   MU = MU - 1
470                END IF
471    80       CONTINUE
472    90    CONTINUE
473       END IF
474 *
475       IF( KU.EQ.0 .AND. KL.GT.0 ) THEN
476 *
477 *        A has been reduced to lower bidiagonal form
478 *
479 *        Transform lower bidiagonal form to upper bidiagonal by applying
480 *        plane rotations from the left, storing diagonal elements in D
481 *        and off-diagonal elements in E
482 *
483          DO 100 I = 1, MIN( M-1, N )
484             CALL SLARTG( AB( 1, I ), AB( 2, I ), RC, RS, RA )
485             D( I ) = RA
486             IF( I.LT.N ) THEN
487                E( I ) = RS*AB( 1, I+1 )
488                AB( 1, I+1 ) = RC*AB( 1, I+1 )
489             END IF
490             IF( WANTQ )
491      $         CALL SROT( M, Q( 1, I ), 1, Q( 1, I+1 ), 1, RC, RS )
492             IF( WANTC )
493      $         CALL SROT( NCC, C( I, 1 ), LDC, C( I+1, 1 ), LDC, RC,
494      $                    RS )
495   100    CONTINUE
496          IF( M.LE.N )
497      $      D( M ) = AB( 1, M )
498       ELSE IF( KU.GT.0 ) THEN
499 *
500 *        A has been reduced to upper bidiagonal form
501 *
502          IF( M.LT.N ) THEN
503 *
504 *           Annihilate a(m,m+1) by applying plane rotations from the
505 *           right, storing diagonal elements in D and off-diagonal
506 *           elements in E
507 *
508             RB = AB( KU, M+1 )
509             DO 110 I = M, 1, -1
510                CALL SLARTG( AB( KU+1, I ), RB, RC, RS, RA )
511                D( I ) = RA
512                IF( I.GT.1 ) THEN
513                   RB = -RS*AB( KU, I )
514                   E( I-1 ) = RC*AB( KU, I )
515                END IF
516                IF( WANTPT )
517      $            CALL SROT( N, PT( I, 1 ), LDPT, PT( M+1, 1 ), LDPT,
518      $                       RC, RS )
519   110       CONTINUE
520          ELSE
521 *
522 *           Copy off-diagonal elements to E and diagonal elements to D
523 *
524             DO 120 I = 1, MINMN - 1
525                E( I ) = AB( KU, I+1 )
526   120       CONTINUE
527             DO 130 I = 1, MINMN
528                D( I ) = AB( KU+1, I )
529   130       CONTINUE
530          END IF
531       ELSE
532 *
533 *        A is diagonal. Set elements of E to zero and copy diagonal
534 *        elements to D.
535 *
536          DO 140 I = 1, MINMN - 1
537             E( I ) = ZERO
538   140    CONTINUE
539          DO 150 I = 1, MINMN
540             D( I ) = AB( 1, I )
541   150    CONTINUE
542       END IF
543       RETURN
544 *
545 *     End of SGBBRD
546 *
547       END