 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ zlasyf_rk()

 subroutine zlasyf_rk ( character UPLO, integer N, integer NB, integer KB, complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( * ) E, integer, dimension( * ) IPIV, complex*16, dimension( ldw, * ) W, integer LDW, integer INFO )

ZLASYF_RK computes a partial factorization of a complex symmetric indefinite matrix using bounded Bunch-Kaufman (rook) diagonal pivoting method.

Purpose:
``` ZLASYF_RK computes a partial factorization of a complex symmetric
matrix A using the bounded Bunch-Kaufman (rook) diagonal
pivoting method. The partial factorization has the form:

A  =  ( I  U12 ) ( A11  0  ) (  I       0    )  if UPLO = 'U', or:
( 0  U22 ) (  0   D  ) ( U12**T U22**T )

A  =  ( L11  0 ) (  D   0  ) ( L11**T L21**T )  if UPLO = 'L',
( L21  I ) (  0  A22 ) (  0       I    )

where the order of D is at most NB. The actual order is returned in
the argument KB, and is either NB or NB-1, or N if N <= NB.

ZLASYF_RK is an auxiliary routine called by ZSYTRF_RK. It uses
blocked code (calling Level 3 BLAS) to update the submatrix
A11 (if UPLO = 'U') or A22 (if UPLO = 'L').```
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the symmetric matrix A is stored: = 'U': Upper triangular = 'L': Lower triangular``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in] NB ``` NB is INTEGER The maximum number of columns of the matrix A that should be factored. NB should be at least 2 to allow for 2-by-2 pivot blocks.``` [out] KB ``` KB is INTEGER The number of columns of A that were actually factored. KB is either NB-1 or NB, or N if N <= NB.``` [in,out] A ``` A is COMPLEX*16 array, dimension (LDA,N) On entry, the symmetric matrix A. If UPLO = 'U': the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L': the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. On exit, contains: a) ONLY diagonal elements of the symmetric block diagonal matrix D on the diagonal of A, i.e. D(k,k) = A(k,k); (superdiagonal (or subdiagonal) elements of D are stored on exit in array E), and b) If UPLO = 'U': factor U in the superdiagonal part of A. If UPLO = 'L': factor L in the subdiagonal part of A.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).``` [out] E ``` E is COMPLEX*16 array, dimension (N) On exit, contains the superdiagonal (or subdiagonal) elements of the symmetric block diagonal matrix D with 1-by-1 or 2-by-2 diagonal blocks, where If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0; If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0. NOTE: For 1-by-1 diagonal block D(k), where 1 <= k <= N, the element E(k) is set to 0 in both UPLO = 'U' or UPLO = 'L' cases.``` [out] IPIV ``` IPIV is INTEGER array, dimension (N) IPIV describes the permutation matrix P in the factorization of matrix A as follows. The absolute value of IPIV(k) represents the index of row and column that were interchanged with the k-th row and column. The value of UPLO describes the order in which the interchanges were applied. Also, the sign of IPIV represents the block structure of the symmetric block diagonal matrix D with 1-by-1 or 2-by-2 diagonal blocks which correspond to 1 or 2 interchanges at each factorization step. If UPLO = 'U', ( in factorization order, k decreases from N to 1 ): a) A single positive entry IPIV(k) > 0 means: D(k,k) is a 1-by-1 diagonal block. If IPIV(k) != k, rows and columns k and IPIV(k) were interchanged in the submatrix A(1:N,N-KB+1:N); If IPIV(k) = k, no interchange occurred. b) A pair of consecutive negative entries IPIV(k) < 0 and IPIV(k-1) < 0 means: D(k-1:k,k-1:k) is a 2-by-2 diagonal block. (NOTE: negative entries in IPIV appear ONLY in pairs). 1) If -IPIV(k) != k, rows and columns k and -IPIV(k) were interchanged in the matrix A(1:N,N-KB+1:N). If -IPIV(k) = k, no interchange occurred. 2) If -IPIV(k-1) != k-1, rows and columns k-1 and -IPIV(k-1) were interchanged in the submatrix A(1:N,N-KB+1:N). If -IPIV(k-1) = k-1, no interchange occurred. c) In both cases a) and b) is always ABS( IPIV(k) ) <= k. d) NOTE: Any entry IPIV(k) is always NONZERO on output. If UPLO = 'L', ( in factorization order, k increases from 1 to N ): a) A single positive entry IPIV(k) > 0 means: D(k,k) is a 1-by-1 diagonal block. If IPIV(k) != k, rows and columns k and IPIV(k) were interchanged in the submatrix A(1:N,1:KB). If IPIV(k) = k, no interchange occurred. b) A pair of consecutive negative entries IPIV(k) < 0 and IPIV(k+1) < 0 means: D(k:k+1,k:k+1) is a 2-by-2 diagonal block. (NOTE: negative entries in IPIV appear ONLY in pairs). 1) If -IPIV(k) != k, rows and columns k and -IPIV(k) were interchanged in the submatrix A(1:N,1:KB). If -IPIV(k) = k, no interchange occurred. 2) If -IPIV(k+1) != k+1, rows and columns k-1 and -IPIV(k-1) were interchanged in the submatrix A(1:N,1:KB). If -IPIV(k+1) = k+1, no interchange occurred. c) In both cases a) and b) is always ABS( IPIV(k) ) >= k. d) NOTE: Any entry IPIV(k) is always NONZERO on output.``` [out] W ` W is COMPLEX*16 array, dimension (LDW,NB)` [in] LDW ``` LDW is INTEGER The leading dimension of the array W. LDW >= max(1,N).``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: If INFO = -k, the k-th argument had an illegal value > 0: If INFO = k, the matrix A is singular, because: If UPLO = 'U': column k in the upper triangular part of A contains all zeros. If UPLO = 'L': column k in the lower triangular part of A contains all zeros. Therefore D(k,k) is exactly zero, and superdiagonal elements of column k of U (or subdiagonal elements of column k of L ) are all zeros. 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. NOTE: INFO only stores the first occurrence of a singularity, any subsequent occurrence of singularity is not stored in INFO even though the factorization always completes.```
Contributors:
```  December 2016,  Igor Kozachenko,
Computer Science Division,
University of California, Berkeley

September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
School of Mathematics,
University of Manchester```

Definition at line 260 of file zlasyf_rk.f.

262*
263* -- LAPACK computational routine --
264* -- LAPACK is a software package provided by Univ. of Tennessee, --
265* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
266*
267* .. Scalar Arguments ..
268 CHARACTER UPLO
269 INTEGER INFO, KB, LDA, LDW, N, NB
270* ..
271* .. Array Arguments ..
272 INTEGER IPIV( * )
273 COMPLEX*16 A( LDA, * ), E( * ), W( LDW, * )
274* ..
275*
276* =====================================================================
277*
278* .. Parameters ..
279 DOUBLE PRECISION ZERO, ONE
280 parameter( zero = 0.0d+0, one = 1.0d+0 )
281 DOUBLE PRECISION EIGHT, SEVTEN
282 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
283 COMPLEX*16 CONE, CZERO
284 parameter( cone = ( 1.0d+0, 0.0d+0 ),
285 \$ czero = ( 0.0d+0, 0.0d+0 ) )
286* ..
287* .. Local Scalars ..
288 LOGICAL DONE
289 INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, K, KK, KW, KKW,
290 \$ KP, KSTEP, P, II
291 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX, SFMIN, DTEMP
292 COMPLEX*16 D11, D12, D21, D22, R1, T, Z
293* ..
294* .. External Functions ..
295 LOGICAL LSAME
296 INTEGER IZAMAX
297 DOUBLE PRECISION DLAMCH
298 EXTERNAL lsame, izamax, dlamch
299* ..
300* .. External Subroutines ..
301 EXTERNAL zcopy, zgemm, zgemv, zscal, zswap
302* ..
303* .. Intrinsic Functions ..
304 INTRINSIC abs, dble, dimag, max, min, sqrt
305* ..
306* .. Statement Functions ..
307 DOUBLE PRECISION CABS1
308* ..
309* .. Statement Function definitions ..
310 cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
311* ..
312* .. Executable Statements ..
313*
314 info = 0
315*
316* Initialize ALPHA for use in choosing pivot block size.
317*
318 alpha = ( one+sqrt( sevten ) ) / eight
319*
320* Compute machine safe minimum
321*
322 sfmin = dlamch( 'S' )
323*
324 IF( lsame( uplo, 'U' ) ) THEN
325*
326* Factorize the trailing columns of A using the upper triangle
327* of A and working backwards, and compute the matrix W = U12*D
328* for use in updating A11
329*
330* Initialize the first entry of array E, where superdiagonal
331* elements of D are stored
332*
333 e( 1 ) = czero
334*
335* K is the main loop index, decreasing from N in steps of 1 or 2
336*
337 k = n
338 10 CONTINUE
339*
340* KW is the column of W which corresponds to column K of A
341*
342 kw = nb + k - n
343*
344* Exit from loop
345*
346 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
347 \$ GO TO 30
348*
349 kstep = 1
350 p = k
351*
352* Copy column K of A to column KW of W and update it
353*
354 CALL zcopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
355 IF( k.LT.n )
356 \$ CALL zgemv( 'No transpose', k, n-k, -cone, a( 1, k+1 ),
357 \$ lda, w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
358*
359* Determine rows and columns to be interchanged and whether
360* a 1-by-1 or 2-by-2 pivot block will be used
361*
362 absakk = cabs1( w( k, kw ) )
363*
364* IMAX is the row-index of the largest off-diagonal element in
365* column K, and COLMAX is its absolute value.
366* Determine both COLMAX and IMAX.
367*
368 IF( k.GT.1 ) THEN
369 imax = izamax( k-1, w( 1, kw ), 1 )
370 colmax = cabs1( w( imax, kw ) )
371 ELSE
372 colmax = zero
373 END IF
374*
375 IF( max( absakk, colmax ).EQ.zero ) THEN
376*
377* Column K is zero or underflow: set INFO and continue
378*
379 IF( info.EQ.0 )
380 \$ info = k
381 kp = k
382 CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
383*
384* Set E( K ) to zero
385*
386 IF( k.GT.1 )
387 \$ e( k ) = czero
388*
389 ELSE
390*
391* ============================================================
392*
393* Test for interchange
394*
395* Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
396* (used to handle NaN and Inf)
397*
398 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
399*
400* no interchange, use 1-by-1 pivot block
401*
402 kp = k
403*
404 ELSE
405*
406 done = .false.
407*
408* Loop until pivot found
409*
410 12 CONTINUE
411*
412* Begin pivot search loop body
413*
414*
415* Copy column IMAX to column KW-1 of W and update it
416*
417 CALL zcopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
418 CALL zcopy( k-imax, a( imax, imax+1 ), lda,
419 \$ w( imax+1, kw-1 ), 1 )
420*
421 IF( k.LT.n )
422 \$ CALL zgemv( 'No transpose', k, n-k, -cone,
423 \$ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
424 \$ cone, w( 1, kw-1 ), 1 )
425*
426* JMAX is the column-index of the largest off-diagonal
427* element in row IMAX, and ROWMAX is its absolute value.
428* Determine both ROWMAX and JMAX.
429*
430 IF( imax.NE.k ) THEN
431 jmax = imax + izamax( k-imax, w( imax+1, kw-1 ),
432 \$ 1 )
433 rowmax = cabs1( w( jmax, kw-1 ) )
434 ELSE
435 rowmax = zero
436 END IF
437*
438 IF( imax.GT.1 ) THEN
439 itemp = izamax( imax-1, w( 1, kw-1 ), 1 )
440 dtemp = cabs1( w( itemp, kw-1 ) )
441 IF( dtemp.GT.rowmax ) THEN
442 rowmax = dtemp
443 jmax = itemp
444 END IF
445 END IF
446*
447* Equivalent to testing for
448* CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX
449* (used to handle NaN and Inf)
450*
451 IF( .NOT.(cabs1( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
452 \$ THEN
453*
454* interchange rows and columns K and IMAX,
455* use 1-by-1 pivot block
456*
457 kp = imax
458*
459* copy column KW-1 of W to column KW of W
460*
461 CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
462*
463 done = .true.
464*
465* Equivalent to testing for ROWMAX.EQ.COLMAX,
466* (used to handle NaN and Inf)
467*
468 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
469 \$ THEN
470*
471* interchange rows and columns K-1 and IMAX,
472* use 2-by-2 pivot block
473*
474 kp = imax
475 kstep = 2
476 done = .true.
477 ELSE
478*
480*
481 p = imax
482 colmax = rowmax
483 imax = jmax
484*
485* Copy updated JMAXth (next IMAXth) column to Kth of W
486*
487 CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
488*
489 END IF
490*
491* End pivot search loop body
492*
493 IF( .NOT. done ) GOTO 12
494*
495 END IF
496*
497* ============================================================
498*
499 kk = k - kstep + 1
500*
501* KKW is the column of W which corresponds to column KK of A
502*
503 kkw = nb + kk - n
504*
505 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
506*
507* Copy non-updated column K to column P
508*
509 CALL zcopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
510 CALL zcopy( p, a( 1, k ), 1, a( 1, p ), 1 )
511*
512* Interchange rows K and P in last N-K+1 columns of A
513* and last N-K+2 columns of W
514*
515 CALL zswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
516 CALL zswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
517 END IF
518*
519* Updated column KP is already stored in column KKW of W
520*
521 IF( kp.NE.kk ) THEN
522*
523* Copy non-updated column KK to column KP
524*
525 a( kp, k ) = a( kk, k )
526 CALL zcopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
527 \$ lda )
528 CALL zcopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
529*
530* Interchange rows KK and KP in last N-KK+1 columns
531* of A and W
532*
533 CALL zswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
534 CALL zswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
535 \$ ldw )
536 END IF
537*
538 IF( kstep.EQ.1 ) THEN
539*
540* 1-by-1 pivot block D(k): column KW of W now holds
541*
542* W(k) = U(k)*D(k)
543*
544* where U(k) is the k-th column of U
545*
546* Store U(k) in column k of A
547*
548 CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
549 IF( k.GT.1 ) THEN
550 IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
551 r1 = cone / a( k, k )
552 CALL zscal( k-1, r1, a( 1, k ), 1 )
553 ELSE IF( a( k, k ).NE.czero ) THEN
554 DO 14 ii = 1, k - 1
555 a( ii, k ) = a( ii, k ) / a( k, k )
556 14 CONTINUE
557 END IF
558*
559* Store the superdiagonal element of D in array E
560*
561 e( k ) = czero
562*
563 END IF
564*
565 ELSE
566*
567* 2-by-2 pivot block D(k): columns KW and KW-1 of W now
568* hold
569*
570* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
571*
572* where U(k) and U(k-1) are the k-th and (k-1)-th columns
573* of U
574*
575 IF( k.GT.2 ) THEN
576*
577* Store U(k) and U(k-1) in columns k and k-1 of A
578*
579 d12 = w( k-1, kw )
580 d11 = w( k, kw ) / d12
581 d22 = w( k-1, kw-1 ) / d12
582 t = cone / ( d11*d22-cone )
583 DO 20 j = 1, k - 2
584 a( j, k-1 ) = t*( (d11*w( j, kw-1 )-w( j, kw ) ) /
585 \$ d12 )
586 a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
587 \$ d12 )
588 20 CONTINUE
589 END IF
590*
591* Copy diagonal elements of D(K) to A,
592* copy superdiagonal element of D(K) to E(K) and
593* ZERO out superdiagonal entry of A
594*
595 a( k-1, k-1 ) = w( k-1, kw-1 )
596 a( k-1, k ) = czero
597 a( k, k ) = w( k, kw )
598 e( k ) = w( k-1, kw )
599 e( k-1 ) = czero
600*
601 END IF
602*
603* End column K is nonsingular
604*
605 END IF
606*
607* Store details of the interchanges in IPIV
608*
609 IF( kstep.EQ.1 ) THEN
610 ipiv( k ) = kp
611 ELSE
612 ipiv( k ) = -p
613 ipiv( k-1 ) = -kp
614 END IF
615*
616* Decrease K and return to the start of the main loop
617*
618 k = k - kstep
619 GO TO 10
620*
621 30 CONTINUE
622*
623* Update the upper triangle of A11 (= A(1:k,1:k)) as
624*
625* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
626*
627* computing blocks of NB columns at a time
628*
629 DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
630 jb = min( nb, k-j+1 )
631*
632* Update the upper triangle of the diagonal block
633*
634 DO 40 jj = j, j + jb - 1
635 CALL zgemv( 'No transpose', jj-j+1, n-k, -cone,
636 \$ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
637 \$ a( j, jj ), 1 )
638 40 CONTINUE
639*
640* Update the rectangular superdiagonal block
641*
642 IF( j.GE.2 )
643 \$ CALL zgemm( 'No transpose', 'Transpose', j-1, jb,
644 \$ n-k, -cone, a( 1, k+1 ), lda, w( j, kw+1 ),
645 \$ ldw, cone, a( 1, j ), lda )
646 50 CONTINUE
647*
648* Set KB to the number of columns factorized
649*
650 kb = n - k
651*
652 ELSE
653*
654* Factorize the leading columns of A using the lower triangle
655* of A and working forwards, and compute the matrix W = L21*D
656* for use in updating A22
657*
658* Initialize the unused last entry of the subdiagonal array E.
659*
660 e( n ) = czero
661*
662* K is the main loop index, increasing from 1 in steps of 1 or 2
663*
664 k = 1
665 70 CONTINUE
666*
667* Exit from loop
668*
669 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
670 \$ GO TO 90
671*
672 kstep = 1
673 p = k
674*
675* Copy column K of A to column K of W and update it
676*
677 CALL zcopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
678 IF( k.GT.1 )
679 \$ CALL zgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
680 \$ lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
681*
682* Determine rows and columns to be interchanged and whether
683* a 1-by-1 or 2-by-2 pivot block will be used
684*
685 absakk = cabs1( w( k, k ) )
686*
687* IMAX is the row-index of the largest off-diagonal element in
688* column K, and COLMAX is its absolute value.
689* Determine both COLMAX and IMAX.
690*
691 IF( k.LT.n ) THEN
692 imax = k + izamax( n-k, w( k+1, k ), 1 )
693 colmax = cabs1( w( imax, k ) )
694 ELSE
695 colmax = zero
696 END IF
697*
698 IF( max( absakk, colmax ).EQ.zero ) THEN
699*
700* Column K is zero or underflow: set INFO and continue
701*
702 IF( info.EQ.0 )
703 \$ info = k
704 kp = k
705 CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
706*
707* Set E( K ) to zero
708*
709 IF( k.LT.n )
710 \$ e( k ) = czero
711*
712 ELSE
713*
714* ============================================================
715*
716* Test for interchange
717*
718* Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
719* (used to handle NaN and Inf)
720*
721 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
722*
723* no interchange, use 1-by-1 pivot block
724*
725 kp = k
726*
727 ELSE
728*
729 done = .false.
730*
731* Loop until pivot found
732*
733 72 CONTINUE
734*
735* Begin pivot search loop body
736*
737*
738* Copy column IMAX to column K+1 of W and update it
739*
740 CALL zcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
741 CALL zcopy( n-imax+1, a( imax, imax ), 1,
742 \$ w( imax, k+1 ), 1 )
743 IF( k.GT.1 )
744 \$ CALL zgemv( 'No transpose', n-k+1, k-1, -cone,
745 \$ a( k, 1 ), lda, w( imax, 1 ), ldw,
746 \$ cone, w( k, k+1 ), 1 )
747*
748* JMAX is the column-index of the largest off-diagonal
749* element in row IMAX, and ROWMAX is its absolute value.
750* Determine both ROWMAX and JMAX.
751*
752 IF( imax.NE.k ) THEN
753 jmax = k - 1 + izamax( imax-k, w( k, k+1 ), 1 )
754 rowmax = cabs1( w( jmax, k+1 ) )
755 ELSE
756 rowmax = zero
757 END IF
758*
759 IF( imax.LT.n ) THEN
760 itemp = imax + izamax( n-imax, w( imax+1, k+1 ), 1)
761 dtemp = cabs1( w( itemp, k+1 ) )
762 IF( dtemp.GT.rowmax ) THEN
763 rowmax = dtemp
764 jmax = itemp
765 END IF
766 END IF
767*
768* Equivalent to testing for
769* CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX
770* (used to handle NaN and Inf)
771*
772 IF( .NOT.( cabs1( w( imax, k+1 ) ).LT.alpha*rowmax ) )
773 \$ THEN
774*
775* interchange rows and columns K and IMAX,
776* use 1-by-1 pivot block
777*
778 kp = imax
779*
780* copy column K+1 of W to column K of W
781*
782 CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
783*
784 done = .true.
785*
786* Equivalent to testing for ROWMAX.EQ.COLMAX,
787* (used to handle NaN and Inf)
788*
789 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
790 \$ THEN
791*
792* interchange rows and columns K+1 and IMAX,
793* use 2-by-2 pivot block
794*
795 kp = imax
796 kstep = 2
797 done = .true.
798 ELSE
799*
801*
802 p = imax
803 colmax = rowmax
804 imax = jmax
805*
806* Copy updated JMAXth (next IMAXth) column to Kth of W
807*
808 CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
809*
810 END IF
811*
812* End pivot search loop body
813*
814 IF( .NOT. done ) GOTO 72
815*
816 END IF
817*
818* ============================================================
819*
820 kk = k + kstep - 1
821*
822 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
823*
824* Copy non-updated column K to column P
825*
826 CALL zcopy( p-k, a( k, k ), 1, a( p, k ), lda )
827 CALL zcopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
828*
829* Interchange rows K and P in first K columns of A
830* and first K+1 columns of W
831*
832 CALL zswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
833 CALL zswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
834 END IF
835*
836* Updated column KP is already stored in column KK of W
837*
838 IF( kp.NE.kk ) THEN
839*
840* Copy non-updated column KK to column KP
841*
842 a( kp, k ) = a( kk, k )
843 CALL zcopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
844 CALL zcopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
845*
846* Interchange rows KK and KP in first KK columns of A and W
847*
848 CALL zswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
849 CALL zswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
850 END IF
851*
852 IF( kstep.EQ.1 ) THEN
853*
854* 1-by-1 pivot block D(k): column k of W now holds
855*
856* W(k) = L(k)*D(k)
857*
858* where L(k) is the k-th column of L
859*
860* Store L(k) in column k of A
861*
862 CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
863 IF( k.LT.n ) THEN
864 IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
865 r1 = cone / a( k, k )
866 CALL zscal( n-k, r1, a( k+1, k ), 1 )
867 ELSE IF( a( k, k ).NE.czero ) THEN
868 DO 74 ii = k + 1, n
869 a( ii, k ) = a( ii, k ) / a( k, k )
870 74 CONTINUE
871 END IF
872*
873* Store the subdiagonal element of D in array E
874*
875 e( k ) = czero
876*
877 END IF
878*
879 ELSE
880*
881* 2-by-2 pivot block D(k): columns k and k+1 of W now hold
882*
883* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
884*
885* where L(k) and L(k+1) are the k-th and (k+1)-th columns
886* of L
887*
888 IF( k.LT.n-1 ) THEN
889*
890* Store L(k) and L(k+1) in columns k and k+1 of A
891*
892 d21 = w( k+1, k )
893 d11 = w( k+1, k+1 ) / d21
894 d22 = w( k, k ) / d21
895 t = cone / ( d11*d22-cone )
896 DO 80 j = k + 2, n
897 a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
898 \$ d21 )
899 a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
900 \$ d21 )
901 80 CONTINUE
902 END IF
903*
904* Copy diagonal elements of D(K) to A,
905* copy subdiagonal element of D(K) to E(K) and
906* ZERO out subdiagonal entry of A
907*
908 a( k, k ) = w( k, k )
909 a( k+1, k ) = czero
910 a( k+1, k+1 ) = w( k+1, k+1 )
911 e( k ) = w( k+1, k )
912 e( k+1 ) = czero
913*
914 END IF
915*
916* End column K is nonsingular
917*
918 END IF
919*
920* Store details of the interchanges in IPIV
921*
922 IF( kstep.EQ.1 ) THEN
923 ipiv( k ) = kp
924 ELSE
925 ipiv( k ) = -p
926 ipiv( k+1 ) = -kp
927 END IF
928*
929* Increase K and return to the start of the main loop
930*
931 k = k + kstep
932 GO TO 70
933*
934 90 CONTINUE
935*
936* Update the lower triangle of A22 (= A(k:n,k:n)) as
937*
938* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
939*
940* computing blocks of NB columns at a time
941*
942 DO 110 j = k, n, nb
943 jb = min( nb, n-j+1 )
944*
945* Update the lower triangle of the diagonal block
946*
947 DO 100 jj = j, j + jb - 1
948 CALL zgemv( 'No transpose', j+jb-jj, k-1, -cone,
949 \$ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
950 \$ a( jj, jj ), 1 )
951 100 CONTINUE
952*
953* Update the rectangular subdiagonal block
954*
955 IF( j+jb.LE.n )
956 \$ CALL zgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
957 \$ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
958 \$ ldw, cone, a( j+jb, j ), lda )
959 110 CONTINUE
960*
961* Set KB to the number of columns factorized
962*
963 kb = k - 1
964*
965 END IF
966*
967 RETURN
968*
969* End of ZLASYF_RK
970*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:71
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
Here is the call graph for this function:
Here is the caller graph for this function: