Skip to content

Use GEMMTR for SY/HE linear updates #1069

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Oct 26, 2024
Merged

Use GEMMTR for SY/HE linear updates #1069

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 10 additions & 45 deletions SRC/clahef.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -200,7 +200,7 @@ SUBROUTINE CLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
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
Expand All @@ -211,7 +211,7 @@ SUBROUTINE CLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
EXTERNAL LSAME, ICAMAX
* ..
* .. External Subroutines ..
EXTERNAL CCOPY, CGEMM, CGEMV, CLACGV, CSSCAL,
EXTERNAL CCOPY, CGEMMTR, CGEMV, CLACGV, CSSCAL,
$ CSWAP
* ..
* .. Intrinsic Functions ..
Expand Down Expand Up @@ -552,28 +552,11 @@ SUBROUTINE CLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
*
* 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
Expand Down Expand Up @@ -916,29 +899,11 @@ SUBROUTINE CLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
*
* 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)
*
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
* (note that conjg(W) is actually stored)
*
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
Expand Down
54 changes: 9 additions & 45 deletions SRC/clahef_rk.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -286,7 +286,7 @@ SUBROUTINE CLAHEF_RK( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
* ..
* .. 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
Expand Down Expand Up @@ -755,29 +755,11 @@ SUBROUTINE CLAHEF_RK( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
*
* 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
*
Expand Down Expand Up @@ -1203,29 +1185,11 @@ SUBROUTINE CLAHEF_RK( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
*
* 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)
*
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
* (note that conjg(W) is actually stored)
*
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
*
Expand Down
49 changes: 8 additions & 41 deletions SRC/clasyf.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -200,7 +200,7 @@ SUBROUTINE CLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
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
Expand All @@ -211,7 +211,7 @@ SUBROUTINE CLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
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
Expand Down Expand Up @@ -482,25 +482,9 @@ SUBROUTINE CLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
*
* 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
Expand Down Expand Up @@ -778,26 +762,9 @@ SUBROUTINE CLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
*
* 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
Expand Down
48 changes: 7 additions & 41 deletions SRC/clasyf_rk.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -298,7 +298,7 @@ SUBROUTINE CLASYF_RK( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
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
Expand Down Expand Up @@ -627,26 +627,9 @@ SUBROUTINE CLASYF_RK( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
*
* 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
*
Expand Down Expand Up @@ -945,26 +928,9 @@ SUBROUTINE CLASYF_RK( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
*
* 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
*
Expand Down
50 changes: 10 additions & 40 deletions SRC/clasyf_rook.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -208,7 +208,7 @@ SUBROUTINE CLASYF_ROOK( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
* ..
* .. 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
Expand All @@ -220,7 +220,7 @@ SUBROUTINE CLASYF_ROOK( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
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
Expand Down Expand Up @@ -525,26 +525,11 @@ SUBROUTINE CLASYF_ROOK( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
*
* 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
Expand Down Expand Up @@ -846,26 +831,11 @@ SUBROUTINE CLASYF_ROOK( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
*
* 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
* (note that conjg(W) is actually stored)
*
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
Expand Down
Loading
Loading