LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
slasyf_rk.f
Go to the documentation of this file.
1 *> \brief \b SLASYF_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 SLASYF_RK + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasyf_rk.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasyf_rk.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasyf_rk.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLASYF_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 * REAL A( LDA, * ), E( * ), W( LDW, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *> SLASYF_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 *> SLASYF_RK is an auxiliary routine called by SSYTRF_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 REAL 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 REAL 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 REAL 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 singleSYcomputational
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 slasyf_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  REAL A( LDA, * ), E( * ), W( LDW, * )
274 * ..
275 *
276 * =====================================================================
277 *
278 * .. Parameters ..
279  REAL ZERO, ONE
280  parameter( zero = 0.0e+0, one = 1.0e+0 )
281  REAL EIGHT, SEVTEN
282  parameter( eight = 8.0e+0, sevten = 17.0e+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  REAL ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
289  $ stemp, r1, rowmax, t, sfmin
290 * ..
291 * .. External Functions ..
292  LOGICAL LSAME
293  INTEGER ISAMAX
294  REAL SLAMCH
295  EXTERNAL lsame, isamax, slamch
296 * ..
297 * .. External Subroutines ..
298  EXTERNAL scopy, sgemm, sgemv, sscal, sswap
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 = slamch( '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 scopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
346  IF( k.LT.n )
347  $ CALL sgemv( '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 = isamax( 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 scopy( 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 scopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
409  CALL scopy( k-imax, a( imax, imax+1 ), lda,
410  $ w( imax+1, kw-1 ), 1 )
411 *
412  IF( k.LT.n )
413  $ CALL sgemv( '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 + isamax( 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 = isamax( imax-1, w( 1, kw-1 ), 1 )
431  stemp = abs( w( itemp, kw-1 ) )
432  IF( stemp.GT.rowmax ) THEN
433  rowmax = stemp
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 scopy( 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 scopy( 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 scopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
501  CALL scopy( 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 sswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
507  CALL sswap( 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 scopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
518  $ lda )
519  CALL scopy( 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 sswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
525  CALL sswap( 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 scopy( 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 sscal( 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 sgemv( '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 sgemm( '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 scopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
669  IF( k.GT.1 )
670  $ CALL sgemv( '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 + isamax( 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 scopy( 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 scopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
732  CALL scopy( n-imax+1, a( imax, imax ), 1,
733  $ w( imax, k+1 ), 1 )
734  IF( k.GT.1 )
735  $ CALL sgemv( '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 + isamax( 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 + isamax( n-imax, w( imax+1, k+1 ), 1)
752  stemp = abs( w( itemp, k+1 ) )
753  IF( stemp.GT.rowmax ) THEN
754  rowmax = stemp
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 scopy( 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 scopy( 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 scopy( p-k, a( k, k ), 1, a( p, k ), lda )
818  CALL scopy( 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 sswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
824  CALL sswap( 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 scopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
835  CALL scopy( 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 sswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
840  CALL sswap( 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 scopy( 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 sscal( 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 sgemv( '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 sgemm( '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 SLASYF_RK
961 *
962  END
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
subroutine slasyf_rk(UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW, INFO)
SLASYF_RK computes a partial factorization of a real symmetric indefinite matrix using bounded Bunch-...
Definition: slasyf_rk.f:262