Merge branch 'Reference-LAPACK:master' into master

This commit is contained in:
Johnathan Rhyne
2024-11-03 19:28:18 -07:00
committed by GitHub
16 changed files with 136 additions and 676 deletions
+10 -45
View File
@@ -200,7 +200,7 @@
PARAMETER ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
* ..
* .. Local Scalars ..
INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
INTEGER IMAX, J, JJ, JMAX, JP, K, KK, KKW, KP,
$ KSTEP, KW
REAL ABSAKK, ALPHA, COLMAX, R1, ROWMAX, T
COMPLEX D11, D21, D22, Z
@@ -211,7 +211,7 @@
EXTERNAL LSAME, ICAMAX
* ..
* .. External Subroutines ..
EXTERNAL CCOPY, CGEMM, CGEMV, CLACGV, CSSCAL,
EXTERNAL CCOPY, CGEMMTR, CGEMV, CLACGV, CSSCAL,
$ CSWAP
* ..
* .. Intrinsic Functions ..
@@ -552,28 +552,11 @@
*
* A11 := A11 - U12*D*U12**H = A11 - U12*W**H
*
* computing blocks of NB columns at a time (note that conjg(W) is
* actually stored)
* (note that conjg(W) is actually stored)
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
A( JJ, JJ ) = REAL( A( JJ, JJ ) )
CALL CGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
$ A( J, JJ ), 1 )
A( JJ, JJ ) = REAL( A( JJ, JJ ) )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
CALL CGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
$ CONE, A( 1, J ), LDA )
50 CONTINUE
CALL CGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ CONE, A( 1, 1 ), LDA )
*
* Put U12 in standard form by partially undoing the interchanges
* in of rows in columns k+1:n looping backwards from k+1 to n
@@ -916,29 +899,11 @@
*
* A22 := A22 - L21*D*L21**H = A22 - L21*W**H
*
* computing blocks of NB columns at a time (note that conjg(W) is
* actually stored)
* (note that conjg(W) is actually stored)
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
A( JJ, JJ ) = REAL( A( JJ, JJ ) )
CALL CGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
$ A( JJ, JJ ), 1 )
A( JJ, JJ ) = REAL( A( JJ, JJ ) )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL CGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
$ LDW, CONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL CGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -CONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ CONE, A( K, K ), LDA )
*
* Put L21 in standard form by partially undoing the interchanges
* of rows in columns 1:k-1 looping backwards from k-1 to 1
+9 -45
View File
@@ -286,7 +286,7 @@
* ..
* .. Local Scalars ..
LOGICAL DONE
INTEGER IMAX, ITEMP, II, J, JB, JJ, JMAX, K, KK, KKW,
INTEGER IMAX, ITEMP, II, J, JMAX, K, KK, KKW,
$ KP, KSTEP, KW, P
REAL ABSAKK, ALPHA, COLMAX, STEMP, R1, ROWMAX, T,
$ SFMIN
@@ -755,29 +755,11 @@
*
* A11 := A11 - U12*D*U12**H = A11 - U12*W**H
*
* computing blocks of NB columns at a time (note that conjg(W) is
* actually stored)
* (note that conjg(W) is actually stored)
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
A( JJ, JJ ) = REAL( A( JJ, JJ ) )
CALL CGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
$ A( J, JJ ), 1 )
A( JJ, JJ ) = REAL( A( JJ, JJ ) )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
IF( J.GE.2 )
$ CALL CGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
$ CONE, A( 1, J ), LDA )
50 CONTINUE
CALL CGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ CONE, A( 1, 1 ), LDA )
*
* Set KB to the number of columns factorized
*
@@ -1203,29 +1185,11 @@
*
* A22 := A22 - L21*D*L21**H = A22 - L21*W**H
*
* computing blocks of NB columns at a time (note that conjg(W) is
* actually stored)
* (note that conjg(W) is actually stored)
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
A( JJ, JJ ) = REAL( A( JJ, JJ ) )
CALL CGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
$ A( JJ, JJ ), 1 )
A( JJ, JJ ) = REAL( A( JJ, JJ ) )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL CGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
$ LDW, CONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL CGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -CONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ CONE, A( K, K ), LDA )
*
* Set KB to the number of columns factorized
*
+8 -41
View File
@@ -200,7 +200,7 @@
PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) )
* ..
* .. Local Scalars ..
INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
INTEGER IMAX, J, JJ, JMAX, JP, K, KK, KKW, KP,
$ KSTEP, KW
REAL ABSAKK, ALPHA, COLMAX, ROWMAX
COMPLEX D11, D21, D22, R1, T, Z
@@ -211,7 +211,7 @@
EXTERNAL LSAME, ICAMAX
* ..
* .. External Subroutines ..
EXTERNAL CCOPY, CGEMM, CGEMV, CSCAL, CSWAP
EXTERNAL CCOPY, CGEMMTR, CGEMV, CSCAL, CSWAP
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, AIMAG, MAX, MIN, REAL, SQRT
@@ -482,25 +482,9 @@
*
* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
*
* computing blocks of NB columns at a time
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
CALL CGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
$ A( J, JJ ), 1 )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
CALL CGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
$ CONE, A( 1, J ), LDA )
50 CONTINUE
CALL CGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ CONE, A( 1, 1 ), LDA )
*
* Put U12 in standard form by partially undoing the interchanges
* in columns k+1:n looping backwards from k+1 to n
@@ -778,26 +762,9 @@
*
* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
*
* computing blocks of NB columns at a time
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
CALL CGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
$ A( JJ, JJ ), 1 )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL CGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
$ LDW, CONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL CGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -CONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ CONE, A( K, K ), LDA )
*
* Put L21 in standard form by partially undoing the interchanges
* of rows in columns 1:k-1 looping backwards from k-1 to 1
+7 -41
View File
@@ -298,7 +298,7 @@
EXTERNAL LSAME, ICAMAX, SLAMCH
* ..
* .. External Subroutines ..
EXTERNAL CCOPY, CGEMM, CGEMV, CSCAL, CSWAP
EXTERNAL CCOPY, CGEMMTR, CGEMV, CSCAL, CSWAP
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, AIMAG, MAX, MIN, REAL, SQRT
@@ -627,26 +627,9 @@
*
* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
*
* computing blocks of NB columns at a time
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
CALL CGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
$ A( J, JJ ), 1 )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
IF( J.GE.2 )
$ CALL CGEMM( 'No transpose', 'Transpose', J-1, JB,
$ N-K, -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ),
$ LDW, CONE, A( 1, J ), LDA )
50 CONTINUE
CALL CGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ CONE, A( 1, 1 ), LDA )
*
* Set KB to the number of columns factorized
*
@@ -945,26 +928,9 @@
*
* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
*
* computing blocks of NB columns at a time
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
CALL CGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
$ A( JJ, JJ ), 1 )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL CGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
$ LDW, CONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL CGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -CONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ CONE, A( K, K ), LDA )
*
* Set KB to the number of columns factorized
*
+10 -40
View File
@@ -208,7 +208,7 @@
* ..
* .. Local Scalars ..
LOGICAL DONE
INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, JP1, JP2, K, KK,
INTEGER IMAX, ITEMP, J, JJ, JMAX, JP1, JP2, K, KK,
$ KW, KKW, KP, KSTEP, P, II
REAL ABSAKK, ALPHA, COLMAX, ROWMAX, STEMP, SFMIN
COMPLEX D11, D12, D21, D22, R1, T, Z
@@ -220,7 +220,7 @@
EXTERNAL LSAME, ICAMAX, SLAMCH
* ..
* .. External Subroutines ..
EXTERNAL CCOPY, CGEMM, CGEMV, CSCAL, CSWAP
EXTERNAL CCOPY, CGEMMTR, CGEMV, CSCAL, CSWAP
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT, AIMAG, REAL
@@ -525,26 +525,11 @@
*
* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
*
* computing blocks of NB columns at a time
* (note that conjg(W) is actually stored)
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
CALL CGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
$ A( J, JJ ), 1 )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
IF( J.GE.2 )
$ CALL CGEMM( 'No transpose', 'Transpose', J-1, JB,
$ N-K, -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
$ CONE, A( 1, J ), LDA )
50 CONTINUE
CALL CGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ CONE, A( 1, 1 ), LDA )
*
* Put U12 in standard form by partially undoing the interchanges
* in columns k+1:n
@@ -846,26 +831,11 @@
*
* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
*
* computing blocks of NB columns at a time
* (note that conjg(W) is actually stored)
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
CALL CGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
$ A( JJ, JJ ), 1 )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL CGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW,
$ CONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL CGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -CONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ CONE, A( K, K ), LDA )
*
* Put L21 in standard form by partially undoing the interchanges
* in columns 1:k-1
+8 -42
View File
@@ -197,7 +197,7 @@
PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
* ..
* .. Local Scalars ..
INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
INTEGER IMAX, J, JJ, JMAX, JP, K, KK, KKW, KP,
$ KSTEP, KW
DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D21, D22, R1,
$ ROWMAX, T
@@ -208,7 +208,7 @@
EXTERNAL LSAME, IDAMAX
* ..
* .. External Subroutines ..
EXTERNAL DCOPY, DGEMM, DGEMV, DSCAL, DSWAP
EXTERNAL DCOPY, DGEMMTR, DGEMV, DSCAL, DSWAP
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT
@@ -474,26 +474,9 @@
*
* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
*
* computing blocks of NB columns at a time
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
CALL DGEMV( 'No transpose', JJ-J+1, N-K, -ONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, ONE,
$ A( J, JJ ), 1 )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
CALL DGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
$ -ONE,
$ A( 1, K+1 ), LDA, W( J, KW+1 ), LDW, ONE,
$ A( 1, J ), LDA )
50 CONTINUE
CALL DGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -ONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ ONE, A( 1, 1 ), LDA )
*
* Put U12 in standard form by partially undoing the interchanges
* in columns k+1:n looping backwards from k+1 to n
@@ -770,26 +753,9 @@
*
* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
*
* computing blocks of NB columns at a time
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
CALL DGEMV( 'No transpose', J+JB-JJ, K-1, -ONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, ONE,
$ A( JJ, JJ ), 1 )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -ONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW,
$ ONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL DGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -ONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ ONE, A( K, K ), LDA )
*
* Put L21 in standard form by partially undoing the interchanges
* of rows in columns 1:k-1 looping backwards from k-1 to 1
+8 -42
View File
@@ -283,7 +283,7 @@
* ..
* .. Local Scalars ..
LOGICAL DONE
INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, K, KK, KW, KKW,
INTEGER IMAX, ITEMP, J, JMAX, K, KK, KW, KKW,
$ KP, KSTEP, P, II
DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
$ DTEMP, R1, ROWMAX, T, SFMIN
@@ -295,7 +295,7 @@
EXTERNAL LSAME, IDAMAX, DLAMCH
* ..
* .. External Subroutines ..
EXTERNAL DCOPY, DGEMM, DGEMV, DSCAL, DSWAP
EXTERNAL DCOPY, DGEMMTR, DGEMV, DSCAL, DSWAP
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT
@@ -618,26 +618,9 @@
*
* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
*
* computing blocks of NB columns at a time
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
CALL DGEMV( 'No transpose', JJ-J+1, N-K, -ONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, ONE,
$ A( J, JJ ), 1 )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
IF( J.GE.2 )
$ CALL DGEMM( 'No transpose', 'Transpose', J-1, JB,
$ N-K, -ONE, A( 1, K+1 ), LDA, W( J, KW+1 ),
$ LDW, ONE, A( 1, J ), LDA )
50 CONTINUE
CALL DGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -ONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ ONE, A( 1, 1 ), LDA )
*
* Set KB to the number of columns factorized
*
@@ -936,26 +919,9 @@
*
* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
*
* computing blocks of NB columns at a time
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
CALL DGEMV( 'No transpose', J+JB-JJ, K-1, -ONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, ONE,
$ A( JJ, JJ ), 1 )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -ONE, A( J+JB, 1 ), LDA, W( J, 1 ),
$ LDW, ONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL DGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -ONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ ONE, A( K, K ), LDA )
*
* Set KB to the number of columns factorized
*
+8 -42
View File
@@ -205,7 +205,7 @@
* ..
* .. Local Scalars ..
LOGICAL DONE
INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, JP1, JP2, K, KK,
INTEGER IMAX, ITEMP, J, JJ, JMAX, JP1, JP2, K, KK,
$ KW, KKW, KP, KSTEP, P, II
DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
@@ -218,7 +218,7 @@
EXTERNAL LSAME, IDAMAX, DLAMCH
* ..
* .. External Subroutines ..
EXTERNAL DCOPY, DGEMM, DGEMV, DSCAL, DSWAP
EXTERNAL DCOPY, DGEMMTR, DGEMV, DSCAL, DSWAP
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT
@@ -517,26 +517,9 @@
*
* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
*
* computing blocks of NB columns at a time
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
CALL DGEMV( 'No transpose', JJ-J+1, N-K, -ONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, ONE,
$ A( J, JJ ), 1 )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
IF( J.GE.2 )
$ CALL DGEMM( 'No transpose', 'Transpose', J-1, JB,
$ N-K, -ONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
$ ONE, A( 1, J ), LDA )
50 CONTINUE
CALL DGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -ONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ ONE, A( 1, 1 ), LDA )
*
* Put U12 in standard form by partially undoing the interchanges
* in columns k+1:n
@@ -838,26 +821,9 @@
*
* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
*
* computing blocks of NB columns at a time
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
CALL DGEMV( 'No transpose', J+JB-JJ, K-1, -ONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, ONE,
$ A( JJ, JJ ), 1 )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -ONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW,
$ ONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL DGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -ONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ ONE, A( K, K ), LDA )
*
* Put L21 in standard form by partially undoing the interchanges
* in columns 1:k-1
+8 -42
View File
@@ -197,7 +197,7 @@
PARAMETER ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
* ..
* .. Local Scalars ..
INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
INTEGER IMAX, J, JJ, JMAX, JP, K, KK, KKW, KP,
$ KSTEP, KW
REAL ABSAKK, ALPHA, COLMAX, D11, D21, D22, R1,
$ ROWMAX, T
@@ -208,7 +208,7 @@
EXTERNAL LSAME, ISAMAX
* ..
* .. External Subroutines ..
EXTERNAL SCOPY, SGEMM, SGEMV, SSCAL, SSWAP
EXTERNAL SCOPY, SGEMMTR, SGEMV, SSCAL, SSWAP
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT
@@ -474,26 +474,9 @@
*
* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
*
* computing blocks of NB columns at a time
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
CALL SGEMV( 'No transpose', JJ-J+1, N-K, -ONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, ONE,
$ A( J, JJ ), 1 )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
CALL SGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
$ -ONE,
$ A( 1, K+1 ), LDA, W( J, KW+1 ), LDW, ONE,
$ A( 1, J ), LDA )
50 CONTINUE
CALL SGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -ONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ ONE, A( 1, 1 ), LDA )
*
* Put U12 in standard form by partially undoing the interchanges
* in columns k+1:n looping backwards from k+1 to n
@@ -770,26 +753,9 @@
*
* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
*
* computing blocks of NB columns at a time
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
CALL SGEMV( 'No transpose', J+JB-JJ, K-1, -ONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, ONE,
$ A( JJ, JJ ), 1 )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL SGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -ONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW,
$ ONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL SGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -ONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ ONE, A( K, K ), LDA )
*
* Put L21 in standard form by partially undoing the interchanges
* of rows in columns 1:k-1 looping backwards from k-1 to 1
+7 -41
View File
@@ -295,7 +295,7 @@
EXTERNAL LSAME, ISAMAX, SLAMCH
* ..
* .. External Subroutines ..
EXTERNAL SCOPY, SGEMM, SGEMV, SSCAL, SSWAP
EXTERNAL SCOPY, SGEMMTR, SGEMV, SSCAL, SSWAP
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT
@@ -618,26 +618,9 @@
*
* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
*
* computing blocks of NB columns at a time
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
CALL SGEMV( 'No transpose', JJ-J+1, N-K, -ONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, ONE,
$ A( J, JJ ), 1 )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
IF( J.GE.2 )
$ CALL SGEMM( 'No transpose', 'Transpose', J-1, JB,
$ N-K, -ONE, A( 1, K+1 ), LDA, W( J, KW+1 ),
$ LDW, ONE, A( 1, J ), LDA )
50 CONTINUE
CALL SGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -ONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ ONE, A( 1, 1 ), LDA )
*
* Set KB to the number of columns factorized
*
@@ -936,26 +919,9 @@
*
* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
*
* computing blocks of NB columns at a time
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
CALL SGEMV( 'No transpose', J+JB-JJ, K-1, -ONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, ONE,
$ A( JJ, JJ ), 1 )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL SGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -ONE, A( J+JB, 1 ), LDA, W( J, 1 ),
$ LDW, ONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL SGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -ONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ ONE, A( K, K ), LDA )
*
* Set KB to the number of columns factorized
*
+8 -42
View File
@@ -205,7 +205,7 @@
* ..
* .. Local Scalars ..
LOGICAL DONE
INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, JP1, JP2, K, KK,
INTEGER IMAX, ITEMP, J, JJ, JMAX, JP1, JP2, K, KK,
$ KW, KKW, KP, KSTEP, P, II
REAL ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
@@ -218,7 +218,7 @@
EXTERNAL LSAME, ISAMAX, SLAMCH
* ..
* .. External Subroutines ..
EXTERNAL SCOPY, SGEMM, SGEMV, SSCAL, SSWAP
EXTERNAL SCOPY, SGEMMTR, SGEMV, SSCAL, SSWAP
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT
@@ -517,26 +517,9 @@
*
* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
*
* computing blocks of NB columns at a time
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
CALL SGEMV( 'No transpose', JJ-J+1, N-K, -ONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, ONE,
$ A( J, JJ ), 1 )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
IF( J.GE.2 )
$ CALL SGEMM( 'No transpose', 'Transpose', J-1, JB,
$ N-K, -ONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
$ ONE, A( 1, J ), LDA )
50 CONTINUE
CALL SGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -ONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ ONE, A( 1, 1 ), LDA )
*
* Put U12 in standard form by partially undoing the interchanges
* in columns k+1:n
@@ -838,26 +821,9 @@
*
* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
*
* computing blocks of NB columns at a time
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
CALL SGEMV( 'No transpose', J+JB-JJ, K-1, -ONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, ONE,
$ A( JJ, JJ ), 1 )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL SGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -ONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW,
$ ONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL SGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -ONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ ONE, A( K, K ), LDA )
*
* Put L21 in standard form by partially undoing the interchanges
* in columns 1:k-1
+10 -45
View File
@@ -200,7 +200,7 @@
PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
* ..
* .. Local Scalars ..
INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
INTEGER IMAX, J, JJ, JMAX, JP, K, KK, KKW, KP,
$ KSTEP, KW
DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, R1, ROWMAX, T
COMPLEX*16 D11, D21, D22, Z
@@ -211,7 +211,7 @@
EXTERNAL LSAME, IZAMAX
* ..
* .. External Subroutines ..
EXTERNAL ZCOPY, ZDSCAL, ZGEMM, ZGEMV, ZLACGV,
EXTERNAL ZCOPY, ZDSCAL, ZGEMMTR, ZGEMV, ZLACGV,
$ ZSWAP
* ..
* .. Intrinsic Functions ..
@@ -551,28 +551,11 @@
*
* A11 := A11 - U12*D*U12**H = A11 - U12*W**H
*
* computing blocks of NB columns at a time (note that conjg(W) is
* actually stored)
* (note that conjg(W) is actually stored)
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
CALL ZGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
$ A( J, JJ ), 1 )
A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
CALL ZGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
$ CONE, A( 1, J ), LDA )
50 CONTINUE
CALL ZGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ CONE, A( 1, 1 ), LDA )
*
* Put U12 in standard form by partially undoing the interchanges
* in columns k+1:n looping backwards from k+1 to n
@@ -915,29 +898,11 @@
*
* A22 := A22 - L21*D*L21**H = A22 - L21*W**H
*
* computing blocks of NB columns at a time (note that conjg(W) is
* actually stored)
* (note that conjg(W) is actually stored)
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
CALL ZGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
$ A( JJ, JJ ), 1 )
A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL ZGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
$ LDW, CONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL ZGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -CONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ CONE, A( K, K ), LDA )
*
* Put L21 in standard form by partially undoing the interchanges
* of rows in columns 1:k-1 looping backwards from k-1 to 1
+10 -46
View File
@@ -287,7 +287,7 @@
* ..
* .. Local Scalars ..
LOGICAL DONE
INTEGER IMAX, ITEMP, II, J, JB, JJ, JMAX, K, KK, KKW,
INTEGER IMAX, ITEMP, II, J, JMAX, K, KK, KKW,
$ KP, KSTEP, KW, P
DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, DTEMP, R1, ROWMAX, T,
$ SFMIN
@@ -300,7 +300,7 @@
EXTERNAL LSAME, IZAMAX, DLAMCH
* ..
* .. External Subroutines ..
EXTERNAL ZCOPY, ZDSCAL, ZGEMM, ZGEMV, ZLACGV,
EXTERNAL ZCOPY, ZDSCAL, ZGEMMTR, ZGEMV, ZLACGV,
$ ZSWAP
* ..
* .. Intrinsic Functions ..
@@ -755,29 +755,11 @@
*
* A11 := A11 - U12*D*U12**H = A11 - U12*W**H
*
* computing blocks of NB columns at a time (note that conjg(W) is
* actually stored)
* (note that conjg(W) is actually stored)
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
CALL ZGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
$ A( J, JJ ), 1 )
A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
IF( J.GE.2 )
$ CALL ZGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
$ CONE, A( 1, J ), LDA )
50 CONTINUE
CALL ZGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ CONE, A( 1, 1 ), LDA )
*
* Set KB to the number of columns factorized
*
@@ -1203,29 +1185,11 @@
*
* A22 := A22 - L21*D*L21**H = A22 - L21*W**H
*
* computing blocks of NB columns at a time (note that conjg(W) is
* actually stored)
* (note that conjg(W) is actually stored)
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
CALL ZGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
$ A( JJ, JJ ), 1 )
A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL ZGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
$ LDW, CONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL ZGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -CONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ CONE, A( K, K ), LDA )
*
* Set KB to the number of columns factorized
*
+8 -41
View File
@@ -200,7 +200,7 @@
PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
* ..
* .. Local Scalars ..
INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
INTEGER IMAX, J, JJ, JMAX, JP, K, KK, KKW, KP,
$ KSTEP, KW
DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX
COMPLEX*16 D11, D21, D22, R1, T, Z
@@ -211,7 +211,7 @@
EXTERNAL LSAME, IZAMAX
* ..
* .. External Subroutines ..
EXTERNAL ZCOPY, ZGEMM, ZGEMV, ZSCAL, ZSWAP
EXTERNAL ZCOPY, ZGEMMTR, ZGEMV, ZSCAL, ZSWAP
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, DIMAG, MAX, MIN, SQRT
@@ -481,25 +481,9 @@
*
* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
*
* computing blocks of NB columns at a time
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
CALL ZGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
$ A( J, JJ ), 1 )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
CALL ZGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
$ CONE, A( 1, J ), LDA )
50 CONTINUE
CALL ZGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ CONE, A( 1, 1 ), LDA )
*
* Put U12 in standard form by partially undoing the interchanges
* in columns k+1:n looping backwards from k+1 to n
@@ -776,26 +760,9 @@
*
* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
*
* computing blocks of NB columns at a time
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
CALL ZGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
$ A( JJ, JJ ), 1 )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL ZGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
$ LDW, CONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL ZGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -CONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ CONE, A( K, K ), LDA )
*
* Put L21 in standard form by partially undoing the interchanges
* of rows in columns 1:k-1 looping backwards from k-1 to 1
+7 -41
View File
@@ -298,7 +298,7 @@
EXTERNAL LSAME, IZAMAX, DLAMCH
* ..
* .. External Subroutines ..
EXTERNAL ZCOPY, ZGEMM, ZGEMV, ZSCAL, ZSWAP
EXTERNAL ZCOPY, ZGEMMTR, ZGEMV, ZSCAL, ZSWAP
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, DIMAG, MAX, MIN, SQRT
@@ -627,26 +627,9 @@
*
* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
*
* computing blocks of NB columns at a time
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
CALL ZGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
$ A( J, JJ ), 1 )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
IF( J.GE.2 )
$ CALL ZGEMM( 'No transpose', 'Transpose', J-1, JB,
$ N-K, -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ),
$ LDW, CONE, A( 1, J ), LDA )
50 CONTINUE
CALL ZGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ CONE, A( 1, 1 ), LDA )
*
* Set KB to the number of columns factorized
*
@@ -945,26 +928,9 @@
*
* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
*
* computing blocks of NB columns at a time
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
CALL ZGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
$ A( JJ, JJ ), 1 )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL ZGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
$ LDW, CONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL ZGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -CONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ CONE, A( K, K ), LDA )
*
* Set KB to the number of columns factorized
*
+10 -40
View File
@@ -208,7 +208,7 @@
* ..
* .. Local Scalars ..
LOGICAL DONE
INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, JP1, JP2, K, KK,
INTEGER IMAX, ITEMP, J, JJ, JMAX, JP1, JP2, K, KK,
$ KW, KKW, KP, KSTEP, P, II
DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX, DTEMP, SFMIN
COMPLEX*16 D11, D12, D21, D22, R1, T, Z
@@ -220,7 +220,7 @@
EXTERNAL LSAME, IZAMAX, DLAMCH
* ..
* .. External Subroutines ..
EXTERNAL ZCOPY, ZGEMM, ZGEMV, ZSCAL, ZSWAP
EXTERNAL ZCOPY, ZGEMMTR, ZGEMV, ZSCAL, ZSWAP
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT, DIMAG, DBLE
@@ -525,26 +525,11 @@
*
* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
*
* computing blocks of NB columns at a time
* (note that conjg(W) is actually stored)
*
DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
JB = MIN( NB, K-J+1 )
*
* Update the upper triangle of the diagonal block
*
DO 40 JJ = J, J + JB - 1
CALL ZGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
$ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
$ A( J, JJ ), 1 )
40 CONTINUE
*
* Update the rectangular superdiagonal block
*
IF( J.GE.2 )
$ CALL ZGEMM( 'No transpose', 'Transpose', J-1, JB,
$ N-K, -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
$ CONE, A( 1, J ), LDA )
50 CONTINUE
CALL ZGEMMTR( 'Upper', 'No transpose', 'Transpose', K, N-K,
$ -CONE, A( 1, K+1 ), LDA, W( 1, KW+1 ), LDW,
$ CONE, A( 1, 1 ), LDA )
*
* Put U12 in standard form by partially undoing the interchanges
* in columns k+1:n
@@ -846,26 +831,11 @@
*
* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
*
* computing blocks of NB columns at a time
* (note that conjg(W) is actually stored)
*
DO 110 J = K, N, NB
JB = MIN( NB, N-J+1 )
*
* Update the lower triangle of the diagonal block
*
DO 100 JJ = J, J + JB - 1
CALL ZGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
$ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
$ A( JJ, JJ ), 1 )
100 CONTINUE
*
* Update the rectangular subdiagonal block
*
IF( J+JB.LE.N )
$ CALL ZGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
$ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW,
$ CONE, A( J+JB, J ), LDA )
110 CONTINUE
CALL ZGEMMTR( 'Lower', 'No transpose', 'Transpose', N-K+1,
$ K-1, -CONE, A( K, 1 ), LDA, W( K, 1 ), LDW,
$ CONE, A( K, K ), LDA )
*
* Put L21 in standard form by partially undoing the interchanges
* in columns 1:k-1