LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine chpr(uplo, n, alpha, x, incx, ap)
CHPR
Definition chpr.f:130
integer function icamax(n, cx, incx)
ICAMAX
Definition icamax.f:71
real function slapy2(x, y)
SLAPY2 returns sqrt(x2+y2).
Definition slapy2.f:63
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: