LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ chptrf()

subroutine chptrf ( character  UPLO,
integer  N,
complex, dimension( * )  AP,
integer, dimension( * )  IPIV,
integer  INFO 
)

CHPTRF

Download CHPTRF + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 CHPTRF computes the factorization of a complex Hermitian packed
 matrix A using the Bunch-Kaufman diagonal pivoting method:

    A = U*D*U**H  or  A = L*D*L**H

 where U (or L) is a product of permutation and unit upper (lower)
 triangular matrices, and D is Hermitian and block diagonal with
 1-by-1 and 2-by-2 diagonal blocks.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]AP
          AP is COMPLEX array, dimension (N*(N+1)/2)
          On entry, the upper or lower triangle of the Hermitian matrix
          A, packed columnwise in a linear array.  The j-th column of A
          is stored in the array AP as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

          On exit, the block diagonal matrix D and the multipliers used
          to obtain the factor U or L, stored as a packed triangular
          matrix overwriting A (see below for further details).
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D.
          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
          interchanged and D(k,k) is a 1-by-1 diagonal block.
          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          > 0: if INFO = i, D(i,i) is exactly zero.  The factorization
               has been completed, but the block diagonal matrix D is
               exactly singular, and division by zero will occur if it
               is used to solve a system of equations.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  If UPLO = 'U', then A = U*D*U**H, where
     U = P(n)*U(n)* ... *P(k)U(k)* ...,
  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

             (   I    v    0   )   k-s
     U(k) =  (   0    I    0   )   s
             (   0    0    I   )   n-k
                k-s   s   n-k

  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
  and A(k,k), and v overwrites A(1:k-2,k-1:k).

  If UPLO = 'L', then A = L*D*L**H, where
     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

             (   I    0     0   )  k-1
     L(k) =  (   0    I     0   )  s
             (   0    v     I   )  n-k-s+1
                k-1   s  n-k-s+1

  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
Contributors:
J. Lewis, Boeing Computer Services Company

Definition at line 158 of file chptrf.f.

159 *
160 * -- LAPACK computational routine --
161 * -- LAPACK is a software package provided by Univ. of Tennessee, --
162 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
163 *
164 * .. Scalar Arguments ..
165  CHARACTER UPLO
166  INTEGER INFO, N
167 * ..
168 * .. Array Arguments ..
169  INTEGER IPIV( * )
170  COMPLEX AP( * )
171 * ..
172 *
173 * =====================================================================
174 *
175 * .. Parameters ..
176  REAL ZERO, ONE
177  parameter( zero = 0.0e+0, one = 1.0e+0 )
178  REAL EIGHT, SEVTEN
179  parameter( eight = 8.0e+0, sevten = 17.0e+0 )
180 * ..
181 * .. Local Scalars ..
182  LOGICAL UPPER
183  INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
184  $ KSTEP, KX, NPP
185  REAL ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
186  $ TT
187  COMPLEX D12, D21, T, WK, WKM1, WKP1, ZDUM
188 * ..
189 * .. External Functions ..
190  LOGICAL LSAME
191  INTEGER ICAMAX
192  REAL SLAPY2
193  EXTERNAL lsame, icamax, slapy2
194 * ..
195 * .. External Subroutines ..
196  EXTERNAL chpr, csscal, cswap, xerbla
197 * ..
198 * .. Intrinsic Functions ..
199  INTRINSIC abs, aimag, cmplx, conjg, max, real, sqrt
200 * ..
201 * .. Statement Functions ..
202  REAL CABS1
203 * ..
204 * .. Statement Function definitions ..
205  cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
206 * ..
207 * .. Executable Statements ..
208 *
209 * Test the input parameters.
210 *
211  info = 0
212  upper = lsame( uplo, 'U' )
213  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
214  info = -1
215  ELSE IF( n.LT.0 ) THEN
216  info = -2
217  END IF
218  IF( info.NE.0 ) THEN
219  CALL xerbla( 'CHPTRF', -info )
220  RETURN
221  END IF
222 *
223 * Initialize ALPHA for use in choosing pivot block size.
224 *
225  alpha = ( one+sqrt( sevten ) ) / eight
226 *
227  IF( upper ) THEN
228 *
229 * Factorize A as U*D*U**H using the upper triangle of A
230 *
231 * K is the main loop index, decreasing from N to 1 in steps of
232 * 1 or 2
233 *
234  k = n
235  kc = ( n-1 )*n / 2 + 1
236  10 CONTINUE
237  knc = kc
238 *
239 * If K < 1, exit from loop
240 *
241  IF( k.LT.1 )
242  $ GO TO 110
243  kstep = 1
244 *
245 * Determine rows and columns to be interchanged and whether
246 * a 1-by-1 or 2-by-2 pivot block will be used
247 *
248  absakk = abs( real( ap( kc+k-1 ) ) )
249 *
250 * IMAX is the row-index of the largest off-diagonal element in
251 * column K, and COLMAX is its absolute value
252 *
253  IF( k.GT.1 ) THEN
254  imax = icamax( k-1, ap( kc ), 1 )
255  colmax = cabs1( ap( kc+imax-1 ) )
256  ELSE
257  colmax = zero
258  END IF
259 *
260  IF( max( absakk, colmax ).EQ.zero ) THEN
261 *
262 * Column K is zero: set INFO and continue
263 *
264  IF( info.EQ.0 )
265  $ info = k
266  kp = k
267  ap( kc+k-1 ) = real( ap( kc+k-1 ) )
268  ELSE
269  IF( absakk.GE.alpha*colmax ) THEN
270 *
271 * no interchange, use 1-by-1 pivot block
272 *
273  kp = k
274  ELSE
275 *
276 * JMAX is the column-index of the largest off-diagonal
277 * element in row IMAX, and ROWMAX is its absolute value
278 *
279  rowmax = zero
280  jmax = imax
281  kx = imax*( imax+1 ) / 2 + imax
282  DO 20 j = imax + 1, k
283  IF( cabs1( ap( kx ) ).GT.rowmax ) THEN
284  rowmax = cabs1( ap( kx ) )
285  jmax = j
286  END IF
287  kx = kx + j
288  20 CONTINUE
289  kpc = ( imax-1 )*imax / 2 + 1
290  IF( imax.GT.1 ) THEN
291  jmax = icamax( imax-1, ap( kpc ), 1 )
292  rowmax = max( rowmax, cabs1( ap( kpc+jmax-1 ) ) )
293  END IF
294 *
295  IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
296 *
297 * no interchange, use 1-by-1 pivot block
298 *
299  kp = k
300  ELSE IF( abs( real( ap( kpc+imax-1 ) ) ).GE.alpha*
301  $ rowmax ) THEN
302 *
303 * interchange rows and columns K and IMAX, use 1-by-1
304 * pivot block
305 *
306  kp = imax
307  ELSE
308 *
309 * interchange rows and columns K-1 and IMAX, use 2-by-2
310 * pivot block
311 *
312  kp = imax
313  kstep = 2
314  END IF
315  END IF
316 *
317  kk = k - kstep + 1
318  IF( kstep.EQ.2 )
319  $ knc = knc - k + 1
320  IF( kp.NE.kk ) THEN
321 *
322 * Interchange rows and columns KK and KP in the leading
323 * submatrix A(1:k,1:k)
324 *
325  CALL cswap( kp-1, ap( knc ), 1, ap( kpc ), 1 )
326  kx = kpc + kp - 1
327  DO 30 j = kp + 1, kk - 1
328  kx = kx + j - 1
329  t = conjg( ap( knc+j-1 ) )
330  ap( knc+j-1 ) = conjg( ap( kx ) )
331  ap( kx ) = t
332  30 CONTINUE
333  ap( kx+kk-1 ) = conjg( ap( kx+kk-1 ) )
334  r1 = real( ap( knc+kk-1 ) )
335  ap( knc+kk-1 ) = real( ap( kpc+kp-1 ) )
336  ap( kpc+kp-1 ) = r1
337  IF( kstep.EQ.2 ) THEN
338  ap( kc+k-1 ) = real( ap( kc+k-1 ) )
339  t = ap( kc+k-2 )
340  ap( kc+k-2 ) = ap( kc+kp-1 )
341  ap( kc+kp-1 ) = t
342  END IF
343  ELSE
344  ap( kc+k-1 ) = real( ap( kc+k-1 ) )
345  IF( kstep.EQ.2 )
346  $ ap( kc-1 ) = real( ap( kc-1 ) )
347  END IF
348 *
349 * Update the leading submatrix
350 *
351  IF( kstep.EQ.1 ) THEN
352 *
353 * 1-by-1 pivot block D(k): column k now holds
354 *
355 * W(k) = U(k)*D(k)
356 *
357 * where U(k) is the k-th column of U
358 *
359 * Perform a rank-1 update of A(1:k-1,1:k-1) as
360 *
361 * A := A - U(k)*D(k)*U(k)**H = A - W(k)*1/D(k)*W(k)**H
362 *
363  r1 = one / real( ap( kc+k-1 ) )
364  CALL chpr( uplo, k-1, -r1, ap( kc ), 1, ap )
365 *
366 * Store U(k) in column k
367 *
368  CALL csscal( k-1, r1, ap( kc ), 1 )
369  ELSE
370 *
371 * 2-by-2 pivot block D(k): columns k and k-1 now hold
372 *
373 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
374 *
375 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
376 * of U
377 *
378 * Perform a rank-2 update of A(1:k-2,1:k-2) as
379 *
380 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**H
381 * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**H
382 *
383  IF( k.GT.2 ) THEN
384 *
385  d = slapy2( real( ap( k-1+( k-1 )*k / 2 ) ),
386  $ aimag( ap( k-1+( k-1 )*k / 2 ) ) )
387  d22 = real( ap( k-1+( k-2 )*( k-1 ) / 2 ) ) / d
388  d11 = real( ap( k+( k-1 )*k / 2 ) ) / d
389  tt = one / ( d11*d22-one )
390  d12 = ap( k-1+( k-1 )*k / 2 ) / d
391  d = tt / d
392 *
393  DO 50 j = k - 2, 1, -1
394  wkm1 = d*( d11*ap( j+( k-2 )*( k-1 ) / 2 )-
395  $ conjg( d12 )*ap( j+( k-1 )*k / 2 ) )
396  wk = d*( d22*ap( j+( k-1 )*k / 2 )-d12*
397  $ ap( j+( k-2 )*( k-1 ) / 2 ) )
398  DO 40 i = j, 1, -1
399  ap( i+( j-1 )*j / 2 ) = ap( i+( j-1 )*j / 2 ) -
400  $ ap( i+( k-1 )*k / 2 )*conjg( wk ) -
401  $ ap( i+( k-2 )*( k-1 ) / 2 )*conjg( wkm1 )
402  40 CONTINUE
403  ap( j+( k-1 )*k / 2 ) = wk
404  ap( j+( k-2 )*( k-1 ) / 2 ) = wkm1
405  ap( j+( j-1 )*j / 2 ) = cmplx( real( ap( j+( j-1 )*
406  $ j / 2 ) ), 0.0e+0 )
407  50 CONTINUE
408 *
409  END IF
410 *
411  END IF
412  END IF
413 *
414 * Store details of the interchanges in IPIV
415 *
416  IF( kstep.EQ.1 ) THEN
417  ipiv( k ) = kp
418  ELSE
419  ipiv( k ) = -kp
420  ipiv( k-1 ) = -kp
421  END IF
422 *
423 * Decrease K and return to the start of the main loop
424 *
425  k = k - kstep
426  kc = knc - k
427  GO TO 10
428 *
429  ELSE
430 *
431 * Factorize A as L*D*L**H using the lower triangle of A
432 *
433 * K is the main loop index, increasing from 1 to N in steps of
434 * 1 or 2
435 *
436  k = 1
437  kc = 1
438  npp = n*( n+1 ) / 2
439  60 CONTINUE
440  knc = kc
441 *
442 * If K > N, exit from loop
443 *
444  IF( k.GT.n )
445  $ GO TO 110
446  kstep = 1
447 *
448 * Determine rows and columns to be interchanged and whether
449 * a 1-by-1 or 2-by-2 pivot block will be used
450 *
451  absakk = abs( real( ap( kc ) ) )
452 *
453 * IMAX is the row-index of the largest off-diagonal element in
454 * column K, and COLMAX is its absolute value
455 *
456  IF( k.LT.n ) THEN
457  imax = k + icamax( n-k, ap( kc+1 ), 1 )
458  colmax = cabs1( ap( kc+imax-k ) )
459  ELSE
460  colmax = zero
461  END IF
462 *
463  IF( max( absakk, colmax ).EQ.zero ) THEN
464 *
465 * Column K is zero: set INFO and continue
466 *
467  IF( info.EQ.0 )
468  $ info = k
469  kp = k
470  ap( kc ) = real( ap( kc ) )
471  ELSE
472  IF( absakk.GE.alpha*colmax ) THEN
473 *
474 * no interchange, use 1-by-1 pivot block
475 *
476  kp = k
477  ELSE
478 *
479 * JMAX is the column-index of the largest off-diagonal
480 * element in row IMAX, and ROWMAX is its absolute value
481 *
482  rowmax = zero
483  kx = kc + imax - k
484  DO 70 j = k, imax - 1
485  IF( cabs1( ap( kx ) ).GT.rowmax ) THEN
486  rowmax = cabs1( ap( kx ) )
487  jmax = j
488  END IF
489  kx = kx + n - j
490  70 CONTINUE
491  kpc = npp - ( n-imax+1 )*( n-imax+2 ) / 2 + 1
492  IF( imax.LT.n ) THEN
493  jmax = imax + icamax( n-imax, ap( kpc+1 ), 1 )
494  rowmax = max( rowmax, cabs1( ap( kpc+jmax-imax ) ) )
495  END IF
496 *
497  IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
498 *
499 * no interchange, use 1-by-1 pivot block
500 *
501  kp = k
502  ELSE IF( abs( real( ap( kpc ) ) ).GE.alpha*rowmax ) THEN
503 *
504 * interchange rows and columns K and IMAX, use 1-by-1
505 * pivot block
506 *
507  kp = imax
508  ELSE
509 *
510 * interchange rows and columns K+1 and IMAX, use 2-by-2
511 * pivot block
512 *
513  kp = imax
514  kstep = 2
515  END IF
516  END IF
517 *
518  kk = k + kstep - 1
519  IF( kstep.EQ.2 )
520  $ knc = knc + n - k + 1
521  IF( kp.NE.kk ) THEN
522 *
523 * Interchange rows and columns KK and KP in the trailing
524 * submatrix A(k:n,k:n)
525 *
526  IF( kp.LT.n )
527  $ CALL cswap( n-kp, ap( knc+kp-kk+1 ), 1, ap( kpc+1 ),
528  $ 1 )
529  kx = knc + kp - kk
530  DO 80 j = kk + 1, kp - 1
531  kx = kx + n - j + 1
532  t = conjg( ap( knc+j-kk ) )
533  ap( knc+j-kk ) = conjg( ap( kx ) )
534  ap( kx ) = t
535  80 CONTINUE
536  ap( knc+kp-kk ) = conjg( ap( knc+kp-kk ) )
537  r1 = real( ap( knc ) )
538  ap( knc ) = real( ap( kpc ) )
539  ap( kpc ) = r1
540  IF( kstep.EQ.2 ) THEN
541  ap( kc ) = real( ap( kc ) )
542  t = ap( kc+1 )
543  ap( kc+1 ) = ap( kc+kp-k )
544  ap( kc+kp-k ) = t
545  END IF
546  ELSE
547  ap( kc ) = real( ap( kc ) )
548  IF( kstep.EQ.2 )
549  $ ap( knc ) = real( ap( knc ) )
550  END IF
551 *
552 * Update the trailing submatrix
553 *
554  IF( kstep.EQ.1 ) THEN
555 *
556 * 1-by-1 pivot block D(k): column k now holds
557 *
558 * W(k) = L(k)*D(k)
559 *
560 * where L(k) is the k-th column of L
561 *
562  IF( k.LT.n ) THEN
563 *
564 * Perform a rank-1 update of A(k+1:n,k+1:n) as
565 *
566 * A := A - L(k)*D(k)*L(k)**H = A - W(k)*(1/D(k))*W(k)**H
567 *
568  r1 = one / real( ap( kc ) )
569  CALL chpr( uplo, n-k, -r1, ap( kc+1 ), 1,
570  $ ap( kc+n-k+1 ) )
571 *
572 * Store L(k) in column K
573 *
574  CALL csscal( n-k, r1, ap( kc+1 ), 1 )
575  END IF
576  ELSE
577 *
578 * 2-by-2 pivot block D(k): columns K and K+1 now hold
579 *
580 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
581 *
582 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
583 * of L
584 *
585  IF( k.LT.n-1 ) THEN
586 *
587 * Perform a rank-2 update of A(k+2:n,k+2:n) as
588 *
589 * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**H
590 * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**H
591 *
592 * where L(k) and L(k+1) are the k-th and (k+1)-th
593 * columns of L
594 *
595  d = slapy2( real( ap( k+1+( k-1 )*( 2*n-k ) / 2 ) ),
596  $ aimag( ap( k+1+( k-1 )*( 2*n-k ) / 2 ) ) )
597  d11 = real( ap( k+1+k*( 2*n-k-1 ) / 2 ) ) / d
598  d22 = real( ap( k+( k-1 )*( 2*n-k ) / 2 ) ) / d
599  tt = one / ( d11*d22-one )
600  d21 = ap( k+1+( k-1 )*( 2*n-k ) / 2 ) / d
601  d = tt / d
602 *
603  DO 100 j = k + 2, n
604  wk = d*( d11*ap( j+( k-1 )*( 2*n-k ) / 2 )-d21*
605  $ ap( j+k*( 2*n-k-1 ) / 2 ) )
606  wkp1 = d*( d22*ap( j+k*( 2*n-k-1 ) / 2 )-
607  $ conjg( d21 )*ap( j+( k-1 )*( 2*n-k ) / 2 ) )
608  DO 90 i = j, n
609  ap( i+( j-1 )*( 2*n-j ) / 2 ) = ap( i+( j-1 )*
610  $ ( 2*n-j ) / 2 ) - ap( i+( k-1 )*( 2*n-k ) /
611  $ 2 )*conjg( wk ) - ap( i+k*( 2*n-k-1 ) / 2 )*
612  $ conjg( wkp1 )
613  90 CONTINUE
614  ap( j+( k-1 )*( 2*n-k ) / 2 ) = wk
615  ap( j+k*( 2*n-k-1 ) / 2 ) = wkp1
616  ap( j+( j-1 )*( 2*n-j ) / 2 )
617  $ = cmplx( real( ap( j+( j-1 )*( 2*n-j ) / 2 ) ),
618  $ 0.0e+0 )
619  100 CONTINUE
620  END IF
621  END IF
622  END IF
623 *
624 * Store details of the interchanges in IPIV
625 *
626  IF( kstep.EQ.1 ) THEN
627  ipiv( k ) = kp
628  ELSE
629  ipiv( k ) = -kp
630  ipiv( k+1 ) = -kp
631  END IF
632 *
633 * Increase K and return to the start of the main loop
634 *
635  k = k + kstep
636  kc = knc + n - k + 2
637  GO TO 60
638 *
639  END IF
640 *
641  110 CONTINUE
642  RETURN
643 *
644 * End of CHPTRF
645 *
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:63
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:71
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:78
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine chpr(UPLO, N, ALPHA, X, INCX, AP)
CHPR
Definition: chpr.f:130
Here is the call graph for this function:
Here is the caller graph for this function: