LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
dsytf2_rk.f
Go to the documentation of this file.
1*> \brief \b DSYTF2_RK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Kaufman (rook) diagonal pivoting method (BLAS2 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_rk.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytf2_rk.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytf2_rk.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DSYTF2_RK( UPLO, N, A, LDA, E, 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, * ), E ( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*> DSYTF2_RK computes the factorization of a real symmetric matrix A
38*> using the bounded Bunch-Kaufman (rook) diagonal pivoting method:
39*>
40*> A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T),
41*>
42*> where U (or L) is unit upper (or lower) triangular matrix,
43*> U**T (or L**T) is the transpose of U (or L), P is a permutation
44*> matrix, P**T is the transpose of P, and D is symmetric and block
45*> 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.
49*> \endverbatim
50*
51* Arguments:
52* ==========
53*
54*> \param[in] UPLO
55*> \verbatim
56*> UPLO is CHARACTER*1
57*> Specifies whether the upper or lower triangular part of the
58*> symmetric matrix A is stored:
59*> = 'U': Upper triangular
60*> = 'L': Lower triangular
61*> \endverbatim
62*>
63*> \param[in] N
64*> \verbatim
65*> N is INTEGER
66*> The order of the matrix A. N >= 0.
67*> \endverbatim
68*>
69*> \param[in,out] A
70*> \verbatim
71*> A is DOUBLE PRECISION array, dimension (LDA,N)
72*> On entry, the symmetric matrix A.
73*> If UPLO = 'U': the leading N-by-N upper triangular part
74*> of A contains the upper triangular part of the matrix A,
75*> and the strictly lower triangular part of A is not
76*> referenced.
77*>
78*> If UPLO = 'L': the leading N-by-N lower triangular part
79*> of A contains the lower triangular part of the matrix A,
80*> and the strictly upper triangular part of A is not
81*> referenced.
82*>
83*> On exit, contains:
84*> a) ONLY diagonal elements of the symmetric block diagonal
85*> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
86*> (superdiagonal (or subdiagonal) elements of D
87*> are stored on exit in array E), and
88*> b) If UPLO = 'U': factor U in the superdiagonal part of A.
89*> If UPLO = 'L': factor L in the subdiagonal part of A.
90*> \endverbatim
91*>
92*> \param[in] LDA
93*> \verbatim
94*> LDA is INTEGER
95*> The leading dimension of the array A. LDA >= max(1,N).
96*> \endverbatim
97*>
98*> \param[out] E
99*> \verbatim
100*> E is DOUBLE PRECISION array, dimension (N)
101*> On exit, contains the superdiagonal (or subdiagonal)
102*> elements of the symmetric block diagonal matrix D
103*> with 1-by-1 or 2-by-2 diagonal blocks, where
104*> If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
105*> If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.
106*>
107*> NOTE: For 1-by-1 diagonal block D(k), where
108*> 1 <= k <= N, the element E(k) is set to 0 in both
109*> UPLO = 'U' or UPLO = 'L' cases.
110*> \endverbatim
111*>
112*> \param[out] IPIV
113*> \verbatim
114*> IPIV is INTEGER array, dimension (N)
115*> IPIV describes the permutation matrix P in the factorization
116*> of matrix A as follows. The absolute value of IPIV(k)
117*> represents the index of row and column that were
118*> interchanged with the k-th row and column. The value of UPLO
119*> describes the order in which the interchanges were applied.
120*> Also, the sign of IPIV represents the block structure of
121*> the symmetric block diagonal matrix D with 1-by-1 or 2-by-2
122*> diagonal blocks which correspond to 1 or 2 interchanges
124*> Details section.
125*>
126*> If UPLO = 'U',
127*> ( in factorization order, k decreases from N to 1 ):
128*> a) A single positive entry IPIV(k) > 0 means:
129*> D(k,k) is a 1-by-1 diagonal block.
130*> If IPIV(k) != k, rows and columns k and IPIV(k) were
131*> interchanged in the matrix A(1:N,1:N);
132*> If IPIV(k) = k, no interchange occurred.
133*>
134*> b) A pair of consecutive negative entries
135*> IPIV(k) < 0 and IPIV(k-1) < 0 means:
136*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
137*> (NOTE: negative entries in IPIV appear ONLY in pairs).
138*> 1) If -IPIV(k) != k, rows and columns
139*> k and -IPIV(k) were interchanged
140*> in the matrix A(1:N,1:N).
141*> If -IPIV(k) = k, no interchange occurred.
142*> 2) If -IPIV(k-1) != k-1, rows and columns
143*> k-1 and -IPIV(k-1) were interchanged
144*> in the matrix A(1:N,1:N).
145*> If -IPIV(k-1) = k-1, no interchange occurred.
146*>
147*> c) In both cases a) and b), always ABS( IPIV(k) ) <= k.
148*>
149*> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
150*>
151*> If UPLO = 'L',
152*> ( in factorization order, k increases from 1 to N ):
153*> a) A single positive entry IPIV(k) > 0 means:
154*> D(k,k) is a 1-by-1 diagonal block.
155*> If IPIV(k) != k, rows and columns k and IPIV(k) were
156*> interchanged in the matrix A(1:N,1:N).
157*> If IPIV(k) = k, no interchange occurred.
158*>
159*> b) A pair of consecutive negative entries
160*> IPIV(k) < 0 and IPIV(k+1) < 0 means:
161*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
162*> (NOTE: negative entries in IPIV appear ONLY in pairs).
163*> 1) If -IPIV(k) != k, rows and columns
164*> k and -IPIV(k) were interchanged
165*> in the matrix A(1:N,1:N).
166*> If -IPIV(k) = k, no interchange occurred.
167*> 2) If -IPIV(k+1) != k+1, rows and columns
168*> k-1 and -IPIV(k-1) were interchanged
169*> in the matrix A(1:N,1:N).
170*> If -IPIV(k+1) = k+1, no interchange occurred.
171*>
172*> c) In both cases a) and b), always ABS( IPIV(k) ) >= k.
173*>
174*> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
175*> \endverbatim
176*>
177*> \param[out] INFO
178*> \verbatim
179*> INFO is INTEGER
180*> = 0: successful exit
181*>
182*> < 0: If INFO = -k, the k-th argument had an illegal value
183*>
184*> > 0: If INFO = k, the matrix A is singular, because:
185*> If UPLO = 'U': column k in the upper
186*> triangular part of A contains all zeros.
187*> If UPLO = 'L': column k in the lower
188*> triangular part of A contains all zeros.
189*>
190*> Therefore D(k,k) is exactly zero, and superdiagonal
191*> elements of column k of U (or subdiagonal elements of
192*> column k of L ) are all zeros. The factorization has
193*> been completed, but the block diagonal matrix D is
194*> exactly singular, and division by zero will occur if
195*> it is used to solve a system of equations.
196*>
197*> NOTE: INFO only stores the first occurrence of
198*> a singularity, any subsequent occurrence of singularity
199*> is not stored in INFO even though the factorization
200*> always completes.
201*> \endverbatim
202*
203* Authors:
204* ========
205*
206*> \author Univ. of Tennessee
207*> \author Univ. of California Berkeley
208*> \author Univ. of Colorado Denver
209*> \author NAG Ltd.
210*
211*> \ingroup hetf2_rk
212*
213*> \par Further Details:
214* =====================
215*>
216*> \verbatim
217*> TODO: put further details
218*> \endverbatim
219*
220*> \par Contributors:
221* ==================
222*>
223*> \verbatim
224*>
225*> December 2016, Igor Kozachenko,
226*> Computer Science Division,
227*> University of California, Berkeley
228*>
229*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
230*> School of Mathematics,
231*> University of Manchester
232*>
233*> 01-01-96 - Based on modifications by
234*> J. Lewis, Boeing Computer Services Company
235*> A. Petitet, Computer Science Dept.,
236*> Univ. of Tenn., Knoxville abd , USA
237*> \endverbatim
238*
239* =====================================================================
240 SUBROUTINE dsytf2_rk( UPLO, N, A, LDA, E, IPIV, INFO )
241*
242* -- LAPACK computational routine --
243* -- LAPACK is a software package provided by Univ. of Tennessee, --
244* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
245*
246* .. Scalar Arguments ..
247 CHARACTER UPLO
248 INTEGER INFO, LDA, N
249* ..
250* .. Array Arguments ..
251 INTEGER IPIV( * )
252 DOUBLE PRECISION A( LDA, * ), E( * )
253* ..
254*
255* =====================================================================
256*
257* .. Parameters ..
258 DOUBLE PRECISION ZERO, ONE
259 parameter( zero = 0.0d+0, one = 1.0d+0 )
260 DOUBLE PRECISION EIGHT, SEVTEN
261 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
262* ..
263* .. Local Scalars ..
264 LOGICAL UPPER, DONE
265 INTEGER I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP,
266 \$ P, II
267 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
268 \$ ROWMAX, DTEMP, T, WK, WKM1, WKP1, SFMIN
269* ..
270* .. External Functions ..
271 LOGICAL LSAME
272 INTEGER IDAMAX
273 DOUBLE PRECISION DLAMCH
274 EXTERNAL lsame, idamax, dlamch
275* ..
276* .. External Subroutines ..
277 EXTERNAL dscal, dswap, dsyr, xerbla
278* ..
279* .. Intrinsic Functions ..
280 INTRINSIC abs, max, sqrt
281* ..
282* .. Executable Statements ..
283*
284* Test the input parameters.
285*
286 info = 0
287 upper = lsame( uplo, 'U' )
288 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
289 info = -1
290 ELSE IF( n.LT.0 ) THEN
291 info = -2
292 ELSE IF( lda.LT.max( 1, n ) ) THEN
293 info = -4
294 END IF
295 IF( info.NE.0 ) THEN
296 CALL xerbla( 'DSYTF2_RK', -info )
297 RETURN
298 END IF
299*
300* Initialize ALPHA for use in choosing pivot block size.
301*
302 alpha = ( one+sqrt( sevten ) ) / eight
303*
304* Compute machine safe minimum
305*
306 sfmin = dlamch( 'S' )
307*
308 IF( upper ) THEN
309*
310* Factorize A as U*D*U**T using the upper triangle of A
311*
312* Initialize the first entry of array E, where superdiagonal
313* elements of D are stored
314*
315 e( 1 ) = zero
316*
317* K is the main loop index, decreasing from N to 1 in steps of
318* 1 or 2
319*
320 k = n
321 10 CONTINUE
322*
323* If K < 1, exit from loop
324*
325 IF( k.LT.1 )
326 \$ GO TO 34
327 kstep = 1
328 p = k
329*
330* Determine rows and columns to be interchanged and whether
331* a 1-by-1 or 2-by-2 pivot block will be used
332*
333 absakk = abs( a( k, k ) )
334*
335* IMAX is the row-index of the largest off-diagonal element in
336* column K, and COLMAX is its absolute value.
337* Determine both COLMAX and IMAX.
338*
339 IF( k.GT.1 ) THEN
340 imax = idamax( k-1, a( 1, k ), 1 )
341 colmax = abs( a( imax, k ) )
342 ELSE
343 colmax = zero
344 END IF
345*
346 IF( (max( absakk, colmax ).EQ.zero) ) THEN
347*
348* Column K is zero or underflow: set INFO and continue
349*
350 IF( info.EQ.0 )
351 \$ info = k
352 kp = k
353*
354* Set E( K ) to zero
355*
356 IF( k.GT.1 )
357 \$ e( k ) = zero
358*
359 ELSE
360*
361* Test for interchange
362*
363* Equivalent to testing for (used to handle NaN and Inf)
364* ABSAKK.GE.ALPHA*COLMAX
365*
366 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
367*
368* no interchange,
369* use 1-by-1 pivot block
370*
371 kp = k
372 ELSE
373*
374 done = .false.
375*
376* Loop until pivot found
377*
378 12 CONTINUE
379*
380* Begin pivot search loop body
381*
382* JMAX is the column-index of the largest off-diagonal
383* element in row IMAX, and ROWMAX is its absolute value.
384* Determine both ROWMAX and JMAX.
385*
386 IF( imax.NE.k ) THEN
387 jmax = imax + idamax( k-imax, a( imax, imax+1 ),
388 \$ lda )
389 rowmax = abs( a( imax, jmax ) )
390 ELSE
391 rowmax = zero
392 END IF
393*
394 IF( imax.GT.1 ) THEN
395 itemp = idamax( imax-1, a( 1, imax ), 1 )
396 dtemp = abs( a( itemp, imax ) )
397 IF( dtemp.GT.rowmax ) THEN
398 rowmax = dtemp
399 jmax = itemp
400 END IF
401 END IF
402*
403* Equivalent to testing for (used to handle NaN and Inf)
404* ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
405*
406 IF( .NOT.( abs( a( imax, imax ) ).LT.alpha*rowmax ) )
407 \$ THEN
408*
409* interchange rows and columns K and IMAX,
410* use 1-by-1 pivot block
411*
412 kp = imax
413 done = .true.
414*
415* Equivalent to testing for ROWMAX .EQ. COLMAX,
416* used to handle NaN and Inf
417*
418 ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
419*
420* interchange rows and columns K+1 and IMAX,
421* use 2-by-2 pivot block
422*
423 kp = imax
424 kstep = 2
425 done = .true.
426 ELSE
427*
429*
430 p = imax
431 colmax = rowmax
432 imax = jmax
433 END IF
434*
435* End pivot search loop body
436*
437 IF( .NOT. done ) GOTO 12
438*
439 END IF
440*
441* Swap TWO rows and TWO columns
442*
443* First swap
444*
445 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
446*
447* Interchange rows and column K and P in the leading
448* submatrix A(1:k,1:k) if we have a 2-by-2 pivot
449*
450 IF( p.GT.1 )
451 \$ CALL dswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
452 IF( p.LT.(k-1) )
453 \$ CALL dswap( k-p-1, a( p+1, k ), 1, a( p, p+1 ),
454 \$ lda )
455 t = a( k, k )
456 a( k, k ) = a( p, p )
457 a( p, p ) = t
458*
459* Convert upper triangle of A into U form by applying
460* the interchanges in columns k+1:N.
461*
462 IF( k.LT.n )
463 \$ CALL dswap( n-k, a( k, k+1 ), lda, a( p, k+1 ), lda )
464*
465 END IF
466*
467* Second swap
468*
469 kk = k - kstep + 1
470 IF( kp.NE.kk ) THEN
471*
472* Interchange rows and columns KK and KP in the leading
473* submatrix A(1:k,1:k)
474*
475 IF( kp.GT.1 )
476 \$ CALL dswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
477 IF( ( kk.GT.1 ) .AND. ( kp.LT.(kk-1) ) )
478 \$ CALL dswap( kk-kp-1, a( kp+1, kk ), 1, a( kp, kp+1 ),
479 \$ lda )
480 t = a( kk, kk )
481 a( kk, kk ) = a( kp, kp )
482 a( kp, kp ) = t
483 IF( kstep.EQ.2 ) THEN
484 t = a( k-1, k )
485 a( k-1, k ) = a( kp, k )
486 a( kp, k ) = t
487 END IF
488*
489* Convert upper triangle of A into U form by applying
490* the interchanges in columns k+1:N.
491*
492 IF( k.LT.n )
493 \$ CALL dswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
494 \$ lda )
495*
496 END IF
497*
499*
500 IF( kstep.EQ.1 ) THEN
501*
502* 1-by-1 pivot block D(k): column k now holds
503*
504* W(k) = U(k)*D(k)
505*
506* where U(k) is the k-th column of U
507*
508 IF( k.GT.1 ) THEN
509*
510* Perform a rank-1 update of A(1:k-1,1:k-1) and
511* store U(k) in column k
512*
513 IF( abs( a( k, k ) ).GE.sfmin ) THEN
514*
515* Perform a rank-1 update of A(1:k-1,1:k-1) as
516* A := A - U(k)*D(k)*U(k)**T
517* = A - W(k)*1/D(k)*W(k)**T
518*
519 d11 = one / a( k, k )
520 CALL dsyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
521*
522* Store U(k) in column k
523*
524 CALL dscal( k-1, d11, a( 1, k ), 1 )
525 ELSE
526*
527* Store L(k) in column K
528*
529 d11 = a( k, k )
530 DO 16 ii = 1, k - 1
531 a( ii, k ) = a( ii, k ) / d11
532 16 CONTINUE
533*
534* Perform a rank-1 update of A(k+1:n,k+1:n) as
535* A := A - U(k)*D(k)*U(k)**T
536* = A - W(k)*(1/D(k))*W(k)**T
537* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
538*
539 CALL dsyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
540 END IF
541*
542* Store the superdiagonal element of D in array E
543*
544 e( k ) = zero
545*
546 END IF
547*
548 ELSE
549*
550* 2-by-2 pivot block D(k): columns k and k-1 now hold
551*
552* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
553*
554* where U(k) and U(k-1) are the k-th and (k-1)-th columns
555* of U
556*
557* Perform a rank-2 update of A(1:k-2,1:k-2) as
558*
559* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
560* = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
561*
562* and store L(k) and L(k+1) in columns k and k+1
563*
564 IF( k.GT.2 ) THEN
565*
566 d12 = a( k-1, k )
567 d22 = a( k-1, k-1 ) / d12
568 d11 = a( k, k ) / d12
569 t = one / ( d11*d22-one )
570*
571 DO 30 j = k - 2, 1, -1
572*
573 wkm1 = t*( d11*a( j, k-1 )-a( j, k ) )
574 wk = t*( d22*a( j, k )-a( j, k-1 ) )
575*
576 DO 20 i = j, 1, -1
577 a( i, j ) = a( i, j ) - (a( i, k ) / d12 )*wk -
578 \$ ( a( i, k-1 ) / d12 )*wkm1
579 20 CONTINUE
580*
581* Store U(k) and U(k-1) in cols k and k-1 for row J
582*
583 a( j, k ) = wk / d12
584 a( j, k-1 ) = wkm1 / d12
585*
586 30 CONTINUE
587*
588 END IF
589*
590* Copy superdiagonal elements of D(K) to E(K) and
591* ZERO out superdiagonal entry of A
592*
593 e( k ) = a( k-1, k )
594 e( k-1 ) = zero
595 a( k-1, k ) = zero
596*
597 END IF
598*
599* End column K is nonsingular
600*
601 END IF
602*
603* Store details of the interchanges in IPIV
604*
605 IF( kstep.EQ.1 ) THEN
606 ipiv( k ) = kp
607 ELSE
608 ipiv( k ) = -p
609 ipiv( k-1 ) = -kp
610 END IF
611*
612* Decrease K and return to the start of the main loop
613*
614 k = k - kstep
615 GO TO 10
616*
617 34 CONTINUE
618*
619 ELSE
620*
621* Factorize A as L*D*L**T using the lower triangle of A
622*
623* Initialize the unused last entry of the subdiagonal array E.
624*
625 e( n ) = zero
626*
627* K is the main loop index, increasing from 1 to N in steps of
628* 1 or 2
629*
630 k = 1
631 40 CONTINUE
632*
633* If K > N, exit from loop
634*
635 IF( k.GT.n )
636 \$ GO TO 64
637 kstep = 1
638 p = k
639*
640* Determine rows and columns to be interchanged and whether
641* a 1-by-1 or 2-by-2 pivot block will be used
642*
643 absakk = abs( a( k, k ) )
644*
645* IMAX is the row-index of the largest off-diagonal element in
646* column K, and COLMAX is its absolute value.
647* Determine both COLMAX and IMAX.
648*
649 IF( k.LT.n ) THEN
650 imax = k + idamax( n-k, a( k+1, k ), 1 )
651 colmax = abs( a( imax, k ) )
652 ELSE
653 colmax = zero
654 END IF
655*
656 IF( ( max( absakk, colmax ).EQ.zero ) ) THEN
657*
658* Column K is zero or underflow: set INFO and continue
659*
660 IF( info.EQ.0 )
661 \$ info = k
662 kp = k
663*
664* Set E( K ) to zero
665*
666 IF( k.LT.n )
667 \$ e( k ) = zero
668*
669 ELSE
670*
671* Test for interchange
672*
673* Equivalent to testing for (used to handle NaN and Inf)
674* ABSAKK.GE.ALPHA*COLMAX
675*
676 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
677*
678* no interchange, use 1-by-1 pivot block
679*
680 kp = k
681*
682 ELSE
683*
684 done = .false.
685*
686* Loop until pivot found
687*
688 42 CONTINUE
689*
690* Begin pivot search loop body
691*
692* JMAX is the column-index of the largest off-diagonal
693* element in row IMAX, and ROWMAX is its absolute value.
694* Determine both ROWMAX and JMAX.
695*
696 IF( imax.NE.k ) THEN
697 jmax = k - 1 + idamax( imax-k, a( imax, k ), lda )
698 rowmax = abs( a( imax, jmax ) )
699 ELSE
700 rowmax = zero
701 END IF
702*
703 IF( imax.LT.n ) THEN
704 itemp = imax + idamax( n-imax, a( imax+1, imax ),
705 \$ 1 )
706 dtemp = abs( a( itemp, imax ) )
707 IF( dtemp.GT.rowmax ) THEN
708 rowmax = dtemp
709 jmax = itemp
710 END IF
711 END IF
712*
713* Equivalent to testing for (used to handle NaN and Inf)
714* ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
715*
716 IF( .NOT.( abs( a( imax, imax ) ).LT.alpha*rowmax ) )
717 \$ THEN
718*
719* interchange rows and columns K and IMAX,
720* use 1-by-1 pivot block
721*
722 kp = imax
723 done = .true.
724*
725* Equivalent to testing for ROWMAX .EQ. COLMAX,
726* used to handle NaN and Inf
727*
728 ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
729*
730* interchange rows and columns K+1 and IMAX,
731* use 2-by-2 pivot block
732*
733 kp = imax
734 kstep = 2
735 done = .true.
736 ELSE
737*
739*
740 p = imax
741 colmax = rowmax
742 imax = jmax
743 END IF
744*
745* End pivot search loop body
746*
747 IF( .NOT. done ) GOTO 42
748*
749 END IF
750*
751* Swap TWO rows and TWO columns
752*
753* First swap
754*
755 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
756*
757* Interchange rows and column K and P in the trailing
758* submatrix A(k:n,k:n) if we have a 2-by-2 pivot
759*
760 IF( p.LT.n )
761 \$ CALL dswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
762 IF( p.GT.(k+1) )
763 \$ CALL dswap( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
764 t = a( k, k )
765 a( k, k ) = a( p, p )
766 a( p, p ) = t
767*
768* Convert lower triangle of A into L form by applying
769* the interchanges in columns 1:k-1.
770*
771 IF ( k.GT.1 )
772 \$ CALL dswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
773*
774 END IF
775*
776* Second swap
777*
778 kk = k + kstep - 1
779 IF( kp.NE.kk ) THEN
780*
781* Interchange rows and columns KK and KP in the trailing
782* submatrix A(k:n,k:n)
783*
784 IF( kp.LT.n )
785 \$ CALL dswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
786 IF( ( kk.LT.n ) .AND. ( kp.GT.(kk+1) ) )
787 \$ CALL dswap( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
788 \$ lda )
789 t = a( kk, kk )
790 a( kk, kk ) = a( kp, kp )
791 a( kp, kp ) = t
792 IF( kstep.EQ.2 ) THEN
793 t = a( k+1, k )
794 a( k+1, k ) = a( kp, k )
795 a( kp, k ) = t
796 END IF
797*
798* Convert lower triangle of A into L form by applying
799* the interchanges in columns 1:k-1.
800*
801 IF ( k.GT.1 )
802 \$ CALL dswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
803*
804 END IF
805*
806* Update the trailing submatrix
807*
808 IF( kstep.EQ.1 ) THEN
809*
810* 1-by-1 pivot block D(k): column k now holds
811*
812* W(k) = L(k)*D(k)
813*
814* where L(k) is the k-th column of L
815*
816 IF( k.LT.n ) THEN
817*
818* Perform a rank-1 update of A(k+1:n,k+1:n) and
819* store L(k) in column k
820*
821 IF( abs( a( k, k ) ).GE.sfmin ) THEN
822*
823* Perform a rank-1 update of A(k+1:n,k+1:n) as
824* A := A - L(k)*D(k)*L(k)**T
825* = A - W(k)*(1/D(k))*W(k)**T
826*
827 d11 = one / a( k, k )
828 CALL dsyr( uplo, n-k, -d11, a( k+1, k ), 1,
829 \$ a( k+1, k+1 ), lda )
830*
831* Store L(k) in column k
832*
833 CALL dscal( n-k, d11, a( k+1, k ), 1 )
834 ELSE
835*
836* Store L(k) in column k
837*
838 d11 = a( k, k )
839 DO 46 ii = k + 1, n
840 a( ii, k ) = a( ii, k ) / d11
841 46 CONTINUE
842*
843* Perform a rank-1 update of A(k+1:n,k+1:n) as
844* A := A - L(k)*D(k)*L(k)**T
845* = A - W(k)*(1/D(k))*W(k)**T
846* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
847*
848 CALL dsyr( uplo, n-k, -d11, a( k+1, k ), 1,
849 \$ a( k+1, k+1 ), lda )
850 END IF
851*
852* Store the subdiagonal element of D in array E
853*
854 e( k ) = zero
855*
856 END IF
857*
858 ELSE
859*
860* 2-by-2 pivot block D(k): columns k and k+1 now hold
861*
862* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
863*
864* where L(k) and L(k+1) are the k-th and (k+1)-th columns
865* of L
866*
867*
868* Perform a rank-2 update of A(k+2:n,k+2:n) as
869*
870* A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
871* = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
872*
873* and store L(k) and L(k+1) in columns k and k+1
874*
875 IF( k.LT.n-1 ) THEN
876*
877 d21 = a( k+1, k )
878 d11 = a( k+1, k+1 ) / d21
879 d22 = a( k, k ) / d21
880 t = one / ( d11*d22-one )
881*
882 DO 60 j = k + 2, n
883*
884* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
885*
886 wk = t*( d11*a( j, k )-a( j, k+1 ) )
887 wkp1 = t*( d22*a( j, k+1 )-a( j, k ) )
888*
889* Perform a rank-2 update of A(k+2:n,k+2:n)
890*
891 DO 50 i = j, n
892 a( i, j ) = a( i, j ) - ( a( i, k ) / d21 )*wk -
893 \$ ( a( i, k+1 ) / d21 )*wkp1
894 50 CONTINUE
895*
896* Store L(k) and L(k+1) in cols k and k+1 for row J
897*
898 a( j, k ) = wk / d21
899 a( j, k+1 ) = wkp1 / d21
900*
901 60 CONTINUE
902*
903 END IF
904*
905* Copy subdiagonal elements of D(K) to E(K) and
906* ZERO out subdiagonal entry of A
907*
908 e( k ) = a( k+1, k )
909 e( k+1 ) = zero
910 a( k+1, k ) = zero
911*
912 END IF
913*
914* End column K is nonsingular
915*
916 END IF
917*
918* Store details of the interchanges in IPIV
919*
920 IF( kstep.EQ.1 ) THEN
921 ipiv( k ) = kp
922 ELSE
923 ipiv( k ) = -p
924 ipiv( k+1 ) = -kp
925 END IF
926*
927* Increase K and return to the start of the main loop
928*
929 k = k + kstep
930 GO TO 40
931*
932 64 CONTINUE
933*
934 END IF
935*
936 RETURN
937*
938* End of DSYTF2_RK
939*
940 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dsyr(uplo, n, alpha, x, incx, a, lda)
DSYR
Definition dsyr.f:132
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 dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82