LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ssytrd_sy2sb.f
Go to the documentation of this file.
1*> \brief \b SSYTRD_SY2SB
2*
3* @generated from zhetrd_he2hb.f, fortran z -> s, Wed Dec 7 08:22:40 2016
4*
5* =========== DOCUMENTATION ===========
6*
7* Online html documentation available at
8* http://www.netlib.org/lapack/explore-html/
9*
10*> \htmlonly
11*> Download SSYTRD_SY2SB + dependencies
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssytrd_sy2sb.f">
13*> [TGZ]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssytrd_sy2sb.f">
15*> [ZIP]</a>
16*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssytrd_sy2sb.f">
17*> [TXT]</a>
18*> \endhtmlonly
19*
20* Definition:
21* ===========
22*
23* SUBROUTINE SSYTRD_SY2SB( UPLO, N, KD, A, LDA, AB, LDAB, TAU,
24* WORK, LWORK, INFO )
25*
26* IMPLICIT NONE
27*
28* .. Scalar Arguments ..
29* CHARACTER UPLO
30* INTEGER INFO, LDA, LDAB, LWORK, N, KD
31* ..
32* .. Array Arguments ..
33* REAL A( LDA, * ), AB( LDAB, * ),
34* TAU( * ), WORK( * )
35* ..
36*
37*
38*> \par Purpose:
39* =============
40*>
41*> \verbatim
42*>
43*> SSYTRD_SY2SB reduces a real symmetric matrix A to real symmetric
44*> band-diagonal form AB by a orthogonal similarity transformation:
45*> Q**T * A * Q = AB.
46*> \endverbatim
47*
48* Arguments:
49* ==========
50*
51*> \param[in] UPLO
52*> \verbatim
53*> UPLO is CHARACTER*1
54*> = 'U': Upper triangle of A is stored;
55*> = 'L': Lower triangle of A is stored.
56*> \endverbatim
57*>
58*> \param[in] N
59*> \verbatim
60*> N is INTEGER
61*> The order of the matrix A. N >= 0.
62*> \endverbatim
63*>
64*> \param[in] KD
65*> \verbatim
66*> KD is INTEGER
67*> The number of superdiagonals of the reduced matrix if UPLO = 'U',
68*> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
69*> The reduced matrix is stored in the array AB.
70*> \endverbatim
71*>
72*> \param[in,out] A
73*> \verbatim
74*> A is REAL array, dimension (LDA,N)
75*> On entry, the symmetric matrix A. If UPLO = 'U', the leading
76*> N-by-N upper triangular part of A contains the upper
77*> triangular part of the matrix A, and the strictly lower
78*> triangular part of A is not referenced. If UPLO = 'L', the
79*> leading N-by-N lower triangular part of A contains the lower
80*> triangular part of the matrix A, and the strictly upper
81*> triangular part of A is not referenced.
82*> On exit, if UPLO = 'U', the diagonal and first superdiagonal
83*> of A are overwritten by the corresponding elements of the
84*> tridiagonal matrix T, and the elements above the first
85*> superdiagonal, with the array TAU, represent the orthogonal
86*> matrix Q as a product of elementary reflectors; if UPLO
87*> = 'L', the diagonal and first subdiagonal of A are over-
88*> written by the corresponding elements of the tridiagonal
89*> matrix T, and the elements below the first subdiagonal, with
90*> the array TAU, represent the orthogonal matrix Q as a product
91*> of elementary reflectors. See Further Details.
92*> \endverbatim
93*>
94*> \param[in] LDA
95*> \verbatim
96*> LDA is INTEGER
97*> The leading dimension of the array A. LDA >= max(1,N).
98*> \endverbatim
99*>
100*> \param[out] AB
101*> \verbatim
102*> AB is REAL array, dimension (LDAB,N)
103*> On exit, the upper or lower triangle of the symmetric band
104*> matrix A, stored in the first KD+1 rows of the array. The
105*> j-th column of A is stored in the j-th column of the array AB
106*> as follows:
107*> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
108*> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
109*> \endverbatim
110*>
111*> \param[in] LDAB
112*> \verbatim
113*> LDAB is INTEGER
114*> The leading dimension of the array AB. LDAB >= KD+1.
115*> \endverbatim
116*>
117*> \param[out] TAU
118*> \verbatim
119*> TAU is REAL array, dimension (N-KD)
120*> The scalar factors of the elementary reflectors (see Further
121*> Details).
122*> \endverbatim
123*>
124*> \param[out] WORK
125*> \verbatim
126*> WORK is REAL array, dimension (LWORK)
127*> On exit, if INFO = 0, or if LWORK=-1,
128*> WORK(1) returns the size of LWORK.
129*> \endverbatim
130*>
131*> \param[in] LWORK
132*> \verbatim
133*> LWORK is INTEGER
134*> The dimension of the array WORK which should be calculated
135*> by a workspace query. LWORK = MAX(1, LWORK_QUERY)
136*> If LWORK = -1, then a workspace query is assumed; the routine
137*> only calculates the optimal size of the WORK array, returns
138*> this value as the first entry of the WORK array, and no error
139*> message related to LWORK is issued by XERBLA.
140*> LWORK_QUERY = N*KD + N*max(KD,FACTOPTNB) + 2*KD*KD
141*> where FACTOPTNB is the blocking used by the QR or LQ
142*> algorithm, usually FACTOPTNB=128 is a good choice otherwise
143*> putting LWORK=-1 will provide the size of WORK.
144*> \endverbatim
145*>
146*> \param[out] INFO
147*> \verbatim
148*> INFO is INTEGER
149*> = 0: successful exit
150*> < 0: if INFO = -i, the i-th argument had an illegal value
151*> \endverbatim
152*
153* Authors:
154* ========
155*
156*> \author Univ. of Tennessee
157*> \author Univ. of California Berkeley
158*> \author Univ. of Colorado Denver
159*> \author NAG Ltd.
160*
161*> \ingroup hetrd_he2hb
162*
163*> \par Further Details:
164* =====================
165*>
166*> \verbatim
167*>
168*> Implemented by Azzam Haidar.
169*>
170*> All details are available on technical report, SC11, SC13 papers.
171*>
172*> Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
173*> Parallel reduction to condensed forms for symmetric eigenvalue problems
174*> using aggregated fine-grained and memory-aware kernels. In Proceedings
175*> of 2011 International Conference for High Performance Computing,
176*> Networking, Storage and Analysis (SC '11), New York, NY, USA,
177*> Article 8 , 11 pages.
178*> http://doi.acm.org/10.1145/2063384.2063394
179*>
180*> A. Haidar, J. Kurzak, P. Luszczek, 2013.
181*> An improved parallel singular value algorithm and its implementation
182*> for multicore hardware, In Proceedings of 2013 International Conference
183*> for High Performance Computing, Networking, Storage and Analysis (SC '13).
184*> Denver, Colorado, USA, 2013.
185*> Article 90, 12 pages.
186*> http://doi.acm.org/10.1145/2503210.2503292
187*>
188*> A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
189*> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
190*> calculations based on fine-grained memory aware tasks.
191*> International Journal of High Performance Computing Applications.
192*> Volume 28 Issue 2, Pages 196-209, May 2014.
193*> http://hpc.sagepub.com/content/28/2/196
194*>
195*> \endverbatim
196*>
197*> \verbatim
198*>
199*> If UPLO = 'U', the matrix Q is represented as a product of elementary
200*> reflectors
201*>
202*> Q = H(k)**T . . . H(2)**T H(1)**T, where k = n-kd.
203*>
204*> Each H(i) has the form
205*>
206*> H(i) = I - tau * v * v**T
207*>
208*> where tau is a real scalar, and v is a real vector with
209*> v(1:i+kd-1) = 0 and v(i+kd) = 1; conjg(v(i+kd+1:n)) is stored on exit in
210*> A(i,i+kd+1:n), and tau in TAU(i).
211*>
212*> If UPLO = 'L', the matrix Q is represented as a product of elementary
213*> reflectors
214*>
215*> Q = H(1) H(2) . . . H(k), where k = n-kd.
216*>
217*> Each H(i) has the form
218*>
219*> H(i) = I - tau * v * v**T
220*>
221*> where tau is a real scalar, and v is a real vector with
222*> v(kd+1:i) = 0 and v(i+kd+1) = 1; v(i+kd+2:n) is stored on exit in
223*> A(i+kd+2:n,i), and tau in TAU(i).
224*>
225*> The contents of A on exit are illustrated by the following examples
226*> with n = 5:
227*>
228*> if UPLO = 'U': if UPLO = 'L':
229*>
230*> ( ab ab/v1 v1 v1 v1 ) ( ab )
231*> ( ab ab/v2 v2 v2 ) ( ab/v1 ab )
232*> ( ab ab/v3 v3 ) ( v1 ab/v2 ab )
233*> ( ab ab/v4 ) ( v1 v2 ab/v3 ab )
234*> ( ab ) ( v1 v2 v3 ab/v4 ab )
235*>
236*> where d and e denote diagonal and off-diagonal elements of T, and vi
237*> denotes an element of the vector defining H(i).
238*> \endverbatim
239*>
240* =====================================================================
241 SUBROUTINE ssytrd_sy2sb( UPLO, N, KD, A, LDA, AB, LDAB, TAU,
242 $ WORK, LWORK, INFO )
243*
244 IMPLICIT NONE
245*
246* -- LAPACK computational routine --
247* -- LAPACK is a software package provided by Univ. of Tennessee, --
248* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
249*
250* .. Scalar Arguments ..
251 CHARACTER UPLO
252 INTEGER INFO, LDA, LDAB, LWORK, N, KD
253* ..
254* .. Array Arguments ..
255 REAL A( LDA, * ), AB( LDAB, * ),
256 $ tau( * ), work( * )
257* ..
258*
259* =====================================================================
260*
261* .. Parameters ..
262 REAL RONE
263 REAL ZERO, ONE, HALF
264 parameter( rone = 1.0e+0,
265 $ zero = 0.0e+0,
266 $ one = 1.0e+0,
267 $ half = 0.5e+0 )
268* ..
269* .. Local Scalars ..
270 LOGICAL LQUERY, UPPER
271 INTEGER I, J, IINFO, LWMIN, PN, PK, LK,
272 $ ldt, ldw, lds2, lds1,
273 $ ls2, ls1, lw, lt,
274 $ tpos, wpos, s2pos, s1pos
275* ..
276* .. External Subroutines ..
277 EXTERNAL xerbla, ssyr2k, ssymm, sgemm, scopy,
279* ..
280* .. Intrinsic Functions ..
281 INTRINSIC min, max
282* ..
283* .. External Functions ..
284 LOGICAL LSAME
285 INTEGER ILAENV2STAGE
286 REAL SROUNDUP_LWORK
287 EXTERNAL lsame, ilaenv2stage, sroundup_lwork
288* ..
289* .. Executable Statements ..
290*
291* Determine the minimal workspace size required
292* and test the input parameters
293*
294 info = 0
295 upper = lsame( uplo, 'U' )
296 lquery = ( lwork.EQ.-1 )
297 lwmin = ilaenv2stage( 4, 'SSYTRD_SY2SB', '', n, kd, -1, -1 )
298
299 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
300 info = -1
301 ELSE IF( n.LT.0 ) THEN
302 info = -2
303 ELSE IF( kd.LT.0 ) THEN
304 info = -3
305 ELSE IF( lda.LT.max( 1, n ) ) THEN
306 info = -5
307 ELSE IF( ldab.LT.max( 1, kd+1 ) ) THEN
308 info = -7
309 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
310 info = -10
311 END IF
312*
313 IF( info.NE.0 ) THEN
314 CALL xerbla( 'SSYTRD_SY2SB', -info )
315 RETURN
316 ELSE IF( lquery ) THEN
317 work( 1 ) = sroundup_lwork(lwmin)
318 RETURN
319 END IF
320*
321* Quick return if possible
322* Copy the upper/lower portion of A into AB
323*
324 IF( n.LE.kd+1 ) THEN
325 IF( upper ) THEN
326 DO 100 i = 1, n
327 lk = min( kd+1, i )
328 CALL scopy( lk, a( i-lk+1, i ), 1,
329 $ ab( kd+1-lk+1, i ), 1 )
330 100 CONTINUE
331 ELSE
332 DO 110 i = 1, n
333 lk = min( kd+1, n-i+1 )
334 CALL scopy( lk, a( i, i ), 1, ab( 1, i ), 1 )
335 110 CONTINUE
336 ENDIF
337 work( 1 ) = 1
338 RETURN
339 END IF
340*
341* Determine the pointer position for the workspace
342*
343 ldt = kd
344 lds1 = kd
345 lt = ldt*kd
346 lw = n*kd
347 ls1 = lds1*kd
348 ls2 = lwmin - lt - lw - ls1
349* LS2 = N*MAX(KD,FACTOPTNB)
350 tpos = 1
351 wpos = tpos + lt
352 s1pos = wpos + lw
353 s2pos = s1pos + ls1
354 IF( upper ) THEN
355 ldw = kd
356 lds2 = kd
357 ELSE
358 ldw = n
359 lds2 = n
360 ENDIF
361*
362*
363* Set the workspace of the triangular matrix T to zero once such a
364* way every time T is generated the upper/lower portion will be always zero
365*
366 CALL slaset( "A", ldt, kd, zero, zero, work( tpos ), ldt )
367*
368 IF( upper ) THEN
369 DO 10 i = 1, n - kd, kd
370 pn = n-i-kd+1
371 pk = min( n-i-kd+1, kd )
372*
373* Compute the LQ factorization of the current block
374*
375 CALL sgelqf( kd, pn, a( i, i+kd ), lda,
376 $ tau( i ), work( s2pos ), ls2, iinfo )
377*
378* Copy the upper portion of A into AB
379*
380 DO 20 j = i, i+pk-1
381 lk = min( kd, n-j ) + 1
382 CALL scopy( lk, a( j, j ), lda, ab( kd+1, j ), ldab-1 )
383 20 CONTINUE
384*
385 CALL slaset( 'Lower', pk, pk, zero, one,
386 $ a( i, i+kd ), lda )
387*
388* Form the matrix T
389*
390 CALL slarft( 'Forward', 'Rowwise', pn, pk,
391 $ a( i, i+kd ), lda, tau( i ),
392 $ work( tpos ), ldt )
393*
394* Compute W:
395*
396 CALL sgemm( 'Conjugate', 'No transpose', pk, pn, pk,
397 $ one, work( tpos ), ldt,
398 $ a( i, i+kd ), lda,
399 $ zero, work( s2pos ), lds2 )
400*
401 CALL ssymm( 'Right', uplo, pk, pn,
402 $ one, a( i+kd, i+kd ), lda,
403 $ work( s2pos ), lds2,
404 $ zero, work( wpos ), ldw )
405*
406 CALL sgemm( 'No transpose', 'Conjugate', pk, pk, pn,
407 $ one, work( wpos ), ldw,
408 $ work( s2pos ), lds2,
409 $ zero, work( s1pos ), lds1 )
410*
411 CALL sgemm( 'No transpose', 'No transpose', pk, pn, pk,
412 $ -half, work( s1pos ), lds1,
413 $ a( i, i+kd ), lda,
414 $ one, work( wpos ), ldw )
415*
416*
417* Update the unreduced submatrix A(i+kd:n,i+kd:n), using
418* an update of the form: A := A - V'*W - W'*V
419*
420 CALL ssyr2k( uplo, 'Conjugate', pn, pk,
421 $ -one, a( i, i+kd ), lda,
422 $ work( wpos ), ldw,
423 $ rone, a( i+kd, i+kd ), lda )
424 10 CONTINUE
425*
426* Copy the upper band to AB which is the band storage matrix
427*
428 DO 30 j = n-kd+1, n
429 lk = min(kd, n-j) + 1
430 CALL scopy( lk, a( j, j ), lda, ab( kd+1, j ), ldab-1 )
431 30 CONTINUE
432*
433 ELSE
434*
435* Reduce the lower triangle of A to lower band matrix
436*
437 DO 40 i = 1, n - kd, kd
438 pn = n-i-kd+1
439 pk = min( n-i-kd+1, kd )
440*
441* Compute the QR factorization of the current block
442*
443 CALL sgeqrf( pn, kd, a( i+kd, i ), lda,
444 $ tau( i ), work( s2pos ), ls2, iinfo )
445*
446* Copy the upper portion of A into AB
447*
448 DO 50 j = i, i+pk-1
449 lk = min( kd, n-j ) + 1
450 CALL scopy( lk, a( j, j ), 1, ab( 1, j ), 1 )
451 50 CONTINUE
452*
453 CALL slaset( 'Upper', pk, pk, zero, one,
454 $ a( i+kd, i ), lda )
455*
456* Form the matrix T
457*
458 CALL slarft( 'Forward', 'Columnwise', pn, pk,
459 $ a( i+kd, i ), lda, tau( i ),
460 $ work( tpos ), ldt )
461*
462* Compute W:
463*
464 CALL sgemm( 'No transpose', 'No transpose', pn, pk, pk,
465 $ one, a( i+kd, i ), lda,
466 $ work( tpos ), ldt,
467 $ zero, work( s2pos ), lds2 )
468*
469 CALL ssymm( 'Left', uplo, pn, pk,
470 $ one, a( i+kd, i+kd ), lda,
471 $ work( s2pos ), lds2,
472 $ zero, work( wpos ), ldw )
473*
474 CALL sgemm( 'Conjugate', 'No transpose', pk, pk, pn,
475 $ one, work( s2pos ), lds2,
476 $ work( wpos ), ldw,
477 $ zero, work( s1pos ), lds1 )
478*
479 CALL sgemm( 'No transpose', 'No transpose', pn, pk, pk,
480 $ -half, a( i+kd, i ), lda,
481 $ work( s1pos ), lds1,
482 $ one, work( wpos ), ldw )
483*
484*
485* Update the unreduced submatrix A(i+kd:n,i+kd:n), using
486* an update of the form: A := A - V*W' - W*V'
487*
488 CALL ssyr2k( uplo, 'No transpose', pn, pk,
489 $ -one, a( i+kd, i ), lda,
490 $ work( wpos ), ldw,
491 $ rone, a( i+kd, i+kd ), lda )
492* ==================================================================
493* RESTORE A FOR COMPARISON AND CHECKING TO BE REMOVED
494* DO 45 J = I, I+PK-1
495* LK = MIN( KD, N-J ) + 1
496* CALL SCOPY( LK, AB( 1, J ), 1, A( J, J ), 1 )
497* 45 CONTINUE
498* ==================================================================
499 40 CONTINUE
500*
501* Copy the lower band to AB which is the band storage matrix
502*
503 DO 60 j = n-kd+1, n
504 lk = min(kd, n-j) + 1
505 CALL scopy( lk, a( j, j ), 1, ab( 1, j ), 1 )
506 60 CONTINUE
507
508 END IF
509*
510 work( 1 ) = sroundup_lwork(lwmin)
511 RETURN
512*
513* End of SSYTRD_SY2SB
514*
515 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine sgelqf(m, n, a, lda, tau, work, lwork, info)
SGELQF
Definition sgelqf.f:143
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine sgeqrf(m, n, a, lda, tau, work, lwork, info)
SGEQRF
Definition sgeqrf.f:146
subroutine ssymm(side, uplo, m, n, alpha, a, lda, b, ldb, beta, c, ldc)
SSYMM
Definition ssymm.f:189
subroutine ssyr2k(uplo, trans, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SSYR2K
Definition ssyr2k.f:192
subroutine ssytrd_sy2sb(uplo, n, kd, a, lda, ab, ldab, tau, work, lwork, info)
SSYTRD_SY2SB
subroutine slarft(direct, storev, n, k, v, ldv, tau, t, ldt)
SLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition slarft.f:163
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