Lots of trailing whitespaces in the files of Syd. Cleaning this. No big deal.
[platform/upstream/lapack.git] / SRC / zsytf2.f
1 *> \brief \b ZSYTF2 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 ZSYTF2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytf2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytf2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytf2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 *  Definition:
19 *  ===========
20 *
21 *       SUBROUTINE ZSYTF2( 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*16         A( LDA, * )
30 *       ..
31 *
32 *
33 *> \par Purpose:
34 *  =============
35 *>
36 *> \verbatim
37 *>
38 *> ZSYTF2 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*16 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 complex16SYcomputational
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. DISNAN(ABSAKK) ) THEN
186 *>
187 *>  1-96 - Based on modifications by J. Lewis, Boeing Computer Services
188 *>         Company
189 *> \endverbatim
190 *
191 *  =====================================================================
192       SUBROUTINE ZSYTF2( 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*16         A( LDA, * )
206 *     ..
207 *
208 *  =====================================================================
209 *
210 *     .. Parameters ..
211       DOUBLE PRECISION   ZERO, ONE
212       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
213       DOUBLE PRECISION   EIGHT, SEVTEN
214       PARAMETER          ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
215       COMPLEX*16         CONE
216       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
217 *     ..
218 *     .. Local Scalars ..
219       LOGICAL            UPPER
220       INTEGER            I, IMAX, J, JMAX, K, KK, KP, KSTEP
221       DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, ROWMAX
222       COMPLEX*16         D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, Z
223 *     ..
224 *     .. External Functions ..
225       LOGICAL            DISNAN, LSAME
226       INTEGER            IZAMAX
227       EXTERNAL           DISNAN, LSAME, IZAMAX
228 *     ..
229 *     .. External Subroutines ..
230       EXTERNAL           XERBLA, ZSCAL, ZSWAP, ZSYR
231 *     ..
232 *     .. Intrinsic Functions ..
233       INTRINSIC          ABS, DBLE, DIMAG, MAX, SQRT
234 *     ..
235 *     .. Statement Functions ..
236       DOUBLE PRECISION   CABS1
237 *     ..
238 *     .. Statement Function definitions ..
239       CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( 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( 'ZSYTF2', -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 = IZAMAX( 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. DISNAN(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 + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
315                ROWMAX = CABS1( A( IMAX, JMAX ) )
316                IF( IMAX.GT.1 ) THEN
317                   JMAX = IZAMAX( 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 ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
349                CALL ZSWAP( 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 ZSYR( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
377 *
378 *              Store U(k) in column k
379 *
380                CALL ZSCAL( 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 + IZAMAX( 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. DISNAN(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 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
485                ROWMAX = CABS1( A( IMAX, JMAX ) )
486                IF( IMAX.LT.N ) THEN
487                   JMAX = IMAX + IZAMAX( 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 ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
520                CALL ZSWAP( 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 ZSYR( 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 ZSCAL( 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 ZSYTF2
610 *
611       END