LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
chetf2.f
Go to the documentation of this file.
1 *> \brief \b CHETF2 computes the factorization of a complex Hermitian 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 CHETF2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetf2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetf2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetf2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CHETF2( 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 *> CHETF2 computes the factorization of a complex Hermitian matrix A
39 *> using the Bunch-Kaufman diagonal pivoting method:
40 *>
41 *> A = U*D*U**H or A = L*D*L**H
42 *>
43 *> where U (or L) is a product of permutation and unit upper (lower)
44 *> triangular matrices, U**H is the conjugate transpose of U, and D is
45 *> Hermitian and 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 *> Hermitian 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 Hermitian 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 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
94 *> interchanged and D(k,k) is a 1-by-1 diagonal block.
95 *> If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
96 *> columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
97 *> is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
98 *> IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
99 *> interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
100 *> \endverbatim
101 *>
102 *> \param[out] INFO
103 *> \verbatim
104 *> INFO is INTEGER
105 *> = 0: successful exit
106 *> < 0: if INFO = -k, the k-th argument had an illegal value
107 *> > 0: if INFO = k, D(k,k) is exactly zero. The factorization
108 *> has been completed, but the block diagonal matrix D is
109 *> exactly singular, and division by zero will occur if it
110 *> is used to solve a system of equations.
111 *> \endverbatim
112 *
113 * Authors:
114 * ========
115 *
116 *> \author Univ. of Tennessee
117 *> \author Univ. of California Berkeley
118 *> \author Univ. of Colorado Denver
119 *> \author NAG Ltd.
120 *
121 *> \date September 2012
122 *
123 *> \ingroup complexHEcomputational
124 *
125 *> \par Further Details:
126 * =====================
127 *>
128 *> \verbatim
129 *>
130 *> 09-29-06 - patch from
131 *> Bobby Cheng, MathWorks
132 *>
133 *> Replace l.210 and l.392
134 *> IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
135 *> by
136 *> IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN
137 *>
138 *> 01-01-96 - Based on modifications by
139 *> J. Lewis, Boeing Computer Services Company
140 *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
141 *>
142 *> If UPLO = 'U', then A = U*D*U**H, where
143 *> U = P(n)*U(n)* ... *P(k)U(k)* ...,
144 *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to
145 *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
146 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
147 *> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
148 *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
149 *>
150 *> ( I v 0 ) k-s
151 *> U(k) = ( 0 I 0 ) s
152 *> ( 0 0 I ) n-k
153 *> k-s s n-k
154 *>
155 *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
156 *> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
157 *> and A(k,k), and v overwrites A(1:k-2,k-1:k).
158 *>
159 *> If UPLO = 'L', then A = L*D*L**H, where
160 *> L = P(1)*L(1)* ... *P(k)*L(k)* ...,
161 *> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
162 *> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
163 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
164 *> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
165 *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
166 *>
167 *> ( I 0 0 ) k-1
168 *> L(k) = ( 0 I 0 ) s
169 *> ( 0 v I ) n-k-s+1
170 *> k-1 s n-k-s+1
171 *>
172 *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
173 *> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
174 *> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
175 *> \endverbatim
176 *>
177 * =====================================================================
178  SUBROUTINE chetf2( UPLO, N, A, LDA, IPIV, INFO )
179 *
180 * -- LAPACK computational routine (version 3.4.2) --
181 * -- LAPACK is a software package provided by Univ. of Tennessee, --
182 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
183 * September 2012
184 *
185 * .. Scalar Arguments ..
186  CHARACTER uplo
187  INTEGER info, lda, n
188 * ..
189 * .. Array Arguments ..
190  INTEGER ipiv( * )
191  COMPLEX a( lda, * )
192 * ..
193 *
194 * =====================================================================
195 *
196 * .. Parameters ..
197  REAL zero, one
198  parameter( zero = 0.0e+0, one = 1.0e+0 )
199  REAL eight, sevten
200  parameter( eight = 8.0e+0, sevten = 17.0e+0 )
201 * ..
202 * .. Local Scalars ..
203  LOGICAL upper
204  INTEGER i, imax, j, jmax, k, kk, kp, kstep
205  REAL absakk, alpha, colmax, d, d11, d22, r1, rowmax,
206  $ tt
207  COMPLEX d12, d21, t, wk, wkm1, wkp1, zdum
208 * ..
209 * .. External Functions ..
210  LOGICAL lsame, sisnan
211  INTEGER icamax
212  REAL slapy2
213  EXTERNAL lsame, icamax, slapy2, sisnan
214 * ..
215 * .. External Subroutines ..
216  EXTERNAL cher, csscal, cswap, xerbla
217 * ..
218 * .. Intrinsic Functions ..
219  INTRINSIC abs, aimag, cmplx, conjg, max, REAL, sqrt
220 * ..
221 * .. Statement Functions ..
222  REAL cabs1
223 * ..
224 * .. Statement Function definitions ..
225  cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( aimag( zdum ) )
226 * ..
227 * .. Executable Statements ..
228 *
229 * Test the input parameters.
230 *
231  info = 0
232  upper = lsame( uplo, 'U' )
233  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
234  info = -1
235  ELSE IF( n.LT.0 ) THEN
236  info = -2
237  ELSE IF( lda.LT.max( 1, n ) ) THEN
238  info = -4
239  END IF
240  IF( info.NE.0 ) THEN
241  CALL xerbla( 'CHETF2', -info )
242  return
243  END IF
244 *
245 * Initialize ALPHA for use in choosing pivot block size.
246 *
247  alpha = ( one+sqrt( sevten ) ) / eight
248 *
249  IF( upper ) THEN
250 *
251 * Factorize A as U*D*U**H using the upper triangle of A
252 *
253 * K is the main loop index, decreasing from N to 1 in steps of
254 * 1 or 2
255 *
256  k = n
257  10 continue
258 *
259 * If K < 1, exit from loop
260 *
261  IF( k.LT.1 )
262  $ go to 90
263  kstep = 1
264 *
265 * Determine rows and columns to be interchanged and whether
266 * a 1-by-1 or 2-by-2 pivot block will be used
267 *
268  absakk = abs( REAL( A( K, K ) ) )
269 *
270 * IMAX is the row-index of the largest off-diagonal element in
271 * column K, and COLMAX is its absolute value
272 *
273  IF( k.GT.1 ) THEN
274  imax = icamax( k-1, a( 1, k ), 1 )
275  colmax = cabs1( a( imax, k ) )
276  ELSE
277  colmax = zero
278  END IF
279 *
280  IF( (max( absakk, colmax ).EQ.zero) .OR. sisnan(absakk) ) THEN
281 *
282 * Column K is zero or contains a NaN: set INFO and continue
283 *
284  IF( info.EQ.0 )
285  $ info = k
286  kp = k
287  a( k, k ) = REAL( A( K, K ) )
288  ELSE
289  IF( absakk.GE.alpha*colmax ) THEN
290 *
291 * no interchange, use 1-by-1 pivot block
292 *
293  kp = k
294  ELSE
295 *
296 * JMAX is the column-index of the largest off-diagonal
297 * element in row IMAX, and ROWMAX is its absolute value
298 *
299  jmax = imax + icamax( k-imax, a( imax, imax+1 ), lda )
300  rowmax = cabs1( a( imax, jmax ) )
301  IF( imax.GT.1 ) THEN
302  jmax = icamax( imax-1, a( 1, imax ), 1 )
303  rowmax = max( rowmax, cabs1( a( jmax, imax ) ) )
304  END IF
305 *
306  IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
307 *
308 * no interchange, use 1-by-1 pivot block
309 *
310  kp = k
311  ELSE IF( abs( REAL( A( IMAX, IMAX ) ) ).GE.alpha*rowmax )
312  $ THEN
313 *
314 * interchange rows and columns K and IMAX, use 1-by-1
315 * pivot block
316 *
317  kp = imax
318  ELSE
319 *
320 * interchange rows and columns K-1 and IMAX, use 2-by-2
321 * pivot block
322 *
323  kp = imax
324  kstep = 2
325  END IF
326  END IF
327 *
328  kk = k - kstep + 1
329  IF( kp.NE.kk ) THEN
330 *
331 * Interchange rows and columns KK and KP in the leading
332 * submatrix A(1:k,1:k)
333 *
334  CALL cswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
335  DO 20 j = kp + 1, kk - 1
336  t = conjg( a( j, kk ) )
337  a( j, kk ) = conjg( a( kp, j ) )
338  a( kp, j ) = t
339  20 continue
340  a( kp, kk ) = conjg( a( kp, kk ) )
341  r1 = REAL( A( KK, KK ) )
342  a( kk, kk ) = REAL( A( KP, KP ) )
343  a( kp, kp ) = r1
344  IF( kstep.EQ.2 ) THEN
345  a( k, k ) = REAL( A( K, K ) )
346  t = a( k-1, k )
347  a( k-1, k ) = a( kp, k )
348  a( kp, k ) = t
349  END IF
350  ELSE
351  a( k, k ) = REAL( A( K, K ) )
352  IF( kstep.EQ.2 )
353  $ a( k-1, k-1 ) = REAL( A( K-1, K-1 ) )
354  END IF
355 *
356 * Update the leading submatrix
357 *
358  IF( kstep.EQ.1 ) THEN
359 *
360 * 1-by-1 pivot block D(k): column k now holds
361 *
362 * W(k) = U(k)*D(k)
363 *
364 * where U(k) is the k-th column of U
365 *
366 * Perform a rank-1 update of A(1:k-1,1:k-1) as
367 *
368 * A := A - U(k)*D(k)*U(k)**H = A - W(k)*1/D(k)*W(k)**H
369 *
370  r1 = one / REAL( A( K, K ) )
371  CALL cher( uplo, k-1, -r1, a( 1, k ), 1, a, lda )
372 *
373 * Store U(k) in column k
374 *
375  CALL csscal( k-1, r1, a( 1, k ), 1 )
376  ELSE
377 *
378 * 2-by-2 pivot block D(k): columns k and k-1 now hold
379 *
380 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
381 *
382 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
383 * of U
384 *
385 * Perform a rank-2 update of A(1:k-2,1:k-2) as
386 *
387 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**H
388 * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**H
389 *
390  IF( k.GT.2 ) THEN
391 *
392  d = slapy2( REAL( A( K-1, K ) ),
393  $ aimag( a( k-1, k ) ) )
394  d22 = REAL( A( K-1, K-1 ) ) / d
395  d11 = REAL( A( K, K ) ) / d
396  tt = one / ( d11*d22-one )
397  d12 = a( k-1, k ) / d
398  d = tt / d
399 *
400  DO 40 j = k - 2, 1, -1
401  wkm1 = d*( d11*a( j, k-1 )-conjg( d12 )*a( j, k ) )
402  wk = d*( d22*a( j, k )-d12*a( j, k-1 ) )
403  DO 30 i = j, 1, -1
404  a( i, j ) = a( i, j ) - a( i, k )*conjg( wk ) -
405  $ a( i, k-1 )*conjg( wkm1 )
406  30 continue
407  a( j, k ) = wk
408  a( j, k-1 ) = wkm1
409  a( j, j ) = cmplx( REAL( A( J, J ) ), 0.0e+0 )
410  40 continue
411 *
412  END IF
413 *
414  END IF
415  END IF
416 *
417 * Store details of the interchanges in IPIV
418 *
419  IF( kstep.EQ.1 ) THEN
420  ipiv( k ) = kp
421  ELSE
422  ipiv( k ) = -kp
423  ipiv( k-1 ) = -kp
424  END IF
425 *
426 * Decrease K and return to the start of the main loop
427 *
428  k = k - kstep
429  go to 10
430 *
431  ELSE
432 *
433 * Factorize A as L*D*L**H using the lower triangle of A
434 *
435 * K is the main loop index, increasing from 1 to N in steps of
436 * 1 or 2
437 *
438  k = 1
439  50 continue
440 *
441 * If K > N, exit from loop
442 *
443  IF( k.GT.n )
444  $ go to 90
445  kstep = 1
446 *
447 * Determine rows and columns to be interchanged and whether
448 * a 1-by-1 or 2-by-2 pivot block will be used
449 *
450  absakk = abs( REAL( A( K, K ) ) )
451 *
452 * IMAX is the row-index of the largest off-diagonal element in
453 * column K, and COLMAX is its absolute value
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 contains a NaN: set INFO and continue
465 *
466  IF( info.EQ.0 )
467  $ info = k
468  kp = k
469  a( k, k ) = REAL( A( K, 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( abs( REAL( A( IMAX, IMAX ) ) ).GE.alpha*rowmax )
494  $ THEN
495 *
496 * interchange rows and columns K and IMAX, use 1-by-1
497 * pivot block
498 *
499  kp = imax
500  ELSE
501 *
502 * interchange rows and columns K+1 and IMAX, use 2-by-2
503 * pivot block
504 *
505  kp = imax
506  kstep = 2
507  END IF
508  END IF
509 *
510  kk = k + kstep - 1
511  IF( kp.NE.kk ) THEN
512 *
513 * Interchange rows and columns KK and KP in the trailing
514 * submatrix A(k:n,k:n)
515 *
516  IF( kp.LT.n )
517  $ CALL cswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
518  DO 60 j = kk + 1, kp - 1
519  t = conjg( a( j, kk ) )
520  a( j, kk ) = conjg( a( kp, j ) )
521  a( kp, j ) = t
522  60 continue
523  a( kp, kk ) = conjg( a( kp, kk ) )
524  r1 = REAL( A( KK, KK ) )
525  a( kk, kk ) = REAL( A( KP, KP ) )
526  a( kp, kp ) = r1
527  IF( kstep.EQ.2 ) THEN
528  a( k, k ) = REAL( A( K, K ) )
529  t = a( k+1, k )
530  a( k+1, k ) = a( kp, k )
531  a( kp, k ) = t
532  END IF
533  ELSE
534  a( k, k ) = REAL( A( K, K ) )
535  IF( kstep.EQ.2 )
536  $ a( k+1, k+1 ) = REAL( A( K+1, K+1 ) )
537  END IF
538 *
539 * Update the trailing submatrix
540 *
541  IF( kstep.EQ.1 ) THEN
542 *
543 * 1-by-1 pivot block D(k): column k now holds
544 *
545 * W(k) = L(k)*D(k)
546 *
547 * where L(k) is the k-th column of L
548 *
549  IF( k.LT.n ) THEN
550 *
551 * Perform a rank-1 update of A(k+1:n,k+1:n) as
552 *
553 * A := A - L(k)*D(k)*L(k)**H = A - W(k)*(1/D(k))*W(k)**H
554 *
555  r1 = one / REAL( A( K, K ) )
556  CALL cher( uplo, n-k, -r1, a( k+1, k ), 1,
557  $ a( k+1, k+1 ), lda )
558 *
559 * Store L(k) in column K
560 *
561  CALL csscal( n-k, r1, a( k+1, k ), 1 )
562  END IF
563  ELSE
564 *
565 * 2-by-2 pivot block D(k)
566 *
567  IF( k.LT.n-1 ) THEN
568 *
569 * Perform a rank-2 update of A(k+2:n,k+2:n) as
570 *
571 * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**H
572 * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**H
573 *
574 * where L(k) and L(k+1) are the k-th and (k+1)-th
575 * columns of L
576 *
577  d = slapy2( REAL( A( K+1, K ) ),
578  $ aimag( a( k+1, k ) ) )
579  d11 = REAL( A( K+1, K+1 ) ) / d
580  d22 = REAL( A( K, K ) ) / d
581  tt = one / ( d11*d22-one )
582  d21 = a( k+1, k ) / d
583  d = tt / d
584 *
585  DO 80 j = k + 2, n
586  wk = d*( d11*a( j, k )-d21*a( j, k+1 ) )
587  wkp1 = d*( d22*a( j, k+1 )-conjg( d21 )*a( j, k ) )
588  DO 70 i = j, n
589  a( i, j ) = a( i, j ) - a( i, k )*conjg( wk ) -
590  $ a( i, k+1 )*conjg( wkp1 )
591  70 continue
592  a( j, k ) = wk
593  a( j, k+1 ) = wkp1
594  a( j, j ) = cmplx( REAL( A( J, J ) ), 0.0e+0 )
595  80 continue
596  END IF
597  END IF
598  END IF
599 *
600 * Store details of the interchanges in IPIV
601 *
602  IF( kstep.EQ.1 ) THEN
603  ipiv( k ) = kp
604  ELSE
605  ipiv( k ) = -kp
606  ipiv( k+1 ) = -kp
607  END IF
608 *
609 * Increase K and return to the start of the main loop
610 *
611  k = k + kstep
612  go to 50
613 *
614  END IF
615 *
616  90 continue
617  return
618 *
619 * End of CHETF2
620 *
621  END