LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
sgghd3.f
Go to the documentation of this file.
1 *> \brief \b SGGHD3
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SGGHD3 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgghd3.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgghd3.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgghd3.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SGGHD3( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
22 * LDQ, Z, LDZ, WORK, LWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER COMPQ, COMPZ
26 * INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N, LWORK
27 * ..
28 * .. Array Arguments ..
29 * REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30 * $ Z( LDZ, * ), WORK( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> SGGHD3 reduces a pair of real matrices (A,B) to generalized upper
40 *> Hessenberg form using orthogonal transformations, where A is a
41 *> general matrix and B is upper triangular. The form of the
42 *> generalized eigenvalue problem is
43 *> A*x = lambda*B*x,
44 *> and B is typically made upper triangular by computing its QR
45 *> factorization and moving the orthogonal matrix Q to the left side
46 *> of the equation.
47 *>
48 *> This subroutine simultaneously reduces A to a Hessenberg matrix H:
49 *> Q**T*A*Z = H
50 *> and transforms B to another upper triangular matrix T:
51 *> Q**T*B*Z = T
52 *> in order to reduce the problem to its standard form
53 *> H*y = lambda*T*y
54 *> where y = Z**T*x.
55 *>
56 *> The orthogonal matrices Q and Z are determined as products of Givens
57 *> rotations. They may either be formed explicitly, or they may be
58 *> postmultiplied into input matrices Q1 and Z1, so that
59 *>
60 *> Q1 * A * Z1**T = (Q1*Q) * H * (Z1*Z)**T
61 *>
62 *> Q1 * B * Z1**T = (Q1*Q) * T * (Z1*Z)**T
63 *>
64 *> If Q1 is the orthogonal matrix from the QR factorization of B in the
65 *> original equation A*x = lambda*B*x, then SGGHD3 reduces the original
66 *> problem to generalized Hessenberg form.
67 *>
68 *> This is a blocked variant of SGGHRD, using matrix-matrix
69 *> multiplications for parts of the computation to enhance performance.
70 *> \endverbatim
71 *
72 * Arguments:
73 * ==========
74 *
75 *> \param[in] COMPQ
76 *> \verbatim
77 *> COMPQ is CHARACTER*1
78 *> = 'N': do not compute Q;
79 *> = 'I': Q is initialized to the unit matrix, and the
80 *> orthogonal matrix Q is returned;
81 *> = 'V': Q must contain an orthogonal matrix Q1 on entry,
82 *> and the product Q1*Q is returned.
83 *> \endverbatim
84 *>
85 *> \param[in] COMPZ
86 *> \verbatim
87 *> COMPZ is CHARACTER*1
88 *> = 'N': do not compute Z;
89 *> = 'I': Z is initialized to the unit matrix, and the
90 *> orthogonal matrix Z is returned;
91 *> = 'V': Z must contain an orthogonal matrix Z1 on entry,
92 *> and the product Z1*Z is returned.
93 *> \endverbatim
94 *>
95 *> \param[in] N
96 *> \verbatim
97 *> N is INTEGER
98 *> The order of the matrices A and B. N >= 0.
99 *> \endverbatim
100 *>
101 *> \param[in] ILO
102 *> \verbatim
103 *> ILO is INTEGER
104 *> \endverbatim
105 *>
106 *> \param[in] IHI
107 *> \verbatim
108 *> IHI is INTEGER
109 *>
110 *> ILO and IHI mark the rows and columns of A which are to be
111 *> reduced. It is assumed that A is already upper triangular
112 *> in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are
113 *> normally set by a previous call to SGGBAL; otherwise they
114 *> should be set to 1 and N respectively.
115 *> 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
116 *> \endverbatim
117 *>
118 *> \param[in,out] A
119 *> \verbatim
120 *> A is REAL array, dimension (LDA, N)
121 *> On entry, the N-by-N general matrix to be reduced.
122 *> On exit, the upper triangle and the first subdiagonal of A
123 *> are overwritten with the upper Hessenberg matrix H, and the
124 *> rest is set to zero.
125 *> \endverbatim
126 *>
127 *> \param[in] LDA
128 *> \verbatim
129 *> LDA is INTEGER
130 *> The leading dimension of the array A. LDA >= max(1,N).
131 *> \endverbatim
132 *>
133 *> \param[in,out] B
134 *> \verbatim
135 *> B is REAL array, dimension (LDB, N)
136 *> On entry, the N-by-N upper triangular matrix B.
137 *> On exit, the upper triangular matrix T = Q**T B Z. The
138 *> elements below the diagonal are set to zero.
139 *> \endverbatim
140 *>
141 *> \param[in] LDB
142 *> \verbatim
143 *> LDB is INTEGER
144 *> The leading dimension of the array B. LDB >= max(1,N).
145 *> \endverbatim
146 *>
147 *> \param[in,out] Q
148 *> \verbatim
149 *> Q is REAL array, dimension (LDQ, N)
150 *> On entry, if COMPQ = 'V', the orthogonal matrix Q1,
151 *> typically from the QR factorization of B.
152 *> On exit, if COMPQ='I', the orthogonal matrix Q, and if
153 *> COMPQ = 'V', the product Q1*Q.
154 *> Not referenced if COMPQ='N'.
155 *> \endverbatim
156 *>
157 *> \param[in] LDQ
158 *> \verbatim
159 *> LDQ is INTEGER
160 *> The leading dimension of the array Q.
161 *> LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
162 *> \endverbatim
163 *>
164 *> \param[in,out] Z
165 *> \verbatim
166 *> Z is REAL array, dimension (LDZ, N)
167 *> On entry, if COMPZ = 'V', the orthogonal matrix Z1.
168 *> On exit, if COMPZ='I', the orthogonal matrix Z, and if
169 *> COMPZ = 'V', the product Z1*Z.
170 *> Not referenced if COMPZ='N'.
171 *> \endverbatim
172 *>
173 *> \param[in] LDZ
174 *> \verbatim
175 *> LDZ is INTEGER
176 *> The leading dimension of the array Z.
177 *> LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
178 *> \endverbatim
179 *>
180 *> \param[out] WORK
181 *> \verbatim
182 *> WORK is REAL array, dimension (LWORK)
183 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
184 *> \endverbatim
185 *>
186 *> \param[in] LWORK
187 *> \verbatim
188 *> LWORK is INTEGER
189 *> The length of the array WORK. LWORK >= 1.
190 *> For optimum performance LWORK >= 6*N*NB, where NB is the
191 *> optimal blocksize.
192 *>
193 *> If LWORK = -1, then a workspace query is assumed; the routine
194 *> only calculates the optimal size of the WORK array, returns
195 *> this value as the first entry of the WORK array, and no error
196 *> message related to LWORK is issued by XERBLA.
197 *> \endverbatim
198 *>
199 *> \param[out] INFO
200 *> \verbatim
201 *> INFO is INTEGER
202 *> = 0: successful exit.
203 *> < 0: if INFO = -i, the i-th argument had an illegal value.
204 *> \endverbatim
205 *
206 * Authors:
207 * ========
208 *
209 *> \author Univ. of Tennessee
210 *> \author Univ. of California Berkeley
211 *> \author Univ. of Colorado Denver
212 *> \author NAG Ltd.
213 *
214 *> \ingroup realOTHERcomputational
215 *
216 *> \par Further Details:
217 * =====================
218 *>
219 *> \verbatim
220 *>
221 *> This routine reduces A to Hessenberg form and maintains B in triangular form
222 *> using a blocked variant of Moler and Stewart's original algorithm,
223 *> as described by Kagstrom, Kressner, Quintana-Orti, and Quintana-Orti
224 *> (BIT 2008).
225 *> \endverbatim
226 *>
227 * =====================================================================
228  SUBROUTINE sgghd3( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
229  $ LDQ, Z, LDZ, WORK, LWORK, INFO )
230 *
231 * -- LAPACK computational routine --
232 * -- LAPACK is a software package provided by Univ. of Tennessee, --
233 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
234 *
235  IMPLICIT NONE
236 *
237 * .. Scalar Arguments ..
238  CHARACTER COMPQ, COMPZ
239  INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N, LWORK
240 * ..
241 * .. Array Arguments ..
242  REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
243  $ z( ldz, * ), work( * )
244 * ..
245 *
246 * =====================================================================
247 *
248 * .. Parameters ..
249  REAL ZERO, ONE
250  parameter( zero = 0.0e+0, one = 1.0e+0 )
251 * ..
252 * .. Local Scalars ..
253  LOGICAL BLK22, INITQ, INITZ, LQUERY, WANTQ, WANTZ
254  CHARACTER*1 COMPQ2, COMPZ2
255  INTEGER COLA, I, IERR, J, J0, JCOL, JJ, JROW, K,
256  $ kacc22, len, lwkopt, n2nb, nb, nblst, nbmin,
257  $ nh, nnb, nx, ppw, ppwo, pw, top, topq
258  REAL C, C1, C2, S, S1, S2, TEMP, TEMP1, TEMP2, TEMP3
259 * ..
260 * .. External Functions ..
261  LOGICAL LSAME
262  INTEGER ILAENV
263  EXTERNAL ilaenv, lsame
264 * ..
265 * .. External Subroutines ..
266  EXTERNAL sgghrd, slartg, slaset, sorm22, srot, sgemm,
267  $ sgemv, strmv, slacpy, xerbla
268 * ..
269 * .. Intrinsic Functions ..
270  INTRINSIC real, max
271 * ..
272 * .. Executable Statements ..
273 *
274 * Decode and test the input parameters.
275 *
276  info = 0
277  nb = ilaenv( 1, 'SGGHD3', ' ', n, ilo, ihi, -1 )
278  lwkopt = max( 6*n*nb, 1 )
279  work( 1 ) = real( lwkopt )
280  initq = lsame( compq, 'I' )
281  wantq = initq .OR. lsame( compq, 'V' )
282  initz = lsame( compz, 'I' )
283  wantz = initz .OR. lsame( compz, 'V' )
284  lquery = ( lwork.EQ.-1 )
285 *
286  IF( .NOT.lsame( compq, 'N' ) .AND. .NOT.wantq ) THEN
287  info = -1
288  ELSE IF( .NOT.lsame( compz, 'N' ) .AND. .NOT.wantz ) THEN
289  info = -2
290  ELSE IF( n.LT.0 ) THEN
291  info = -3
292  ELSE IF( ilo.LT.1 ) THEN
293  info = -4
294  ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
295  info = -5
296  ELSE IF( lda.LT.max( 1, n ) ) THEN
297  info = -7
298  ELSE IF( ldb.LT.max( 1, n ) ) THEN
299  info = -9
300  ELSE IF( ( wantq .AND. ldq.LT.n ) .OR. ldq.LT.1 ) THEN
301  info = -11
302  ELSE IF( ( wantz .AND. ldz.LT.n ) .OR. ldz.LT.1 ) THEN
303  info = -13
304  ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
305  info = -15
306  END IF
307  IF( info.NE.0 ) THEN
308  CALL xerbla( 'SGGHD3', -info )
309  RETURN
310  ELSE IF( lquery ) THEN
311  RETURN
312  END IF
313 *
314 * Initialize Q and Z if desired.
315 *
316  IF( initq )
317  $ CALL slaset( 'All', n, n, zero, one, q, ldq )
318  IF( initz )
319  $ CALL slaset( 'All', n, n, zero, one, z, ldz )
320 *
321 * Zero out lower triangle of B.
322 *
323  IF( n.GT.1 )
324  $ CALL slaset( 'Lower', n-1, n-1, zero, zero, b(2, 1), ldb )
325 *
326 * Quick return if possible
327 *
328  nh = ihi - ilo + 1
329  IF( nh.LE.1 ) THEN
330  work( 1 ) = one
331  RETURN
332  END IF
333 *
334 * Determine the blocksize.
335 *
336  nbmin = ilaenv( 2, 'SGGHD3', ' ', n, ilo, ihi, -1 )
337  IF( nb.GT.1 .AND. nb.LT.nh ) THEN
338 *
339 * Determine when to use unblocked instead of blocked code.
340 *
341  nx = max( nb, ilaenv( 3, 'SGGHD3', ' ', n, ilo, ihi, -1 ) )
342  IF( nx.LT.nh ) THEN
343 *
344 * Determine if workspace is large enough for blocked code.
345 *
346  IF( lwork.LT.lwkopt ) THEN
347 *
348 * Not enough workspace to use optimal NB: determine the
349 * minimum value of NB, and reduce NB or force use of
350 * unblocked code.
351 *
352  nbmin = max( 2, ilaenv( 2, 'SGGHD3', ' ', n, ilo, ihi,
353  $ -1 ) )
354  IF( lwork.GE.6*n*nbmin ) THEN
355  nb = lwork / ( 6*n )
356  ELSE
357  nb = 1
358  END IF
359  END IF
360  END IF
361  END IF
362 *
363  IF( nb.LT.nbmin .OR. nb.GE.nh ) THEN
364 *
365 * Use unblocked code below
366 *
367  jcol = ilo
368 *
369  ELSE
370 *
371 * Use blocked code
372 *
373  kacc22 = ilaenv( 16, 'SGGHD3', ' ', n, ilo, ihi, -1 )
374  blk22 = kacc22.EQ.2
375  DO jcol = ilo, ihi-2, nb
376  nnb = min( nb, ihi-jcol-1 )
377 *
378 * Initialize small orthogonal factors that will hold the
379 * accumulated Givens rotations in workspace.
380 * N2NB denotes the number of 2*NNB-by-2*NNB factors
381 * NBLST denotes the (possibly smaller) order of the last
382 * factor.
383 *
384  n2nb = ( ihi-jcol-1 ) / nnb - 1
385  nblst = ihi - jcol - n2nb*nnb
386  CALL slaset( 'All', nblst, nblst, zero, one, work, nblst )
387  pw = nblst * nblst + 1
388  DO i = 1, n2nb
389  CALL slaset( 'All', 2*nnb, 2*nnb, zero, one,
390  $ work( pw ), 2*nnb )
391  pw = pw + 4*nnb*nnb
392  END DO
393 *
394 * Reduce columns JCOL:JCOL+NNB-1 of A to Hessenberg form.
395 *
396  DO j = jcol, jcol+nnb-1
397 *
398 * Reduce Jth column of A. Store cosines and sines in Jth
399 * column of A and B, respectively.
400 *
401  DO i = ihi, j+2, -1
402  temp = a( i-1, j )
403  CALL slartg( temp, a( i, j ), c, s, a( i-1, j ) )
404  a( i, j ) = c
405  b( i, j ) = s
406  END DO
407 *
408 * Accumulate Givens rotations into workspace array.
409 *
410  ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
411  len = 2 + j - jcol
412  jrow = j + n2nb*nnb + 2
413  DO i = ihi, jrow, -1
414  c = a( i, j )
415  s = b( i, j )
416  DO jj = ppw, ppw+len-1
417  temp = work( jj + nblst )
418  work( jj + nblst ) = c*temp - s*work( jj )
419  work( jj ) = s*temp + c*work( jj )
420  END DO
421  len = len + 1
422  ppw = ppw - nblst - 1
423  END DO
424 *
425  ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
426  j0 = jrow - nnb
427  DO jrow = j0, j+2, -nnb
428  ppw = ppwo
429  len = 2 + j - jcol
430  DO i = jrow+nnb-1, jrow, -1
431  c = a( i, j )
432  s = b( i, j )
433  DO jj = ppw, ppw+len-1
434  temp = work( jj + 2*nnb )
435  work( jj + 2*nnb ) = c*temp - s*work( jj )
436  work( jj ) = s*temp + c*work( jj )
437  END DO
438  len = len + 1
439  ppw = ppw - 2*nnb - 1
440  END DO
441  ppwo = ppwo + 4*nnb*nnb
442  END DO
443 *
444 * TOP denotes the number of top rows in A and B that will
445 * not be updated during the next steps.
446 *
447  IF( jcol.LE.2 ) THEN
448  top = 0
449  ELSE
450  top = jcol
451  END IF
452 *
453 * Propagate transformations through B and replace stored
454 * left sines/cosines by right sines/cosines.
455 *
456  DO jj = n, j+1, -1
457 *
458 * Update JJth column of B.
459 *
460  DO i = min( jj+1, ihi ), j+2, -1
461  c = a( i, j )
462  s = b( i, j )
463  temp = b( i, jj )
464  b( i, jj ) = c*temp - s*b( i-1, jj )
465  b( i-1, jj ) = s*temp + c*b( i-1, jj )
466  END DO
467 *
468 * Annihilate B( JJ+1, JJ ).
469 *
470  IF( jj.LT.ihi ) THEN
471  temp = b( jj+1, jj+1 )
472  CALL slartg( temp, b( jj+1, jj ), c, s,
473  $ b( jj+1, jj+1 ) )
474  b( jj+1, jj ) = zero
475  CALL srot( jj-top, b( top+1, jj+1 ), 1,
476  $ b( top+1, jj ), 1, c, s )
477  a( jj+1, j ) = c
478  b( jj+1, j ) = -s
479  END IF
480  END DO
481 *
482 * Update A by transformations from right.
483 * Explicit loop unrolling provides better performance
484 * compared to SLASR.
485 * CALL SLASR( 'Right', 'Variable', 'Backward', IHI-TOP,
486 * $ IHI-J, A( J+2, J ), B( J+2, J ),
487 * $ A( TOP+1, J+1 ), LDA )
488 *
489  jj = mod( ihi-j-1, 3 )
490  DO i = ihi-j-3, jj+1, -3
491  c = a( j+1+i, j )
492  s = -b( j+1+i, j )
493  c1 = a( j+2+i, j )
494  s1 = -b( j+2+i, j )
495  c2 = a( j+3+i, j )
496  s2 = -b( j+3+i, j )
497 *
498  DO k = top+1, ihi
499  temp = a( k, j+i )
500  temp1 = a( k, j+i+1 )
501  temp2 = a( k, j+i+2 )
502  temp3 = a( k, j+i+3 )
503  a( k, j+i+3 ) = c2*temp3 + s2*temp2
504  temp2 = -s2*temp3 + c2*temp2
505  a( k, j+i+2 ) = c1*temp2 + s1*temp1
506  temp1 = -s1*temp2 + c1*temp1
507  a( k, j+i+1 ) = c*temp1 + s*temp
508  a( k, j+i ) = -s*temp1 + c*temp
509  END DO
510  END DO
511 *
512  IF( jj.GT.0 ) THEN
513  DO i = jj, 1, -1
514  CALL srot( ihi-top, a( top+1, j+i+1 ), 1,
515  $ a( top+1, j+i ), 1, a( j+1+i, j ),
516  $ -b( j+1+i, j ) )
517  END DO
518  END IF
519 *
520 * Update (J+1)th column of A by transformations from left.
521 *
522  IF ( j .LT. jcol + nnb - 1 ) THEN
523  len = 1 + j - jcol
524 *
525 * Multiply with the trailing accumulated orthogonal
526 * matrix, which takes the form
527 *
528 * [ U11 U12 ]
529 * U = [ ],
530 * [ U21 U22 ]
531 *
532 * where U21 is a LEN-by-LEN matrix and U12 is lower
533 * triangular.
534 *
535  jrow = ihi - nblst + 1
536  CALL sgemv( 'Transpose', nblst, len, one, work,
537  $ nblst, a( jrow, j+1 ), 1, zero,
538  $ work( pw ), 1 )
539  ppw = pw + len
540  DO i = jrow, jrow+nblst-len-1
541  work( ppw ) = a( i, j+1 )
542  ppw = ppw + 1
543  END DO
544  CALL strmv( 'Lower', 'Transpose', 'Non-unit',
545  $ nblst-len, work( len*nblst + 1 ), nblst,
546  $ work( pw+len ), 1 )
547  CALL sgemv( 'Transpose', len, nblst-len, one,
548  $ work( (len+1)*nblst - len + 1 ), nblst,
549  $ a( jrow+nblst-len, j+1 ), 1, one,
550  $ work( pw+len ), 1 )
551  ppw = pw
552  DO i = jrow, jrow+nblst-1
553  a( i, j+1 ) = work( ppw )
554  ppw = ppw + 1
555  END DO
556 *
557 * Multiply with the other accumulated orthogonal
558 * matrices, which take the form
559 *
560 * [ U11 U12 0 ]
561 * [ ]
562 * U = [ U21 U22 0 ],
563 * [ ]
564 * [ 0 0 I ]
565 *
566 * where I denotes the (NNB-LEN)-by-(NNB-LEN) identity
567 * matrix, U21 is a LEN-by-LEN upper triangular matrix
568 * and U12 is an NNB-by-NNB lower triangular matrix.
569 *
570  ppwo = 1 + nblst*nblst
571  j0 = jrow - nnb
572  DO jrow = j0, jcol+1, -nnb
573  ppw = pw + len
574  DO i = jrow, jrow+nnb-1
575  work( ppw ) = a( i, j+1 )
576  ppw = ppw + 1
577  END DO
578  ppw = pw
579  DO i = jrow+nnb, jrow+nnb+len-1
580  work( ppw ) = a( i, j+1 )
581  ppw = ppw + 1
582  END DO
583  CALL strmv( 'Upper', 'Transpose', 'Non-unit', len,
584  $ work( ppwo + nnb ), 2*nnb, work( pw ),
585  $ 1 )
586  CALL strmv( 'Lower', 'Transpose', 'Non-unit', nnb,
587  $ work( ppwo + 2*len*nnb ),
588  $ 2*nnb, work( pw + len ), 1 )
589  CALL sgemv( 'Transpose', nnb, len, one,
590  $ work( ppwo ), 2*nnb, a( jrow, j+1 ), 1,
591  $ one, work( pw ), 1 )
592  CALL sgemv( 'Transpose', len, nnb, one,
593  $ work( ppwo + 2*len*nnb + nnb ), 2*nnb,
594  $ a( jrow+nnb, j+1 ), 1, one,
595  $ work( pw+len ), 1 )
596  ppw = pw
597  DO i = jrow, jrow+len+nnb-1
598  a( i, j+1 ) = work( ppw )
599  ppw = ppw + 1
600  END DO
601  ppwo = ppwo + 4*nnb*nnb
602  END DO
603  END IF
604  END DO
605 *
606 * Apply accumulated orthogonal matrices to A.
607 *
608  cola = n - jcol - nnb + 1
609  j = ihi - nblst + 1
610  CALL sgemm( 'Transpose', 'No Transpose', nblst,
611  $ cola, nblst, one, work, nblst,
612  $ a( j, jcol+nnb ), lda, zero, work( pw ),
613  $ nblst )
614  CALL slacpy( 'All', nblst, cola, work( pw ), nblst,
615  $ a( j, jcol+nnb ), lda )
616  ppwo = nblst*nblst + 1
617  j0 = j - nnb
618  DO j = j0, jcol+1, -nnb
619  IF ( blk22 ) THEN
620 *
621 * Exploit the structure of
622 *
623 * [ U11 U12 ]
624 * U = [ ]
625 * [ U21 U22 ],
626 *
627 * where all blocks are NNB-by-NNB, U21 is upper
628 * triangular and U12 is lower triangular.
629 *
630  CALL sorm22( 'Left', 'Transpose', 2*nnb, cola, nnb,
631  $ nnb, work( ppwo ), 2*nnb,
632  $ a( j, jcol+nnb ), lda, work( pw ),
633  $ lwork-pw+1, ierr )
634  ELSE
635 *
636 * Ignore the structure of U.
637 *
638  CALL sgemm( 'Transpose', 'No Transpose', 2*nnb,
639  $ cola, 2*nnb, one, work( ppwo ), 2*nnb,
640  $ a( j, jcol+nnb ), lda, zero, work( pw ),
641  $ 2*nnb )
642  CALL slacpy( 'All', 2*nnb, cola, work( pw ), 2*nnb,
643  $ a( j, jcol+nnb ), lda )
644  END IF
645  ppwo = ppwo + 4*nnb*nnb
646  END DO
647 *
648 * Apply accumulated orthogonal matrices to Q.
649 *
650  IF( wantq ) THEN
651  j = ihi - nblst + 1
652  IF ( initq ) THEN
653  topq = max( 2, j - jcol + 1 )
654  nh = ihi - topq + 1
655  ELSE
656  topq = 1
657  nh = n
658  END IF
659  CALL sgemm( 'No Transpose', 'No Transpose', nh,
660  $ nblst, nblst, one, q( topq, j ), ldq,
661  $ work, nblst, zero, work( pw ), nh )
662  CALL slacpy( 'All', nh, nblst, work( pw ), nh,
663  $ q( topq, j ), ldq )
664  ppwo = nblst*nblst + 1
665  j0 = j - nnb
666  DO j = j0, jcol+1, -nnb
667  IF ( initq ) THEN
668  topq = max( 2, j - jcol + 1 )
669  nh = ihi - topq + 1
670  END IF
671  IF ( blk22 ) THEN
672 *
673 * Exploit the structure of U.
674 *
675  CALL sorm22( 'Right', 'No Transpose', nh, 2*nnb,
676  $ nnb, nnb, work( ppwo ), 2*nnb,
677  $ q( topq, j ), ldq, work( pw ),
678  $ lwork-pw+1, ierr )
679  ELSE
680 *
681 * Ignore the structure of U.
682 *
683  CALL sgemm( 'No Transpose', 'No Transpose', nh,
684  $ 2*nnb, 2*nnb, one, q( topq, j ), ldq,
685  $ work( ppwo ), 2*nnb, zero, work( pw ),
686  $ nh )
687  CALL slacpy( 'All', nh, 2*nnb, work( pw ), nh,
688  $ q( topq, j ), ldq )
689  END IF
690  ppwo = ppwo + 4*nnb*nnb
691  END DO
692  END IF
693 *
694 * Accumulate right Givens rotations if required.
695 *
696  IF ( wantz .OR. top.GT.0 ) THEN
697 *
698 * Initialize small orthogonal factors that will hold the
699 * accumulated Givens rotations in workspace.
700 *
701  CALL slaset( 'All', nblst, nblst, zero, one, work,
702  $ nblst )
703  pw = nblst * nblst + 1
704  DO i = 1, n2nb
705  CALL slaset( 'All', 2*nnb, 2*nnb, zero, one,
706  $ work( pw ), 2*nnb )
707  pw = pw + 4*nnb*nnb
708  END DO
709 *
710 * Accumulate Givens rotations into workspace array.
711 *
712  DO j = jcol, jcol+nnb-1
713  ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
714  len = 2 + j - jcol
715  jrow = j + n2nb*nnb + 2
716  DO i = ihi, jrow, -1
717  c = a( i, j )
718  a( i, j ) = zero
719  s = b( i, j )
720  b( i, j ) = zero
721  DO jj = ppw, ppw+len-1
722  temp = work( jj + nblst )
723  work( jj + nblst ) = c*temp - s*work( jj )
724  work( jj ) = s*temp + c*work( jj )
725  END DO
726  len = len + 1
727  ppw = ppw - nblst - 1
728  END DO
729 *
730  ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
731  j0 = jrow - nnb
732  DO jrow = j0, j+2, -nnb
733  ppw = ppwo
734  len = 2 + j - jcol
735  DO i = jrow+nnb-1, jrow, -1
736  c = a( i, j )
737  a( i, j ) = zero
738  s = b( i, j )
739  b( i, j ) = zero
740  DO jj = ppw, ppw+len-1
741  temp = work( jj + 2*nnb )
742  work( jj + 2*nnb ) = c*temp - s*work( jj )
743  work( jj ) = s*temp + c*work( jj )
744  END DO
745  len = len + 1
746  ppw = ppw - 2*nnb - 1
747  END DO
748  ppwo = ppwo + 4*nnb*nnb
749  END DO
750  END DO
751  ELSE
752 *
753  CALL slaset( 'Lower', ihi - jcol - 1, nnb, zero, zero,
754  $ a( jcol + 2, jcol ), lda )
755  CALL slaset( 'Lower', ihi - jcol - 1, nnb, zero, zero,
756  $ b( jcol + 2, jcol ), ldb )
757  END IF
758 *
759 * Apply accumulated orthogonal matrices to A and B.
760 *
761  IF ( top.GT.0 ) THEN
762  j = ihi - nblst + 1
763  CALL sgemm( 'No Transpose', 'No Transpose', top,
764  $ nblst, nblst, one, a( 1, j ), lda,
765  $ work, nblst, zero, work( pw ), top )
766  CALL slacpy( 'All', top, nblst, work( pw ), top,
767  $ a( 1, j ), lda )
768  ppwo = nblst*nblst + 1
769  j0 = j - nnb
770  DO j = j0, jcol+1, -nnb
771  IF ( blk22 ) THEN
772 *
773 * Exploit the structure of U.
774 *
775  CALL sorm22( 'Right', 'No Transpose', top, 2*nnb,
776  $ nnb, nnb, work( ppwo ), 2*nnb,
777  $ a( 1, j ), lda, work( pw ),
778  $ lwork-pw+1, ierr )
779  ELSE
780 *
781 * Ignore the structure of U.
782 *
783  CALL sgemm( 'No Transpose', 'No Transpose', top,
784  $ 2*nnb, 2*nnb, one, a( 1, j ), lda,
785  $ work( ppwo ), 2*nnb, zero,
786  $ work( pw ), top )
787  CALL slacpy( 'All', top, 2*nnb, work( pw ), top,
788  $ a( 1, j ), lda )
789  END IF
790  ppwo = ppwo + 4*nnb*nnb
791  END DO
792 *
793  j = ihi - nblst + 1
794  CALL sgemm( 'No Transpose', 'No Transpose', top,
795  $ nblst, nblst, one, b( 1, j ), ldb,
796  $ work, nblst, zero, work( pw ), top )
797  CALL slacpy( 'All', top, nblst, work( pw ), top,
798  $ b( 1, j ), ldb )
799  ppwo = nblst*nblst + 1
800  j0 = j - nnb
801  DO j = j0, jcol+1, -nnb
802  IF ( blk22 ) THEN
803 *
804 * Exploit the structure of U.
805 *
806  CALL sorm22( 'Right', 'No Transpose', top, 2*nnb,
807  $ nnb, nnb, work( ppwo ), 2*nnb,
808  $ b( 1, j ), ldb, work( pw ),
809  $ lwork-pw+1, ierr )
810  ELSE
811 *
812 * Ignore the structure of U.
813 *
814  CALL sgemm( 'No Transpose', 'No Transpose', top,
815  $ 2*nnb, 2*nnb, one, b( 1, j ), ldb,
816  $ work( ppwo ), 2*nnb, zero,
817  $ work( pw ), top )
818  CALL slacpy( 'All', top, 2*nnb, work( pw ), top,
819  $ b( 1, j ), ldb )
820  END IF
821  ppwo = ppwo + 4*nnb*nnb
822  END DO
823  END IF
824 *
825 * Apply accumulated orthogonal matrices to Z.
826 *
827  IF( wantz ) THEN
828  j = ihi - nblst + 1
829  IF ( initq ) THEN
830  topq = max( 2, j - jcol + 1 )
831  nh = ihi - topq + 1
832  ELSE
833  topq = 1
834  nh = n
835  END IF
836  CALL sgemm( 'No Transpose', 'No Transpose', nh,
837  $ nblst, nblst, one, z( topq, j ), ldz,
838  $ work, nblst, zero, work( pw ), nh )
839  CALL slacpy( 'All', nh, nblst, work( pw ), nh,
840  $ z( topq, j ), ldz )
841  ppwo = nblst*nblst + 1
842  j0 = j - nnb
843  DO j = j0, jcol+1, -nnb
844  IF ( initq ) THEN
845  topq = max( 2, j - jcol + 1 )
846  nh = ihi - topq + 1
847  END IF
848  IF ( blk22 ) THEN
849 *
850 * Exploit the structure of U.
851 *
852  CALL sorm22( 'Right', 'No Transpose', nh, 2*nnb,
853  $ nnb, nnb, work( ppwo ), 2*nnb,
854  $ z( topq, j ), ldz, work( pw ),
855  $ lwork-pw+1, ierr )
856  ELSE
857 *
858 * Ignore the structure of U.
859 *
860  CALL sgemm( 'No Transpose', 'No Transpose', nh,
861  $ 2*nnb, 2*nnb, one, z( topq, j ), ldz,
862  $ work( ppwo ), 2*nnb, zero, work( pw ),
863  $ nh )
864  CALL slacpy( 'All', nh, 2*nnb, work( pw ), nh,
865  $ z( topq, j ), ldz )
866  END IF
867  ppwo = ppwo + 4*nnb*nnb
868  END DO
869  END IF
870  END DO
871  END IF
872 *
873 * Use unblocked code to reduce the rest of the matrix
874 * Avoid re-initialization of modified Q and Z.
875 *
876  compq2 = compq
877  compz2 = compz
878  IF ( jcol.NE.ilo ) THEN
879  IF ( wantq )
880  $ compq2 = 'V'
881  IF ( wantz )
882  $ compz2 = 'V'
883  END IF
884 *
885  IF ( jcol.LT.ihi )
886  $ CALL sgghrd( compq2, compz2, n, jcol, ihi, a, lda, b, ldb, q,
887  $ ldq, z, ldz, ierr )
888  work( 1 ) = real( lwkopt )
889 *
890  RETURN
891 *
892 * End of SGGHD3
893 *
894  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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
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 sorm22(SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC, WORK, LWORK, INFO)
SORM22 multiplies a general matrix by a banded orthogonal matrix.
Definition: sorm22.f:163
subroutine sgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SGGHD3
Definition: sgghd3.f:230
subroutine sgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
SGGHRD
Definition: sgghrd.f:207
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:92
subroutine strmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
STRMV
Definition: strmv.f:147
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187