fixed comments in LAPACK routines (c,z)lahef.f and (c,z)lahef_rook.f

1 files changed, 62 insertions, 19 deletions

diff --git a/SRC/zlahef_rook.f b/SRC/zlahef_rook.f
index 01713a58..9b1e05e7 100644
--- a/SRC/zlahef_rook.f
+++ b/SRC/zlahef_rook.f
@@ -310,8 +310,9 @@
 *
 *           ============================================================
 *
-*           Test for interchange
+*           BEGIN pivot search
 *
+*           Case(1)
 *           Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
 *           (used to handle NaN and Inf)
 *
@@ -329,7 +330,7 @@
 *
    12          CONTINUE
 *
-*                 Begin pivot search loop body
+*                 BEGIN pivot search loop body
 *
 *
 *                 Copy column IMAX to column KW-1 of W and update it
@@ -369,6 +370,7 @@
                      END IF
                   END IF
 *
+*                 Case(2)
 *                 Equivalent to testing for
 *                 ABS( DBLE( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
 *                 (used to handle NaN and Inf)
@@ -387,6 +389,7 @@
 *
                      DONE = .TRUE.
 *
+*                 Case(3)
 *                 Equivalent to testing for ROWMAX.EQ.COLMAX,
 *                 (used to handle NaN and Inf)
 *
@@ -399,6 +402,8 @@
                      KP = IMAX
                      KSTEP = 2
                      DONE = .TRUE.
+*
+*                 Case(4)
                   ELSE
 *
 *                    Pivot not found: set params and repeat
@@ -413,12 +418,14 @@
 *
                   END IF
 *
-*                 End pivot search loop body
+*                 END pivot search loop body
 *
                IF( .NOT.DONE ) GOTO 12
 *
             END IF
 *
+*           END pivot search
+*
 *           ============================================================
 *
 *           KK is the column of A where pivoting step stopped
@@ -508,6 +515,10 @@
                CALL ZCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
                IF( K.GT.1 ) THEN
 *
+*                 (NOTE: No need to check if A(k,k) is NOT ZERO,
+*                  since that was ensured earlier in pivot search:
+*                  case A(k,k) = 0 falls into 2x2 pivot case(3))
+*
 *                 Handle division by a small number
 *
                   T = DBLE( A( K, K ) )
@@ -515,8 +526,6 @@
                      R1 = ONE / T
                      CALL ZDSCAL( K-1, R1, A( 1, K ), 1 )
                   ELSE
-*                    (NOTE: No need to check if T=D(k,k) is NOT ZERO,
-*                    since that was ensured earlier in pivot search)
                      DO 14 II = 1, K-1
                         A( II, K ) = A( II, K ) / T
    14                CONTINUE
@@ -546,10 +555,10 @@
 *
                IF( K.GT.2 ) THEN
 *
-*                 Compose the columns of the inverse of 2-by-2 pivot
-*                 block D in the following way to reduce the number
-*                 of FLOPS when we multiply panel ( W(kw-1) W(kw) ) by
-*                 this inverse
+*                 Factor out the columns of the inverse of 2-by-2 pivot
+*                 block D, so that each column contains 1, to reduce the
+*                 number of FLOPS when we multiply panel
+*                 ( W(kw-1) W(kw) ) by this inverse, i.e. by D**(-1).
 *
 *                 D**(-1) = ( d11 cj(d21) )**(-1) =
 *                           ( d21    d22 )
@@ -577,7 +586,19 @@
 *                 operations is important)
 *
 *                 = ( T*(( D11 )/conj(D21)) T*((  -1 )/D21 ) )
-*                   (   ((  -1 )          )   (( D22 )     ) )
+*                   (   ((  -1 )          )   (( D22 )     ) ),
+*
+*                 where D11 = d22/d21,
+*                       D22 = d11/conj(d21),
+*                       D21 = d21,
+*                       T = 1/(D22*D11-1).
+*
+*                 (NOTE: No need to check for division by ZERO,
+*                  since that was ensured earlier in pivot search:
+*                  (a) d21 != 0 in 2x2 pivot case(4),
+*                      since |d21| should be larger than |d11| and |d22|;
+*                  (b) (D22*D11 - 1) != 0, since from (a),
+*                      both |D11| < 1, |D22| < 1, hence |D22*D11| << 1.)
 *
                   D21 = W( K-1, KW )
                   D11 = W( K, KW ) / DCONJG( D21 )
@@ -748,8 +769,9 @@
 *
 *           ============================================================
 *
-*           Test for interchange
+*           BEGIN pivot search
 *
+*           Case(1)
 *           Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
 *           (used to handle NaN and Inf)
 *
@@ -767,7 +789,7 @@
 *
    72          CONTINUE
 *
-*                 Begin pivot search loop body
+*                 BEGIN pivot search loop body
 *
 *
 *                 Copy column IMAX to column k+1 of W and update it
@@ -805,6 +827,7 @@
                      END IF
                   END IF
 *
+*                 Case(2)
 *                 Equivalent to testing for
 *                 ABS( DBLE( W( IMAX,K+1 ) ) ).GE.ALPHA*ROWMAX
 *                 (used to handle NaN and Inf)
@@ -823,6 +846,7 @@
 *
                      DONE = .TRUE.
 *
+*                 Case(3)
 *                 Equivalent to testing for ROWMAX.EQ.COLMAX,
 *                 (used to handle NaN and Inf)
 *
@@ -835,6 +859,8 @@
                      KP = IMAX
                      KSTEP = 2
                      DONE = .TRUE.
+*
+*                 Case(4)
                   ELSE
 *
 *                    Pivot not found: set params and repeat
@@ -849,12 +875,15 @@
 *
                   END IF
 *
+*
 *                 End pivot search loop body
 *
                IF( .NOT.DONE ) GOTO 72
 *
             END IF
 *
+*           END pivot search
+*
 *           ============================================================
 *
 *           KK is the column of A where pivoting step stopped
@@ -935,6 +964,10 @@
                CALL ZCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
                IF( K.LT.N ) THEN
 *
+*                 (NOTE: No need to check if A(k,k) is NOT ZERO,
+*                  since that was ensured earlier in pivot search:
+*                  case A(k,k) = 0 falls into 2x2 pivot case(3))
+*
 *                 Handle division by a small number
 *
                   T = DBLE( A( K, K ) )
@@ -942,8 +975,6 @@
                      R1 = ONE / T
                      CALL ZDSCAL( N-K, R1, A( K+1, K ), 1 )
                   ELSE
-*                    (NOTE: No need to check if T=D(k,k) is NOT ZERO,
-*                    since that was ensured earlier in pivot search)
                      DO 74 II = K + 1, N
                         A( II, K ) = A( II, K ) / T
    74                CONTINUE
@@ -973,10 +1004,10 @@
 *
                IF( K.LT.N-1 ) THEN
 *
-*                 Compose the columns of the inverse of 2-by-2 pivot
-*                 block D in the following way to reduce the number
-*                 of FLOPS when we multiply panel ( W(k) W(k+1) ) by
-*                 this inverse
+*                 Factor out the columns of the inverse of 2-by-2 pivot
+*                 block D, so that each column contains 1, to reduce the
+*                 number of FLOPS when we multiply panel
+*                 ( W(kw-1) W(kw) ) by this inverse, i.e. by D**(-1).
 *
 *                 D**(-1) = ( d11 cj(d21) )**(-1) =
 *                           ( d21    d22 )
@@ -1004,7 +1035,19 @@
 *                 operations is important)
 *
 *                 = ( T*(( D11 )/conj(D21)) T*((  -1 )/D21 ) )
-*                   (   ((  -1 )          )   (( D22 )     ) )
+*                   (   ((  -1 )          )   (( D22 )     ) ),
+*
+*                 where D11 = d22/d21,
+*                       D22 = d11/conj(d21),
+*                       D21 = d21,
+*                       T = 1/(D22*D11-1).
+*
+*                 (NOTE: No need to check for division by ZERO,
+*                  since that was ensured earlier in pivot search:
+*                  (a) d21 != 0 in 2x2 pivot case(4),
+*                      since |d21| should be larger than |d11| and |d22|;
+*                  (b) (D22*D11 - 1) != 0, since from (a),
+*                      both |D11| < 1, |D22| < 1, hence |D22*D11| << 1.)
 *
                   D21 = W( K+1, K )
                   D11 = W( K+1, K+1 ) / D21