LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
dsytf2.f
Go to the documentation of this file.
1*> \brief \b DSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting method (unblocked algorithm).
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytf2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytf2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytf2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DSYTF2( UPLO, N, A, LDA, IPIV, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER UPLO
25* INTEGER INFO, LDA, N
26* ..
27* .. Array Arguments ..
28* INTEGER IPIV( * )
29* DOUBLE PRECISION A( LDA, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> DSYTF2 computes the factorization of a real symmetric matrix A using
39*> the Bunch-Kaufman diagonal pivoting method:
40*>
41*> A = U*D*U**T or A = L*D*L**T
42*>
43*> where U (or L) is a product of permutation and unit upper (lower)
44*> triangular matrices, U**T is the transpose of U, and D is symmetric and
45*> block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
46*>
47*> This is the unblocked version of the algorithm, calling Level 2 BLAS.
48*> \endverbatim
49*
50* Arguments:
51* ==========
52*
53*> \param[in] UPLO
54*> \verbatim
55*> UPLO is CHARACTER*1
56*> Specifies whether the upper or lower triangular part of the
57*> symmetric matrix A is stored:
58*> = 'U': Upper triangular
59*> = 'L': Lower triangular
60*> \endverbatim
61*>
62*> \param[in] N
63*> \verbatim
64*> N is INTEGER
65*> The order of the matrix A. N >= 0.
66*> \endverbatim
67*>
68*> \param[in,out] A
69*> \verbatim
70*> A is DOUBLE PRECISION array, dimension (LDA,N)
71*> On entry, the symmetric matrix A. If UPLO = 'U', the leading
72*> n-by-n upper triangular part of A contains the upper
73*> triangular part of the matrix A, and the strictly lower
74*> triangular part of A is not referenced. If UPLO = 'L', the
75*> leading n-by-n lower triangular part of A contains the lower
76*> triangular part of the matrix A, and the strictly upper
77*> triangular part of A is not referenced.
78*>
79*> On exit, the block diagonal matrix D and the multipliers used
80*> to obtain the factor U or L (see below for further details).
81*> \endverbatim
82*>
83*> \param[in] LDA
84*> \verbatim
85*> LDA is INTEGER
86*> The leading dimension of the array A. LDA >= max(1,N).
87*> \endverbatim
88*>
89*> \param[out] IPIV
90*> \verbatim
91*> IPIV is INTEGER array, dimension (N)
92*> Details of the interchanges and the block structure of D.
93*>
94*> If UPLO = 'U':
95*> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
96*> interchanged and D(k,k) is a 1-by-1 diagonal block.
97*>
98*> If IPIV(k) = IPIV(k-1) < 0, then rows and columns
99*> k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
100*> is a 2-by-2 diagonal block.
101*>
102*> If UPLO = 'L':
103*> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
104*> interchanged and D(k,k) is a 1-by-1 diagonal block.
105*>
106*> If IPIV(k) = IPIV(k+1) < 0, then rows and columns
107*> k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
108*> is a 2-by-2 diagonal block.
109*> \endverbatim
110*>
111*> \param[out] INFO
112*> \verbatim
113*> INFO is INTEGER
114*> = 0: successful exit
115*> < 0: if INFO = -k, the k-th argument had an illegal value
116*> > 0: if INFO = k, D(k,k) is exactly zero. The factorization
117*> has been completed, but the block diagonal matrix D is
118*> exactly singular, and division by zero will occur if it
119*> is used to solve a system of equations.
120*> \endverbatim
121*
122* Authors:
123* ========
124*
125*> \author Univ. of Tennessee
126*> \author Univ. of California Berkeley
127*> \author Univ. of Colorado Denver
128*> \author NAG Ltd.
129*
130*> \ingroup doubleSYcomputational
131*
132*> \par Further Details:
133* =====================
134*>
135*> \verbatim
136*>
137*> If UPLO = 'U', then A = U*D*U**T, where
138*> U = P(n)*U(n)* ... *P(k)U(k)* ...,
139*> i.e., U is a product of terms P(k)*U(k), where k decreases from n to
140*> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
141*> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
142*> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
143*> that if the diagonal block D(k) is of order s (s = 1 or 2), then
144*>
145*> ( I v 0 ) k-s
146*> U(k) = ( 0 I 0 ) s
147*> ( 0 0 I ) n-k
148*> k-s s n-k
149*>
150*> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
151*> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
152*> and A(k,k), and v overwrites A(1:k-2,k-1:k).
153*>
154*> If UPLO = 'L', then A = L*D*L**T, where
155*> L = P(1)*L(1)* ... *P(k)*L(k)* ...,
156*> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
157*> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
158*> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
159*> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
160*> that if the diagonal block D(k) is of order s (s = 1 or 2), then
161*>
162*> ( I 0 0 ) k-1
163*> L(k) = ( 0 I 0 ) s
164*> ( 0 v I ) n-k-s+1
165*> k-1 s n-k-s+1
166*>
167*> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
168*> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
169*> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
170*> \endverbatim
171*
172*> \par Contributors:
173* ==================
174*>
175*> \verbatim
176*>
177*> 09-29-06 - patch from
178*> Bobby Cheng, MathWorks
179*>
180*> Replace l.204 and l.372
181*> IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
182*> by
183*> IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
184*>
185*> 01-01-96 - Based on modifications by
186*> J. Lewis, Boeing Computer Services Company
187*> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
188*> 1-96 - Based on modifications by J. Lewis, Boeing Computer Services
189*> Company
190*> \endverbatim
191*
192* =====================================================================
193 SUBROUTINE dsytf2( UPLO, N, A, LDA, IPIV, INFO )
194*
195* -- LAPACK computational routine --
196* -- LAPACK is a software package provided by Univ. of Tennessee, --
197* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
198*
199* .. Scalar Arguments ..
200 CHARACTER UPLO
201 INTEGER INFO, LDA, N
202* ..
203* .. Array Arguments ..
204 INTEGER IPIV( * )
205 DOUBLE PRECISION A( LDA, * )
206* ..
207*
208* =====================================================================
209*
210* .. Parameters ..
211 DOUBLE PRECISION ZERO, ONE
212 parameter( zero = 0.0d+0, one = 1.0d+0 )
213 DOUBLE PRECISION EIGHT, SEVTEN
214 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
215* ..
216* .. Local Scalars ..
217 LOGICAL UPPER
218 INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP
219 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1,
220 \$ ROWMAX, T, WK, WKM1, WKP1
221* ..
222* .. External Functions ..
223 LOGICAL LSAME, DISNAN
224 INTEGER IDAMAX
225 EXTERNAL lsame, idamax, disnan
226* ..
227* .. External Subroutines ..
228 EXTERNAL dscal, dswap, dsyr, xerbla
229* ..
230* .. Intrinsic Functions ..
231 INTRINSIC abs, max, sqrt
232* ..
233* .. Executable Statements ..
234*
235* Test the input parameters.
236*
237 info = 0
238 upper = lsame( uplo, 'U' )
239 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
240 info = -1
241 ELSE IF( n.LT.0 ) THEN
242 info = -2
243 ELSE IF( lda.LT.max( 1, n ) ) THEN
244 info = -4
245 END IF
246 IF( info.NE.0 ) THEN
247 CALL xerbla( 'DSYTF2', -info )
248 RETURN
249 END IF
250*
251* Initialize ALPHA for use in choosing pivot block size.
252*
253 alpha = ( one+sqrt( sevten ) ) / eight
254*
255 IF( upper ) THEN
256*
257* Factorize A as U*D*U**T using the upper triangle of A
258*
259* K is the main loop index, decreasing from N to 1 in steps of
260* 1 or 2
261*
262 k = n
263 10 CONTINUE
264*
265* If K < 1, exit from loop
266*
267 IF( k.LT.1 )
268 \$ GO TO 70
269 kstep = 1
270*
271* Determine rows and columns to be interchanged and whether
272* a 1-by-1 or 2-by-2 pivot block will be used
273*
274 absakk = abs( a( k, k ) )
275*
276* IMAX is the row-index of the largest off-diagonal element in
277* column K, and COLMAX is its absolute value.
278* Determine both COLMAX and IMAX.
279*
280 IF( k.GT.1 ) THEN
281 imax = idamax( k-1, a( 1, k ), 1 )
282 colmax = abs( a( imax, k ) )
283 ELSE
284 colmax = zero
285 END IF
286*
287 IF( (max( absakk, colmax ).EQ.zero) .OR. disnan(absakk) ) THEN
288*
289* Column K is zero or underflow, or contains a NaN:
290* set INFO and continue
291*
292 IF( info.EQ.0 )
293 \$ info = k
294 kp = k
295 ELSE
296 IF( absakk.GE.alpha*colmax ) THEN
297*
298* no interchange, use 1-by-1 pivot block
299*
300 kp = k
301 ELSE
302*
303* JMAX is the column-index of the largest off-diagonal
304* element in row IMAX, and ROWMAX is its absolute value
305*
306 jmax = imax + idamax( k-imax, a( imax, imax+1 ), lda )
307 rowmax = abs( a( imax, jmax ) )
308 IF( imax.GT.1 ) THEN
309 jmax = idamax( imax-1, a( 1, imax ), 1 )
310 rowmax = max( rowmax, abs( a( jmax, imax ) ) )
311 END IF
312*
313 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
314*
315* no interchange, use 1-by-1 pivot block
316*
317 kp = k
318 ELSE IF( abs( a( imax, imax ) ).GE.alpha*rowmax ) THEN
319*
320* interchange rows and columns K and IMAX, use 1-by-1
321* pivot block
322*
323 kp = imax
324 ELSE
325*
326* interchange rows and columns K-1 and IMAX, use 2-by-2
327* pivot block
328*
329 kp = imax
330 kstep = 2
331 END IF
332 END IF
333*
334 kk = k - kstep + 1
335 IF( kp.NE.kk ) THEN
336*
337* Interchange rows and columns KK and KP in the leading
338* submatrix A(1:k,1:k)
339*
340 CALL dswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
341 CALL dswap( kk-kp-1, a( kp+1, kk ), 1, a( kp, kp+1 ),
342 \$ lda )
343 t = a( kk, kk )
344 a( kk, kk ) = a( kp, kp )
345 a( kp, kp ) = t
346 IF( kstep.EQ.2 ) THEN
347 t = a( k-1, k )
348 a( k-1, k ) = a( kp, k )
349 a( kp, k ) = t
350 END IF
351 END IF
352*
354*
355 IF( kstep.EQ.1 ) THEN
356*
357* 1-by-1 pivot block D(k): column k now holds
358*
359* W(k) = U(k)*D(k)
360*
361* where U(k) is the k-th column of U
362*
363* Perform a rank-1 update of A(1:k-1,1:k-1) as
364*
365* A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T
366*
367 r1 = one / a( k, k )
368 CALL dsyr( uplo, k-1, -r1, a( 1, k ), 1, a, lda )
369*
370* Store U(k) in column k
371*
372 CALL dscal( k-1, r1, a( 1, k ), 1 )
373 ELSE
374*
375* 2-by-2 pivot block D(k): columns k and k-1 now hold
376*
377* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
378*
379* where U(k) and U(k-1) are the k-th and (k-1)-th columns
380* of U
381*
382* Perform a rank-2 update of A(1:k-2,1:k-2) as
383*
384* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
385* = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T
386*
387 IF( k.GT.2 ) THEN
388*
389 d12 = a( k-1, k )
390 d22 = a( k-1, k-1 ) / d12
391 d11 = a( k, k ) / d12
392 t = one / ( d11*d22-one )
393 d12 = t / d12
394*
395 DO 30 j = k - 2, 1, -1
396 wkm1 = d12*( d11*a( j, k-1 )-a( j, k ) )
397 wk = d12*( d22*a( j, k )-a( j, k-1 ) )
398 DO 20 i = j, 1, -1
399 a( i, j ) = a( i, j ) - a( i, k )*wk -
400 \$ a( i, k-1 )*wkm1
401 20 CONTINUE
402 a( j, k ) = wk
403 a( j, k-1 ) = wkm1
404 30 CONTINUE
405*
406 END IF
407*
408 END IF
409 END IF
410*
411* Store details of the interchanges in IPIV
412*
413 IF( kstep.EQ.1 ) THEN
414 ipiv( k ) = kp
415 ELSE
416 ipiv( k ) = -kp
417 ipiv( k-1 ) = -kp
418 END IF
419*
420* Decrease K and return to the start of the main loop
421*
422 k = k - kstep
423 GO TO 10
424*
425 ELSE
426*
427* Factorize A as L*D*L**T using the lower triangle of A
428*
429* K is the main loop index, increasing from 1 to N in steps of
430* 1 or 2
431*
432 k = 1
433 40 CONTINUE
434*
435* If K > N, exit from loop
436*
437 IF( k.GT.n )
438 \$ GO TO 70
439 kstep = 1
440*
441* Determine rows and columns to be interchanged and whether
442* a 1-by-1 or 2-by-2 pivot block will be used
443*
444 absakk = abs( a( k, k ) )
445*
446* IMAX is the row-index of the largest off-diagonal element in
447* column K, and COLMAX is its absolute value.
448* Determine both COLMAX and IMAX.
449*
450 IF( k.LT.n ) THEN
451 imax = k + idamax( n-k, a( k+1, k ), 1 )
452 colmax = abs( a( imax, k ) )
453 ELSE
454 colmax = zero
455 END IF
456*
457 IF( (max( absakk, colmax ).EQ.zero) .OR. disnan(absakk) ) THEN
458*
459* Column K is zero or underflow, or contains a NaN:
460* set INFO and continue
461*
462 IF( info.EQ.0 )
463 \$ info = k
464 kp = k
465 ELSE
466 IF( absakk.GE.alpha*colmax ) THEN
467*
468* no interchange, use 1-by-1 pivot block
469*
470 kp = k
471 ELSE
472*
473* JMAX is the column-index of the largest off-diagonal
474* element in row IMAX, and ROWMAX is its absolute value
475*
476 jmax = k - 1 + idamax( imax-k, a( imax, k ), lda )
477 rowmax = abs( a( imax, jmax ) )
478 IF( imax.LT.n ) THEN
479 jmax = imax + idamax( n-imax, a( imax+1, imax ), 1 )
480 rowmax = max( rowmax, abs( a( jmax, imax ) ) )
481 END IF
482*
483 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
484*
485* no interchange, use 1-by-1 pivot block
486*
487 kp = k
488 ELSE IF( abs( a( imax, imax ) ).GE.alpha*rowmax ) THEN
489*
490* interchange rows and columns K and IMAX, use 1-by-1
491* pivot block
492*
493 kp = imax
494 ELSE
495*
496* interchange rows and columns K+1 and IMAX, use 2-by-2
497* pivot block
498*
499 kp = imax
500 kstep = 2
501 END IF
502 END IF
503*
504 kk = k + kstep - 1
505 IF( kp.NE.kk ) THEN
506*
507* Interchange rows and columns KK and KP in the trailing
508* submatrix A(k:n,k:n)
509*
510 IF( kp.LT.n )
511 \$ CALL dswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
512 CALL dswap( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
513 \$ lda )
514 t = a( kk, kk )
515 a( kk, kk ) = a( kp, kp )
516 a( kp, kp ) = t
517 IF( kstep.EQ.2 ) THEN
518 t = a( k+1, k )
519 a( k+1, k ) = a( kp, k )
520 a( kp, k ) = t
521 END IF
522 END IF
523*
524* Update the trailing submatrix
525*
526 IF( kstep.EQ.1 ) THEN
527*
528* 1-by-1 pivot block D(k): column k now holds
529*
530* W(k) = L(k)*D(k)
531*
532* where L(k) is the k-th column of L
533*
534 IF( k.LT.n ) THEN
535*
536* Perform a rank-1 update of A(k+1:n,k+1:n) as
537*
538* A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T
539*
540 d11 = one / a( k, k )
541 CALL dsyr( uplo, n-k, -d11, a( k+1, k ), 1,
542 \$ a( k+1, k+1 ), lda )
543*
544* Store L(k) in column K
545*
546 CALL dscal( n-k, d11, a( k+1, k ), 1 )
547 END IF
548 ELSE
549*
550* 2-by-2 pivot block D(k)
551*
552 IF( k.LT.n-1 ) THEN
553*
554* Perform a rank-2 update of A(k+2:n,k+2:n) as
555*
556* A := A - ( (A(k) A(k+1))*D(k)**(-1) ) * (A(k) A(k+1))**T
557*
558* where L(k) and L(k+1) are the k-th and (k+1)-th
559* columns of L
560*
561 d21 = a( k+1, k )
562 d11 = a( k+1, k+1 ) / d21
563 d22 = a( k, k ) / d21
564 t = one / ( d11*d22-one )
565 d21 = t / d21
566*
567 DO 60 j = k + 2, n
568*
569 wk = d21*( d11*a( j, k )-a( j, k+1 ) )
570 wkp1 = d21*( d22*a( j, k+1 )-a( j, k ) )
571*
572 DO 50 i = j, n
573 a( i, j ) = a( i, j ) - a( i, k )*wk -
574 \$ a( i, k+1 )*wkp1
575 50 CONTINUE
576*
577 a( j, k ) = wk
578 a( j, k+1 ) = wkp1
579*
580 60 CONTINUE
581 END IF
582 END IF
583 END IF
584*
585* Store details of the interchanges in IPIV
586*
587 IF( kstep.EQ.1 ) THEN
588 ipiv( k ) = kp
589 ELSE
590 ipiv( k ) = -kp
591 ipiv( k+1 ) = -kp
592 END IF
593*
594* Increase K and return to the start of the main loop
595*
596 k = k + kstep
597 GO TO 40
598*
599 END IF
600*
601 70 CONTINUE
602*
603 RETURN
604*
605* End of DSYTF2
606*
607 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
subroutine dsyr(UPLO, N, ALPHA, X, INCX, A, LDA)
DSYR
Definition: dsyr.f:132
subroutine dsytf2(UPLO, N, A, LDA, IPIV, INFO)
DSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting ...
Definition: dsytf2.f:194