LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
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 *> \date December 2016
243 *
244 *> \ingroup doubleSYcomputational
245 *
246 *> \par Contributors:
247 * ==================
248 *>
249 *> \verbatim
250 *>
251 *> December 2016, Igor Kozachenko,
252 *> Computer Science Division,
253 *> University of California, Berkeley
254 *>
255 *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
256 *> School of Mathematics,
257 *> University of Manchester
258 *>
259 *> \endverbatim
260 *
261 * =====================================================================
262  SUBROUTINE dlasyf_rk( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
263  $ INFO )
264 *
265 * -- LAPACK computational routine (version 3.7.0) --
266 * -- LAPACK is a software package provided by Univ. of Tennessee, --
267 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
268 * December 2016
269 *
270 * .. Scalar Arguments ..
271  CHARACTER UPLO
272  INTEGER INFO, KB, LDA, LDW, N, NB
273 * ..
274 * .. Array Arguments ..
275  INTEGER IPIV( * )
276  DOUBLE PRECISION A( LDA, * ), E( * ), W( LDW, * )
277 * ..
278 *
279 * =====================================================================
280 *
281 * .. Parameters ..
282  DOUBLE PRECISION ZERO, ONE
283  parameter( zero = 0.0d+0, one = 1.0d+0 )
284  DOUBLE PRECISION EIGHT, SEVTEN
285  parameter( eight = 8.0d+0, sevten = 17.0d+0 )
286 * ..
287 * .. Local Scalars ..
288  LOGICAL DONE
289  INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, K, KK, KW, KKW,
290  $ kp, kstep, p, ii
291  DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
292  $ dtemp, r1, rowmax, t, sfmin
293 * ..
294 * .. External Functions ..
295  LOGICAL LSAME
296  INTEGER IDAMAX
297  DOUBLE PRECISION DLAMCH
298  EXTERNAL lsame, idamax, dlamch
299 * ..
300 * .. External Subroutines ..
301  EXTERNAL dcopy, dgemm, dgemv, dscal, dswap
302 * ..
303 * .. Intrinsic Functions ..
304  INTRINSIC abs, max, min, sqrt
305 * ..
306 * .. Executable Statements ..
307 *
308  info = 0
309 *
310 * Initialize ALPHA for use in choosing pivot block size.
311 *
312  alpha = ( one+sqrt( sevten ) ) / eight
313 *
314 * Compute machine safe minimum
315 *
316  sfmin = dlamch( 'S' )
317 *
318  IF( lsame( uplo, 'U' ) ) THEN
319 *
320 * Factorize the trailing columns of A using the upper triangle
321 * of A and working backwards, and compute the matrix W = U12*D
322 * for use in updating A11
323 *
324 * Initialize the first entry of array E, where superdiagonal
325 * elements of D are stored
326 *
327  e( 1 ) = zero
328 *
329 * K is the main loop index, decreasing from N in steps of 1 or 2
330 *
331  k = n
332  10 CONTINUE
333 *
334 * KW is the column of W which corresponds to column K of A
335 *
336  kw = nb + k - n
337 *
338 * Exit from loop
339 *
340  IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
341  $ GO TO 30
342 *
343  kstep = 1
344  p = k
345 *
346 * Copy column K of A to column KW of W and update it
347 *
348  CALL dcopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
349  IF( k.LT.n )
350  $ CALL dgemv( 'No transpose', k, n-k, -one, a( 1, k+1 ),
351  $ lda, w( k, kw+1 ), ldw, one, w( 1, kw ), 1 )
352 *
353 * Determine rows and columns to be interchanged and whether
354 * a 1-by-1 or 2-by-2 pivot block will be used
355 *
356  absakk = abs( w( k, kw ) )
357 *
358 * IMAX is the row-index of the largest off-diagonal element in
359 * column K, and COLMAX is its absolute value.
360 * Determine both COLMAX and IMAX.
361 *
362  IF( k.GT.1 ) THEN
363  imax = idamax( k-1, w( 1, kw ), 1 )
364  colmax = abs( w( imax, kw ) )
365  ELSE
366  colmax = zero
367  END IF
368 *
369  IF( max( absakk, colmax ).EQ.zero ) THEN
370 *
371 * Column K is zero or underflow: set INFO and continue
372 *
373  IF( info.EQ.0 )
374  $ info = k
375  kp = k
376  CALL dcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
377 *
378 * Set E( K ) to zero
379 *
380  IF( k.GT.1 )
381  $ e( k ) = zero
382 *
383  ELSE
384 *
385 * ============================================================
386 *
387 * Test for interchange
388 *
389 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
390 * (used to handle NaN and Inf)
391 *
392  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
393 *
394 * no interchange, use 1-by-1 pivot block
395 *
396  kp = k
397 *
398  ELSE
399 *
400  done = .false.
401 *
402 * Loop until pivot found
403 *
404  12 CONTINUE
405 *
406 * Begin pivot search loop body
407 *
408 *
409 * Copy column IMAX to column KW-1 of W and update it
410 *
411  CALL dcopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
412  CALL dcopy( k-imax, a( imax, imax+1 ), lda,
413  $ w( imax+1, kw-1 ), 1 )
414 *
415  IF( k.LT.n )
416  $ CALL dgemv( 'No transpose', k, n-k, -one,
417  $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
418  $ one, w( 1, kw-1 ), 1 )
419 *
420 * JMAX is the column-index of the largest off-diagonal
421 * element in row IMAX, and ROWMAX is its absolute value.
422 * Determine both ROWMAX and JMAX.
423 *
424  IF( imax.NE.k ) THEN
425  jmax = imax + idamax( k-imax, w( imax+1, kw-1 ),
426  $ 1 )
427  rowmax = abs( w( jmax, kw-1 ) )
428  ELSE
429  rowmax = zero
430  END IF
431 *
432  IF( imax.GT.1 ) THEN
433  itemp = idamax( imax-1, w( 1, kw-1 ), 1 )
434  dtemp = abs( w( itemp, kw-1 ) )
435  IF( dtemp.GT.rowmax ) THEN
436  rowmax = dtemp
437  jmax = itemp
438  END IF
439  END IF
440 *
441 * Equivalent to testing for
442 * ABS( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX
443 * (used to handle NaN and Inf)
444 *
445  IF( .NOT.(abs( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
446  $ THEN
447 *
448 * interchange rows and columns K and IMAX,
449 * use 1-by-1 pivot block
450 *
451  kp = imax
452 *
453 * copy column KW-1 of W to column KW of W
454 *
455  CALL dcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
456 *
457  done = .true.
458 *
459 * Equivalent to testing for ROWMAX.EQ.COLMAX,
460 * (used to handle NaN and Inf)
461 *
462  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
463  $ THEN
464 *
465 * interchange rows and columns K-1 and IMAX,
466 * use 2-by-2 pivot block
467 *
468  kp = imax
469  kstep = 2
470  done = .true.
471  ELSE
472 *
473 * Pivot not found: set params and repeat
474 *
475  p = imax
476  colmax = rowmax
477  imax = jmax
478 *
479 * Copy updated JMAXth (next IMAXth) column to Kth of W
480 *
481  CALL dcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
482 *
483  END IF
484 *
485 * End pivot search loop body
486 *
487  IF( .NOT. done ) GOTO 12
488 *
489  END IF
490 *
491 * ============================================================
492 *
493  kk = k - kstep + 1
494 *
495 * KKW is the column of W which corresponds to column KK of A
496 *
497  kkw = nb + kk - n
498 *
499  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
500 *
501 * Copy non-updated column K to column P
502 *
503  CALL dcopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
504  CALL dcopy( p, a( 1, k ), 1, a( 1, p ), 1 )
505 *
506 * Interchange rows K and P in last N-K+1 columns of A
507 * and last N-K+2 columns of W
508 *
509  CALL dswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
510  CALL dswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
511  END IF
512 *
513 * Updated column KP is already stored in column KKW of W
514 *
515  IF( kp.NE.kk ) THEN
516 *
517 * Copy non-updated column KK to column KP
518 *
519  a( kp, k ) = a( kk, k )
520  CALL dcopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
521  $ lda )
522  CALL dcopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
523 *
524 * Interchange rows KK and KP in last N-KK+1 columns
525 * of A and W
526 *
527  CALL dswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
528  CALL dswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
529  $ ldw )
530  END IF
531 *
532  IF( kstep.EQ.1 ) THEN
533 *
534 * 1-by-1 pivot block D(k): column KW of W now holds
535 *
536 * W(k) = U(k)*D(k)
537 *
538 * where U(k) is the k-th column of U
539 *
540 * Store U(k) in column k of A
541 *
542  CALL dcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
543  IF( k.GT.1 ) THEN
544  IF( abs( a( k, k ) ).GE.sfmin ) THEN
545  r1 = one / a( k, k )
546  CALL dscal( k-1, r1, a( 1, k ), 1 )
547  ELSE IF( a( k, k ).NE.zero ) THEN
548  DO 14 ii = 1, k - 1
549  a( ii, k ) = a( ii, k ) / a( k, k )
550  14 CONTINUE
551  END IF
552 *
553 * Store the superdiagonal element of D in array E
554 *
555  e( k ) = zero
556 *
557  END IF
558 *
559  ELSE
560 *
561 * 2-by-2 pivot block D(k): columns KW and KW-1 of W now
562 * hold
563 *
564 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
565 *
566 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
567 * of U
568 *
569  IF( k.GT.2 ) THEN
570 *
571 * Store U(k) and U(k-1) in columns k and k-1 of A
572 *
573  d12 = w( k-1, kw )
574  d11 = w( k, kw ) / d12
575  d22 = w( k-1, kw-1 ) / d12
576  t = one / ( d11*d22-one )
577  DO 20 j = 1, k - 2
578  a( j, k-1 ) = t*( (d11*w( j, kw-1 )-w( j, kw ) ) /
579  $ d12 )
580  a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
581  $ d12 )
582  20 CONTINUE
583  END IF
584 *
585 * Copy diagonal elements of D(K) to A,
586 * copy superdiagonal element of D(K) to E(K) and
587 * ZERO out superdiagonal entry of A
588 *
589  a( k-1, k-1 ) = w( k-1, kw-1 )
590  a( k-1, k ) = zero
591  a( k, k ) = w( k, kw )
592  e( k ) = w( k-1, kw )
593  e( k-1 ) = zero
594 *
595  END IF
596 *
597 * End column K is nonsingular
598 *
599  END IF
600 *
601 * Store details of the interchanges in IPIV
602 *
603  IF( kstep.EQ.1 ) THEN
604  ipiv( k ) = kp
605  ELSE
606  ipiv( k ) = -p
607  ipiv( k-1 ) = -kp
608  END IF
609 *
610 * Decrease K and return to the start of the main loop
611 *
612  k = k - kstep
613  GO TO 10
614 *
615  30 CONTINUE
616 *
617 * Update the upper triangle of A11 (= A(1:k,1:k)) as
618 *
619 * A11 := A11 - U12*D*U12**T = A11 - U12*W**T
620 *
621 * computing blocks of NB columns at a time
622 *
623  DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
624  jb = min( nb, k-j+1 )
625 *
626 * Update the upper triangle of the diagonal block
627 *
628  DO 40 jj = j, j + jb - 1
629  CALL dgemv( 'No transpose', jj-j+1, n-k, -one,
630  $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, one,
631  $ a( j, jj ), 1 )
632  40 CONTINUE
633 *
634 * Update the rectangular superdiagonal block
635 *
636  IF( j.GE.2 )
637  $ CALL dgemm( 'No transpose', 'Transpose', j-1, jb,
638  $ n-k, -one, a( 1, k+1 ), lda, w( j, kw+1 ),
639  $ ldw, one, a( 1, j ), lda )
640  50 CONTINUE
641 *
642 * Set KB to the number of columns factorized
643 *
644  kb = n - k
645 *
646  ELSE
647 *
648 * Factorize the leading columns of A using the lower triangle
649 * of A and working forwards, and compute the matrix W = L21*D
650 * for use in updating A22
651 *
652 * Initialize the unused last entry of the subdiagonal array E.
653 *
654  e( n ) = zero
655 *
656 * K is the main loop index, increasing from 1 in steps of 1 or 2
657 *
658  k = 1
659  70 CONTINUE
660 *
661 * Exit from loop
662 *
663  IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
664  $ GO TO 90
665 *
666  kstep = 1
667  p = k
668 *
669 * Copy column K of A to column K of W and update it
670 *
671  CALL dcopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
672  IF( k.GT.1 )
673  $ CALL dgemv( 'No transpose', n-k+1, k-1, -one, a( k, 1 ),
674  $ lda, w( k, 1 ), ldw, one, w( k, k ), 1 )
675 *
676 * Determine rows and columns to be interchanged and whether
677 * a 1-by-1 or 2-by-2 pivot block will be used
678 *
679  absakk = abs( w( k, k ) )
680 *
681 * IMAX is the row-index of the largest off-diagonal element in
682 * column K, and COLMAX is its absolute value.
683 * Determine both COLMAX and IMAX.
684 *
685  IF( k.LT.n ) THEN
686  imax = k + idamax( n-k, w( k+1, k ), 1 )
687  colmax = abs( w( imax, k ) )
688  ELSE
689  colmax = zero
690  END IF
691 *
692  IF( max( absakk, colmax ).EQ.zero ) THEN
693 *
694 * Column K is zero or underflow: set INFO and continue
695 *
696  IF( info.EQ.0 )
697  $ info = k
698  kp = k
699  CALL dcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
700 *
701 * Set E( K ) to zero
702 *
703  IF( k.LT.n )
704  $ e( k ) = zero
705 *
706  ELSE
707 *
708 * ============================================================
709 *
710 * Test for interchange
711 *
712 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
713 * (used to handle NaN and Inf)
714 *
715  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
716 *
717 * no interchange, use 1-by-1 pivot block
718 *
719  kp = k
720 *
721  ELSE
722 *
723  done = .false.
724 *
725 * Loop until pivot found
726 *
727  72 CONTINUE
728 *
729 * Begin pivot search loop body
730 *
731 *
732 * Copy column IMAX to column K+1 of W and update it
733 *
734  CALL dcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
735  CALL dcopy( n-imax+1, a( imax, imax ), 1,
736  $ w( imax, k+1 ), 1 )
737  IF( k.GT.1 )
738  $ CALL dgemv( 'No transpose', n-k+1, k-1, -one,
739  $ a( k, 1 ), lda, w( imax, 1 ), ldw,
740  $ one, w( k, k+1 ), 1 )
741 *
742 * JMAX is the column-index of the largest off-diagonal
743 * element in row IMAX, and ROWMAX is its absolute value.
744 * Determine both ROWMAX and JMAX.
745 *
746  IF( imax.NE.k ) THEN
747  jmax = k - 1 + idamax( imax-k, w( k, k+1 ), 1 )
748  rowmax = abs( w( jmax, k+1 ) )
749  ELSE
750  rowmax = zero
751  END IF
752 *
753  IF( imax.LT.n ) THEN
754  itemp = imax + idamax( n-imax, w( imax+1, k+1 ), 1)
755  dtemp = abs( w( itemp, k+1 ) )
756  IF( dtemp.GT.rowmax ) THEN
757  rowmax = dtemp
758  jmax = itemp
759  END IF
760  END IF
761 *
762 * Equivalent to testing for
763 * ABS( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX
764 * (used to handle NaN and Inf)
765 *
766  IF( .NOT.( abs( w( imax, k+1 ) ).LT.alpha*rowmax ) )
767  $ THEN
768 *
769 * interchange rows and columns K and IMAX,
770 * use 1-by-1 pivot block
771 *
772  kp = imax
773 *
774 * copy column K+1 of W to column K of W
775 *
776  CALL dcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
777 *
778  done = .true.
779 *
780 * Equivalent to testing for ROWMAX.EQ.COLMAX,
781 * (used to handle NaN and Inf)
782 *
783  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
784  $ THEN
785 *
786 * interchange rows and columns K+1 and IMAX,
787 * use 2-by-2 pivot block
788 *
789  kp = imax
790  kstep = 2
791  done = .true.
792  ELSE
793 *
794 * Pivot not found: set params and repeat
795 *
796  p = imax
797  colmax = rowmax
798  imax = jmax
799 *
800 * Copy updated JMAXth (next IMAXth) column to Kth of W
801 *
802  CALL dcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
803 *
804  END IF
805 *
806 * End pivot search loop body
807 *
808  IF( .NOT. done ) GOTO 72
809 *
810  END IF
811 *
812 * ============================================================
813 *
814  kk = k + kstep - 1
815 *
816  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
817 *
818 * Copy non-updated column K to column P
819 *
820  CALL dcopy( p-k, a( k, k ), 1, a( p, k ), lda )
821  CALL dcopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
822 *
823 * Interchange rows K and P in first K columns of A
824 * and first K+1 columns of W
825 *
826  CALL dswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
827  CALL dswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
828  END IF
829 *
830 * Updated column KP is already stored in column KK of W
831 *
832  IF( kp.NE.kk ) THEN
833 *
834 * Copy non-updated column KK to column KP
835 *
836  a( kp, k ) = a( kk, k )
837  CALL dcopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
838  CALL dcopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
839 *
840 * Interchange rows KK and KP in first KK columns of A and W
841 *
842  CALL dswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
843  CALL dswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
844  END IF
845 *
846  IF( kstep.EQ.1 ) THEN
847 *
848 * 1-by-1 pivot block D(k): column k of W now holds
849 *
850 * W(k) = L(k)*D(k)
851 *
852 * where L(k) is the k-th column of L
853 *
854 * Store L(k) in column k of A
855 *
856  CALL dcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
857  IF( k.LT.n ) THEN
858  IF( abs( a( k, k ) ).GE.sfmin ) THEN
859  r1 = one / a( k, k )
860  CALL dscal( n-k, r1, a( k+1, k ), 1 )
861  ELSE IF( a( k, k ).NE.zero ) THEN
862  DO 74 ii = k + 1, n
863  a( ii, k ) = a( ii, k ) / a( k, k )
864  74 CONTINUE
865  END IF
866 *
867 * Store the subdiagonal element of D in array E
868 *
869  e( k ) = zero
870 *
871  END IF
872 *
873  ELSE
874 *
875 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
876 *
877 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
878 *
879 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
880 * of L
881 *
882  IF( k.LT.n-1 ) THEN
883 *
884 * Store L(k) and L(k+1) in columns k and k+1 of A
885 *
886  d21 = w( k+1, k )
887  d11 = w( k+1, k+1 ) / d21
888  d22 = w( k, k ) / d21
889  t = one / ( d11*d22-one )
890  DO 80 j = k + 2, n
891  a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
892  $ d21 )
893  a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
894  $ d21 )
895  80 CONTINUE
896  END IF
897 *
898 * Copy diagonal elements of D(K) to A,
899 * copy subdiagonal element of D(K) to E(K) and
900 * ZERO out subdiagonal entry of A
901 *
902  a( k, k ) = w( k, k )
903  a( k+1, k ) = zero
904  a( k+1, k+1 ) = w( k+1, k+1 )
905  e( k ) = w( k+1, k )
906  e( k+1 ) = zero
907 *
908  END IF
909 *
910 * End column K is nonsingular
911 *
912  END IF
913 *
914 * Store details of the interchanges in IPIV
915 *
916  IF( kstep.EQ.1 ) THEN
917  ipiv( k ) = kp
918  ELSE
919  ipiv( k ) = -p
920  ipiv( k+1 ) = -kp
921  END IF
922 *
923 * Increase K and return to the start of the main loop
924 *
925  k = k + kstep
926  GO TO 70
927 *
928  90 CONTINUE
929 *
930 * Update the lower triangle of A22 (= A(k:n,k:n)) as
931 *
932 * A22 := A22 - L21*D*L21**T = A22 - L21*W**T
933 *
934 * computing blocks of NB columns at a time
935 *
936  DO 110 j = k, n, nb
937  jb = min( nb, n-j+1 )
938 *
939 * Update the lower triangle of the diagonal block
940 *
941  DO 100 jj = j, j + jb - 1
942  CALL dgemv( 'No transpose', j+jb-jj, k-1, -one,
943  $ a( jj, 1 ), lda, w( jj, 1 ), ldw, one,
944  $ a( jj, jj ), 1 )
945  100 CONTINUE
946 *
947 * Update the rectangular subdiagonal block
948 *
949  IF( j+jb.LE.n )
950  $ CALL dgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
951  $ k-1, -one, a( j+jb, 1 ), lda, w( j, 1 ),
952  $ ldw, one, a( j+jb, j ), lda )
953  110 CONTINUE
954 *
955 * Set KB to the number of columns factorized
956 *
957  kb = k - 1
958 *
959  END IF
960 *
961  RETURN
962 *
963 * End of DLASYF_RK
964 *
965  END
dlasyf_rk
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:264
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
dgemv
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
dswap
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:84
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81