LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dsytf2_rook.f
Go to the documentation of this file.
1 *> \brief \b DSYTF2_ROOK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Kaufman ("rook") 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
9 *> Download DSYTF2_ROOK + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytf2_rook.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytf2_rook.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytf2_rook.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DSYTF2_ROOK( 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_ROOK computes the factorization of a real symmetric matrix A
39 *> using the bounded Bunch-Kaufman ("rook") 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)
96 *> were interchanged and D(k,k) is a 1-by-1 diagonal block.
97 *>
98 *> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
99 *> columns k and -IPIV(k) were interchanged and rows and
100 *> columns k-1 and -IPIV(k-1) were inerchaged,
101 *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
102 *>
103 *> If UPLO = 'L':
104 *> If IPIV(k) > 0, then rows and columns k and IPIV(k)
105 *> were interchanged and D(k,k) is a 1-by-1 diagonal block.
106 *>
107 *> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
108 *> columns k and -IPIV(k) were interchanged and rows and
109 *> columns k+1 and -IPIV(k+1) were inerchaged,
110 *> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
111 *> \endverbatim
112 *>
113 *> \param[out] INFO
114 *> \verbatim
115 *> INFO is INTEGER
116 *> = 0: successful exit
117 *> < 0: if INFO = -k, the k-th argument had an illegal value
118 *> > 0: if INFO = k, D(k,k) is exactly zero. The factorization
119 *> has been completed, but the block diagonal matrix D is
120 *> exactly singular, and division by zero will occur if it
121 *> is used to solve a system of equations.
122 *> \endverbatim
123 *
124 * Authors:
125 * ========
126 *
127 *> \author Univ. of Tennessee
128 *> \author Univ. of California Berkeley
129 *> \author Univ. of Colorado Denver
130 *> \author NAG Ltd.
131 *
132 *> \ingroup doubleSYcomputational
133 *
134 *> \par Further Details:
135 * =====================
136 *>
137 *> \verbatim
138 *>
139 *> If UPLO = 'U', then A = U*D*U**T, where
140 *> U = P(n)*U(n)* ... *P(k)U(k)* ...,
141 *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to
142 *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
143 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
144 *> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
145 *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
146 *>
147 *> ( I v 0 ) k-s
148 *> U(k) = ( 0 I 0 ) s
149 *> ( 0 0 I ) n-k
150 *> k-s s n-k
151 *>
152 *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
153 *> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
154 *> and A(k,k), and v overwrites A(1:k-2,k-1:k).
155 *>
156 *> If UPLO = 'L', then A = L*D*L**T, where
157 *> L = P(1)*L(1)* ... *P(k)*L(k)* ...,
158 *> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
159 *> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
160 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
161 *> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
162 *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
163 *>
164 *> ( I 0 0 ) k-1
165 *> L(k) = ( 0 I 0 ) s
166 *> ( 0 v I ) n-k-s+1
167 *> k-1 s n-k-s+1
168 *>
169 *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
170 *> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
171 *> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
172 *> \endverbatim
173 *
174 *> \par Contributors:
175 * ==================
176 *>
177 *> \verbatim
178 *>
179 *> November 2013, Igor Kozachenko,
180 *> Computer Science Division,
181 *> University of California, Berkeley
182 *>
183 *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
184 *> School of Mathematics,
185 *> University of Manchester
186 *>
187 *> 01-01-96 - Based on modifications by
188 *> J. Lewis, Boeing Computer Services Company
189 *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville abd , USA
190 *> \endverbatim
191 *
192 * =====================================================================
193  SUBROUTINE dsytf2_rook( 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, DONE
218  INTEGER I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP,
219  $ P, II
220  DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
221  $ ROWMAX, DTEMP, T, WK, WKM1, WKP1, SFMIN
222 * ..
223 * .. External Functions ..
224  LOGICAL LSAME
225  INTEGER IDAMAX
226  DOUBLE PRECISION DLAMCH
227  EXTERNAL lsame, idamax, dlamch
228 * ..
229 * .. External Subroutines ..
230  EXTERNAL dscal, dswap, dsyr, xerbla
231 * ..
232 * .. Intrinsic Functions ..
233  INTRINSIC abs, max, sqrt
234 * ..
235 * .. Executable Statements ..
236 *
237 * Test the input parameters.
238 *
239  info = 0
240  upper = lsame( uplo, 'U' )
241  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
242  info = -1
243  ELSE IF( n.LT.0 ) THEN
244  info = -2
245  ELSE IF( lda.LT.max( 1, n ) ) THEN
246  info = -4
247  END IF
248  IF( info.NE.0 ) THEN
249  CALL xerbla( 'DSYTF2_ROOK', -info )
250  RETURN
251  END IF
252 *
253 * Initialize ALPHA for use in choosing pivot block size.
254 *
255  alpha = ( one+sqrt( sevten ) ) / eight
256 *
257 * Compute machine safe minimum
258 *
259  sfmin = dlamch( 'S' )
260 *
261  IF( upper ) THEN
262 *
263 * Factorize A as U*D*U**T using the upper triangle of A
264 *
265 * K is the main loop index, decreasing from N to 1 in steps of
266 * 1 or 2
267 *
268  k = n
269  10 CONTINUE
270 *
271 * If K < 1, exit from loop
272 *
273  IF( k.LT.1 )
274  $ GO TO 70
275  kstep = 1
276  p = k
277 *
278 * Determine rows and columns to be interchanged and whether
279 * a 1-by-1 or 2-by-2 pivot block will be used
280 *
281  absakk = abs( a( k, k ) )
282 *
283 * IMAX is the row-index of the largest off-diagonal element in
284 * column K, and COLMAX is its absolute value.
285 * Determine both COLMAX and IMAX.
286 *
287  IF( k.GT.1 ) THEN
288  imax = idamax( k-1, a( 1, k ), 1 )
289  colmax = abs( a( imax, k ) )
290  ELSE
291  colmax = zero
292  END IF
293 *
294  IF( (max( absakk, colmax ).EQ.zero) ) THEN
295 *
296 * Column K is zero or underflow: set INFO and continue
297 *
298  IF( info.EQ.0 )
299  $ info = k
300  kp = k
301  ELSE
302 *
303 * Test for interchange
304 *
305 * Equivalent to testing for (used to handle NaN and Inf)
306 * ABSAKK.GE.ALPHA*COLMAX
307 *
308  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
309 *
310 * no interchange,
311 * use 1-by-1 pivot block
312 *
313  kp = k
314  ELSE
315 *
316  done = .false.
317 *
318 * Loop until pivot found
319 *
320  12 CONTINUE
321 *
322 * Begin pivot search loop body
323 *
324 * JMAX is the column-index of the largest off-diagonal
325 * element in row IMAX, and ROWMAX is its absolute value.
326 * Determine both ROWMAX and JMAX.
327 *
328  IF( imax.NE.k ) THEN
329  jmax = imax + idamax( k-imax, a( imax, imax+1 ),
330  $ lda )
331  rowmax = abs( a( imax, jmax ) )
332  ELSE
333  rowmax = zero
334  END IF
335 *
336  IF( imax.GT.1 ) THEN
337  itemp = idamax( imax-1, a( 1, imax ), 1 )
338  dtemp = abs( a( itemp, imax ) )
339  IF( dtemp.GT.rowmax ) THEN
340  rowmax = dtemp
341  jmax = itemp
342  END IF
343  END IF
344 *
345 * Equivalent to testing for (used to handle NaN and Inf)
346 * ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
347 *
348  IF( .NOT.( abs( a( imax, imax ) ).LT.alpha*rowmax ) )
349  $ THEN
350 *
351 * interchange rows and columns K and IMAX,
352 * use 1-by-1 pivot block
353 *
354  kp = imax
355  done = .true.
356 *
357 * Equivalent to testing for ROWMAX .EQ. COLMAX,
358 * used to handle NaN and Inf
359 *
360  ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
361 *
362 * interchange rows and columns K+1 and IMAX,
363 * use 2-by-2 pivot block
364 *
365  kp = imax
366  kstep = 2
367  done = .true.
368  ELSE
369 *
370 * Pivot NOT found, set variables and repeat
371 *
372  p = imax
373  colmax = rowmax
374  imax = jmax
375  END IF
376 *
377 * End pivot search loop body
378 *
379  IF( .NOT. done ) GOTO 12
380 *
381  END IF
382 *
383 * Swap TWO rows and TWO columns
384 *
385 * First swap
386 *
387  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
388 *
389 * Interchange rows and column K and P in the leading
390 * submatrix A(1:k,1:k) if we have a 2-by-2 pivot
391 *
392  IF( p.GT.1 )
393  $ CALL dswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
394  IF( p.LT.(k-1) )
395  $ CALL dswap( k-p-1, a( p+1, k ), 1, a( p, p+1 ),
396  $ lda )
397  t = a( k, k )
398  a( k, k ) = a( p, p )
399  a( p, p ) = t
400  END IF
401 *
402 * Second swap
403 *
404  kk = k - kstep + 1
405  IF( kp.NE.kk ) THEN
406 *
407 * Interchange rows and columns KK and KP in the leading
408 * submatrix A(1:k,1:k)
409 *
410  IF( kp.GT.1 )
411  $ CALL dswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
412  IF( ( kk.GT.1 ) .AND. ( kp.LT.(kk-1) ) )
413  $ CALL dswap( kk-kp-1, a( kp+1, kk ), 1, a( kp, kp+1 ),
414  $ lda )
415  t = a( kk, kk )
416  a( kk, kk ) = a( kp, kp )
417  a( kp, kp ) = t
418  IF( kstep.EQ.2 ) THEN
419  t = a( k-1, k )
420  a( k-1, k ) = a( kp, k )
421  a( kp, k ) = t
422  END IF
423  END IF
424 *
425 * Update the leading submatrix
426 *
427  IF( kstep.EQ.1 ) THEN
428 *
429 * 1-by-1 pivot block D(k): column k now holds
430 *
431 * W(k) = U(k)*D(k)
432 *
433 * where U(k) is the k-th column of U
434 *
435  IF( k.GT.1 ) THEN
436 *
437 * Perform a rank-1 update of A(1:k-1,1:k-1) and
438 * store U(k) in column k
439 *
440  IF( abs( a( k, k ) ).GE.sfmin ) THEN
441 *
442 * Perform a rank-1 update of A(1:k-1,1:k-1) as
443 * A := A - U(k)*D(k)*U(k)**T
444 * = A - W(k)*1/D(k)*W(k)**T
445 *
446  d11 = one / a( k, k )
447  CALL dsyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
448 *
449 * Store U(k) in column k
450 *
451  CALL dscal( k-1, d11, a( 1, k ), 1 )
452  ELSE
453 *
454 * Store L(k) in column K
455 *
456  d11 = a( k, k )
457  DO 16 ii = 1, k - 1
458  a( ii, k ) = a( ii, k ) / d11
459  16 CONTINUE
460 *
461 * Perform a rank-1 update of A(k+1:n,k+1:n) as
462 * A := A - U(k)*D(k)*U(k)**T
463 * = A - W(k)*(1/D(k))*W(k)**T
464 * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
465 *
466  CALL dsyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
467  END IF
468  END IF
469 *
470  ELSE
471 *
472 * 2-by-2 pivot block D(k): columns k and k-1 now hold
473 *
474 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
475 *
476 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
477 * of U
478 *
479 * Perform a rank-2 update of A(1:k-2,1:k-2) as
480 *
481 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
482 * = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
483 *
484 * and store L(k) and L(k+1) in columns k and k+1
485 *
486  IF( k.GT.2 ) THEN
487 *
488  d12 = a( k-1, k )
489  d22 = a( k-1, k-1 ) / d12
490  d11 = a( k, k ) / d12
491  t = one / ( d11*d22-one )
492 *
493  DO 30 j = k - 2, 1, -1
494 *
495  wkm1 = t*( d11*a( j, k-1 )-a( j, k ) )
496  wk = t*( d22*a( j, k )-a( j, k-1 ) )
497 *
498  DO 20 i = j, 1, -1
499  a( i, j ) = a( i, j ) - (a( i, k ) / d12 )*wk -
500  $ ( a( i, k-1 ) / d12 )*wkm1
501  20 CONTINUE
502 *
503 * Store U(k) and U(k-1) in cols k and k-1 for row J
504 *
505  a( j, k ) = wk / d12
506  a( j, k-1 ) = wkm1 / d12
507 *
508  30 CONTINUE
509 *
510  END IF
511 *
512  END IF
513  END IF
514 *
515 * Store details of the interchanges in IPIV
516 *
517  IF( kstep.EQ.1 ) THEN
518  ipiv( k ) = kp
519  ELSE
520  ipiv( k ) = -p
521  ipiv( k-1 ) = -kp
522  END IF
523 *
524 * Decrease K and return to the start of the main loop
525 *
526  k = k - kstep
527  GO TO 10
528 *
529  ELSE
530 *
531 * Factorize A as L*D*L**T using the lower triangle of A
532 *
533 * K is the main loop index, increasing from 1 to N in steps of
534 * 1 or 2
535 *
536  k = 1
537  40 CONTINUE
538 *
539 * If K > N, exit from loop
540 *
541  IF( k.GT.n )
542  $ GO TO 70
543  kstep = 1
544  p = k
545 *
546 * Determine rows and columns to be interchanged and whether
547 * a 1-by-1 or 2-by-2 pivot block will be used
548 *
549  absakk = abs( a( k, k ) )
550 *
551 * IMAX is the row-index of the largest off-diagonal element in
552 * column K, and COLMAX is its absolute value.
553 * Determine both COLMAX and IMAX.
554 *
555  IF( k.LT.n ) THEN
556  imax = k + idamax( n-k, a( k+1, k ), 1 )
557  colmax = abs( a( imax, k ) )
558  ELSE
559  colmax = zero
560  END IF
561 *
562  IF( ( max( absakk, colmax ).EQ.zero ) ) THEN
563 *
564 * Column K is zero or underflow: set INFO and continue
565 *
566  IF( info.EQ.0 )
567  $ info = k
568  kp = k
569  ELSE
570 *
571 * Test for interchange
572 *
573 * Equivalent to testing for (used to handle NaN and Inf)
574 * ABSAKK.GE.ALPHA*COLMAX
575 *
576  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
577 *
578 * no interchange, use 1-by-1 pivot block
579 *
580  kp = k
581  ELSE
582 *
583  done = .false.
584 *
585 * Loop until pivot found
586 *
587  42 CONTINUE
588 *
589 * Begin pivot search loop body
590 *
591 * JMAX is the column-index of the largest off-diagonal
592 * element in row IMAX, and ROWMAX is its absolute value.
593 * Determine both ROWMAX and JMAX.
594 *
595  IF( imax.NE.k ) THEN
596  jmax = k - 1 + idamax( imax-k, a( imax, k ), lda )
597  rowmax = abs( a( imax, jmax ) )
598  ELSE
599  rowmax = zero
600  END IF
601 *
602  IF( imax.LT.n ) THEN
603  itemp = imax + idamax( n-imax, a( imax+1, imax ),
604  $ 1 )
605  dtemp = abs( a( itemp, imax ) )
606  IF( dtemp.GT.rowmax ) THEN
607  rowmax = dtemp
608  jmax = itemp
609  END IF
610  END IF
611 *
612 * Equivalent to testing for (used to handle NaN and Inf)
613 * ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
614 *
615  IF( .NOT.( abs( a( imax, imax ) ).LT.alpha*rowmax ) )
616  $ THEN
617 *
618 * interchange rows and columns K and IMAX,
619 * use 1-by-1 pivot block
620 *
621  kp = imax
622  done = .true.
623 *
624 * Equivalent to testing for ROWMAX .EQ. COLMAX,
625 * used to handle NaN and Inf
626 *
627  ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
628 *
629 * interchange rows and columns K+1 and IMAX,
630 * use 2-by-2 pivot block
631 *
632  kp = imax
633  kstep = 2
634  done = .true.
635  ELSE
636 *
637 * Pivot NOT found, set variables and repeat
638 *
639  p = imax
640  colmax = rowmax
641  imax = jmax
642  END IF
643 *
644 * End pivot search loop body
645 *
646  IF( .NOT. done ) GOTO 42
647 *
648  END IF
649 *
650 * Swap TWO rows and TWO columns
651 *
652 * First swap
653 *
654  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
655 *
656 * Interchange rows and column K and P in the trailing
657 * submatrix A(k:n,k:n) if we have a 2-by-2 pivot
658 *
659  IF( p.LT.n )
660  $ CALL dswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
661  IF( p.GT.(k+1) )
662  $ CALL dswap( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
663  t = a( k, k )
664  a( k, k ) = a( p, p )
665  a( p, p ) = t
666  END IF
667 *
668 * Second swap
669 *
670  kk = k + kstep - 1
671  IF( kp.NE.kk ) THEN
672 *
673 * Interchange rows and columns KK and KP in the trailing
674 * submatrix A(k:n,k:n)
675 *
676  IF( kp.LT.n )
677  $ CALL dswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
678  IF( ( kk.LT.n ) .AND. ( kp.GT.(kk+1) ) )
679  $ CALL dswap( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
680  $ lda )
681  t = a( kk, kk )
682  a( kk, kk ) = a( kp, kp )
683  a( kp, kp ) = t
684  IF( kstep.EQ.2 ) THEN
685  t = a( k+1, k )
686  a( k+1, k ) = a( kp, k )
687  a( kp, k ) = t
688  END IF
689  END IF
690 *
691 * Update the trailing submatrix
692 *
693  IF( kstep.EQ.1 ) THEN
694 *
695 * 1-by-1 pivot block D(k): column k now holds
696 *
697 * W(k) = L(k)*D(k)
698 *
699 * where L(k) is the k-th column of L
700 *
701  IF( k.LT.n ) THEN
702 *
703 * Perform a rank-1 update of A(k+1:n,k+1:n) and
704 * store L(k) in column k
705 *
706  IF( abs( a( k, k ) ).GE.sfmin ) THEN
707 *
708 * Perform a rank-1 update of A(k+1:n,k+1:n) as
709 * A := A - L(k)*D(k)*L(k)**T
710 * = A - W(k)*(1/D(k))*W(k)**T
711 *
712  d11 = one / a( k, k )
713  CALL dsyr( uplo, n-k, -d11, a( k+1, k ), 1,
714  $ a( k+1, k+1 ), lda )
715 *
716 * Store L(k) in column k
717 *
718  CALL dscal( n-k, d11, a( k+1, k ), 1 )
719  ELSE
720 *
721 * Store L(k) in column k
722 *
723  d11 = a( k, k )
724  DO 46 ii = k + 1, n
725  a( ii, k ) = a( ii, k ) / d11
726  46 CONTINUE
727 *
728 * Perform a rank-1 update of A(k+1:n,k+1:n) as
729 * A := A - L(k)*D(k)*L(k)**T
730 * = A - W(k)*(1/D(k))*W(k)**T
731 * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
732 *
733  CALL dsyr( uplo, n-k, -d11, a( k+1, k ), 1,
734  $ a( k+1, k+1 ), lda )
735  END IF
736  END IF
737 *
738  ELSE
739 *
740 * 2-by-2 pivot block D(k): columns k and k+1 now hold
741 *
742 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
743 *
744 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
745 * of L
746 *
747 *
748 * Perform a rank-2 update of A(k+2:n,k+2:n) as
749 *
750 * A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
751 * = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
752 *
753 * and store L(k) and L(k+1) in columns k and k+1
754 *
755  IF( k.LT.n-1 ) THEN
756 *
757  d21 = a( k+1, k )
758  d11 = a( k+1, k+1 ) / d21
759  d22 = a( k, k ) / d21
760  t = one / ( d11*d22-one )
761 *
762  DO 60 j = k + 2, n
763 *
764 * Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
765 *
766  wk = t*( d11*a( j, k )-a( j, k+1 ) )
767  wkp1 = t*( d22*a( j, k+1 )-a( j, k ) )
768 *
769 * Perform a rank-2 update of A(k+2:n,k+2:n)
770 *
771  DO 50 i = j, n
772  a( i, j ) = a( i, j ) - ( a( i, k ) / d21 )*wk -
773  $ ( a( i, k+1 ) / d21 )*wkp1
774  50 CONTINUE
775 *
776 * Store L(k) and L(k+1) in cols k and k+1 for row J
777 *
778  a( j, k ) = wk / d21
779  a( j, k+1 ) = wkp1 / d21
780 *
781  60 CONTINUE
782 *
783  END IF
784 *
785  END IF
786  END IF
787 *
788 * Store details of the interchanges in IPIV
789 *
790  IF( kstep.EQ.1 ) THEN
791  ipiv( k ) = kp
792  ELSE
793  ipiv( k ) = -p
794  ipiv( k+1 ) = -kp
795  END IF
796 *
797 * Increase K and return to the start of the main loop
798 *
799  k = k + kstep
800  GO TO 40
801 *
802  END IF
803 *
804  70 CONTINUE
805 *
806  RETURN
807 *
808 * End of DSYTF2_ROOK
809 *
810  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_rook(UPLO, N, A, LDA, IPIV, INFO)
DSYTF2_ROOK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-...
Definition: dsytf2_rook.f:194