LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
csytf2.f
Go to the documentation of this file.
1 *> \brief \b CSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting method (unblocked algorithm).
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CSYTF2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytf2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytf2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytf2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CSYTF2( UPLO, N, A, LDA, IPIV, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER UPLO
25 * INTEGER INFO, LDA, N
26 * ..
27 * .. Array Arguments ..
28 * INTEGER IPIV( * )
29 * COMPLEX A( LDA, * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> CSYTF2 computes the factorization of a complex symmetric matrix A
39 *> using the Bunch-Kaufman diagonal pivoting method:
40 *>
41 *> A = U*D*U**T or A = L*D*L**T
42 *>
43 *> where U (or L) is a product of permutation and unit upper (lower)
44 *> triangular matrices, U**T is the transpose of U, and D is symmetric and
45 *> block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
46 *>
47 *> This is the unblocked version of the algorithm, calling Level 2 BLAS.
48 *> \endverbatim
49 *
50 * Arguments:
51 * ==========
52 *
53 *> \param[in] UPLO
54 *> \verbatim
55 *> UPLO is CHARACTER*1
56 *> Specifies whether the upper or lower triangular part of the
57 *> symmetric matrix A is stored:
58 *> = 'U': Upper triangular
59 *> = 'L': Lower triangular
60 *> \endverbatim
61 *>
62 *> \param[in] N
63 *> \verbatim
64 *> N is INTEGER
65 *> The order of the matrix A. N >= 0.
66 *> \endverbatim
67 *>
68 *> \param[in,out] A
69 *> \verbatim
70 *> A is COMPLEX array, dimension (LDA,N)
71 *> On entry, the symmetric matrix A. If UPLO = 'U', the leading
72 *> n-by-n upper triangular part of A contains the upper
73 *> triangular part of the matrix A, and the strictly lower
74 *> triangular part of A is not referenced. If UPLO = 'L', the
75 *> leading n-by-n lower triangular part of A contains the lower
76 *> triangular part of the matrix A, and the strictly upper
77 *> triangular part of A is not referenced.
78 *>
79 *> On exit, the block diagonal matrix D and the multipliers used
80 *> to obtain the factor U or L (see below for further details).
81 *> \endverbatim
82 *>
83 *> \param[in] LDA
84 *> \verbatim
85 *> LDA is INTEGER
86 *> The leading dimension of the array A. LDA >= max(1,N).
87 *> \endverbatim
88 *>
89 *> \param[out] IPIV
90 *> \verbatim
91 *> IPIV is INTEGER array, dimension (N)
92 *> Details of the interchanges and the block structure of D.
93 *>
94 *> If UPLO = 'U':
95 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
96 *> interchanged and D(k,k) is a 1-by-1 diagonal block.
97 *>
98 *> If IPIV(k) = IPIV(k-1) < 0, then rows and columns
99 *> k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
100 *> is a 2-by-2 diagonal block.
101 *>
102 *> If UPLO = 'L':
103 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
104 *> interchanged and D(k,k) is a 1-by-1 diagonal block.
105 *>
106 *> If IPIV(k) = IPIV(k+1) < 0, then rows and columns
107 *> k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
108 *> is a 2-by-2 diagonal block.
109 *> \endverbatim
110 *>
111 *> \param[out] INFO
112 *> \verbatim
113 *> INFO is INTEGER
114 *> = 0: successful exit
115 *> < 0: if INFO = -k, the k-th argument had an illegal value
116 *> > 0: if INFO = k, D(k,k) is exactly zero. The factorization
117 *> has been completed, but the block diagonal matrix D is
118 *> exactly singular, and division by zero will occur if it
119 *> is used to solve a system of equations.
120 *> \endverbatim
121 *
122 * Authors:
123 * ========
124 *
125 *> \author Univ. of Tennessee
126 *> \author Univ. of California Berkeley
127 *> \author Univ. of Colorado Denver
128 *> \author NAG Ltd.
129 *
130 *> \ingroup complexSYcomputational
131 *
132 *> \par Further Details:
133 * =====================
134 *>
135 *> \verbatim
136 *>
137 *> If UPLO = 'U', then A = U*D*U**T, where
138 *> U = P(n)*U(n)* ... *P(k)U(k)* ...,
139 *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to
140 *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
141 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
142 *> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
143 *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
144 *>
145 *> ( I v 0 ) k-s
146 *> U(k) = ( 0 I 0 ) s
147 *> ( 0 0 I ) n-k
148 *> k-s s n-k
149 *>
150 *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
151 *> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
152 *> and A(k,k), and v overwrites A(1:k-2,k-1:k).
153 *>
154 *> If UPLO = 'L', then A = L*D*L**T, where
155 *> L = P(1)*L(1)* ... *P(k)*L(k)* ...,
156 *> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
157 *> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
158 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
159 *> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
160 *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
161 *>
162 *> ( I 0 0 ) k-1
163 *> L(k) = ( 0 I 0 ) s
164 *> ( 0 v I ) n-k-s+1
165 *> k-1 s n-k-s+1
166 *>
167 *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
168 *> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
169 *> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
170 *> \endverbatim
171 *
172 *> \par Contributors:
173 * ==================
174 *>
175 *> \verbatim
176 *>
177 *> 09-29-06 - patch from
178 *> Bobby Cheng, MathWorks
179 *>
180 *> Replace l.209 and l.377
181 *> IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
182 *> by
183 *> IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN
184 *>
185 *> 1-96 - Based on modifications by J. Lewis, Boeing Computer Services
186 *> Company
187 *> \endverbatim
188 *
189 * =====================================================================
190  SUBROUTINE csytf2( UPLO, N, A, LDA, IPIV, INFO )
191 *
192 * -- LAPACK computational routine --
193 * -- LAPACK is a software package provided by Univ. of Tennessee, --
194 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
195 *
196 * .. Scalar Arguments ..
197  CHARACTER UPLO
198  INTEGER INFO, LDA, N
199 * ..
200 * .. Array Arguments ..
201  INTEGER IPIV( * )
202  COMPLEX A( LDA, * )
203 * ..
204 *
205 * =====================================================================
206 *
207 * .. Parameters ..
208  REAL ZERO, ONE
209  parameter( zero = 0.0e+0, one = 1.0e+0 )
210  REAL EIGHT, SEVTEN
211  parameter( eight = 8.0e+0, sevten = 17.0e+0 )
212  COMPLEX CONE
213  parameter( cone = ( 1.0e+0, 0.0e+0 ) )
214 * ..
215 * .. Local Scalars ..
216  LOGICAL UPPER
217  INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP
218  REAL ABSAKK, ALPHA, COLMAX, ROWMAX
219  COMPLEX D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, Z
220 * ..
221 * .. External Functions ..
222  LOGICAL LSAME, SISNAN
223  INTEGER ICAMAX
224  EXTERNAL lsame, icamax, sisnan
225 * ..
226 * .. External Subroutines ..
227  EXTERNAL cscal, cswap, csyr, xerbla
228 * ..
229 * .. Intrinsic Functions ..
230  INTRINSIC abs, aimag, max, real, sqrt
231 * ..
232 * .. Statement Functions ..
233  REAL CABS1
234 * ..
235 * .. Statement Function definitions ..
236  cabs1( z ) = abs( real( z ) ) + abs( aimag( z ) )
237 * ..
238 * .. Executable Statements ..
239 *
240 * Test the input parameters.
241 *
242  info = 0
243  upper = lsame( uplo, 'U' )
244  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
245  info = -1
246  ELSE IF( n.LT.0 ) THEN
247  info = -2
248  ELSE IF( lda.LT.max( 1, n ) ) THEN
249  info = -4
250  END IF
251  IF( info.NE.0 ) THEN
252  CALL xerbla( 'CSYTF2', -info )
253  RETURN
254  END IF
255 *
256 * Initialize ALPHA for use in choosing pivot block size.
257 *
258  alpha = ( one+sqrt( sevten ) ) / eight
259 *
260  IF( upper ) THEN
261 *
262 * Factorize A as U*D*U**T using the upper triangle of A
263 *
264 * K is the main loop index, decreasing from N to 1 in steps of
265 * 1 or 2
266 *
267  k = n
268  10 CONTINUE
269 *
270 * If K < 1, exit from loop
271 *
272  IF( k.LT.1 )
273  $ GO TO 70
274  kstep = 1
275 *
276 * Determine rows and columns to be interchanged and whether
277 * a 1-by-1 or 2-by-2 pivot block will be used
278 *
279  absakk = cabs1( a( k, k ) )
280 *
281 * IMAX is the row-index of the largest off-diagonal element in
282 * column K, and COLMAX is its absolute value.
283 * Determine both COLMAX and IMAX.
284 *
285  IF( k.GT.1 ) THEN
286  imax = icamax( k-1, a( 1, k ), 1 )
287  colmax = cabs1( a( imax, k ) )
288  ELSE
289  colmax = zero
290  END IF
291 *
292  IF( max( absakk, colmax ).EQ.zero .OR. sisnan(absakk) ) THEN
293 *
294 * Column K is zero or underflow, or contains a NaN:
295 * set INFO and continue
296 *
297  IF( info.EQ.0 )
298  $ info = k
299  kp = k
300  ELSE
301  IF( absakk.GE.alpha*colmax ) THEN
302 *
303 * no interchange, use 1-by-1 pivot block
304 *
305  kp = k
306  ELSE
307 *
308 * JMAX is the column-index of the largest off-diagonal
309 * element in row IMAX, and ROWMAX is its absolute value
310 *
311  jmax = imax + icamax( k-imax, a( imax, imax+1 ), lda )
312  rowmax = cabs1( a( imax, jmax ) )
313  IF( imax.GT.1 ) THEN
314  jmax = icamax( imax-1, a( 1, imax ), 1 )
315  rowmax = max( rowmax, cabs1( a( jmax, imax ) ) )
316  END IF
317 *
318  IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
319 *
320 * no interchange, use 1-by-1 pivot block
321 *
322  kp = k
323  ELSE IF( cabs1( a( imax, imax ) ).GE.alpha*rowmax ) THEN
324 *
325 * interchange rows and columns K and IMAX, use 1-by-1
326 * pivot block
327 *
328  kp = imax
329  ELSE
330 *
331 * interchange rows and columns K-1 and IMAX, use 2-by-2
332 * pivot block
333 *
334  kp = imax
335  kstep = 2
336  END IF
337  END IF
338 *
339  kk = k - kstep + 1
340  IF( kp.NE.kk ) THEN
341 *
342 * Interchange rows and columns KK and KP in the leading
343 * submatrix A(1:k,1:k)
344 *
345  CALL cswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
346  CALL cswap( kk-kp-1, a( kp+1, kk ), 1, a( kp, kp+1 ),
347  $ lda )
348  t = a( kk, kk )
349  a( kk, kk ) = a( kp, kp )
350  a( kp, kp ) = t
351  IF( kstep.EQ.2 ) THEN
352  t = a( k-1, k )
353  a( k-1, k ) = a( kp, k )
354  a( kp, k ) = t
355  END IF
356  END IF
357 *
358 * Update the leading submatrix
359 *
360  IF( kstep.EQ.1 ) THEN
361 *
362 * 1-by-1 pivot block D(k): column k now holds
363 *
364 * W(k) = U(k)*D(k)
365 *
366 * where U(k) is the k-th column of U
367 *
368 * Perform a rank-1 update of A(1:k-1,1:k-1) as
369 *
370 * A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T
371 *
372  r1 = cone / a( k, k )
373  CALL csyr( uplo, k-1, -r1, a( 1, k ), 1, a, lda )
374 *
375 * Store U(k) in column k
376 *
377  CALL cscal( k-1, r1, a( 1, k ), 1 )
378  ELSE
379 *
380 * 2-by-2 pivot block D(k): columns k and k-1 now hold
381 *
382 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
383 *
384 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
385 * of U
386 *
387 * Perform a rank-2 update of A(1:k-2,1:k-2) as
388 *
389 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
390 * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T
391 *
392  IF( k.GT.2 ) THEN
393 *
394  d12 = a( k-1, k )
395  d22 = a( k-1, k-1 ) / d12
396  d11 = a( k, k ) / d12
397  t = cone / ( d11*d22-cone )
398  d12 = t / d12
399 *
400  DO 30 j = k - 2, 1, -1
401  wkm1 = d12*( d11*a( j, k-1 )-a( j, k ) )
402  wk = d12*( d22*a( j, k )-a( j, k-1 ) )
403  DO 20 i = j, 1, -1
404  a( i, j ) = a( i, j ) - a( i, k )*wk -
405  $ a( i, k-1 )*wkm1
406  20 CONTINUE
407  a( j, k ) = wk
408  a( j, k-1 ) = wkm1
409  30 CONTINUE
410 *
411  END IF
412 *
413  END IF
414  END IF
415 *
416 * Store details of the interchanges in IPIV
417 *
418  IF( kstep.EQ.1 ) THEN
419  ipiv( k ) = kp
420  ELSE
421  ipiv( k ) = -kp
422  ipiv( k-1 ) = -kp
423  END IF
424 *
425 * Decrease K and return to the start of the main loop
426 *
427  k = k - kstep
428  GO TO 10
429 *
430  ELSE
431 *
432 * Factorize A as L*D*L**T using the lower triangle of A
433 *
434 * K is the main loop index, increasing from 1 to N in steps of
435 * 1 or 2
436 *
437  k = 1
438  40 CONTINUE
439 *
440 * If K > N, exit from loop
441 *
442  IF( k.GT.n )
443  $ GO TO 70
444  kstep = 1
445 *
446 * Determine rows and columns to be interchanged and whether
447 * a 1-by-1 or 2-by-2 pivot block will be used
448 *
449  absakk = cabs1( a( k, k ) )
450 *
451 * IMAX is the row-index of the largest off-diagonal element in
452 * column K, and COLMAX is its absolute value.
453 * Determine both COLMAX and IMAX.
454 *
455  IF( k.LT.n ) THEN
456  imax = k + icamax( n-k, a( k+1, k ), 1 )
457  colmax = cabs1( a( imax, k ) )
458  ELSE
459  colmax = zero
460  END IF
461 *
462  IF( max( absakk, colmax ).EQ.zero .OR. sisnan(absakk) ) THEN
463 *
464 * Column K is zero or underflow, or contains a NaN:
465 * set INFO and continue
466 *
467  IF( info.EQ.0 )
468  $ info = k
469  kp = k
470  ELSE
471  IF( absakk.GE.alpha*colmax ) THEN
472 *
473 * no interchange, use 1-by-1 pivot block
474 *
475  kp = k
476  ELSE
477 *
478 * JMAX is the column-index of the largest off-diagonal
479 * element in row IMAX, and ROWMAX is its absolute value
480 *
481  jmax = k - 1 + icamax( imax-k, a( imax, k ), lda )
482  rowmax = cabs1( a( imax, jmax ) )
483  IF( imax.LT.n ) THEN
484  jmax = imax + icamax( n-imax, a( imax+1, imax ), 1 )
485  rowmax = max( rowmax, cabs1( a( jmax, imax ) ) )
486  END IF
487 *
488  IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
489 *
490 * no interchange, use 1-by-1 pivot block
491 *
492  kp = k
493  ELSE IF( cabs1( a( imax, imax ) ).GE.alpha*rowmax ) THEN
494 *
495 * interchange rows and columns K and IMAX, use 1-by-1
496 * pivot block
497 *
498  kp = imax
499  ELSE
500 *
501 * interchange rows and columns K+1 and IMAX, use 2-by-2
502 * pivot block
503 *
504  kp = imax
505  kstep = 2
506  END IF
507  END IF
508 *
509  kk = k + kstep - 1
510  IF( kp.NE.kk ) THEN
511 *
512 * Interchange rows and columns KK and KP in the trailing
513 * submatrix A(k:n,k:n)
514 *
515  IF( kp.LT.n )
516  $ CALL cswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
517  CALL cswap( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
518  $ lda )
519  t = a( kk, kk )
520  a( kk, kk ) = a( kp, kp )
521  a( kp, kp ) = t
522  IF( kstep.EQ.2 ) THEN
523  t = a( k+1, k )
524  a( k+1, k ) = a( kp, k )
525  a( kp, k ) = t
526  END IF
527  END IF
528 *
529 * Update the trailing submatrix
530 *
531  IF( kstep.EQ.1 ) THEN
532 *
533 * 1-by-1 pivot block D(k): column k now holds
534 *
535 * W(k) = L(k)*D(k)
536 *
537 * where L(k) is the k-th column of L
538 *
539  IF( k.LT.n ) THEN
540 *
541 * Perform a rank-1 update of A(k+1:n,k+1:n) as
542 *
543 * A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T
544 *
545  r1 = cone / a( k, k )
546  CALL csyr( uplo, n-k, -r1, a( k+1, k ), 1,
547  $ a( k+1, k+1 ), lda )
548 *
549 * Store L(k) in column K
550 *
551  CALL cscal( n-k, r1, a( k+1, k ), 1 )
552  END IF
553  ELSE
554 *
555 * 2-by-2 pivot block D(k)
556 *
557  IF( k.LT.n-1 ) THEN
558 *
559 * Perform a rank-2 update of A(k+2:n,k+2:n) as
560 *
561 * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**T
562 * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**T
563 *
564 * where L(k) and L(k+1) are the k-th and (k+1)-th
565 * columns of L
566 *
567  d21 = a( k+1, k )
568  d11 = a( k+1, k+1 ) / d21
569  d22 = a( k, k ) / d21
570  t = cone / ( d11*d22-cone )
571  d21 = t / d21
572 *
573  DO 60 j = k + 2, n
574  wk = d21*( d11*a( j, k )-a( j, k+1 ) )
575  wkp1 = d21*( d22*a( j, k+1 )-a( j, k ) )
576  DO 50 i = j, n
577  a( i, j ) = a( i, j ) - a( i, k )*wk -
578  $ a( i, k+1 )*wkp1
579  50 CONTINUE
580  a( j, k ) = wk
581  a( j, k+1 ) = wkp1
582  60 CONTINUE
583  END IF
584  END IF
585  END IF
586 *
587 * Store details of the interchanges in IPIV
588 *
589  IF( kstep.EQ.1 ) THEN
590  ipiv( k ) = kp
591  ELSE
592  ipiv( k ) = -kp
593  ipiv( k+1 ) = -kp
594  END IF
595 *
596 * Increase K and return to the start of the main loop
597 *
598  k = k + kstep
599  GO TO 40
600 *
601  END IF
602 *
603  70 CONTINUE
604  RETURN
605 *
606 * End of CSYTF2
607 *
608  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine csyr(UPLO, N, ALPHA, X, INCX, A, LDA)
CSYR performs the symmetric rank-1 update of a complex symmetric matrix.
Definition: csyr.f:135
subroutine csytf2(UPLO, N, A, LDA, IPIV, INFO)
CSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting ...
Definition: csytf2.f:191