LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
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 *> \date November 2013
131 *
132 *> \ingroup complexSYcomputational
133 *
134 *> \par Further Details:
135 * =====================
136 *>
137 *> \verbatim
138 *>
139 *> If UPLO = 'U', then A = U*D*U**T, where
140 *> U = P(n)*U(n)* ... *P(k)U(k)* ...,
141 *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to
142 *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
143 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
144 *> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
145 *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
146 *>
147 *> ( I v 0 ) k-s
148 *> U(k) = ( 0 I 0 ) s
149 *> ( 0 0 I ) n-k
150 *> k-s s n-k
151 *>
152 *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
153 *> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
154 *> and A(k,k), and v overwrites A(1:k-2,k-1:k).
155 *>
156 *> If UPLO = 'L', then A = L*D*L**T, where
157 *> L = P(1)*L(1)* ... *P(k)*L(k)* ...,
158 *> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
159 *> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
160 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
161 *> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
162 *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
163 *>
164 *> ( I 0 0 ) k-1
165 *> L(k) = ( 0 I 0 ) s
166 *> ( 0 v I ) n-k-s+1
167 *> k-1 s n-k-s+1
168 *>
169 *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
170 *> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
171 *> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
172 *> \endverbatim
173 *
174 *> \par Contributors:
175 * ==================
176 *>
177 *> \verbatim
178 *>
179 *> 09-29-06 - patch from
180 *> Bobby Cheng, MathWorks
181 *>
182 *> Replace l.209 and l.377
183 *> IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
184 *> by
185 *> IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN
186 *>
187 *> 1-96 - Based on modifications by J. Lewis, Boeing Computer Services
188 *> Company
189 *> \endverbatim
190 *
191 * =====================================================================
192  SUBROUTINE csytf2( UPLO, N, A, LDA, IPIV, INFO )
193 *
194 * -- LAPACK computational routine (version 3.5.0) --
195 * -- LAPACK is a software package provided by Univ. of Tennessee, --
196 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
197 * November 2013
198 *
199 * .. Scalar Arguments ..
200  CHARACTER uplo
201  INTEGER info, lda, n
202 * ..
203 * .. Array Arguments ..
204  INTEGER ipiv( * )
205  COMPLEX a( lda, * )
206 * ..
207 *
208 * =====================================================================
209 *
210 * .. Parameters ..
211  REAL zero, one
212  parameter( zero = 0.0e+0, one = 1.0e+0 )
213  REAL eight, sevten
214  parameter( eight = 8.0e+0, sevten = 17.0e+0 )
215  COMPLEX cone
216  parameter( cone = ( 1.0e+0, 0.0e+0 ) )
217 * ..
218 * .. Local Scalars ..
219  LOGICAL upper
220  INTEGER i, imax, j, jmax, k, kk, kp, kstep
221  REAL absakk, alpha, colmax, rowmax
222  COMPLEX d11, d12, d21, d22, r1, t, wk, wkm1, wkp1, z
223 * ..
224 * .. External Functions ..
225  LOGICAL lsame, sisnan
226  INTEGER icamax
227  EXTERNAL lsame, icamax, sisnan
228 * ..
229 * .. External Subroutines ..
230  EXTERNAL cscal, cswap, csyr, xerbla
231 * ..
232 * .. Intrinsic Functions ..
233  INTRINSIC abs, aimag, max, REAL, sqrt
234 * ..
235 * .. Statement Functions ..
236  REAL cabs1
237 * ..
238 * .. Statement Function definitions ..
239  cabs1( z ) = abs( REAL( Z ) ) + abs( aimag( z ) )
240 * ..
241 * .. Executable Statements ..
242 *
243 * Test the input parameters.
244 *
245  info = 0
246  upper = lsame( uplo, 'U' )
247  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
248  info = -1
249  ELSE IF( n.LT.0 ) THEN
250  info = -2
251  ELSE IF( lda.LT.max( 1, n ) ) THEN
252  info = -4
253  END IF
254  IF( info.NE.0 ) THEN
255  CALL xerbla( 'CSYTF2', -info )
256  RETURN
257  END IF
258 *
259 * Initialize ALPHA for use in choosing pivot block size.
260 *
261  alpha = ( one+sqrt( sevten ) ) / eight
262 *
263  IF( upper ) THEN
264 *
265 * Factorize A as U*D*U**T using the upper triangle of A
266 *
267 * K is the main loop index, decreasing from N to 1 in steps of
268 * 1 or 2
269 *
270  k = n
271  10 CONTINUE
272 *
273 * If K < 1, exit from loop
274 *
275  IF( k.LT.1 )
276  $ go to 70
277  kstep = 1
278 *
279 * Determine rows and columns to be interchanged and whether
280 * a 1-by-1 or 2-by-2 pivot block will be used
281 *
282  absakk = cabs1( a( k, k ) )
283 *
284 * IMAX is the row-index of the largest off-diagonal element in
285 * column K, and COLMAX is its absolute value.
286 * Determine both COLMAX and IMAX.
287 *
288  IF( k.GT.1 ) THEN
289  imax = icamax( k-1, a( 1, k ), 1 )
290  colmax = cabs1( a( imax, k ) )
291  ELSE
292  colmax = zero
293  END IF
294 *
295  IF( max( absakk, colmax ).EQ.zero .OR. sisnan(absakk) ) THEN
296 *
297 * Column K is zero or underflow, or contains a NaN:
298 * set INFO and continue
299 *
300  IF( info.EQ.0 )
301  $ info = k
302  kp = k
303  ELSE
304  IF( absakk.GE.alpha*colmax ) THEN
305 *
306 * no interchange, use 1-by-1 pivot block
307 *
308  kp = k
309  ELSE
310 *
311 * JMAX is the column-index of the largest off-diagonal
312 * element in row IMAX, and ROWMAX is its absolute value
313 *
314  jmax = imax + icamax( k-imax, a( imax, imax+1 ), lda )
315  rowmax = cabs1( a( imax, jmax ) )
316  IF( imax.GT.1 ) THEN
317  jmax = icamax( imax-1, a( 1, imax ), 1 )
318  rowmax = max( rowmax, cabs1( a( jmax, imax ) ) )
319  END IF
320 *
321  IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
322 *
323 * no interchange, use 1-by-1 pivot block
324 *
325  kp = k
326  ELSE IF( cabs1( a( imax, imax ) ).GE.alpha*rowmax ) THEN
327 *
328 * interchange rows and columns K and IMAX, use 1-by-1
329 * pivot block
330 *
331  kp = imax
332  ELSE
333 *
334 * interchange rows and columns K-1 and IMAX, use 2-by-2
335 * pivot block
336 *
337  kp = imax
338  kstep = 2
339  END IF
340  END IF
341 *
342  kk = k - kstep + 1
343  IF( kp.NE.kk ) THEN
344 *
345 * Interchange rows and columns KK and KP in the leading
346 * submatrix A(1:k,1:k)
347 *
348  CALL cswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
349  CALL cswap( kk-kp-1, a( kp+1, kk ), 1, a( kp, kp+1 ),
350  $ lda )
351  t = a( kk, kk )
352  a( kk, kk ) = a( kp, kp )
353  a( kp, kp ) = t
354  IF( kstep.EQ.2 ) THEN
355  t = a( k-1, k )
356  a( k-1, k ) = a( kp, k )
357  a( kp, k ) = t
358  END IF
359  END IF
360 *
361 * Update the leading submatrix
362 *
363  IF( kstep.EQ.1 ) THEN
364 *
365 * 1-by-1 pivot block D(k): column k now holds
366 *
367 * W(k) = U(k)*D(k)
368 *
369 * where U(k) is the k-th column of U
370 *
371 * Perform a rank-1 update of A(1:k-1,1:k-1) as
372 *
373 * A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T
374 *
375  r1 = cone / a( k, k )
376  CALL csyr( uplo, k-1, -r1, a( 1, k ), 1, a, lda )
377 *
378 * Store U(k) in column k
379 *
380  CALL cscal( k-1, r1, a( 1, k ), 1 )
381  ELSE
382 *
383 * 2-by-2 pivot block D(k): columns k and k-1 now hold
384 *
385 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
386 *
387 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
388 * of U
389 *
390 * Perform a rank-2 update of A(1:k-2,1:k-2) as
391 *
392 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
393 * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T
394 *
395  IF( k.GT.2 ) THEN
396 *
397  d12 = a( k-1, k )
398  d22 = a( k-1, k-1 ) / d12
399  d11 = a( k, k ) / d12
400  t = cone / ( d11*d22-cone )
401  d12 = t / d12
402 *
403  DO 30 j = k - 2, 1, -1
404  wkm1 = d12*( d11*a( j, k-1 )-a( j, k ) )
405  wk = d12*( d22*a( j, k )-a( j, k-1 ) )
406  DO 20 i = j, 1, -1
407  a( i, j ) = a( i, j ) - a( i, k )*wk -
408  $ a( i, k-1 )*wkm1
409  20 CONTINUE
410  a( j, k ) = wk
411  a( j, k-1 ) = wkm1
412  30 CONTINUE
413 *
414  END IF
415 *
416  END IF
417  END IF
418 *
419 * Store details of the interchanges in IPIV
420 *
421  IF( kstep.EQ.1 ) THEN
422  ipiv( k ) = kp
423  ELSE
424  ipiv( k ) = -kp
425  ipiv( k-1 ) = -kp
426  END IF
427 *
428 * Decrease K and return to the start of the main loop
429 *
430  k = k - kstep
431  go to 10
432 *
433  ELSE
434 *
435 * Factorize A as L*D*L**T using the lower triangle of A
436 *
437 * K is the main loop index, increasing from 1 to N in steps of
438 * 1 or 2
439 *
440  k = 1
441  40 CONTINUE
442 *
443 * If K > N, exit from loop
444 *
445  IF( k.GT.n )
446  $ go to 70
447  kstep = 1
448 *
449 * Determine rows and columns to be interchanged and whether
450 * a 1-by-1 or 2-by-2 pivot block will be used
451 *
452  absakk = cabs1( a( k, k ) )
453 *
454 * IMAX is the row-index of the largest off-diagonal element in
455 * column K, and COLMAX is its absolute value.
456 * Determine both COLMAX and IMAX.
457 *
458  IF( k.LT.n ) THEN
459  imax = k + icamax( n-k, a( k+1, k ), 1 )
460  colmax = cabs1( a( imax, k ) )
461  ELSE
462  colmax = zero
463  END IF
464 *
465  IF( max( absakk, colmax ).EQ.zero .OR. sisnan(absakk) ) THEN
466 *
467 * Column K is zero or underflow, or contains a NaN:
468 * set INFO and continue
469 *
470  IF( info.EQ.0 )
471  $ info = k
472  kp = k
473  ELSE
474  IF( absakk.GE.alpha*colmax ) THEN
475 *
476 * no interchange, use 1-by-1 pivot block
477 *
478  kp = k
479  ELSE
480 *
481 * JMAX is the column-index of the largest off-diagonal
482 * element in row IMAX, and ROWMAX is its absolute value
483 *
484  jmax = k - 1 + icamax( imax-k, a( imax, k ), lda )
485  rowmax = cabs1( a( imax, jmax ) )
486  IF( imax.LT.n ) THEN
487  jmax = imax + icamax( n-imax, a( imax+1, imax ), 1 )
488  rowmax = max( rowmax, cabs1( a( jmax, imax ) ) )
489  END IF
490 *
491  IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
492 *
493 * no interchange, use 1-by-1 pivot block
494 *
495  kp = k
496  ELSE IF( cabs1( a( imax, imax ) ).GE.alpha*rowmax ) THEN
497 *
498 * interchange rows and columns K and IMAX, use 1-by-1
499 * pivot block
500 *
501  kp = imax
502  ELSE
503 *
504 * interchange rows and columns K+1 and IMAX, use 2-by-2
505 * pivot block
506 *
507  kp = imax
508  kstep = 2
509  END IF
510  END IF
511 *
512  kk = k + kstep - 1
513  IF( kp.NE.kk ) THEN
514 *
515 * Interchange rows and columns KK and KP in the trailing
516 * submatrix A(k:n,k:n)
517 *
518  IF( kp.LT.n )
519  $ CALL cswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
520  CALL cswap( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
521  $ lda )
522  t = a( kk, kk )
523  a( kk, kk ) = a( kp, kp )
524  a( kp, kp ) = t
525  IF( kstep.EQ.2 ) THEN
526  t = a( k+1, k )
527  a( k+1, k ) = a( kp, k )
528  a( kp, k ) = t
529  END IF
530  END IF
531 *
532 * Update the trailing submatrix
533 *
534  IF( kstep.EQ.1 ) THEN
535 *
536 * 1-by-1 pivot block D(k): column k now holds
537 *
538 * W(k) = L(k)*D(k)
539 *
540 * where L(k) is the k-th column of L
541 *
542  IF( k.LT.n ) THEN
543 *
544 * Perform a rank-1 update of A(k+1:n,k+1:n) as
545 *
546 * A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T
547 *
548  r1 = cone / a( k, k )
549  CALL csyr( uplo, n-k, -r1, a( k+1, k ), 1,
550  $ a( k+1, k+1 ), lda )
551 *
552 * Store L(k) in column K
553 *
554  CALL cscal( n-k, r1, a( k+1, k ), 1 )
555  END IF
556  ELSE
557 *
558 * 2-by-2 pivot block D(k)
559 *
560  IF( k.LT.n-1 ) THEN
561 *
562 * Perform a rank-2 update of A(k+2:n,k+2:n) as
563 *
564 * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**T
565 * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**T
566 *
567 * where L(k) and L(k+1) are the k-th and (k+1)-th
568 * columns of L
569 *
570  d21 = a( k+1, k )
571  d11 = a( k+1, k+1 ) / d21
572  d22 = a( k, k ) / d21
573  t = cone / ( d11*d22-cone )
574  d21 = t / d21
575 *
576  DO 60 j = k + 2, n
577  wk = d21*( d11*a( j, k )-a( j, k+1 ) )
578  wkp1 = d21*( d22*a( j, k+1 )-a( j, k ) )
579  DO 50 i = j, n
580  a( i, j ) = a( i, j ) - a( i, k )*wk -
581  $ a( i, k+1 )*wkp1
582  50 CONTINUE
583  a( j, k ) = wk
584  a( j, k+1 ) = wkp1
585  60 CONTINUE
586  END IF
587  END IF
588  END IF
589 *
590 * Store details of the interchanges in IPIV
591 *
592  IF( kstep.EQ.1 ) THEN
593  ipiv( k ) = kp
594  ELSE
595  ipiv( k ) = -kp
596  ipiv( k+1 ) = -kp
597  END IF
598 *
599 * Increase K and return to the start of the main loop
600 *
601  k = k + kstep
602  go to 40
603 *
604  END IF
605 *
606  70 CONTINUE
607  RETURN
608 *
609 * End of CSYTF2
610 *
611  END