LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ssytrf_rk.f
Go to the documentation of this file.
1*> \brief \b SSYTRF_RK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Kaufman (rook) diagonal pivoting method (BLAS3 blocked algorithm).
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SSYTRF_RK + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssytrf_rk.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssytrf_rk.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssytrf_rk.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SSYTRF_RK( UPLO, N, A, LDA, E, IPIV, WORK, LWORK,
22* INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER UPLO
26* INTEGER INFO, LDA, LWORK, N
27* ..
28* .. Array Arguments ..
29* INTEGER IPIV( * )
30* REAL A( LDA, * ), E ( * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*> SSYTRF_RK computes the factorization of a real symmetric matrix A
39*> using the bounded Bunch-Kaufman (rook) diagonal pivoting method:
40*>
41*> A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T),
42*>
43*> where U (or L) is unit upper (or lower) triangular matrix,
44*> U**T (or L**T) is the transpose of U (or L), P is a permutation
45*> matrix, P**T is the transpose of P, and D is symmetric and block
46*> diagonal with 1-by-1 and 2-by-2 diagonal blocks.
47*>
48*> This is the blocked version of the algorithm, calling Level 3 BLAS.
49*> For more information see Further Details section.
50*> \endverbatim
51*
52* Arguments:
53* ==========
54*
55*> \param[in] UPLO
56*> \verbatim
57*> UPLO is CHARACTER*1
58*> Specifies whether the upper or lower triangular part of the
59*> symmetric matrix A is stored:
60*> = 'U': Upper triangular
61*> = 'L': Lower triangular
62*> \endverbatim
63*>
64*> \param[in] N
65*> \verbatim
66*> N is INTEGER
67*> The order of the matrix A. N >= 0.
68*> \endverbatim
69*>
70*> \param[in,out] A
71*> \verbatim
72*> A is REAL array, dimension (LDA,N)
73*> On entry, the symmetric matrix A.
74*> If UPLO = 'U': the leading N-by-N upper triangular part
75*> of A contains the upper triangular part of the matrix A,
76*> and the strictly lower triangular part of A is not
77*> referenced.
78*>
79*> If UPLO = 'L': the leading N-by-N lower triangular part
80*> of A contains the lower triangular part of the matrix A,
81*> and the strictly upper triangular part of A is not
82*> referenced.
83*>
84*> On exit, contains:
85*> a) ONLY diagonal elements of the symmetric block diagonal
86*> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
87*> (superdiagonal (or subdiagonal) elements of D
88*> are stored on exit in array E), and
89*> b) If UPLO = 'U': factor U in the superdiagonal part of A.
90*> If UPLO = 'L': factor L in the subdiagonal part of A.
91*> \endverbatim
92*>
93*> \param[in] LDA
94*> \verbatim
95*> LDA is INTEGER
96*> The leading dimension of the array A. LDA >= max(1,N).
97*> \endverbatim
98*>
99*> \param[out] E
100*> \verbatim
101*> E is REAL array, dimension (N)
102*> On exit, contains the superdiagonal (or subdiagonal)
103*> elements of the symmetric block diagonal matrix D
104*> with 1-by-1 or 2-by-2 diagonal blocks, where
105*> If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
106*> If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.
107*>
108*> NOTE: For 1-by-1 diagonal block D(k), where
109*> 1 <= k <= N, the element E(k) is set to 0 in both
110*> UPLO = 'U' or UPLO = 'L' cases.
111*> \endverbatim
112*>
113*> \param[out] IPIV
114*> \verbatim
115*> IPIV is INTEGER array, dimension (N)
116*> IPIV describes the permutation matrix P in the factorization
117*> of matrix A as follows. The absolute value of IPIV(k)
118*> represents the index of row and column that were
119*> interchanged with the k-th row and column. The value of UPLO
120*> describes the order in which the interchanges were applied.
121*> Also, the sign of IPIV represents the block structure of
122*> the symmetric block diagonal matrix D with 1-by-1 or 2-by-2
123*> diagonal blocks which correspond to 1 or 2 interchanges
124*> at each factorization step. For more info see Further
125*> Details section.
126*>
127*> If UPLO = 'U',
128*> ( in factorization order, k decreases from N to 1 ):
129*> a) A single positive entry IPIV(k) > 0 means:
130*> D(k,k) is a 1-by-1 diagonal block.
131*> If IPIV(k) != k, rows and columns k and IPIV(k) were
132*> interchanged in the matrix A(1:N,1:N);
133*> If IPIV(k) = k, no interchange occurred.
134*>
135*> b) A pair of consecutive negative entries
136*> IPIV(k) < 0 and IPIV(k-1) < 0 means:
137*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
138*> (NOTE: negative entries in IPIV appear ONLY in pairs).
139*> 1) If -IPIV(k) != k, rows and columns
140*> k and -IPIV(k) were interchanged
141*> in the matrix A(1:N,1:N).
142*> If -IPIV(k) = k, no interchange occurred.
143*> 2) If -IPIV(k-1) != k-1, rows and columns
144*> k-1 and -IPIV(k-1) were interchanged
145*> in the matrix A(1:N,1:N).
146*> If -IPIV(k-1) = k-1, no interchange occurred.
147*>
148*> c) In both cases a) and b), always ABS( IPIV(k) ) <= k.
149*>
150*> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
151*>
152*> If UPLO = 'L',
153*> ( in factorization order, k increases from 1 to N ):
154*> a) A single positive entry IPIV(k) > 0 means:
155*> D(k,k) is a 1-by-1 diagonal block.
156*> If IPIV(k) != k, rows and columns k and IPIV(k) were
157*> interchanged in the matrix A(1:N,1:N).
158*> If IPIV(k) = k, no interchange occurred.
159*>
160*> b) A pair of consecutive negative entries
161*> IPIV(k) < 0 and IPIV(k+1) < 0 means:
162*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
163*> (NOTE: negative entries in IPIV appear ONLY in pairs).
164*> 1) If -IPIV(k) != k, rows and columns
165*> k and -IPIV(k) were interchanged
166*> in the matrix A(1:N,1:N).
167*> If -IPIV(k) = k, no interchange occurred.
168*> 2) If -IPIV(k+1) != k+1, rows and columns
169*> k-1 and -IPIV(k-1) were interchanged
170*> in the matrix A(1:N,1:N).
171*> If -IPIV(k+1) = k+1, no interchange occurred.
172*>
173*> c) In both cases a) and b), always ABS( IPIV(k) ) >= k.
174*>
175*> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
176*> \endverbatim
177*>
178*> \param[out] WORK
179*> \verbatim
180*> WORK is REAL array, dimension ( MAX(1,LWORK) ).
181*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
182*> \endverbatim
183*>
184*> \param[in] LWORK
185*> \verbatim
186*> LWORK is INTEGER
187*> The length of WORK. LWORK >=1. For best performance
188*> LWORK >= N*NB, where NB is the block size returned
189*> by ILAENV.
190*>
191*> If LWORK = -1, then a workspace query is assumed;
192*> the routine only calculates the optimal size of the WORK
193*> array, returns this value as the first entry of the WORK
194*> array, and no error message related to LWORK is issued
195*> by XERBLA.
196*> \endverbatim
197*>
198*> \param[out] INFO
199*> \verbatim
200*> INFO is INTEGER
201*> = 0: successful exit
202*>
203*> < 0: If INFO = -k, the k-th argument had an illegal value
204*>
205*> > 0: If INFO = k, the matrix A is singular, because:
206*> If UPLO = 'U': column k in the upper
207*> triangular part of A contains all zeros.
208*> If UPLO = 'L': column k in the lower
209*> triangular part of A contains all zeros.
210*>
211*> Therefore D(k,k) is exactly zero, and superdiagonal
212*> elements of column k of U (or subdiagonal elements of
213*> column k of L ) are all zeros. The factorization has
214*> been completed, but the block diagonal matrix D is
215*> exactly singular, and division by zero will occur if
216*> it is used to solve a system of equations.
217*>
218*> NOTE: INFO only stores the first occurrence of
219*> a singularity, any subsequent occurrence of singularity
220*> is not stored in INFO even though the factorization
221*> always completes.
222*> \endverbatim
223*
224* Authors:
225* ========
226*
227*> \author Univ. of Tennessee
228*> \author Univ. of California Berkeley
229*> \author Univ. of Colorado Denver
230*> \author NAG Ltd.
231*
232*> \ingroup hetrf_rk
233*
234*> \par Further Details:
235* =====================
236*>
237*> \verbatim
238*> TODO: put correct description
239*> \endverbatim
240*
241*> \par Contributors:
242* ==================
243*>
244*> \verbatim
245*>
246*> December 2016, Igor Kozachenko,
247*> Computer Science Division,
248*> University of California, Berkeley
249*>
250*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
251*> School of Mathematics,
252*> University of Manchester
253*>
254*> \endverbatim
255*
256* =====================================================================
257 SUBROUTINE ssytrf_rk( UPLO, N, A, LDA, E, IPIV, WORK, LWORK,
258 $ INFO )
259*
260* -- LAPACK computational routine --
261* -- LAPACK is a software package provided by Univ. of Tennessee, --
262* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
263*
264* .. Scalar Arguments ..
265 CHARACTER UPLO
266 INTEGER INFO, LDA, LWORK, N
267* ..
268* .. Array Arguments ..
269 INTEGER IPIV( * )
270 REAL A( LDA, * ), E( * ), WORK( * )
271* ..
272*
273* =====================================================================
274*
275* .. Local Scalars ..
276 LOGICAL LQUERY, UPPER
277 INTEGER I, IINFO, IP, IWS, K, KB, LDWORK, LWKOPT,
278 $ nb, nbmin
279* ..
280* .. External Functions ..
281 LOGICAL LSAME
282 INTEGER ILAENV
283 REAL SROUNDUP_LWORK
284 EXTERNAL lsame, ilaenv, sroundup_lwork
285* ..
286* .. External Subroutines ..
287 EXTERNAL slasyf_rk, ssytf2_rk, sswap, xerbla
288* ..
289* .. Intrinsic Functions ..
290 INTRINSIC abs, max
291* ..
292* .. Executable Statements ..
293*
294* Test the input parameters.
295*
296 info = 0
297 upper = lsame( uplo, 'U' )
298 lquery = ( lwork.EQ.-1 )
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( lda.LT.max( 1, n ) ) THEN
304 info = -4
305 ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
306 info = -8
307 END IF
308*
309 IF( info.EQ.0 ) THEN
310*
311* Determine the block size
312*
313 nb = ilaenv( 1, 'SSYTRF_RK', uplo, n, -1, -1, -1 )
314 lwkopt = max( 1, n*nb )
315 work( 1 ) = sroundup_lwork(lwkopt)
316 END IF
317*
318 IF( info.NE.0 ) THEN
319 CALL xerbla( 'SSYTRF_RK', -info )
320 RETURN
321 ELSE IF( lquery ) THEN
322 RETURN
323 END IF
324*
325 nbmin = 2
326 ldwork = n
327 IF( nb.GT.1 .AND. nb.LT.n ) THEN
328 iws = ldwork*nb
329 IF( lwork.LT.iws ) THEN
330 nb = max( lwork / ldwork, 1 )
331 nbmin = max( 2, ilaenv( 2, 'SSYTRF_RK',
332 $ uplo, n, -1, -1, -1 ) )
333 END IF
334 ELSE
335 iws = 1
336 END IF
337 IF( nb.LT.nbmin )
338 $ nb = n
339*
340 IF( upper ) THEN
341*
342* Factorize A as U*D*U**T using the upper triangle of A
343*
344* K is the main loop index, decreasing from N to 1 in steps of
345* KB, where KB is the number of columns factorized by SLASYF_RK;
346* KB is either NB or NB-1, or K for the last block
347*
348 k = n
349 10 CONTINUE
350*
351* If K < 1, exit from loop
352*
353 IF( k.LT.1 )
354 $ GO TO 15
355*
356 IF( k.GT.nb ) THEN
357*
358* Factorize columns k-kb+1:k of A and use blocked code to
359* update columns 1:k-kb
360*
361 CALL slasyf_rk( uplo, k, nb, kb, a, lda, e,
362 $ ipiv, work, ldwork, iinfo )
363 ELSE
364*
365* Use unblocked code to factorize columns 1:k of A
366*
367 CALL ssytf2_rk( uplo, k, a, lda, e, ipiv, iinfo )
368 kb = k
369 END IF
370*
371* Set INFO on the first occurrence of a zero pivot
372*
373 IF( info.EQ.0 .AND. iinfo.GT.0 )
374 $ info = iinfo
375*
376* No need to adjust IPIV
377*
378*
379* Apply permutations to the leading panel 1:k-1
380*
381* Read IPIV from the last block factored, i.e.
382* indices k-kb+1:k and apply row permutations to the
383* last k+1 colunms k+1:N after that block
384* (We can do the simple loop over IPIV with decrement -1,
385* since the ABS value of IPIV( I ) represents the row index
386* of the interchange with row i in both 1x1 and 2x2 pivot cases)
387*
388 IF( k.LT.n ) THEN
389 DO i = k, ( k - kb + 1 ), -1
390 ip = abs( ipiv( i ) )
391 IF( ip.NE.i ) THEN
392 CALL sswap( n-k, a( i, k+1 ), lda,
393 $ a( ip, k+1 ), lda )
394 END IF
395 END DO
396 END IF
397*
398* Decrease K and return to the start of the main loop
399*
400 k = k - kb
401 GO TO 10
402*
403* This label is the exit from main loop over K decreasing
404* from N to 1 in steps of KB
405*
406 15 CONTINUE
407*
408 ELSE
409*
410* Factorize A as L*D*L**T using the lower triangle of A
411*
412* K is the main loop index, increasing from 1 to N in steps of
413* KB, where KB is the number of columns factorized by SLASYF_RK;
414* KB is either NB or NB-1, or N-K+1 for the last block
415*
416 k = 1
417 20 CONTINUE
418*
419* If K > N, exit from loop
420*
421 IF( k.GT.n )
422 $ GO TO 35
423*
424 IF( k.LE.n-nb ) THEN
425*
426* Factorize columns k:k+kb-1 of A and use blocked code to
427* update columns k+kb:n
428*
429 CALL slasyf_rk( uplo, n-k+1, nb, kb, a( k, k ), lda, e( k ),
430 $ ipiv( k ), work, ldwork, iinfo )
431
432
433 ELSE
434*
435* Use unblocked code to factorize columns k:n of A
436*
437 CALL ssytf2_rk( uplo, n-k+1, a( k, k ), lda, e( k ),
438 $ ipiv( k ), iinfo )
439 kb = n - k + 1
440*
441 END IF
442*
443* Set INFO on the first occurrence of a zero pivot
444*
445 IF( info.EQ.0 .AND. iinfo.GT.0 )
446 $ info = iinfo + k - 1
447*
448* Adjust IPIV
449*
450 DO i = k, k + kb - 1
451 IF( ipiv( i ).GT.0 ) THEN
452 ipiv( i ) = ipiv( i ) + k - 1
453 ELSE
454 ipiv( i ) = ipiv( i ) - k + 1
455 END IF
456 END DO
457*
458* Apply permutations to the leading panel 1:k-1
459*
460* Read IPIV from the last block factored, i.e.
461* indices k:k+kb-1 and apply row permutations to the
462* first k-1 colunms 1:k-1 before that block
463* (We can do the simple loop over IPIV with increment 1,
464* since the ABS value of IPIV( I ) represents the row index
465* of the interchange with row i in both 1x1 and 2x2 pivot cases)
466*
467 IF( k.GT.1 ) THEN
468 DO i = k, ( k + kb - 1 ), 1
469 ip = abs( ipiv( i ) )
470 IF( ip.NE.i ) THEN
471 CALL sswap( k-1, a( i, 1 ), lda,
472 $ a( ip, 1 ), lda )
473 END IF
474 END DO
475 END IF
476*
477* Increase K and return to the start of the main loop
478*
479 k = k + kb
480 GO TO 20
481*
482* This label is the exit from main loop over K increasing
483* from 1 to N in steps of KB
484*
485 35 CONTINUE
486*
487* End Lower
488*
489 END IF
490*
491 work( 1 ) = sroundup_lwork(lwkopt)
492 RETURN
493*
494* End of SSYTRF_RK
495*
496 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ssytf2_rk(uplo, n, a, lda, e, ipiv, info)
SSYTF2_RK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Ka...
Definition ssytf2_rk.f:241
subroutine ssytrf_rk(uplo, n, a, lda, e, ipiv, work, lwork, info)
SSYTRF_RK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Ka...
Definition ssytrf_rk.f:259
subroutine slasyf_rk(uplo, n, nb, kb, a, lda, e, ipiv, w, ldw, info)
SLASYF_RK computes a partial factorization of a real symmetric indefinite matrix using bounded Bunch-...
Definition slasyf_rk.f:262
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82