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