Skip to content

Add recursive implementation of ?LAUUM (Reference-LAPACK PR 1402) #6162

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
martin-frbg wants to merge 2 commits into
base: develop
Choose a base branch
from
+7,665 −1,493
Open

Add recursive implementation of ?LAUUM (Reference-LAPACK PR 1402) #6162

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
6fce7af
Add recursive implementation of ?LAUUM (Reference-LAPACK PR 1402)
martin-frbg Oct 23, 2026
93b64fe
fix tabbing
martin-frbg Oct 23, 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
34 changes: 20 additions & 14 deletions lapack-netlib/SRC/Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -185,7 +185,7 @@ SLASRC_O = \
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 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 @@ -235,7 +235,8 @@ SLASRC_O = \
ssbev_2stage.o ssbevx_2stage.o ssbevd_2stage.o ssygv_2stage.o \
sgesvdq.o slatrs3.o strsyl3.o sgelst.o sgedmd.o sgedmdq.o sgecxx.o \
sbdsvdmr3.o sbdsvdmr3_work.o sbdtgk.o sgesvdmr3.o slar1v_tgk.o \
slarrb_tgk.o slarrf_tgk.o slarrv_tgk.o
slarrb_tgk.o slarrf_tgk.o slarrv_tgk.o slaed0_br.o slaed0_merge_pair.o \
slaed4_br.o slaed7_br.o slaed8_br.o slaed9_br.o sstebr.o ssteqr_boundary_br.o
endif

ifneq "$(or $(BUILD_SINGLE),$(BUILD_DOUBLE))" ""
Expand Down Expand Up @@ -303,7 +304,7 @@ CLASRC_O = \
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 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 @@ -397,7 +398,7 @@ DLASRC_O = \
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 dlatrz.o dlauu2.o \
dlauum.o dopgtr.o dopmtr.o dorg2l.o dorg2r.o \
dlauum.o dlauum_blocked.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 @@ -448,7 +449,8 @@ DLASRC_O = \
dsbev_2stage.o dsbevx_2stage.o dsbevd_2stage.o dsygv_2stage.o \
dgesvdq.o dlatrs3.o dtrsyl3.o dgelst.o dgedmd.o dgedmdq.o dgecxx.o \
dbdsvdmr3.o dbdsvdmr3_work.o dbdtgk.o dgesvdmr3.o dlar1v_tgk.o \
dlarrb_tgk.o dlarrf_tgk.o dlarrv_tgk.o
dlarrb_tgk.o dlarrf_tgk.o dlarrv_tgk.o dlaed0_br.o dlaed0_merge_pair.o \
dlaed4_br.o dlaed7_br.o dlaed8_br.o dlaed9_br.o dstebr.o dsteqr_boundary_br.o
endif

ifdef USEXBLAS
Expand Down Expand Up @@ -514,7 +516,7 @@ ZLASRC_O = \
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 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 Expand Up @@ -610,26 +612,30 @@ DZL_AUX_OBJS = dlaed3.o
SLAPACKOBJS = \
sgetrf.o sgetrs.o spotrf.o sgetf2.o \
spotf2.o slaswp.o sgesv.o slauu2.o \
slauum.o strti2.o strtri.o strtrs.o \
ssymv.o ssyr.o sspmv.o sspr.o
slauum.o slauum_blocked.o slauum_recursive.o \
strti2.o strtri.o strtrs.o \
ssymv.o ssyr.o sspmv.o sspr.o

DLAPACKOBJS = \
dgetrf.o dgetrs.o dpotrf.o dgetf2.o \
dpotf2.o dlaswp.o dgesv.o dlauu2.o \
dlauum.o dtrti2.o dtrtri.o dtrtrs.o \
dsymv.o dsyr.o dspmv.o dspr.o
dlauum.o dlauum_blocked.o dlauum_recursive.o \
dtrti2.o dtrtri.o dtrtrs.o \
dsymv.o dsyr.o dspmv.o dspr.o

CLAPACKOBJS = \
cgetrf.o cgetrs.o cpotrf.o cgetf2.o \
cpotf2.o claswp.o cgesv.o clauu2.o \
clauum.o ctrti2.o ctrtri.o ctrtrs.o \
csymv.o csyr.o cspmv.o cspr.o
clauum.o clauum_blocked.o clauum_recursive.o \
ctrti2.o ctrtri.o ctrtrs.o \
csymv.o csyr.o cspmv.o cspr.o

ZLAPACKOBJS = \
zgetrf.o zgetrs.o zpotrf.o zgetf2.o \
zpotf2.o zlaswp.o zgesv.o zlauu2.o \
zlauum.o ztrti2.o ztrtri.o ztrtrs.o \
zsymv.o zsyr.o zspmv.o zspr.o
zlauum.o zlauum_blocked.o zlauum_recursive.o \
ztrti2.o ztrtri.o ztrtrs.o \
zsymv.o zsyr.o zspmv.o zspr.o

ALLAUX = $(filter-out $(ALL_AUX_OBJS),$(ALLAUX_O))
SCLAUX = $(filter-out $(SCL_AUX_OBJS),$(SCLAUX_O))
Expand Down
125 changes: 58 additions & 67 deletions lapack-netlib/SRC/clauu2.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -411,9 +411,6 @@ static inline void zdotu_(doublecomplex *z, integer *n_, doublecomplex *x, integ
-lf2c -lm (in that order)
*/




/* Table of constant values */

static complex c_b1 = {1.f,0.f};
Expand All @@ -427,7 +424,6 @@ nblocked algorithm). */
/* Online html documentation available at */
/* http://www.netlib.org/lapack/explore-html/ */

/* > \htmlonly */
/* > Download CLAUU2 + dependencies */
/* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clauu2.
f"> */
Expand All @@ -438,7 +434,6 @@ f"> */
/* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clauu2.
f"> */
/* > [TXT]</a> */
/* > \endhtmlonly */

/* Definition: */
/* =========== */
Expand Down Expand Up @@ -492,7 +487,9 @@ f"> */
/* > 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. */
/* > the lower triangle of the product L**H * L. Regardless, */
/* > the diagonal elements are explicitly set to be real by casting */
/* . to REAL before exit. */
/* > \endverbatim */
/* > */
/* > \param[in] LDA */
Expand All @@ -516,38 +513,32 @@ f"> */
/* > \author Univ. of Colorado Denver */
/* > \author NAG Ltd. */

/* > \date December 2016 */

/* > \ingroup complexOTHERauxiliary */
/* > \ingroup lauu2 */

/* ===================================================================== */
/* Subroutine */ void clauu2_(char *uplo, integer *n, complex *a, integer *lda,
integer *info)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3;
real r__1=0.;
complex q__1={0.,0.};
real r__1;
complex q__1;

/* Local variables */
integer i__;
extern /* Complex */ VOID cdotc_(complex *, integer *, complex *, integer
*, complex *, integer *);
integer i__, j;
complex aii, ajj;
extern /* Subroutine */ void cscal_(integer *, complex *, complex *,
integer *), cgemm_(char *, char *, integer *, integer *, integer *
, complex *, complex *, integer *, complex *, integer *, complex *
, complex *, integer *);
extern logical lsame_(char *, char *);
extern /* Subroutine */ void cgemv_(char *, integer *, integer *, complex *
, complex *, integer *, complex *, integer *, complex *, complex *
, integer *);
logical upper;
extern /* Subroutine */ void clacgv_(integer *, complex *, integer *),
csscal_(integer *, real *, complex *, integer *);
extern void xerbla_(char *, integer *, ftnlen);
real aii;
extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen);


/* -- LAPACK auxiliary routine (version 3.7.0) -- */
/* -- LAPACK auxiliary routine -- */
/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
/* December 2016 */


/* ===================================================================== */
Expand All @@ -557,7 +548,7 @@ f"> */

/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1 * 1;
a_offset = 1 + a_dim1;
a -= a_offset;

/* Function Body */
Expand Down Expand Up @@ -587,30 +578,30 @@ f"> */
/* Compute the product U * U**H. */

i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = i__ + i__ * a_dim1;
aii = a[i__2].r;
if (i__ < *n) {
i__2 = i__ + i__ * a_dim1;
i__3 = *n - i__;
cdotc_(&q__1, &i__3, &a[i__ + (i__ + 1) * a_dim1], lda, &a[
i__ + (i__ + 1) * a_dim1], lda);
r__1 = aii * aii + q__1.r;
a[i__2].r = r__1, a[i__2].i = 0.f;
i__2 = *n - i__;
clacgv_(&i__2, &a[i__ + (i__ + 1) * a_dim1], lda);
i__2 = i__ - 1;
i__3 = *n - i__;
q__1.r = aii, q__1.i = 0.f;
cgemv_("No transpose", &i__2, &i__3, &c_b1, &a[(i__ + 1) *
a_dim1 + 1], lda, &a[i__ + (i__ + 1) * a_dim1], lda, &
q__1, &a[i__ * a_dim1 + 1], &c__1);
i__2 = *n - i__;
clacgv_(&i__2, &a[i__ + (i__ + 1) * a_dim1], lda);
for (j = 1; j <= i__1; ++j) {
r_cnjg(&q__1, &a[j + j * a_dim1]);
ajj.r = q__1.r, ajj.i = q__1.i;
if (j == *n) {

/* We have this case separate just to ensure we never */
/* access A out of bounds when we do A(J,J+1) = A(N,N+1) */

cscal_(n, &ajj, &a[*n * a_dim1 + 1], &c__1);
} else {
csscal_(&i__, &aii, &a[i__ * a_dim1 + 1], &c__1);
i__2 = *n - j;
cgemm_("No Transpose", "Conjugate Transpose", &j, &c__1, &
i__2, &c_b1, &a[(j + 1) * a_dim1 + 1], lda, &a[j + (j
+ 1) * a_dim1], lda, &ajj, &a[j * a_dim1 + 1], lda);
}
/* L10: */

/* Since in exact arithmetic the diagonal will be real on output */
/* we explicitly cast the diagonal to real (zero out the */
/* imaginary component) */

i__2 = j + j * a_dim1;
i__3 = j + j * a_dim1;
r__1 = a[i__3].r;
a[i__2].r = r__1, a[i__2].i = 0.f;
}

} else {
Expand All @@ -619,29 +610,29 @@ f"> */

i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = i__ + i__ * a_dim1;
aii = a[i__2].r;
if (i__ < *n) {
i__2 = i__ + i__ * a_dim1;
i__3 = *n - i__;
cdotc_(&q__1, &i__3, &a[i__ + 1 + i__ * a_dim1], &c__1, &a[
i__ + 1 + i__ * a_dim1], &c__1);
r__1 = aii * aii + q__1.r;
a[i__2].r = r__1, a[i__2].i = 0.f;
i__2 = i__ - 1;
clacgv_(&i__2, &a[i__ + a_dim1], lda);
i__2 = *n - i__;
i__3 = i__ - 1;
q__1.r = aii, q__1.i = 0.f;
cgemv_("Conjugate transpose", &i__2, &i__3, &c_b1, &a[i__ + 1
+ a_dim1], lda, &a[i__ + 1 + i__ * a_dim1], &c__1, &
q__1, &a[i__ + a_dim1], lda);
i__2 = i__ - 1;
clacgv_(&i__2, &a[i__ + a_dim1], lda);
r_cnjg(&q__1, &a[i__ + i__ * a_dim1]);
aii.r = q__1.r, aii.i = q__1.i;
if (i__ == *n) {

/* We have this case separate just to ensure we never */
/* access A out of bounds when we do A(I+1,I) = A(N+1,N) */

cscal_(n, &aii, &a[*n + a_dim1], lda);
} else {
csscal_(&i__, &aii, &a[i__ + a_dim1], lda);
i__2 = *n - i__;
cgemm_("Conjugate Transpose", "No Transpose", &c__1, &i__, &
i__2, &c_b1, &a[i__ + 1 + i__ * a_dim1], lda, &a[i__
+ 1 + a_dim1], lda, &aii, &a[i__ + a_dim1], lda);
}
/* L20: */

/* Since in exact arithmetic the diagonal will be real on output */
/* we explicitly cast the diagonal to real (zero out the */
/* imaginary component) */

i__2 = i__ + i__ * a_dim1;
i__3 = i__ + i__ * a_dim1;
r__1 = a[i__3].r;
a[i__2].r = r__1, a[i__2].i = 0.f;
}
}

Expand Down
Loading
Loading