Skip to content

Recursive LAUUM #1402

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

Open
jprhyne wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Recursive LAUUM #1402

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8808d9a
adding a recursiving lauum while keeping the blocked variant around. …
Sep 10, 2026
5e9701f
updating ILAENV to have a nx computation for lauum_recursive. Due to …
Sep 10, 2026
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
8 changes: 4 additions & 4 deletions SRC/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,7 @@ set(SLASRC
slarz.f slarzb.f slarzt.f slasy2.f
slasyf.f slasyf_rook.f slasyf_rk.f slasyf_aa.f
slatbs.f slatdf.f slatps.f slatrd.f slatrs.f slatrs3.f slatrz.f
slauu2.f slauum.f sopgtr.f sopmtr.f sorg2l.f sorg2r.f
slauu2.f slauum.f slauum_blocked.f slauum_recursive.f sopgtr.f sopmtr.f sorg2l.f sorg2r.f
sorgbr.f sorghr.f sorgl2.f sorglq.f sorgql.f sorgqr.f sorgr2.f
sorgrq.f sorgtr.f sorgtsqr.f sorgtsqr_row.f sorm2l.f sorm2r.f sorm22.f
sormbr.f sormhr.f sorml2.f sormlq.f sormql.f sormqr.f sormr2.f
Expand Down Expand Up @@ -225,7 +225,7 @@ set(CLASRC
clarz.f clarzb.f clarzt.f clascl.f claset.f clasr.f classq.f90
claswp.f clasyf.f clasyf_rook.f clasyf_rk.f clasyf_aa.f
clatbs.f clatdf.f clatps.f clatrd.f clatrs.f clatrs3.f clatrz.f
clauu2.f clauum.f cpbcon.f cpbequ.f cpbrfs.f cpbstf.f cpbsv.f
clauu2.f clauum.f clauum_blocked.f clauum_recursive.f cpbcon.f cpbequ.f cpbrfs.f cpbstf.f cpbsv.f
cpbsvx.f cpbtf2.f cpbtrf.f cpbtrs.f cpocon.f cpoequ.f cporfs.f
cposv.f cposvx.f cpotf2.f cpotrf2.f cpotri.f cpstrf.f cpstf2.f
cppcon.f cppequ.f cpprfs.f cppsv.f cppsvx.f cpptrf.f cpptri.f cpptrs.f
Expand Down Expand Up @@ -314,7 +314,7 @@ set(DLASRC
dlarz.f dlarzb.f dlarzt.f dlaswp.f dlasy2.f
dlasyf.f dlasyf_rook.f dlasyf_rk.f dlasyf_aa.f
dlatbs.f dlatdf.f dlatps.f dlatrd.f dlatrs.f dlatrs3.f dlatrz.f dlauu2.f
dlauum.f dopgtr.f dopmtr.f dorg2l.f dorg2r.f
dlauum.f dlauum_blocked.f dlauum_recursive.f dopgtr.f dopmtr.f dorg2l.f dorg2r.f
dorgbr.f dorghr.f dorgl2.f dorglq.f dorgql.f dorgqr.f dorgr2.f
dorgrq.f dorgtr.f dorgtsqr.f dorgtsqr_row.f dorm2l.f dorm2r.f dorm22.f
dormbr.f dormhr.f dorml2.f dormlq.f dormql.f dormqr.f dormr2.f
Expand Down Expand Up @@ -426,7 +426,7 @@ set(ZLASRC
zlarz.f zlarzb.f zlarzt.f zlascl.f zlaset.f zlasr.f
zlassq.f90 zlaswp.f zlasyf.f zlasyf_rook.f zlasyf_rk.f zlasyf_aa.f
zlatbs.f zlatdf.f zlatps.f zlatrd.f zlatrs.f zlatrs3.f zlatrz.f zlauu2.f
zlauum.f zpbcon.f zpbequ.f zpbrfs.f zpbstf.f zpbsv.f
zlauum.f zlauum_blocked.f zlauum_recursive.f zpbcon.f zpbequ.f zpbrfs.f zpbstf.f zpbsv.f
zpbsvx.f zpbtf2.f zpbtrf.f zpbtrs.f zpocon.f zpoequ.f zporfs.f
zposv.f zposvx.f zpotf2.f zpotrf.f zpotrf2.f zpotri.f zpotrs.f zpstrf.f zpstf2.f
zppcon.f zppequ.f zpprfs.f zppsv.f zppsvx.f zpptrf.f zpptri.f zpptrs.f
Expand Down
8 changes: 4 additions & 4 deletions SRC/Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,7 @@ SLASRC = \
slarz.o slarzb.o slarzt.o slaswp.o slasy2.o slasyf.o slasyf_rook.o \
slasyf_rk.o \
slatbs.o slatdf.o slatps.o slatrd.o slatrs.o slatrs3.o slatrz.o \
slauu2.o slauum.o sopgtr.o sopmtr.o sorg2l.o sorg2r.o \
slauu2.o slauum.o slauum_blocked.o slauum_recursive.o sopgtr.o sopmtr.o sorg2l.o sorg2r.o \
sorgbr.o sorghr.o sorgl2.o sorglq.o sorgql.o sorgqr.o sorgr2.o \
sorgrq.o sorgtr.o sorgtsqr.o sorgtsqr_row.o sorm2l.o sorm2r.o sorm22.o \
sormbr.o sormhr.o sorml2.o sormlq.o sormql.o sormqr.o sormr2.o \
Expand Down Expand Up @@ -254,7 +254,7 @@ CLASRC = \
clarz.o clarzb.o clarzt.o clascl.o claset.o clasr.o classq.o \
claswp.o clasyf.o clasyf_rook.o clasyf_rk.o clasyf_aa.o \
clatbs.o clatdf.o clatps.o clatrd.o clatrs.o clatrs3.o clatrz.o \
clauu2.o clauum.o cpbcon.o cpbequ.o cpbrfs.o cpbstf.o cpbsv.o \
clauu2.o clauum.o clauum_blocked.o clauum_recursive.o cpbcon.o cpbequ.o cpbrfs.o cpbstf.o cpbsv.o \
cpbsvx.o cpbtf2.o cpbtrf.o cpbtrs.o cpocon.o cpoequ.o cporfs.o \
cposv.o cposvx.o cpotf2.o cpotri.o cpstrf.o cpstf2.o \
cppcon.o cppequ.o cpprfs.o cppsv.o cppsvx.o cpptrf.o cpptri.o cpptrs.o \
Expand Down Expand Up @@ -344,7 +344,7 @@ DLASRC = \
dlarz.o dlarzb.o dlarzt.o dlaswp.o dlasy2.o \
dlasyf.o dlasyf_rook.o dlasyf_rk.o \
dlatbs.o dlatdf.o dlatps.o dlatrd.o dlatrs.o dlatrs3.o dlatrz.o dlauu2.o \
dlauum.o dopgtr.o dopmtr.o dorg2l.o dorg2r.o \
dlauum.o dlauum_blocked.o dlauum_recursive.o dlauum_recursive.o dopgtr.o dopmtr.o dorg2l.o dorg2r.o \
dorgbr.o dorghr.o dorgl2.o dorglq.o dorgql.o dorgqr.o dorgr2.o \
dorgrq.o dorgtr.o dorgtsqr.o dorgtsqr_row.o dorm2l.o dorm2r.o dorm22.o \
dormbr.o dormhr.o dorml2.o dormlq.o dormql.o dormqr.o dormr2.o \
Expand Down Expand Up @@ -459,7 +459,7 @@ ZLASRC = \
zlarz.o zlarzb.o zlarzt.o zlascl.o zlaset.o zlasr.o \
zlassq.o zlaswp.o zlasyf.o zlasyf_rook.o zlasyf_rk.o zlasyf_aa.o \
zlatbs.o zlatdf.o zlatps.o zlatrd.o zlatrs.o zlatrs3.o zlatrz.o zlauu2.o \
zlauum.o zpbcon.o zpbequ.o zpbrfs.o zpbstf.o zpbsv.o \
zlauum.o zlauum_blocked.o zlauum_recursive.o zpbcon.o zpbequ.o zpbrfs.o zpbstf.o zpbsv.o \
zpbsvx.o zpbtf2.o zpbtrf.o zpbtrs.o zpocon.o zpoequ.o zporfs.o \
zposv.o zposvx.o zpotf2.o zpotrf.o zpotri.o zpotrs.o zpstrf.o zpstf2.o \
zppcon.o zppequ.o zpprfs.o zppsv.o zppsvx.o zpptrf.o zpptri.o zpptrs.o \
Expand Down
61 changes: 6 additions & 55 deletions SRC/clauum.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,62 +157,13 @@ SUBROUTINE CLAUUM( UPLO, N, A, LDA, INFO )
IF( N.EQ.0 )
$ RETURN
*
* Determine the block size for this environment.
*
NB = ILAENV( 1, 'CLAUUM', UPLO, N, -1, -1, -1 )
*
IF( NB.LE.1 .OR. NB.GE.N ) THEN
*
* Use unblocked code
*
CALL CLAUU2( UPLO, N, A, LDA, INFO )
ELSE
*
* Use blocked code
*
IF( UPPER ) THEN
*
* Compute the product U * U**H.
*
DO 10 I = 1, N, NB
IB = MIN( NB, N-I+1 )
CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
$ 'Non-unit', I-1, IB, CONE, A( I, I ), LDA,
$ A( 1, I ), LDA )
CALL CLAUU2( 'Upper', IB, A( I, I ), LDA, INFO )
IF( I+IB.LE.N ) THEN
CALL CGEMM( 'No transpose', 'Conjugate transpose',
$ I-1, IB, N-I-IB+1, CONE, A( 1, I+IB ),
$ LDA, A( I, I+IB ), LDA, CONE, A( 1, I ),
$ LDA )
CALL CHERK( 'Upper', 'No transpose', IB, N-I-IB+1,
$ ONE, A( I, I+IB ), LDA, ONE, A( I, I ),
$ LDA )
END IF
10 CONTINUE
ELSE
*
* Compute the product L**H * L.
*
DO 20 I = 1, N, NB
IB = MIN( NB, N-I+1 )
CALL CTRMM( 'Left', 'Lower', 'Conjugate transpose',
$ 'Non-unit', IB, I-1, CONE, A( I, I ), LDA,
$ A( I, 1 ), LDA )
CALL CLAUU2( 'Lower', IB, A( I, I ), LDA, INFO )
IF( I+IB.LE.N ) THEN
CALL CGEMM( 'Conjugate transpose', 'No transpose',
$ IB,
$ I-1, N-I-IB+1, CONE, A( I+IB, I ), LDA,
$ A( I+IB, 1 ), LDA, CONE, A( I, 1 ), LDA )
CALL CHERK( 'Lower', 'Conjugate transpose', IB,
$ N-I-IB+1, ONE, A( I+IB, I ), LDA, ONE,
$ A( I, I ), LDA )
END IF
20 CONTINUE
END IF
END IF
* Here we dispatch to whatever is more efficient in a particular environment
* We are defaulting to recursive, but if you want to use the blocked variant
* Comment out the line starting with `CALL CLAUUM_RECURSIVE...`
* and uncomment the line starting with `CALL CLAUUM_BLOCKED...`
*
CALL CLAUUM_RECURSIVE(UPLO, N, A, LDA, INFO)
* CALL CLAUUM_BLOCKED(UPLO, N, A, LDA, INFO)
RETURN
*
* End of CLAUUM
Expand Down
220 changes: 220 additions & 0 deletions SRC/clauum_blocked.f
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,220 @@
*> \brief \b CLAUUM_BLOCKED computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked algorithm).
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> Download CLAUUM_BLOCKED + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clauum.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clauum.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clauum.f">
*> [TXT]</a>
*
* Definition:
* ===========
*
* SUBROUTINE CLAUUM_BLOCKED( UPLO, N, A, LDA, INFO )
*
* .. Scalar Arguments ..
* CHARACTER UPLO
* INTEGER INFO, LDA, N
* ..
* .. Array Arguments ..
* COMPLEX A( LDA, * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> CLAUUM_BLOCKED computes the product U * U**H or L**H * L, where the triangular
*> factor U or L is stored in the upper or lower triangular part of
*> the array A.
*>
*> If UPLO = 'U' or 'u' then the upper triangle of the result is stored,
*> overwriting the factor U in A.
*> If UPLO = 'L' or 'l' then the lower triangle of the result is stored,
*> overwriting the factor L in A.
*>
*> This is the blocked form of the algorithm, calling Level 3 BLAS.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is CHARACTER*1
*> Specifies whether the triangular factor stored in the array A
*> is upper or lower triangular:
*> = 'U': Upper triangular
*> = 'L': Lower triangular
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The order of the triangular factor U or L. N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*> A is COMPLEX array, dimension (LDA,N)
*> On entry, the triangular factor U or L.
*> On exit, if UPLO = 'U', the upper triangle of A is
*> overwritten with the upper triangle of the product U * U**H;
*> if UPLO = 'L', the lower triangle of A is overwritten with
*> the lower triangle of the product L**H * L.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> The leading dimension of the array A. LDA >= max(1,N).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
*> = 0: successful exit
*> < 0: if INFO = -k, the k-th argument had an illegal value
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup lauum
*
* =====================================================================
SUBROUTINE CLAUUM_BLOCKED( UPLO, N, A, LDA, INFO )
IMPLICIT NONE
*
* -- LAPACK auxiliary routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
CHARACTER UPLO
INTEGER INFO, LDA, N
* ..
* .. Array Arguments ..
COMPLEX A( LDA, * )
* ..
*
* =====================================================================
*
* .. Parameters ..
REAL ONE
PARAMETER ( ONE = 1.0E+0 )
COMPLEX CONE
PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) )
* ..
* .. Local Scalars ..
LOGICAL UPPER
INTEGER I, IB, NB
* ..
* .. External Functions ..
LOGICAL LSAME
INTEGER ILAENV
EXTERNAL LSAME, ILAENV
* ..
* .. External Subroutines ..
EXTERNAL CGEMM, CHERK, CLAUU2, CTRMM, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
UPPER = LSAME( UPLO, 'U' )
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
INFO = -4
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'CLAUUM_BLOCKED', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( N.EQ.0 )
$ RETURN
*
* Determine the block size for this environment.
*
NB = ILAENV( 1, 'CLAUUM_BLOCKED', UPLO, N, -1, -1, -1 )
*
IF( NB.LE.1 .OR. NB.GE.N ) THEN
*
* Use unblocked code
*
CALL CLAUU2( UPLO, N, A, LDA, INFO )
ELSE
*
* Use blocked code
*
IF( UPPER ) THEN
*
* Compute the product U * U**H.
*
DO 10 I = 1, N, NB
IB = MIN( NB, N-I+1 )
CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
$ 'Non-unit', I-1, IB, CONE, A( I, I ), LDA,
$ A( 1, I ), LDA )
CALL CLAUU2( 'Upper', IB, A( I, I ), LDA, INFO )
IF( I+IB.LE.N ) THEN
CALL CGEMM( 'No transpose', 'Conjugate transpose',
$ I-1, IB, N-I-IB+1, CONE, A( 1, I+IB ),
$ LDA, A( I, I+IB ), LDA, CONE, A( 1, I ),
$ LDA )
CALL CHERK( 'Upper', 'No transpose', IB, N-I-IB+1,
$ ONE, A( I, I+IB ), LDA, ONE, A( I, I ),
$ LDA )
END IF
10 CONTINUE
ELSE
*
* Compute the product L**H * L.
*
DO 20 I = 1, N, NB
IB = MIN( NB, N-I+1 )
CALL CTRMM( 'Left', 'Lower', 'Conjugate transpose',
$ 'Non-unit', IB, I-1, CONE, A( I, I ), LDA,
$ A( I, 1 ), LDA )
CALL CLAUU2( 'Lower', IB, A( I, I ), LDA, INFO )
IF( I+IB.LE.N ) THEN
CALL CGEMM( 'Conjugate transpose', 'No transpose',
$ IB,
$ I-1, N-I-IB+1, CONE, A( I+IB, I ), LDA,
$ A( I+IB, 1 ), LDA, CONE, A( I, 1 ), LDA )
CALL CHERK( 'Lower', 'Conjugate transpose', IB,
$ N-I-IB+1, ONE, A( I+IB, I ), LDA, ONE,
$ A( I, I ), LDA )
END IF
20 CONTINUE
END IF
END IF
*
RETURN
*
* End of CLAUUM_BLOCKED
*
END
Loading
Loading