LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dsytrf_rk.f
Go to the documentation of this file.
1*> \brief \b DSYTRF_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 DSYTRF_RK + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytrf_rk.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytrf_rk.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytrf_rk.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DSYTRF_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* DOUBLE PRECISION A( LDA, * ), E ( * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*> DSYTRF_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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 dsytrf_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 DOUBLE PRECISION 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 EXTERNAL lsame, ilaenv
284* ..
285* .. External Subroutines ..
286 EXTERNAL dlasyf_rk, dsytf2_rk, dswap, xerbla
287* ..
288* .. Intrinsic Functions ..
289 INTRINSIC abs, max
290* ..
291* .. Executable Statements ..
292*
293* Test the input parameters.
294*
295 info = 0
296 upper = lsame( uplo, 'U' )
297 lquery = ( lwork.EQ.-1 )
298 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
299 info = -1
300 ELSE IF( n.LT.0 ) THEN
301 info = -2
302 ELSE IF( lda.LT.max( 1, n ) ) THEN
303 info = -4
304 ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
305 info = -8
306 END IF
307*
308 IF( info.EQ.0 ) THEN
309*
310* Determine the block size
311*
312 nb = ilaenv( 1, 'DSYTRF_RK', uplo, n, -1, -1, -1 )
313 lwkopt = max( 1, n*nb )
314 work( 1 ) = lwkopt
315 END IF
316*
317 IF( info.NE.0 ) THEN
318 CALL xerbla( 'DSYTRF_RK', -info )
319 RETURN
320 ELSE IF( lquery ) THEN
321 RETURN
322 END IF
323*
324 nbmin = 2
325 ldwork = n
326 IF( nb.GT.1 .AND. nb.LT.n ) THEN
327 iws = ldwork*nb
328 IF( lwork.LT.iws ) THEN
329 nb = max( lwork / ldwork, 1 )
330 nbmin = max( 2, ilaenv( 2, 'DSYTRF_RK',
331 $ uplo, n, -1, -1, -1 ) )
332 END IF
333 ELSE
334 iws = 1
335 END IF
336 IF( nb.LT.nbmin )
337 $ nb = n
338*
339 IF( upper ) THEN
340*
341* Factorize A as U*D*U**T using the upper triangle of A
342*
343* K is the main loop index, decreasing from N to 1 in steps of
344* KB, where KB is the number of columns factorized by DLASYF_RK;
345* KB is either NB or NB-1, or K for the last block
346*
347 k = n
348 10 CONTINUE
349*
350* If K < 1, exit from loop
351*
352 IF( k.LT.1 )
353 $ GO TO 15
354*
355 IF( k.GT.nb ) THEN
356*
357* Factorize columns k-kb+1:k of A and use blocked code to
358* update columns 1:k-kb
359*
360 CALL dlasyf_rk( uplo, k, nb, kb, a, lda, e,
361 $ ipiv, work, ldwork, iinfo )
362 ELSE
363*
364* Use unblocked code to factorize columns 1:k of A
365*
366 CALL dsytf2_rk( uplo, k, a, lda, e, ipiv, iinfo )
367 kb = k
368 END IF
369*
370* Set INFO on the first occurrence of a zero pivot
371*
372 IF( info.EQ.0 .AND. iinfo.GT.0 )
373 $ info = iinfo
374*
375* No need to adjust IPIV
376*
377*
378* Apply permutations to the leading panel 1:k-1
379*
380* Read IPIV from the last block factored, i.e.
381* indices k-kb+1:k and apply row permutations to the
382* last k+1 colunms k+1:N after that block
383* (We can do the simple loop over IPIV with decrement -1,
384* since the ABS value of IPIV( I ) represents the row index
385* of the interchange with row i in both 1x1 and 2x2 pivot cases)
386*
387 IF( k.LT.n ) THEN
388 DO i = k, ( k - kb + 1 ), -1
389 ip = abs( ipiv( i ) )
390 IF( ip.NE.i ) THEN
391 CALL dswap( n-k, a( i, k+1 ), lda,
392 $ a( ip, k+1 ), lda )
393 END IF
394 END DO
395 END IF
396*
397* Decrease K and return to the start of the main loop
398*
399 k = k - kb
400 GO TO 10
401*
402* This label is the exit from main loop over K decreasing
403* from N to 1 in steps of KB
404*
405 15 CONTINUE
406*
407 ELSE
408*
409* Factorize A as L*D*L**T using the lower triangle of A
410*
411* K is the main loop index, increasing from 1 to N in steps of
412* KB, where KB is the number of columns factorized by DLASYF_RK;
413* KB is either NB or NB-1, or N-K+1 for the last block
414*
415 k = 1
416 20 CONTINUE
417*
418* If K > N, exit from loop
419*
420 IF( k.GT.n )
421 $ GO TO 35
422*
423 IF( k.LE.n-nb ) THEN
424*
425* Factorize columns k:k+kb-1 of A and use blocked code to
426* update columns k+kb:n
427*
428 CALL dlasyf_rk( uplo, n-k+1, nb, kb, a( k, k ), lda, e( k ),
429 $ ipiv( k ), work, ldwork, iinfo )
430
431
432 ELSE
433*
434* Use unblocked code to factorize columns k:n of A
435*
436 CALL dsytf2_rk( uplo, n-k+1, a( k, k ), lda, e( k ),
437 $ ipiv( k ), iinfo )
438 kb = n - k + 1
439*
440 END IF
441*
442* Set INFO on the first occurrence of a zero pivot
443*
444 IF( info.EQ.0 .AND. iinfo.GT.0 )
445 $ info = iinfo + k - 1
446*
447* Adjust IPIV
448*
449 DO i = k, k + kb - 1
450 IF( ipiv( i ).GT.0 ) THEN
451 ipiv( i ) = ipiv( i ) + k - 1
452 ELSE
453 ipiv( i ) = ipiv( i ) - k + 1
454 END IF
455 END DO
456*
457* Apply permutations to the leading panel 1:k-1
458*
459* Read IPIV from the last block factored, i.e.
460* indices k:k+kb-1 and apply row permutations to the
461* first k-1 colunms 1:k-1 before that block
462* (We can do the simple loop over IPIV with increment 1,
463* since the ABS value of IPIV( I ) represents the row index
464* of the interchange with row i in both 1x1 and 2x2 pivot cases)
465*
466 IF( k.GT.1 ) THEN
467 DO i = k, ( k + kb - 1 ), 1
468 ip = abs( ipiv( i ) )
469 IF( ip.NE.i ) THEN
470 CALL dswap( k-1, a( i, 1 ), lda,
471 $ a( ip, 1 ), lda )
472 END IF
473 END DO
474 END IF
475*
476* Increase K and return to the start of the main loop
477*
478 k = k + kb
479 GO TO 20
480*
481* This label is the exit from main loop over K increasing
482* from 1 to N in steps of KB
483*
484 35 CONTINUE
485*
486* End Lower
487*
488 END IF
489*
490 work( 1 ) = lwkopt
491 RETURN
492*
493* End of DSYTRF_RK
494*
495 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dsytf2_rk(uplo, n, a, lda, e, ipiv, info)
DSYTF2_RK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Ka...
Definition dsytf2_rk.f:241
subroutine dsytrf_rk(uplo, n, a, lda, e, ipiv, work, lwork, info)
DSYTRF_RK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Ka...
Definition dsytrf_rk.f:259
subroutine dlasyf_rk(uplo, n, nb, kb, a, lda, e, ipiv, w, ldw, info)
DLASYF_RK computes a partial factorization of a real symmetric indefinite matrix using bounded Bunch-...
Definition dlasyf_rk.f:262
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82