LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
ssytrf_rk.f
Go to the documentation of this file.
1 *> \brief \b SSYTRF_RK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Kaufman (rook) diagonal pivoting method (BLAS3 blocked algorithm).
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SSYTRF_RK + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssytrf_rk.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssytrf_rk.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssytrf_rk.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SSYTRF_RK( UPLO, N, A, LDA, E, IPIV, WORK, LWORK,
22 * INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER UPLO
26 * INTEGER INFO, LDA, LWORK, N
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IPIV( * )
30 * REAL A( LDA, * ), E ( * ), WORK( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *> SSYTRF_RK computes the factorization of a real symmetric matrix A
39 *> using the bounded Bunch-Kaufman (rook) diagonal pivoting method:
40 *>
41 *> A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T),
42 *>
43 *> where U (or L) is unit upper (or lower) triangular matrix,
44 *> U**T (or L**T) is the transpose of U (or L), P is a permutation
45 *> matrix, P**T is the transpose of P, and D is symmetric and block
46 *> diagonal with 1-by-1 and 2-by-2 diagonal blocks.
47 *>
48 *> This is the blocked version of the algorithm, calling Level 3 BLAS.
49 *> For more information see Further Details section.
50 *> \endverbatim
51 *
52 * Arguments:
53 * ==========
54 *
55 *> \param[in] UPLO
56 *> \verbatim
57 *> UPLO is CHARACTER*1
58 *> Specifies whether the upper or lower triangular part of the
59 *> symmetric matrix A is stored:
60 *> = 'U': Upper triangular
61 *> = 'L': Lower triangular
62 *> \endverbatim
63 *>
64 *> \param[in] N
65 *> \verbatim
66 *> N is INTEGER
67 *> The order of the matrix A. N >= 0.
68 *> \endverbatim
69 *>
70 *> \param[in,out] A
71 *> \verbatim
72 *> A is REAL array, dimension (LDA,N)
73 *> On entry, the symmetric matrix A.
74 *> If UPLO = 'U': the leading N-by-N upper triangular part
75 *> of A contains the upper triangular part of the matrix A,
76 *> and the strictly lower triangular part of A is not
77 *> referenced.
78 *>
79 *> If UPLO = 'L': the leading N-by-N lower triangular part
80 *> of A contains the lower triangular part of the matrix A,
81 *> and the strictly upper triangular part of A is not
82 *> referenced.
83 *>
84 *> On exit, contains:
85 *> a) ONLY diagonal elements of the symmetric block diagonal
86 *> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
87 *> (superdiagonal (or subdiagonal) elements of D
88 *> are stored on exit in array E), and
89 *> b) If UPLO = 'U': factor U in the superdiagonal part of A.
90 *> If UPLO = 'L': factor L in the subdiagonal part of A.
91 *> \endverbatim
92 *>
93 *> \param[in] LDA
94 *> \verbatim
95 *> LDA is INTEGER
96 *> The leading dimension of the array A. LDA >= max(1,N).
97 *> \endverbatim
98 *>
99 *> \param[out] E
100 *> \verbatim
101 *> E is REAL array, dimension (N)
102 *> On exit, contains the superdiagonal (or subdiagonal)
103 *> elements of the symmetric block diagonal matrix D
104 *> with 1-by-1 or 2-by-2 diagonal blocks, where
105 *> If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
106 *> If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.
107 *>
108 *> NOTE: For 1-by-1 diagonal block D(k), where
109 *> 1 <= k <= N, the element E(k) is set to 0 in both
110 *> UPLO = 'U' or UPLO = 'L' cases.
111 *> \endverbatim
112 *>
113 *> \param[out] IPIV
114 *> \verbatim
115 *> IPIV is INTEGER array, dimension (N)
116 *> IPIV describes the permutation matrix P in the factorization
117 *> of matrix A as follows. The absolute value of IPIV(k)
118 *> represents the index of row and column that were
119 *> interchanged with the k-th row and column. The value of UPLO
120 *> describes the order in which the interchanges were applied.
121 *> Also, the sign of IPIV represents the block structure of
122 *> the symmetric block diagonal matrix D with 1-by-1 or 2-by-2
123 *> diagonal blocks which correspond to 1 or 2 interchanges
124 *> at each factorization step. For more info see Further
125 *> Details section.
126 *>
127 *> If UPLO = 'U',
128 *> ( in factorization order, k decreases from N to 1 ):
129 *> a) A single positive entry IPIV(k) > 0 means:
130 *> D(k,k) is a 1-by-1 diagonal block.
131 *> If IPIV(k) != k, rows and columns k and IPIV(k) were
132 *> interchanged in the matrix A(1:N,1:N);
133 *> If IPIV(k) = k, no interchange occurred.
134 *>
135 *> b) A pair of consecutive negative entries
136 *> IPIV(k) < 0 and IPIV(k-1) < 0 means:
137 *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
138 *> (NOTE: negative entries in IPIV appear ONLY in pairs).
139 *> 1) If -IPIV(k) != k, rows and columns
140 *> k and -IPIV(k) were interchanged
141 *> in the matrix A(1:N,1:N).
142 *> If -IPIV(k) = k, no interchange occurred.
143 *> 2) If -IPIV(k-1) != k-1, rows and columns
144 *> k-1 and -IPIV(k-1) were interchanged
145 *> in the matrix A(1:N,1:N).
146 *> If -IPIV(k-1) = k-1, no interchange occurred.
147 *>
148 *> c) In both cases a) and b), always ABS( IPIV(k) ) <= k.
149 *>
150 *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
151 *>
152 *> If UPLO = 'L',
153 *> ( in factorization order, k increases from 1 to N ):
154 *> a) A single positive entry IPIV(k) > 0 means:
155 *> D(k,k) is a 1-by-1 diagonal block.
156 *> If IPIV(k) != k, rows and columns k and IPIV(k) were
157 *> interchanged in the matrix A(1:N,1:N).
158 *> If IPIV(k) = k, no interchange occurred.
159 *>
160 *> b) A pair of consecutive negative entries
161 *> IPIV(k) < 0 and IPIV(k+1) < 0 means:
162 *> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
163 *> (NOTE: negative entries in IPIV appear ONLY in pairs).
164 *> 1) If -IPIV(k) != k, rows and columns
165 *> k and -IPIV(k) were interchanged
166 *> in the matrix A(1:N,1:N).
167 *> If -IPIV(k) = k, no interchange occurred.
168 *> 2) If -IPIV(k+1) != k+1, rows and columns
169 *> k-1 and -IPIV(k-1) were interchanged
170 *> in the matrix A(1:N,1:N).
171 *> If -IPIV(k+1) = k+1, no interchange occurred.
172 *>
173 *> c) In both cases a) and b), always ABS( IPIV(k) ) >= k.
174 *>
175 *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
176 *> \endverbatim
177 *>
178 *> \param[out] WORK
179 *> \verbatim
180 *> WORK is REAL array, dimension ( MAX(1,LWORK) ).
181 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
182 *> \endverbatim
183 *>
184 *> \param[in] LWORK
185 *> \verbatim
186 *> LWORK is INTEGER
187 *> The length of WORK. LWORK >=1. For best performance
188 *> LWORK >= N*NB, where NB is the block size returned
189 *> by ILAENV.
190 *>
191 *> If LWORK = -1, then a workspace query is assumed;
192 *> the routine only calculates the optimal size of the WORK
193 *> array, returns this value as the first entry of the WORK
194 *> array, and no error message related to LWORK is issued
195 *> by XERBLA.
196 *> \endverbatim
197 *>
198 *> \param[out] INFO
199 *> \verbatim
200 *> INFO is INTEGER
201 *> = 0: successful exit
202 *>
203 *> < 0: If INFO = -k, the k-th argument had an illegal value
204 *>
205 *> > 0: If INFO = k, the matrix A is singular, because:
206 *> If UPLO = 'U': column k in the upper
207 *> triangular part of A contains all zeros.
208 *> If UPLO = 'L': column k in the lower
209 *> triangular part of A contains all zeros.
210 *>
211 *> Therefore D(k,k) is exactly zero, and superdiagonal
212 *> elements of column k of U (or subdiagonal elements of
213 *> column k of L ) are all zeros. The factorization has
214 *> been completed, but the block diagonal matrix D is
215 *> exactly singular, and division by zero will occur if
216 *> it is used to solve a system of equations.
217 *>
218 *> NOTE: INFO only stores the first occurrence of
219 *> a singularity, any subsequent occurrence of singularity
220 *> is not stored in INFO even though the factorization
221 *> always completes.
222 *> \endverbatim
223 *
224 * Authors:
225 * ========
226 *
227 *> \author Univ. of Tennessee
228 *> \author Univ. of California Berkeley
229 *> \author Univ. of Colorado Denver
230 *> \author NAG Ltd.
231 *
232 *> \date December 2016
233 *
234 *> \ingroup singleSYcomputational
235 *
236 *> \par Further Details:
237 * =====================
238 *>
239 *> \verbatim
240 *> TODO: put correct description
241 *> \endverbatim
242 *
243 *> \par Contributors:
244 * ==================
245 *>
246 *> \verbatim
247 *>
248 *> December 2016, Igor Kozachenko,
249 *> Computer Science Division,
250 *> University of California, Berkeley
251 *>
252 *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
253 *> School of Mathematics,
254 *> University of Manchester
255 *>
256 *> \endverbatim
257 *
258 * =====================================================================
259  SUBROUTINE ssytrf_rk( UPLO, N, A, LDA, E, IPIV, WORK, LWORK,
260  $ INFO )
261 *
262 * -- LAPACK computational routine (version 3.7.0) --
263 * -- LAPACK is a software package provided by Univ. of Tennessee, --
264 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
265 * December 2016
266 *
267 * .. Scalar Arguments ..
268  CHARACTER UPLO
269  INTEGER INFO, LDA, LWORK, N
270 * ..
271 * .. Array Arguments ..
272  INTEGER IPIV( * )
273  REAL A( lda, * ), E( * ), WORK( * )
274 * ..
275 *
276 * =====================================================================
277 *
278 * .. Local Scalars ..
279  LOGICAL LQUERY, UPPER
280  INTEGER I, IINFO, IP, IWS, K, KB, LDWORK, LWKOPT,
281  $ nb, nbmin
282 * ..
283 * .. External Functions ..
284  LOGICAL LSAME
285  INTEGER ILAENV
286  EXTERNAL lsame, ilaenv
287 * ..
288 * .. External Subroutines ..
289  EXTERNAL slasyf_rk, ssytf2_rk, sswap, xerbla
290 * ..
291 * .. Intrinsic Functions ..
292  INTRINSIC abs, max
293 * ..
294 * .. Executable Statements ..
295 *
296 * Test the input parameters.
297 *
298  info = 0
299  upper = lsame( uplo, 'U' )
300  lquery = ( lwork.EQ.-1 )
301  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
302  info = -1
303  ELSE IF( n.LT.0 ) THEN
304  info = -2
305  ELSE IF( lda.LT.max( 1, n ) ) THEN
306  info = -4
307  ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
308  info = -8
309  END IF
310 *
311  IF( info.EQ.0 ) THEN
312 *
313 * Determine the block size
314 *
315  nb = ilaenv( 1, 'SSYTRF_RK', uplo, n, -1, -1, -1 )
316  lwkopt = n*nb
317  work( 1 ) = lwkopt
318  END IF
319 *
320  IF( info.NE.0 ) THEN
321  CALL xerbla( 'SSYTRF_RK', -info )
322  RETURN
323  ELSE IF( lquery ) THEN
324  RETURN
325  END IF
326 *
327  nbmin = 2
328  ldwork = n
329  IF( nb.GT.1 .AND. nb.LT.n ) THEN
330  iws = ldwork*nb
331  IF( lwork.LT.iws ) THEN
332  nb = max( lwork / ldwork, 1 )
333  nbmin = max( 2, ilaenv( 2, 'SSYTRF_RK',
334  $ uplo, n, -1, -1, -1 ) )
335  END IF
336  ELSE
337  iws = 1
338  END IF
339  IF( nb.LT.nbmin )
340  $ nb = n
341 *
342  IF( upper ) THEN
343 *
344 * Factorize A as U*D*U**T using the upper triangle of A
345 *
346 * K is the main loop index, decreasing from N to 1 in steps of
347 * KB, where KB is the number of columns factorized by SLASYF_RK;
348 * KB is either NB or NB-1, or K for the last block
349 *
350  k = n
351  10 CONTINUE
352 *
353 * If K < 1, exit from loop
354 *
355  IF( k.LT.1 )
356  $ GO TO 15
357 *
358  IF( k.GT.nb ) THEN
359 *
360 * Factorize columns k-kb+1:k of A and use blocked code to
361 * update columns 1:k-kb
362 *
363  CALL slasyf_rk( uplo, k, nb, kb, a, lda, e,
364  $ ipiv, work, ldwork, iinfo )
365  ELSE
366 *
367 * Use unblocked code to factorize columns 1:k of A
368 *
369  CALL ssytf2_rk( uplo, k, a, lda, e, ipiv, iinfo )
370  kb = k
371  END IF
372 *
373 * Set INFO on the first occurrence of a zero pivot
374 *
375  IF( info.EQ.0 .AND. iinfo.GT.0 )
376  $ info = iinfo
377 *
378 * No need to adjust IPIV
379 *
380 *
381 * Apply permutations to the leading panel 1:k-1
382 *
383 * Read IPIV from the last block factored, i.e.
384 * indices k-kb+1:k and apply row permutations to the
385 * last k+1 colunms k+1:N after that block
386 * (We can do the simple loop over IPIV with decrement -1,
387 * since the ABS value of IPIV( I ) represents the row index
388 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
389 *
390  IF( k.LT.n ) THEN
391  DO i = k, ( k - kb + 1 ), -1
392  ip = abs( ipiv( i ) )
393  IF( ip.NE.i ) THEN
394  CALL sswap( n-k, a( i, k+1 ), lda,
395  $ a( ip, k+1 ), lda )
396  END IF
397  END DO
398  END IF
399 *
400 * Decrease K and return to the start of the main loop
401 *
402  k = k - kb
403  GO TO 10
404 *
405 * This label is the exit from main loop over K decreasing
406 * from N to 1 in steps of KB
407 *
408  15 CONTINUE
409 *
410  ELSE
411 *
412 * Factorize A as L*D*L**T using the lower triangle of A
413 *
414 * K is the main loop index, increasing from 1 to N in steps of
415 * KB, where KB is the number of columns factorized by SLASYF_RK;
416 * KB is either NB or NB-1, or N-K+1 for the last block
417 *
418  k = 1
419  20 CONTINUE
420 *
421 * If K > N, exit from loop
422 *
423  IF( k.GT.n )
424  $ GO TO 35
425 *
426  IF( k.LE.n-nb ) THEN
427 *
428 * Factorize columns k:k+kb-1 of A and use blocked code to
429 * update columns k+kb:n
430 *
431  CALL slasyf_rk( uplo, n-k+1, nb, kb, a( k, k ), lda, e( k ),
432  $ ipiv( k ), work, ldwork, iinfo )
433 
434 
435  ELSE
436 *
437 * Use unblocked code to factorize columns k:n of A
438 *
439  CALL ssytf2_rk( uplo, n-k+1, a( k, k ), lda, e( k ),
440  $ ipiv( k ), iinfo )
441  kb = n - k + 1
442 *
443  END IF
444 *
445 * Set INFO on the first occurrence of a zero pivot
446 *
447  IF( info.EQ.0 .AND. iinfo.GT.0 )
448  $ info = iinfo + k - 1
449 *
450 * Adjust IPIV
451 *
452  DO i = k, k + kb - 1
453  IF( ipiv( i ).GT.0 ) THEN
454  ipiv( i ) = ipiv( i ) + k - 1
455  ELSE
456  ipiv( i ) = ipiv( i ) - k + 1
457  END IF
458  END DO
459 *
460 * Apply permutations to the leading panel 1:k-1
461 *
462 * Read IPIV from the last block factored, i.e.
463 * indices k:k+kb-1 and apply row permutations to the
464 * first k-1 colunms 1:k-1 before that block
465 * (We can do the simple loop over IPIV with increment 1,
466 * since the ABS value of IPIV( I ) represents the row index
467 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
468 *
469  IF( k.GT.1 ) THEN
470  DO i = k, ( k + kb - 1 ), 1
471  ip = abs( ipiv( i ) )
472  IF( ip.NE.i ) THEN
473  CALL sswap( k-1, a( i, 1 ), lda,
474  $ a( ip, 1 ), lda )
475  END IF
476  END DO
477  END IF
478 *
479 * Increase K and return to the start of the main loop
480 *
481  k = k + kb
482  GO TO 20
483 *
484 * This label is the exit from main loop over K increasing
485 * from 1 to N in steps of KB
486 *
487  35 CONTINUE
488 *
489 * End Lower
490 *
491  END IF
492 *
493  work( 1 ) = lwkopt
494  RETURN
495 *
496 * End of SSYTRF_RK
497 *
498  END
subroutine ssytf2_rk(UPLO, N, A, LDA, E, IPIV, INFO)
SSYTF2_RK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Ka...
Definition: ssytf2_rk.f:243
subroutine ssytrf_rk(UPLO, N, A, LDA, E, IPIV, WORK, LWORK, INFO)
SSYTRF_RK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Ka...
Definition: ssytrf_rk.f:261
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:84
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:264