From 5a837cb48b031bc5c74e88208be66941dc9f3557 Mon Sep 17 00:00:00 2001 From: julie Date: Fri, 18 Jun 2010 22:48:07 +0000 Subject: [PATCH] Step 2 of xlarfp: add new routines and add same test code plus check R(i,i) is nonnegative and real --- SRC/Makefile | 30 +++--- SRC/cgeqr2p.f | 122 ++++++++++++++++++++++ SRC/cgeqrfp.f | 197 ++++++++++++++++++++++++++++++++++++ SRC/{clarfp.f => clarfgp.f} | 6 +- SRC/dgeqr2p.f | 122 ++++++++++++++++++++++ SRC/dgeqrfp.f | 197 ++++++++++++++++++++++++++++++++++++ SRC/{dlarfp.f => dlarfgp.f} | 6 +- SRC/sgeqr2p.f | 122 ++++++++++++++++++++++ SRC/sgeqrfp.f | 197 ++++++++++++++++++++++++++++++++++++ SRC/{slarfp.f => slarfgp.f} | 6 +- SRC/zgeqr2p.f | 122 ++++++++++++++++++++++ SRC/zgeqrfp.f | 197 ++++++++++++++++++++++++++++++++++++ SRC/{zlarfp.f => zlarfgp.f} | 6 +- TESTING/LIN/Makefile | 8 +- TESTING/LIN/alahd.f | 5 + TESTING/LIN/cchkqr.f | 21 +++- TESTING/LIN/cerrqr.f | 33 +++++- TESTING/LIN/cqrt01p.f | 155 ++++++++++++++++++++++++++++ TESTING/LIN/dchkqr.f | 24 ++++- TESTING/LIN/derrqr.f | 33 +++++- TESTING/LIN/dqrt01p.f | 155 ++++++++++++++++++++++++++++ TESTING/LIN/schkqr.f | 19 +++- TESTING/LIN/serrqr.f | 33 +++++- TESTING/LIN/sqrt01p.f | 155 ++++++++++++++++++++++++++++ TESTING/LIN/zchkqr.f | 17 +++- TESTING/LIN/zerrqr.f | 33 +++++- TESTING/LIN/zqrt01p.f | 157 ++++++++++++++++++++++++++++ 27 files changed, 2125 insertions(+), 53 deletions(-) create mode 100644 SRC/cgeqr2p.f create mode 100644 SRC/cgeqrfp.f rename SRC/{clarfp.f => clarfgp.f} (97%) create mode 100644 SRC/dgeqr2p.f create mode 100644 SRC/dgeqrfp.f rename SRC/{dlarfp.f => dlarfgp.f} (97%) create mode 100644 SRC/sgeqr2p.f create mode 100644 SRC/sgeqrfp.f rename SRC/{slarfp.f => slarfgp.f} (97%) create mode 100644 SRC/zgeqr2p.f create mode 100644 SRC/zgeqrfp.f rename SRC/{zlarfp.f => zlarfgp.f} (97%) create mode 100644 TESTING/LIN/cqrt01p.f create mode 100644 TESTING/LIN/dqrt01p.f create mode 100644 TESTING/LIN/sqrt01p.f create mode 100644 TESTING/LIN/zqrt01p.f diff --git a/SRC/Makefile b/SRC/Makefile index 67315be6..a114d038 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -94,9 +94,9 @@ SLASRC = \ sgebrd.o sgecon.o sgeequ.o sgees.o sgeesx.o sgeev.o sgeevx.o \ sgegs.o sgegv.o sgehd2.o sgehrd.o sgelq2.o sgelqf.o \ sgels.o sgelsd.o sgelss.o sgelsx.o sgelsy.o sgeql2.o sgeqlf.o \ - sgeqp3.o sgeqpf.o sgeqr2.o sgeqrf.o sgerfs.o sgerq2.o sgerqf.o \ - sgesc2.o sgesdd.o sgesv.o sgesvd.o sgesvx.o sgetc2.o sgetf2.o \ - sgetrf.o sgetri.o \ + sgeqp3.o sgeqpf.o sgeqr2.o sgeqr2p.o sgeqrf.o sgeqrfp.o sgerfs.o \ + sgerq2.o sgerqf.o sgesc2.o sgesdd.o sgesv.o sgesvd.o sgesvx.o \ + sgetc2.o sgetf2.o sgetrf.o sgetri.o \ sgetrs.o sggbak.o sggbal.o sgges.o sggesx.o sggev.o sggevx.o \ sggglm.o sgghrd.o sgglse.o sggqrf.o \ sggrqf.o sggsvd.o sggsvp.o sgtcon.o sgtrfs.o sgtsv.o \ @@ -110,7 +110,7 @@ SLASRC = \ slaqgb.o slaqge.o slaqp2.o slaqps.o slaqsb.o slaqsp.o slaqsy.o \ slaqr0.o slaqr1.o slaqr2.o slaqr3.o slaqr4.o slaqr5.o \ slaqtr.o slar1v.o slar2v.o ilaslr.o ilaslc.o \ - slarf.o slarfb.o slarfg.o slarft.o slarfx.o slargv.o \ + slarf.o slarfb.o slarfgp.o slarft.o slarfx.o slargv.o \ slarrv.o slartv.o slarfp.o \ slarz.o slarzb.o slarzt.o slaswp.o slasy2.o slasyf.o \ slatbs.o slatdf.o slatps.o slatrd.o slatrs.o slatrz.o slatzm.o \ @@ -156,9 +156,9 @@ CLASRC = \ cgecon.o cgeequ.o cgees.o cgeesx.o cgeev.o cgeevx.o \ cgegs.o cgegv.o cgehd2.o cgehrd.o cgelq2.o cgelqf.o \ cgels.o cgelsd.o cgelss.o cgelsx.o cgelsy.o cgeql2.o cgeqlf.o cgeqp3.o \ - cgeqpf.o cgeqr2.o cgeqrf.o cgerfs.o cgerq2.o cgerqf.o \ - cgesc2.o cgesdd.o cgesv.o cgesvd.o cgesvx.o cgetc2.o cgetf2.o cgetrf.o \ - cgetri.o cgetrs.o \ + cgeqpf.o cgeqr2.o cgeqr2p.o cgeqrf.o cgeqrfp.o cgerfs.o \ + cgerq2.o cgerqf.o cgesc2.o cgesdd.o cgesv.o cgesvd.o \ + cgesvx.o cgetc2.o cgetf2.o cgetrf.o cgetri.o cgetrs.o \ cggbak.o cggbal.o cgges.o cggesx.o cggev.o cggevx.o cggglm.o \ cgghrd.o cgglse.o cggqrf.o cggrqf.o \ cggsvd.o cggsvp.o \ @@ -182,7 +182,7 @@ CLASRC = \ claqhb.o claqhe.o claqhp.o claqp2.o claqps.o claqsb.o \ claqr0.o claqr1.o claqr2.o claqr3.o claqr4.o claqr5.o \ claqsp.o claqsy.o clar1v.o clar2v.o ilaclr.o ilaclc.o \ - clarf.o clarfb.o clarfg.o clarft.o clarfp.o \ + clarf.o clarfb.o clarfg.o clarft.o clarfgp.o \ clarfx.o clargv.o clarnv.o clarrv.o clartg.o clartv.o \ clarz.o clarzb.o clarzt.o clascl.o claset.o clasr.o classq.o \ claswp.o clasyf.o clatbs.o clatdf.o clatps.o clatrd.o clatrs.o clatrz.o \ @@ -226,9 +226,9 @@ DLASRC = \ dgebrd.o dgecon.o dgeequ.o dgees.o dgeesx.o dgeev.o dgeevx.o \ dgegs.o dgegv.o dgehd2.o dgehrd.o dgelq2.o dgelqf.o \ dgels.o dgelsd.o dgelss.o dgelsx.o dgelsy.o dgeql2.o dgeqlf.o \ - dgeqp3.o dgeqpf.o dgeqr2.o dgeqrf.o dgerfs.o dgerq2.o dgerqf.o \ - dgesc2.o dgesdd.o dgesv.o dgesvd.o dgesvx.o dgetc2.o dgetf2.o \ - dgetrf.o dgetri.o \ + dgeqp3.o dgeqpf.o dgeqr2.o dgeqr2p.o dgeqrf.o dgeqrfp.o dgerfs.o \ + dgerq2.o dgerqf.o dgesc2.o dgesdd.o dgesv.o dgesvd.o dgesvx.o \ + dgetc2.o dgetf2.o dgetrf.o dgetri.o \ dgetrs.o dggbak.o dggbal.o dgges.o dggesx.o dggev.o dggevx.o \ dggglm.o dgghrd.o dgglse.o dggqrf.o \ dggrqf.o dggsvd.o dggsvp.o dgtcon.o dgtrfs.o dgtsv.o \ @@ -242,8 +242,8 @@ DLASRC = \ dlaqgb.o dlaqge.o dlaqp2.o dlaqps.o dlaqsb.o dlaqsp.o dlaqsy.o \ dlaqr0.o dlaqr1.o dlaqr2.o dlaqr3.o dlaqr4.o dlaqr5.o \ dlaqtr.o dlar1v.o dlar2v.o iladlr.o iladlc.o \ - dlarf.o dlarfb.o dlarfg.o dlarft.o dlarfx.o dlargv.o \ - dlarrv.o dlartv.o dlarfp.o \ + dlarf.o dlarfb.o dlarfg.o dlarfgp.o dlarft.o dlarfx.o \ + dlargv.o dlarrv.o dlartv.o \ dlarz.o dlarzb.o dlarzt.o dlaswp.o dlasy2.o dlasyf.o \ dlatbs.o dlatdf.o dlatps.o dlatrd.o dlatrs.o dlatrz.o dlatzm.o dlauu2.o \ dlauum.o dopgtr.o dopmtr.o dorg2l.o dorg2r.o \ @@ -290,7 +290,7 @@ ZLASRC = \ zgecon.o zgeequ.o zgees.o zgeesx.o zgeev.o zgeevx.o \ zgegs.o zgegv.o zgehd2.o zgehrd.o zgelq2.o zgelqf.o \ zgels.o zgelsd.o zgelss.o zgelsx.o zgelsy.o zgeql2.o zgeqlf.o zgeqp3.o \ - zgeqpf.o zgeqr2.o zgeqrf.o zgerfs.o zgerq2.o zgerqf.o \ + zgeqpf.o zgeqr2.o zgeqr2p.o zgeqrf.o zgeqrfp.o zgerfs.o zgerq2.o zgerqf.o \ zgesc2.o zgesdd.o zgesv.o zgesvd.o zgesvx.o zgetc2.o zgetf2.o zgetrf.o \ zgetri.o zgetrs.o \ zggbak.o zggbal.o zgges.o zggesx.o zggev.o zggevx.o zggglm.o \ @@ -318,7 +318,7 @@ ZLASRC = \ zlaqr0.o zlaqr1.o zlaqr2.o zlaqr3.o zlaqr4.o zlaqr5.o \ zlaqsp.o zlaqsy.o zlar1v.o zlar2v.o ilazlr.o ilazlc.o \ zlarcm.o zlarf.o zlarfb.o \ - zlarfg.o zlarft.o zlarfp.o \ + zlarfg.o zlarft.o zlarfgp.o \ zlarfx.o zlargv.o zlarnv.o zlarrv.o zlartg.o zlartv.o \ zlarz.o zlarzb.o zlarzt.o zlascl.o zlaset.o zlasr.o \ zlassq.o zlaswp.o zlasyf.o \ diff --git a/SRC/cgeqr2p.f b/SRC/cgeqr2p.f new file mode 100644 index 00000000..3ec2389d --- /dev/null +++ b/SRC/cgeqr2p.f @@ -0,0 +1,122 @@ + SUBROUTINE CGEQR2P( M, N, A, LDA, TAU, WORK, INFO ) +* +* -- LAPACK routine (version 3.2) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2006 +* +* .. Scalar Arguments .. + INTEGER INFO, LDA, M, N +* .. +* .. Array Arguments .. + COMPLEX A( LDA, * ), TAU( * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* CGEQR2P computes a QR factorization of a complex m by n matrix A: +* A = Q * R. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0. +* +* A (input/output) COMPLEX array, dimension (LDA,N) +* On entry, the m by n matrix A. +* On exit, the elements on and above the diagonal of the array +* contain the min(m,n) by n upper trapezoidal matrix R (R is +* upper triangular if m >= n); the elements below the diagonal, +* with the array TAU, represent the unitary matrix Q as a +* product of elementary reflectors (see Further Details). +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,M). +* +* TAU (output) COMPLEX array, dimension (min(M,N)) +* The scalar factors of the elementary reflectors (see Further +* Details). +* +* WORK (workspace) COMPLEX array, dimension (N) +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* Further Details +* =============== +* +* The matrix Q is represented as a product of elementary reflectors +* +* Q = H(1) H(2) . . . H(k), where k = min(m,n). +* +* Each H(i) has the form +* +* H(i) = I - tau * v * v' +* +* where tau is a complex scalar, and v is a complex vector with +* v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), +* and tau in TAU(i). +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX ONE + PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + INTEGER I, K + COMPLEX ALPHA +* .. +* .. External Subroutines .. + EXTERNAL CLARF, CLARFGP, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC CONJG, MAX, MIN +* .. +* .. Executable Statements .. +* +* Test the input arguments +* + INFO = 0 + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -4 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CGEQR2P', -INFO ) + RETURN + END IF +* + K = MIN( M, N ) +* + DO 10 I = 1, K +* +* Generate elementary reflector H(i) to annihilate A(i+1:m,i) +* + CALL CLARFGP( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1, + $ TAU( I ) ) + IF( I.LT.N ) THEN +* +* Apply H(i)' to A(i:m,i+1:n) from the left +* + ALPHA = A( I, I ) + A( I, I ) = ONE + CALL CLARF( 'Left', M-I+1, N-I, A( I, I ), 1, + $ CONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK ) + A( I, I ) = ALPHA + END IF + 10 CONTINUE + RETURN +* +* End of CGEQR2P +* + END diff --git a/SRC/cgeqrfp.f b/SRC/cgeqrfp.f new file mode 100644 index 00000000..6a857777 --- /dev/null +++ b/SRC/cgeqrfp.f @@ -0,0 +1,197 @@ + SUBROUTINE CGEQRFP( M, N, A, LDA, TAU, WORK, LWORK, INFO ) +* +* -- LAPACK routine (version 3.2) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2006 +* +* .. Scalar Arguments .. + INTEGER INFO, LDA, LWORK, M, N +* .. +* .. Array Arguments .. + COMPLEX A( LDA, * ), TAU( * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* CGEQRFP computes a QR factorization of a complex M-by-N matrix A: +* A = Q * R. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0. +* +* A (input/output) COMPLEX array, dimension (LDA,N) +* On entry, the M-by-N matrix A. +* On exit, the elements on and above the diagonal of the array +* contain the min(M,N)-by-N upper trapezoidal matrix R (R is +* upper triangular if m >= n); the elements below the diagonal, +* with the array TAU, represent the unitary matrix Q as a +* product of min(m,n) elementary reflectors (see Further +* Details). +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,M). +* +* TAU (output) COMPLEX array, dimension (min(M,N)) +* The scalar factors of the elementary reflectors (see Further +* Details). +* +* WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) +* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +* +* LWORK (input) INTEGER +* The dimension of the array WORK. LWORK >= max(1,N). +* For optimum performance LWORK >= N*NB, where NB is +* the optimal blocksize. +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* Further Details +* =============== +* +* The matrix Q is represented as a product of elementary reflectors +* +* Q = H(1) H(2) . . . H(k), where k = min(m,n). +* +* Each H(i) has the form +* +* H(i) = I - tau * v * v' +* +* where tau is a complex scalar, and v is a complex vector with +* v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), +* and tau in TAU(i). +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL LQUERY + INTEGER I, IB, IINFO, IWS, K, LDWORK, LWKOPT, NB, + $ NBMIN, NX +* .. +* .. External Subroutines .. + EXTERNAL CGEQR2P, CLARFB, CLARFT, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. External Functions .. + INTEGER ILAENV + EXTERNAL ILAENV +* .. +* .. Executable Statements .. +* +* Test the input arguments +* + INFO = 0 + NB = ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) + LWKOPT = N*NB + WORK( 1 ) = LWKOPT + LQUERY = ( LWORK.EQ.-1 ) + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -4 + ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN + INFO = -7 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CGEQRFP', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Quick return if possible +* + K = MIN( M, N ) + IF( K.EQ.0 ) THEN + WORK( 1 ) = 1 + RETURN + END IF +* + NBMIN = 2 + NX = 0 + IWS = N + IF( NB.GT.1 .AND. NB.LT.K ) THEN +* +* Determine when to cross over from blocked to unblocked code. +* + NX = MAX( 0, ILAENV( 3, 'CGEQRF', ' ', M, N, -1, -1 ) ) + IF( NX.LT.K ) THEN +* +* Determine if workspace is large enough for blocked code. +* + LDWORK = N + IWS = LDWORK*NB + IF( LWORK.LT.IWS ) THEN +* +* Not enough workspace to use optimal NB: reduce NB and +* determine the minimum value of NB. +* + NB = LWORK / LDWORK + NBMIN = MAX( 2, ILAENV( 2, 'CGEQRF', ' ', M, N, -1, + $ -1 ) ) + END IF + END IF + END IF +* + IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN +* +* Use blocked code initially +* + DO 10 I = 1, K - NX, NB + IB = MIN( K-I+1, NB ) +* +* Compute the QR factorization of the current block +* A(i:m,i:i+ib-1) +* + CALL CGEQR2P( M-I+1, IB, A( I, I ), LDA, TAU( I ), WORK, + $ IINFO ) + IF( I+IB.LE.N ) THEN +* +* Form the triangular factor of the block reflector +* H = H(i) H(i+1) . . . H(i+ib-1) +* + CALL CLARFT( 'Forward', 'Columnwise', M-I+1, IB, + $ A( I, I ), LDA, TAU( I ), WORK, LDWORK ) +* +* Apply H' to A(i:m,i+ib:n) from the left +* + CALL CLARFB( 'Left', 'Conjugate transpose', 'Forward', + $ 'Columnwise', M-I+1, N-I-IB+1, IB, + $ A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ), + $ LDA, WORK( IB+1 ), LDWORK ) + END IF + 10 CONTINUE + ELSE + I = 1 + END IF +* +* Use unblocked code to factor the last or only block. +* + IF( I.LE.K ) + $ CALL CGEQR2P( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ), WORK, + $ IINFO ) +* + WORK( 1 ) = IWS + RETURN +* +* End of CGEQRFP +* + END diff --git a/SRC/clarfp.f b/SRC/clarfgp.f similarity index 97% rename from SRC/clarfp.f rename to SRC/clarfgp.f index 357773bf..9465ed22 100644 --- a/SRC/clarfp.f +++ b/SRC/clarfgp.f @@ -1,4 +1,4 @@ - SUBROUTINE CLARFP( N, ALPHA, X, INCX, TAU ) + SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU ) * * -- LAPACK auxiliary routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -16,7 +16,7 @@ * Purpose * ======= * -* CLARFP generates a complex elementary reflector H of order n, such +* CLARFGP generates a complex elementary reflector H of order n, such * that * * H' * ( alpha ) = ( beta ), H' * H = I. @@ -205,6 +205,6 @@ * RETURN * -* End of CLARFP +* End of CLARFGP * END diff --git a/SRC/dgeqr2p.f b/SRC/dgeqr2p.f new file mode 100644 index 00000000..0cc58949 --- /dev/null +++ b/SRC/dgeqr2p.f @@ -0,0 +1,122 @@ + SUBROUTINE DGEQR2P( M, N, A, LDA, TAU, WORK, INFO ) +* +* -- LAPACK routine (version 3.2) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2006 +* +* .. Scalar Arguments .. + INTEGER INFO, LDA, M, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* DGEQR2 computes a QR factorization of a real m by n matrix A: +* A = Q * R. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0. +* +* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +* On entry, the m by n matrix A. +* On exit, the elements on and above the diagonal of the array +* contain the min(m,n) by n upper trapezoidal matrix R (R is +* upper triangular if m >= n); the elements below the diagonal, +* with the array TAU, represent the orthogonal matrix Q as a +* product of elementary reflectors (see Further Details). +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,M). +* +* TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) +* The scalar factors of the elementary reflectors (see Further +* Details). +* +* WORK (workspace) DOUBLE PRECISION array, dimension (N) +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* Further Details +* =============== +* +* The matrix Q is represented as a product of elementary reflectors +* +* Q = H(1) H(2) . . . H(k), where k = min(m,n). +* +* Each H(i) has the form +* +* H(i) = I - tau * v * v' +* +* where tau is a real scalar, and v is a real vector with +* v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), +* and tau in TAU(i). +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + INTEGER I, K + DOUBLE PRECISION AII +* .. +* .. External Subroutines .. + EXTERNAL DLARF, DLARFGP, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* Test the input arguments +* + INFO = 0 + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -4 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DGEQR2P', -INFO ) + RETURN + END IF +* + K = MIN( M, N ) +* + DO 10 I = 1, K +* +* Generate elementary reflector H(i) to annihilate A(i+1:m,i) +* + CALL DLARFGP( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1, + $ TAU( I ) ) + IF( I.LT.N ) THEN +* +* Apply H(i) to A(i:m,i+1:n) from the left +* + AII = A( I, I ) + A( I, I ) = ONE + CALL DLARF( 'Left', M-I+1, N-I, A( I, I ), 1, TAU( I ), + $ A( I, I+1 ), LDA, WORK ) + A( I, I ) = AII + END IF + 10 CONTINUE + RETURN +* +* End of DGEQR2P +* + END diff --git a/SRC/dgeqrfp.f b/SRC/dgeqrfp.f new file mode 100644 index 00000000..99cb3d9b --- /dev/null +++ b/SRC/dgeqrfp.f @@ -0,0 +1,197 @@ + SUBROUTINE DGEQRFP( M, N, A, LDA, TAU, WORK, LWORK, INFO ) +* +* -- LAPACK routine (version 3.2) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2006 +* +* .. Scalar Arguments .. + INTEGER INFO, LDA, LWORK, M, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* DGEQRFP computes a QR factorization of a real M-by-N matrix A: +* A = Q * R. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0. +* +* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +* On entry, the M-by-N matrix A. +* On exit, the elements on and above the diagonal of the array +* contain the min(M,N)-by-N upper trapezoidal matrix R (R is +* upper triangular if m >= n); the elements below the diagonal, +* with the array TAU, represent the orthogonal matrix Q as a +* product of min(m,n) elementary reflectors (see Further +* Details). +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,M). +* +* TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) +* The scalar factors of the elementary reflectors (see Further +* Details). +* +* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) +* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +* +* LWORK (input) INTEGER +* The dimension of the array WORK. LWORK >= max(1,N). +* For optimum performance LWORK >= N*NB, where NB is +* the optimal blocksize. +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* Further Details +* =============== +* +* The matrix Q is represented as a product of elementary reflectors +* +* Q = H(1) H(2) . . . H(k), where k = min(m,n). +* +* Each H(i) has the form +* +* H(i) = I - tau * v * v' +* +* where tau is a real scalar, and v is a real vector with +* v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), +* and tau in TAU(i). +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL LQUERY + INTEGER I, IB, IINFO, IWS, K, LDWORK, LWKOPT, NB, + $ NBMIN, NX +* .. +* .. External Subroutines .. + EXTERNAL DGEQR2P, DLARFB, DLARFT, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. External Functions .. + INTEGER ILAENV + EXTERNAL ILAENV +* .. +* .. Executable Statements .. +* +* Test the input arguments +* + INFO = 0 + NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) + LWKOPT = N*NB + WORK( 1 ) = LWKOPT + LQUERY = ( LWORK.EQ.-1 ) + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -4 + ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN + INFO = -7 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DGEQRFP', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Quick return if possible +* + K = MIN( M, N ) + IF( K.EQ.0 ) THEN + WORK( 1 ) = 1 + RETURN + END IF +* + NBMIN = 2 + NX = 0 + IWS = N + IF( NB.GT.1 .AND. NB.LT.K ) THEN +* +* Determine when to cross over from blocked to unblocked code. +* + NX = MAX( 0, ILAENV( 3, 'DGEQRF', ' ', M, N, -1, -1 ) ) + IF( NX.LT.K ) THEN +* +* Determine if workspace is large enough for blocked code. +* + LDWORK = N + IWS = LDWORK*NB + IF( LWORK.LT.IWS ) THEN +* +* Not enough workspace to use optimal NB: reduce NB and +* determine the minimum value of NB. +* + NB = LWORK / LDWORK + NBMIN = MAX( 2, ILAENV( 2, 'DGEQRF', ' ', M, N, -1, + $ -1 ) ) + END IF + END IF + END IF +* + IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN +* +* Use blocked code initially +* + DO 10 I = 1, K - NX, NB + IB = MIN( K-I+1, NB ) +* +* Compute the QR factorization of the current block +* A(i:m,i:i+ib-1) +* + CALL DGEQR2P( M-I+1, IB, A( I, I ), LDA, TAU( I ), WORK, + $ IINFO ) + IF( I+IB.LE.N ) THEN +* +* Form the triangular factor of the block reflector +* H = H(i) H(i+1) . . . H(i+ib-1) +* + CALL DLARFT( 'Forward', 'Columnwise', M-I+1, IB, + $ A( I, I ), LDA, TAU( I ), WORK, LDWORK ) +* +* Apply H' to A(i:m,i+ib:n) from the left +* + CALL DLARFB( 'Left', 'Transpose', 'Forward', + $ 'Columnwise', M-I+1, N-I-IB+1, IB, + $ A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ), + $ LDA, WORK( IB+1 ), LDWORK ) + END IF + 10 CONTINUE + ELSE + I = 1 + END IF +* +* Use unblocked code to factor the last or only block. +* + IF( I.LE.K ) + $ CALL DGEQR2P( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ), WORK, + $ IINFO ) +* + WORK( 1 ) = IWS + RETURN +* +* End of DGEQRFP +* + END diff --git a/SRC/dlarfp.f b/SRC/dlarfgp.f similarity index 97% rename from SRC/dlarfp.f rename to SRC/dlarfgp.f index 06f7fae4..b0abe084 100644 --- a/SRC/dlarfp.f +++ b/SRC/dlarfgp.f @@ -1,4 +1,4 @@ - SUBROUTINE DLARFP( N, ALPHA, X, INCX, TAU ) + SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU ) * * -- LAPACK auxiliary routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -16,7 +16,7 @@ * Purpose * ======= * -* DLARFP generates a real elementary reflector H of order n, such +* DLARFGP generates a real elementary reflector H of order n, such * that * * H * ( alpha ) = ( beta ), H' * H = I. @@ -175,6 +175,6 @@ * RETURN * -* End of DLARFP +* End of DLARFGP * END diff --git a/SRC/sgeqr2p.f b/SRC/sgeqr2p.f new file mode 100644 index 00000000..dbe9136d --- /dev/null +++ b/SRC/sgeqr2p.f @@ -0,0 +1,122 @@ + SUBROUTINE SGEQR2P( M, N, A, LDA, TAU, WORK, INFO ) +* +* -- LAPACK routine (version 3.2) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2006 +* +* .. Scalar Arguments .. + INTEGER INFO, LDA, M, N +* .. +* .. Array Arguments .. + REAL A( LDA, * ), TAU( * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* SGEQR2P computes a QR factorization of a real m by n matrix A: +* A = Q * R. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0. +* +* A (input/output) REAL array, dimension (LDA,N) +* On entry, the m by n matrix A. +* On exit, the elements on and above the diagonal of the array +* contain the min(m,n) by n upper trapezoidal matrix R (R is +* upper triangular if m >= n); the elements below the diagonal, +* with the array TAU, represent the orthogonal matrix Q as a +* product of elementary reflectors (see Further Details). +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,M). +* +* TAU (output) REAL array, dimension (min(M,N)) +* The scalar factors of the elementary reflectors (see Further +* Details). +* +* WORK (workspace) REAL array, dimension (N) +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* Further Details +* =============== +* +* The matrix Q is represented as a product of elementary reflectors +* +* Q = H(1) H(2) . . . H(k), where k = min(m,n). +* +* Each H(i) has the form +* +* H(i) = I - tau * v * v' +* +* where tau is a real scalar, and v is a real vector with +* v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), +* and tau in TAU(i). +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE + PARAMETER ( ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + INTEGER I, K + REAL AII +* .. +* .. External Subroutines .. + EXTERNAL SLARF, SLARFGP, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* Test the input arguments +* + INFO = 0 + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -4 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SGEQR2P', -INFO ) + RETURN + END IF +* + K = MIN( M, N ) +* + DO 10 I = 1, K +* +* Generate elementary reflector H(i) to annihilate A(i+1:m,i) +* + CALL SLARFGP( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1, + $ TAU( I ) ) + IF( I.LT.N ) THEN +* +* Apply H(i) to A(i:m,i+1:n) from the left +* + AII = A( I, I ) + A( I, I ) = ONE + CALL SLARF( 'Left', M-I+1, N-I, A( I, I ), 1, TAU( I ), + $ A( I, I+1 ), LDA, WORK ) + A( I, I ) = AII + END IF + 10 CONTINUE + RETURN +* +* End of SGEQR2P +* + END diff --git a/SRC/sgeqrfp.f b/SRC/sgeqrfp.f new file mode 100644 index 00000000..d88ecb7a --- /dev/null +++ b/SRC/sgeqrfp.f @@ -0,0 +1,197 @@ + SUBROUTINE SGEQRFP( M, N, A, LDA, TAU, WORK, LWORK, INFO ) +* +* -- LAPACK routine (version 3.2) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2006 +* +* .. Scalar Arguments .. + INTEGER INFO, LDA, LWORK, M, N +* .. +* .. Array Arguments .. + REAL A( LDA, * ), TAU( * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* SGEQRFP computes a QR factorization of a real M-by-N matrix A: +* A = Q * R. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0. +* +* A (input/output) REAL array, dimension (LDA,N) +* On entry, the M-by-N matrix A. +* On exit, the elements on and above the diagonal of the array +* contain the min(M,N)-by-N upper trapezoidal matrix R (R is +* upper triangular if m >= n); the elements below the diagonal, +* with the array TAU, represent the orthogonal matrix Q as a +* product of min(m,n) elementary reflectors (see Further +* Details). +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,M). +* +* TAU (output) REAL array, dimension (min(M,N)) +* The scalar factors of the elementary reflectors (see Further +* Details). +* +* WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) +* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +* +* LWORK (input) INTEGER +* The dimension of the array WORK. LWORK >= max(1,N). +* For optimum performance LWORK >= N*NB, where NB is +* the optimal blocksize. +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* Further Details +* =============== +* +* The matrix Q is represented as a product of elementary reflectors +* +* Q = H(1) H(2) . . . H(k), where k = min(m,n). +* +* Each H(i) has the form +* +* H(i) = I - tau * v * v' +* +* where tau is a real scalar, and v is a real vector with +* v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), +* and tau in TAU(i). +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL LQUERY + INTEGER I, IB, IINFO, IWS, K, LDWORK, LWKOPT, NB, + $ NBMIN, NX +* .. +* .. External Subroutines .. + EXTERNAL SGEQR2P, SLARFB, SLARFT, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. External Functions .. + INTEGER ILAENV + EXTERNAL ILAENV +* .. +* .. Executable Statements .. +* +* Test the input arguments +* + INFO = 0 + NB = ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) + LWKOPT = N*NB + WORK( 1 ) = LWKOPT + LQUERY = ( LWORK.EQ.-1 ) + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -4 + ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN + INFO = -7 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SGEQRFP', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Quick return if possible +* + K = MIN( M, N ) + IF( K.EQ.0 ) THEN + WORK( 1 ) = 1 + RETURN + END IF +* + NBMIN = 2 + NX = 0 + IWS = N + IF( NB.GT.1 .AND. NB.LT.K ) THEN +* +* Determine when to cross over from blocked to unblocked code. +* + NX = MAX( 0, ILAENV( 3, 'SGEQRF', ' ', M, N, -1, -1 ) ) + IF( NX.LT.K ) THEN +* +* Determine if workspace is large enough for blocked code. +* + LDWORK = N + IWS = LDWORK*NB + IF( LWORK.LT.IWS ) THEN +* +* Not enough workspace to use optimal NB: reduce NB and +* determine the minimum value of NB. +* + NB = LWORK / LDWORK + NBMIN = MAX( 2, ILAENV( 2, 'SGEQRF', ' ', M, N, -1, + $ -1 ) ) + END IF + END IF + END IF +* + IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN +* +* Use blocked code initially +* + DO 10 I = 1, K - NX, NB + IB = MIN( K-I+1, NB ) +* +* Compute the QR factorization of the current block +* A(i:m,i:i+ib-1) +* + CALL SGEQR2P( M-I+1, IB, A( I, I ), LDA, TAU( I ), WORK, + $ IINFO ) + IF( I+IB.LE.N ) THEN +* +* Form the triangular factor of the block reflector +* H = H(i) H(i+1) . . . H(i+ib-1) +* + CALL SLARFT( 'Forward', 'Columnwise', M-I+1, IB, + $ A( I, I ), LDA, TAU( I ), WORK, LDWORK ) +* +* Apply H' to A(i:m,i+ib:n) from the left +* + CALL SLARFB( 'Left', 'Transpose', 'Forward', + $ 'Columnwise', M-I+1, N-I-IB+1, IB, + $ A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ), + $ LDA, WORK( IB+1 ), LDWORK ) + END IF + 10 CONTINUE + ELSE + I = 1 + END IF +* +* Use unblocked code to factor the last or only block. +* + IF( I.LE.K ) + $ CALL SGEQR2P( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ), WORK, + $ IINFO ) +* + WORK( 1 ) = IWS + RETURN +* +* End of SGEQRFP +* + END diff --git a/SRC/slarfp.f b/SRC/slarfgp.f similarity index 97% rename from SRC/slarfp.f rename to SRC/slarfgp.f index 9a0b4bd8..c74c79a8 100644 --- a/SRC/slarfp.f +++ b/SRC/slarfgp.f @@ -1,4 +1,4 @@ - SUBROUTINE SLARFP( N, ALPHA, X, INCX, TAU ) + SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU ) * * -- LAPACK auxiliary routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -16,7 +16,7 @@ * Purpose * ======= * -* SLARFP generates a real elementary reflector H of order n, such +* SLARFGP generates a real elementary reflector H of order n, such * that * * H * ( alpha ) = ( beta ), H' * H = I. @@ -175,6 +175,6 @@ * RETURN * -* End of SLARFP +* End of SLARFGP * END diff --git a/SRC/zgeqr2p.f b/SRC/zgeqr2p.f new file mode 100644 index 00000000..ba908b83 --- /dev/null +++ b/SRC/zgeqr2p.f @@ -0,0 +1,122 @@ + SUBROUTINE ZGEQR2P( M, N, A, LDA, TAU, WORK, INFO ) +* +* -- LAPACK routine (version 3.2) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2006 +* +* .. Scalar Arguments .. + INTEGER INFO, LDA, M, N +* .. +* .. Array Arguments .. + COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* ZGEQR2P computes a QR factorization of a complex m by n matrix A: +* A = Q * R. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0. +* +* A (input/output) COMPLEX*16 array, dimension (LDA,N) +* On entry, the m by n matrix A. +* On exit, the elements on and above the diagonal of the array +* contain the min(m,n) by n upper trapezoidal matrix R (R is +* upper triangular if m >= n); the elements below the diagonal, +* with the array TAU, represent the unitary matrix Q as a +* product of elementary reflectors (see Further Details). +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,M). +* +* TAU (output) COMPLEX*16 array, dimension (min(M,N)) +* The scalar factors of the elementary reflectors (see Further +* Details). +* +* WORK (workspace) COMPLEX*16 array, dimension (N) +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* Further Details +* =============== +* +* The matrix Q is represented as a product of elementary reflectors +* +* Q = H(1) H(2) . . . H(k), where k = min(m,n). +* +* Each H(i) has the form +* +* H(i) = I - tau * v * v' +* +* where tau is a complex scalar, and v is a complex vector with +* v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), +* and tau in TAU(i). +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX*16 ONE + PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + INTEGER I, K + COMPLEX*16 ALPHA +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZLARF, ZLARFGP +* .. +* .. Intrinsic Functions .. + INTRINSIC DCONJG, MAX, MIN +* .. +* .. Executable Statements .. +* +* Test the input arguments +* + INFO = 0 + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -4 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZGEQR2P', -INFO ) + RETURN + END IF +* + K = MIN( M, N ) +* + DO 10 I = 1, K +* +* Generate elementary reflector H(i) to annihilate A(i+1:m,i) +* + CALL ZLARFGP( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1, + $ TAU( I ) ) + IF( I.LT.N ) THEN +* +* Apply H(i)' to A(i:m,i+1:n) from the left +* + ALPHA = A( I, I ) + A( I, I ) = ONE + CALL ZLARF( 'Left', M-I+1, N-I, A( I, I ), 1, + $ DCONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK ) + A( I, I ) = ALPHA + END IF + 10 CONTINUE + RETURN +* +* End of ZGEQR2P +* + END diff --git a/SRC/zgeqrfp.f b/SRC/zgeqrfp.f new file mode 100644 index 00000000..23378e69 --- /dev/null +++ b/SRC/zgeqrfp.f @@ -0,0 +1,197 @@ + SUBROUTINE ZGEQRFP( M, N, A, LDA, TAU, WORK, LWORK, INFO ) +* +* -- LAPACK routine (version 3.2) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2006 +* +* .. Scalar Arguments .. + INTEGER INFO, LDA, LWORK, M, N +* .. +* .. Array Arguments .. + COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* ZGEQRFP computes a QR factorization of a complex M-by-N matrix A: +* A = Q * R. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0. +* +* A (input/output) COMPLEX*16 array, dimension (LDA,N) +* On entry, the M-by-N matrix A. +* On exit, the elements on and above the diagonal of the array +* contain the min(M,N)-by-N upper trapezoidal matrix R (R is +* upper triangular if m >= n); the elements below the diagonal, +* with the array TAU, represent the unitary matrix Q as a +* product of min(m,n) elementary reflectors (see Further +* Details). +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,M). +* +* TAU (output) COMPLEX*16 array, dimension (min(M,N)) +* The scalar factors of the elementary reflectors (see Further +* Details). +* +* WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) +* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +* +* LWORK (input) INTEGER +* The dimension of the array WORK. LWORK >= max(1,N). +* For optimum performance LWORK >= N*NB, where NB is +* the optimal blocksize. +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* Further Details +* =============== +* +* The matrix Q is represented as a product of elementary reflectors +* +* Q = H(1) H(2) . . . H(k), where k = min(m,n). +* +* Each H(i) has the form +* +* H(i) = I - tau * v * v' +* +* where tau is a complex scalar, and v is a complex vector with +* v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), +* and tau in TAU(i). +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL LQUERY + INTEGER I, IB, IINFO, IWS, K, LDWORK, LWKOPT, NB, + $ NBMIN, NX +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZGEQR2P, ZLARFB, ZLARFT +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. External Functions .. + INTEGER ILAENV + EXTERNAL ILAENV +* .. +* .. Executable Statements .. +* +* Test the input arguments +* + INFO = 0 + NB = ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) + LWKOPT = N*NB + WORK( 1 ) = LWKOPT + LQUERY = ( LWORK.EQ.-1 ) + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -4 + ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN + INFO = -7 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZGEQRFP', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Quick return if possible +* + K = MIN( M, N ) + IF( K.EQ.0 ) THEN + WORK( 1 ) = 1 + RETURN + END IF +* + NBMIN = 2 + NX = 0 + IWS = N + IF( NB.GT.1 .AND. NB.LT.K ) THEN +* +* Determine when to cross over from blocked to unblocked code. +* + NX = MAX( 0, ILAENV( 3, 'ZGEQRF', ' ', M, N, -1, -1 ) ) + IF( NX.LT.K ) THEN +* +* Determine if workspace is large enough for blocked code. +* + LDWORK = N + IWS = LDWORK*NB + IF( LWORK.LT.IWS ) THEN +* +* Not enough workspace to use optimal NB: reduce NB and +* determine the minimum value of NB. +* + NB = LWORK / LDWORK + NBMIN = MAX( 2, ILAENV( 2, 'ZGEQRF', ' ', M, N, -1, + $ -1 ) ) + END IF + END IF + END IF +* + IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN +* +* Use blocked code initially +* + DO 10 I = 1, K - NX, NB + IB = MIN( K-I+1, NB ) +* +* Compute the QR factorization of the current block +* A(i:m,i:i+ib-1) +* + CALL ZGEQR2P( M-I+1, IB, A( I, I ), LDA, TAU( I ), WORK, + $ IINFO ) + IF( I+IB.LE.N ) THEN +* +* Form the triangular factor of the block reflector +* H = H(i) H(i+1) . . . H(i+ib-1) +* + CALL ZLARFT( 'Forward', 'Columnwise', M-I+1, IB, + $ A( I, I ), LDA, TAU( I ), WORK, LDWORK ) +* +* Apply H' to A(i:m,i+ib:n) from the left +* + CALL ZLARFB( 'Left', 'Conjugate transpose', 'Forward', + $ 'Columnwise', M-I+1, N-I-IB+1, IB, + $ A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ), + $ LDA, WORK( IB+1 ), LDWORK ) + END IF + 10 CONTINUE + ELSE + I = 1 + END IF +* +* Use unblocked code to factor the last or only block. +* + IF( I.LE.K ) + $ CALL ZGEQR2P( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ), WORK, + $ IINFO ) +* + WORK( 1 ) = IWS + RETURN +* +* End of ZGEQRFP +* + END diff --git a/SRC/zlarfp.f b/SRC/zlarfgp.f similarity index 97% rename from SRC/zlarfp.f rename to SRC/zlarfgp.f index 72e83273..8c0508eb 100644 --- a/SRC/zlarfp.f +++ b/SRC/zlarfgp.f @@ -1,4 +1,4 @@ - SUBROUTINE ZLARFP( N, ALPHA, X, INCX, TAU ) + SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU ) * * -- LAPACK auxiliary routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -16,7 +16,7 @@ * Purpose * ======= * -* ZLARFP generates a complex elementary reflector H of order n, such +* ZLARFGP generates a complex elementary reflector H of order n, such * that * * H' * ( alpha ) = ( beta ), H' * H = I. @@ -205,6 +205,6 @@ * RETURN * -* End of ZLARFP +* End of ZLARFGP * END diff --git a/TESTING/LIN/Makefile b/TESTING/LIN/Makefile index b79c3ab6..38b522f4 100644 --- a/TESTING/LIN/Makefile +++ b/TESTING/LIN/Makefile @@ -67,7 +67,7 @@ SLINTST = schkaa.o \ spot02.o spot03.o spot05.o spst01.o sppt01.o \ sppt02.o sppt03.o sppt05.o sptt01.o sptt02.o \ sptt05.o sqlt01.o sqlt02.o sqlt03.o sqpt01.o \ - sqrt01.o sqrt02.o sqrt03.o sqrt11.o sqrt12.o \ + sqrt01.o sqrt01p.o sqrt02.o sqrt03.o sqrt11.o sqrt12.o \ sqrt13.o sqrt14.o sqrt15.o sqrt16.o sqrt17.o \ srqt01.o srqt02.o srqt03.o srzt01.o srzt02.o \ sspt01.o ssyt01.o \ @@ -106,7 +106,7 @@ CLINTST = cchkaa.o \ cpot01.o cpot02.o cpot03.o cpot05.o cpst01.o \ cppt01.o cppt02.o cppt03.o cppt05.o cptt01.o \ cptt02.o cptt05.o cqlt01.o cqlt02.o cqlt03.o \ - cqpt01.o cqrt01.o cqrt02.o cqrt03.o cqrt11.o \ + cqpt01.o cqrt01.o cqrt01p.o cqrt02.o cqrt03.o cqrt11.o \ cqrt12.o cqrt13.o cqrt14.o cqrt15.o cqrt16.o \ cqrt17.o crqt01.o crqt02.o crqt03.o crzt01.o crzt02.o \ csbmv.o cspt01.o \ @@ -144,7 +144,7 @@ DLINTST = dchkaa.o \ dpot02.o dpot03.o dpot05.o dpst01.o dppt01.o \ dppt02.o dppt03.o dppt05.o dptt01.o dptt02.o \ dptt05.o dqlt01.o dqlt02.o dqlt03.o dqpt01.o \ - dqrt01.o dqrt02.o dqrt03.o dqrt11.o dqrt12.o \ + dqrt01.o dqrt01p.o dqrt02.o dqrt03.o dqrt11.o dqrt12.o \ dqrt13.o dqrt14.o dqrt15.o dqrt16.o dqrt17.o \ drqt01.o drqt02.o drqt03.o drzt01.o drzt02.o \ dspt01.o dsyt01.o \ @@ -183,7 +183,7 @@ ZLINTST = zchkaa.o \ zpot01.o zpot02.o zpot03.o zpot05.o zpst01.o \ zppt01.o zppt02.o zppt03.o zppt05.o zptt01.o \ zptt02.o zptt05.o zqlt01.o zqlt02.o zqlt03.o \ - zqpt01.o zqrt01.o zqrt02.o zqrt03.o zqrt11.o \ + zqpt01.o zqrt01.o zqrt01p.o zqrt02.o zqrt03.o zqrt11.o \ zqrt12.o zqrt13.o zqrt14.o zqrt15.o zqrt16.o \ zqrt17.o zrqt01.o zrqt02.o zrqt03.o zrzt01.o zrzt02.o \ zsbmv.o zspt01.o \ diff --git a/TESTING/LIN/alahd.f b/TESTING/LIN/alahd.f index 8a5637ea..b9621be2 100644 --- a/TESTING/LIN/alahd.f +++ b/TESTING/LIN/alahd.f @@ -351,12 +351,14 @@ WRITE( IOUNIT, FMT = 9970 ) WRITE( IOUNIT, FMT = '( '' Test ratios:'' )' ) WRITE( IOUNIT, FMT = 9950 )1 + WRITE( IOUNIT, FMT = 6950 )8 WRITE( IOUNIT, FMT = 9946 )2 WRITE( IOUNIT, FMT = 9944 )3, 'M' WRITE( IOUNIT, FMT = 9943 )4, 'M' WRITE( IOUNIT, FMT = 9942 )5, 'M' WRITE( IOUNIT, FMT = 9941 )6, 'M' WRITE( IOUNIT, FMT = 9960 )7 + WRITE( IOUNIT, FMT = 6660 )9 WRITE( IOUNIT, FMT = '( '' Messages:'' )' ) * ELSE IF( LSAMEN( 2, P2, 'LQ' ) ) THEN @@ -726,6 +728,7 @@ $ '( N * norm(A) * norm(AINV) * EPS )' ) 9960 FORMAT( 3X, I2, ': norm( B - A * X ) / ', $ '( norm(A) * norm(X) * EPS )' ) + 6660 FORMAT( 3X, I2, ': diagonal is not non-negative') 9959 FORMAT( 3X, I2, ': norm( X - XACT ) / ', $ '( norm(XACT) * CNDNUM * EPS )' ) 9958 FORMAT( 3X, I2, ': norm( X - XACT ) / ', @@ -750,6 +753,8 @@ 9951 FORMAT( ' Test ratio for ', A, ':', / 3X, I2, $ ': norm( s*b - A*x ) / ( norm(A) * norm(x) * EPS )' ) 9950 FORMAT( 3X, I2, ': norm( R - Q'' * A ) / ( M * norm(A) * EPS )' ) + 6950 FORMAT( 3X, I2, ': norm( R - Q'' * A ) / ( M * norm(A) * EPS ) + $ [RFPG]' ) 9949 FORMAT( 3X, I2, ': norm( L - A * Q'' ) / ( N * norm(A) * EPS )' ) 9948 FORMAT( 3X, I2, ': norm( L - Q'' * A ) / ( M * norm(A) * EPS )' ) 9947 FORMAT( 3X, I2, ': norm( R - A * Q'' ) / ( N * norm(A) * EPS )' ) diff --git a/TESTING/LIN/cchkqr.f b/TESTING/LIN/cchkqr.f index 0a3f4acd..a7f2748f 100644 --- a/TESTING/LIN/cchkqr.f +++ b/TESTING/LIN/cchkqr.f @@ -103,7 +103,7 @@ * * .. Parameters .. INTEGER NTESTS - PARAMETER ( NTESTS = 7 ) + PARAMETER ( NTESTS = 9 ) INTEGER NTYPES PARAMETER ( NTYPES = 8 ) REAL ZERO @@ -121,10 +121,14 @@ INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) REAL RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL CGENND + EXTERNAL CGENND +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, CERRQR, CGEQRS, CGET02, - $ CLACPY, CLARHS, CLATB4, CLATMS, CQRT01, CQRT02, - $ CQRT03, XLAENV + $ CLACPY, CLARHS, CLATB4, CLATMS, CQRT01, CQRT01P, + $ CQRT02, CQRT03, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -240,7 +244,16 @@ * CALL CQRT01( M, N, A, AF, AQ, AR, LDA, TAU, $ WORK, LWORK, RWORK, RESULT( 1 ) ) - ELSE IF( M.GE.N ) THEN +* +* Test CGEQRFP +* + CALL CQRT01P( M, N, A, AF, AQ, AR, LDA, TAU, + $ WORK, LWORK, RWORK, RESULT( 8 ) ) + + IF( .NOT. CGENND( M, N, AF, LDA ) ) + $ RESULT( 9 ) = 2*THRESH + NT = NT + 1 + ELSE IF( M.GE.N ) THEN * * Test CUNGQR, using factorization * returned by CQRT01 diff --git a/TESTING/LIN/cerrqr.f b/TESTING/LIN/cerrqr.f index d2dbcae9..abd789b0 100644 --- a/TESTING/LIN/cerrqr.f +++ b/TESTING/LIN/cerrqr.f @@ -38,8 +38,8 @@ $ W( NMAX ), X( NMAX ) * .. * .. External Subroutines .. - EXTERNAL ALAESM, CGEQR2, CGEQRF, CGEQRS, CHKXER, CUNG2R, - $ CUNGQR, CUNM2R, CUNMQR + EXTERNAL ALAESM, CGEQR2, CGEQR2P, CGEQRF, CGEQRFP, CGEQRS, + $ CHKXER, CUNG2R, CUNGQR, CUNM2R, CUNMQR * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -89,6 +89,22 @@ CALL CGEQRF( 1, 2, A, 1, B, W, 1, INFO ) CALL CHKXER( 'CGEQRF', INFOT, NOUT, LERR, OK ) * +* CGEQRFP +* + SRNAMT = 'CGEQRFP' + INFOT = 1 + CALL CGEQRFP( -1, 0, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'CGEQRFP', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL CGEQRFP( 0, -1, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'CGEQRFP', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL CGEQRFP( 2, 1, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'CGEQRFP', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL CGEQRFP( 1, 2, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'CGEQRFP', INFOT, NOUT, LERR, OK ) +* * CGEQR2 * SRNAMT = 'CGEQR2' @@ -102,6 +118,19 @@ CALL CGEQR2( 2, 1, A, 1, B, W, INFO ) CALL CHKXER( 'CGEQR2', INFOT, NOUT, LERR, OK ) * +* CGEQR2P +* + SRNAMT = 'CGEQR2P' + INFOT = 1 + CALL CGEQR2P( -1, 0, A, 1, B, W, INFO ) + CALL CHKXER( 'CGEQR2P', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL CGEQR2P( 0, -1, A, 1, B, W, INFO ) + CALL CHKXER( 'CGEQR2P', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL CGEQR2P( 2, 1, A, 1, B, W, INFO ) + CALL CHKXER( 'CGEQR2P', INFOT, NOUT, LERR, OK ) +* * CGEQRS * SRNAMT = 'CGEQRS' diff --git a/TESTING/LIN/cqrt01p.f b/TESTING/LIN/cqrt01p.f new file mode 100644 index 00000000..397d494c --- /dev/null +++ b/TESTING/LIN/cqrt01p.f @@ -0,0 +1,155 @@ + SUBROUTINE CQRT01P( M, N, A, AF, Q, R, LDA, TAU, WORK, LWORK, + $ RWORK, RESULT ) +* +* -- LAPACK test routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER LDA, LWORK, M, N +* .. +* .. Array Arguments .. + REAL RESULT( * ), RWORK( * ) + COMPLEX A( LDA, * ), AF( LDA, * ), Q( LDA, * ), + $ R( LDA, * ), TAU( * ), WORK( LWORK ) +* .. +* +* Purpose +* ======= +* +* CQRT01P tests CGEQRFP, which computes the QR factorization of an m-by-n +* matrix A, and partially tests CUNGQR which forms the m-by-m +* orthogonal matrix Q. +* +* CQRT01P compares R with Q'*A, and checks that Q is orthogonal. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0. +* +* A (input) COMPLEX array, dimension (LDA,N) +* The m-by-n matrix A. +* +* AF (output) COMPLEX array, dimension (LDA,N) +* Details of the QR factorization of A, as returned by CGEQRFP. +* See CGEQRFP for further details. +* +* Q (output) COMPLEX array, dimension (LDA,M) +* The m-by-m orthogonal matrix Q. +* +* R (workspace) COMPLEX array, dimension (LDA,max(M,N)) +* +* LDA (input) INTEGER +* The leading dimension of the arrays A, AF, Q and R. +* LDA >= max(M,N). +* +* TAU (output) COMPLEX array, dimension (min(M,N)) +* The scalar factors of the elementary reflectors, as returned +* by CGEQRFP. +* +* WORK (workspace) COMPLEX array, dimension (LWORK) +* +* LWORK (input) INTEGER +* The dimension of the array WORK. +* +* RWORK (workspace) REAL array, dimension (M) +* +* RESULT (output) REAL array, dimension (2) +* The test ratios: +* RESULT(1) = norm( R - Q'*A ) / ( M * norm(A) * EPS ) +* RESULT(2) = norm( I - Q'*Q ) / ( M * EPS ) +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) + COMPLEX ROGUE + PARAMETER ( ROGUE = ( -1.0E+10, -1.0E+10 ) ) +* .. +* .. Local Scalars .. + INTEGER INFO, MINMN + REAL ANORM, EPS, RESID +* .. +* .. External Functions .. + REAL CLANGE, CLANSY, SLAMCH + EXTERNAL CLANGE, CLANSY, SLAMCH +* .. +* .. External Subroutines .. + EXTERNAL CGEMM, CGEQRFP, CHERK, CLACPY, CLASET, CUNGQR +* .. +* .. Intrinsic Functions .. + INTRINSIC CMPLX, MAX, MIN, REAL +* .. +* .. Scalars in Common .. + CHARACTER*32 SRNAMT +* .. +* .. Common blocks .. + COMMON / SRNAMC / SRNAMT +* .. +* .. Executable Statements .. +* + MINMN = MIN( M, N ) + EPS = SLAMCH( 'Epsilon' ) +* +* Copy the matrix A to the array AF. +* + CALL CLACPY( 'Full', M, N, A, LDA, AF, LDA ) +* +* Factorize the matrix A in the array AF. +* + SRNAMT = 'CGEQRFP' + CALL CGEQRFP( M, N, AF, LDA, TAU, WORK, LWORK, INFO ) +* +* Copy details of Q +* + CALL CLASET( 'Full', M, M, ROGUE, ROGUE, Q, LDA ) + CALL CLACPY( 'Lower', M-1, N, AF( 2, 1 ), LDA, Q( 2, 1 ), LDA ) +* +* Generate the m-by-m matrix Q +* + SRNAMT = 'CUNGQR' + CALL CUNGQR( M, M, MINMN, Q, LDA, TAU, WORK, LWORK, INFO ) +* +* Copy R +* + CALL CLASET( 'Full', M, N, CMPLX( ZERO ), CMPLX( ZERO ), R, LDA ) + CALL CLACPY( 'Upper', M, N, AF, LDA, R, LDA ) +* +* Compute R - Q'*A +* + CALL CGEMM( 'Conjugate transpose', 'No transpose', M, N, M, + $ CMPLX( -ONE ), Q, LDA, A, LDA, CMPLX( ONE ), R, LDA ) +* +* Compute norm( R - Q'*A ) / ( M * norm(A) * EPS ) . +* + ANORM = CLANGE( '1', M, N, A, LDA, RWORK ) + RESID = CLANGE( '1', M, N, R, LDA, RWORK ) + IF( ANORM.GT.ZERO ) THEN + RESULT( 1 ) = ( ( RESID / REAL( MAX( 1, M ) ) ) / ANORM ) / EPS + ELSE + RESULT( 1 ) = ZERO + END IF +* +* Compute I - Q'*Q +* + CALL CLASET( 'Full', M, M, CMPLX( ZERO ), CMPLX( ONE ), R, LDA ) + CALL CHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, Q, LDA, + $ ONE, R, LDA ) +* +* Compute norm( I - Q'*Q ) / ( M * EPS ) . +* + RESID = CLANSY( '1', 'Upper', M, R, LDA, RWORK ) +* + RESULT( 2 ) = ( RESID / REAL( MAX( 1, M ) ) ) / EPS +* + RETURN +* +* End of CQRT01P +* + END diff --git a/TESTING/LIN/dchkqr.f b/TESTING/LIN/dchkqr.f index bab6359e..10030b6a 100644 --- a/TESTING/LIN/dchkqr.f +++ b/TESTING/LIN/dchkqr.f @@ -103,7 +103,7 @@ * * .. Parameters .. INTEGER NTESTS - PARAMETER ( NTESTS = 7 ) + PARAMETER ( NTESTS = 9 ) INTEGER NTYPES PARAMETER ( NTYPES = 8 ) DOUBLE PRECISION ZERO @@ -121,10 +121,14 @@ INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) DOUBLE PRECISION RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL DGENND + EXTERNAL DGENND +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, DERRQR, DGEQRS, DGET02, - $ DLACPY, DLARHS, DLATB4, DLATMS, DQRT01, DQRT02, - $ DQRT03, XLAENV + $ DLACPY, DLARHS, DLATB4, DLATMS, DQRT01, DQRT01P, + $ DQRT02, DQRT03, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -240,7 +244,17 @@ * CALL DQRT01( M, N, A, AF, AQ, AR, LDA, TAU, $ WORK, LWORK, RWORK, RESULT( 1 ) ) - ELSE IF( M.GE.N ) THEN + +* +* Test DGEQRFP +* + CALL DQRT01P( M, N, A, AF, AQ, AR, LDA, TAU, + $ WORK, LWORK, RWORK, RESULT( 8 ) ) + + IF( .NOT. DGENND( M, N, AF, LDA ) ) + $ RESULT( 9 ) = 2*THRESH + NT = NT + 1 + ELSE IF( M.GE.N ) THEN * * Test DORGQR, using factorization * returned by DQRT01 @@ -295,7 +309,7 @@ * Print information about the tests that did not * pass the threshold. * - DO 20 I = 1, NT + DO 20 I = 1, NTESTS IF( RESULT( I ).GE.THRESH ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) diff --git a/TESTING/LIN/derrqr.f b/TESTING/LIN/derrqr.f index 9503f992..60d62864 100644 --- a/TESTING/LIN/derrqr.f +++ b/TESTING/LIN/derrqr.f @@ -38,8 +38,8 @@ $ W( NMAX ), X( NMAX ) * .. * .. External Subroutines .. - EXTERNAL ALAESM, CHKXER, DGEQR2, DGEQRF, DGEQRS, DORG2R, - $ DORGQR, DORM2R, DORMQR + EXTERNAL ALAESM, CHKXER, DGEQR2, DGEQR2P, DGEQRF, DGEQRFP, + $ DGEQRS, DORG2R, DORGQR, DORM2R, DORMQR * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -89,6 +89,22 @@ CALL DGEQRF( 1, 2, A, 1, B, W, 1, INFO ) CALL CHKXER( 'DGEQRF', INFOT, NOUT, LERR, OK ) * +* DGEQRFP +* + SRNAMT = 'DGEQRFP' + INFOT = 1 + CALL DGEQRFP( -1, 0, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'DGEQRFP', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL DGEQRFP( 0, -1, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'DGEQRFP', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL DGEQRFP( 2, 1, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'DGEQRFP', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL DGEQRFP( 1, 2, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'DGEQRFP', INFOT, NOUT, LERR, OK ) +* * DGEQR2 * SRNAMT = 'DGEQR2' @@ -102,6 +118,19 @@ CALL DGEQR2( 2, 1, A, 1, B, W, INFO ) CALL CHKXER( 'DGEQR2', INFOT, NOUT, LERR, OK ) * +* DGEQR2P +* + SRNAMT = 'DGEQR2P' + INFOT = 1 + CALL DGEQR2P( -1, 0, A, 1, B, W, INFO ) + CALL CHKXER( 'DGEQR2P', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL DGEQR2P( 0, -1, A, 1, B, W, INFO ) + CALL CHKXER( 'DGEQR2P', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL DGEQR2P( 2, 1, A, 1, B, W, INFO ) + CALL CHKXER( 'DGEQR2P', INFOT, NOUT, LERR, OK ) +* * DGEQRS * SRNAMT = 'DGEQRS' diff --git a/TESTING/LIN/dqrt01p.f b/TESTING/LIN/dqrt01p.f new file mode 100644 index 00000000..97dc68c4 --- /dev/null +++ b/TESTING/LIN/dqrt01p.f @@ -0,0 +1,155 @@ + SUBROUTINE DQRT01P( M, N, A, AF, Q, R, LDA, TAU, WORK, LWORK, + $ RWORK, RESULT ) +* +* -- LAPACK test routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER LDA, LWORK, M, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), AF( LDA, * ), Q( LDA, * ), + $ R( LDA, * ), RESULT( * ), RWORK( * ), TAU( * ), + $ WORK( LWORK ) +* .. +* +* Purpose +* ======= +* +* DQRT01P tests DGEQRFP, which computes the QR factorization of an m-by-n +* matrix A, and partially tests DORGQR which forms the m-by-m +* orthogonal matrix Q. +* +* DQRT01P compares R with Q'*A, and checks that Q is orthogonal. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0. +* +* A (input) DOUBLE PRECISION array, dimension (LDA,N) +* The m-by-n matrix A. +* +* AF (output) DOUBLE PRECISION array, dimension (LDA,N) +* Details of the QR factorization of A, as returned by DGEQRFP. +* See DGEQRFP for further details. +* +* Q (output) DOUBLE PRECISION array, dimension (LDA,M) +* The m-by-m orthogonal matrix Q. +* +* R (workspace) DOUBLE PRECISION array, dimension (LDA,max(M,N)) +* +* LDA (input) INTEGER +* The leading dimension of the arrays A, AF, Q and R. +* LDA >= max(M,N). +* +* TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) +* The scalar factors of the elementary reflectors, as returned +* by DGEQRFP. +* +* WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) +* +* LWORK (input) INTEGER +* The dimension of the array WORK. +* +* RWORK (workspace) DOUBLE PRECISION array, dimension (M) +* +* RESULT (output) DOUBLE PRECISION array, dimension (2) +* The test ratios: +* RESULT(1) = norm( R - Q'*A ) / ( M * norm(A) * EPS ) +* RESULT(2) = norm( I - Q'*Q ) / ( M * EPS ) +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) + DOUBLE PRECISION ROGUE + PARAMETER ( ROGUE = -1.0D+10 ) +* .. +* .. Local Scalars .. + INTEGER INFO, MINMN + DOUBLE PRECISION ANORM, EPS, RESID +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH, DLANGE, DLANSY + EXTERNAL DLAMCH, DLANGE, DLANSY +* .. +* .. External Subroutines .. + EXTERNAL DGEMM, DGEQRFP, DLACPY, DLASET, DORGQR, DSYRK +* .. +* .. Intrinsic Functions .. + INTRINSIC DBLE, MAX, MIN +* .. +* .. Scalars in Common .. + CHARACTER*32 SRNAMT +* .. +* .. Common blocks .. + COMMON / SRNAMC / SRNAMT +* .. +* .. Executable Statements .. +* + MINMN = MIN( M, N ) + EPS = DLAMCH( 'Epsilon' ) +* +* Copy the matrix A to the array AF. +* + CALL DLACPY( 'Full', M, N, A, LDA, AF, LDA ) +* +* Factorize the matrix A in the array AF. +* + SRNAMT = 'DGEQRFP' + CALL DGEQRFP( M, N, AF, LDA, TAU, WORK, LWORK, INFO ) +* +* Copy details of Q +* + CALL DLASET( 'Full', M, M, ROGUE, ROGUE, Q, LDA ) + CALL DLACPY( 'Lower', M-1, N, AF( 2, 1 ), LDA, Q( 2, 1 ), LDA ) +* +* Generate the m-by-m matrix Q +* + SRNAMT = 'DORGQR' + CALL DORGQR( M, M, MINMN, Q, LDA, TAU, WORK, LWORK, INFO ) +* +* Copy R +* + CALL DLASET( 'Full', M, N, ZERO, ZERO, R, LDA ) + CALL DLACPY( 'Upper', M, N, AF, LDA, R, LDA ) +* +* Compute R - Q'*A +* + CALL DGEMM( 'Transpose', 'No transpose', M, N, M, -ONE, Q, LDA, A, + $ LDA, ONE, R, LDA ) +* +* Compute norm( R - Q'*A ) / ( M * norm(A) * EPS ) . +* + ANORM = DLANGE( '1', M, N, A, LDA, RWORK ) + RESID = DLANGE( '1', M, N, R, LDA, RWORK ) + IF( ANORM.GT.ZERO ) THEN + RESULT( 1 ) = ( ( RESID / DBLE( MAX( 1, M ) ) ) / ANORM ) / EPS + ELSE + RESULT( 1 ) = ZERO + END IF +* +* Compute I - Q'*Q +* + CALL DLASET( 'Full', M, M, ZERO, ONE, R, LDA ) + CALL DSYRK( 'Upper', 'Transpose', M, M, -ONE, Q, LDA, ONE, R, + $ LDA ) +* +* Compute norm( I - Q'*Q ) / ( M * EPS ) . +* + RESID = DLANSY( '1', 'Upper', M, R, LDA, RWORK ) +* + RESULT( 2 ) = ( RESID / DBLE( MAX( 1, M ) ) ) / EPS +* + RETURN +* +* End of DQRT01P +* + END diff --git a/TESTING/LIN/schkqr.f b/TESTING/LIN/schkqr.f index c97ea392..36f08a84 100644 --- a/TESTING/LIN/schkqr.f +++ b/TESTING/LIN/schkqr.f @@ -103,7 +103,7 @@ * * .. Parameters .. INTEGER NTESTS - PARAMETER ( NTESTS = 7 ) + PARAMETER ( NTESTS = 9 ) INTEGER NTYPES PARAMETER ( NTYPES = 8 ) REAL ZERO @@ -121,10 +121,14 @@ INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) REAL RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL SGENND + EXTERNAL SGENND +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, SERRQR, SGEQRS, SGET02, - $ SLACPY, SLARHS, SLATB4, SLATMS, SQRT01, SQRT02, - $ SQRT03, XLAENV + $ SLACPY, SLARHS, SLATB4, SLATMS, SQRT01, SQRT01P, + $ SQRT02, SQRT03, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -240,6 +244,15 @@ * CALL SQRT01( M, N, A, AF, AQ, AR, LDA, TAU, $ WORK, LWORK, RWORK, RESULT( 1 ) ) +* +* Test SGEQRFP +* + CALL SQRT01P( M, N, A, AF, AQ, AR, LDA, TAU, + $ WORK, LWORK, RWORK, RESULT( 8 ) ) + + IF( .NOT. SGENND( M, N, AF, LDA ) ) + $ RESULT( 9 ) = 2*THRESH + NT = NT + 1 ELSE IF( M.GE.N ) THEN * * Test SORGQR, using factorization diff --git a/TESTING/LIN/serrqr.f b/TESTING/LIN/serrqr.f index b1c7eaba..1ae93a8e 100644 --- a/TESTING/LIN/serrqr.f +++ b/TESTING/LIN/serrqr.f @@ -38,8 +38,8 @@ $ W( NMAX ), X( NMAX ) * .. * .. External Subroutines .. - EXTERNAL ALAESM, CHKXER, SGEQR2, SGEQRF, SGEQRS, SORG2R, - $ SORGQR, SORM2R, SORMQR + EXTERNAL ALAESM, CHKXER, SGEQR2, SGEQR2P, SGEQRF, SGEQRFP, + $ SGEQRS, SORG2R, SORGQR, SORM2R, SORMQR * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -89,6 +89,22 @@ CALL SGEQRF( 1, 2, A, 1, B, W, 1, INFO ) CALL CHKXER( 'SGEQRF', INFOT, NOUT, LERR, OK ) * +* SGEQRFP +* + SRNAMT = 'SGEQRFP' + INFOT = 1 + CALL SGEQRFP( -1, 0, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'SGEQRFP', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL SGEQRFP( 0, -1, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'SGEQRFP', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL SGEQRFP( 2, 1, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'SGEQRFP', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL SGEQRFP( 1, 2, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'SGEQRFP', INFOT, NOUT, LERR, OK ) +* * SGEQR2 * SRNAMT = 'SGEQR2' @@ -102,6 +118,19 @@ CALL SGEQR2( 2, 1, A, 1, B, W, INFO ) CALL CHKXER( 'SGEQR2', INFOT, NOUT, LERR, OK ) * +* SGEQR2P +* + SRNAMT = 'SGEQR2P' + INFOT = 1 + CALL SGEQR2P( -1, 0, A, 1, B, W, INFO ) + CALL CHKXER( 'SGEQR2P', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL SGEQR2P( 0, -1, A, 1, B, W, INFO ) + CALL CHKXER( 'SGEQR2P', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL SGEQR2P( 2, 1, A, 1, B, W, INFO ) + CALL CHKXER( 'SGEQR2P', INFOT, NOUT, LERR, OK ) +* * SGEQRS * SRNAMT = 'SGEQRS' diff --git a/TESTING/LIN/sqrt01p.f b/TESTING/LIN/sqrt01p.f new file mode 100644 index 00000000..23faa58c --- /dev/null +++ b/TESTING/LIN/sqrt01p.f @@ -0,0 +1,155 @@ + SUBROUTINE SQRT01P( M, N, A, AF, Q, R, LDA, TAU, WORK, LWORK, + $ RWORK, RESULT ) +* +* -- LAPACK test routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER LDA, LWORK, M, N +* .. +* .. Array Arguments .. + REAL A( LDA, * ), AF( LDA, * ), Q( LDA, * ), + $ R( LDA, * ), RESULT( * ), RWORK( * ), TAU( * ), + $ WORK( LWORK ) +* .. +* +* Purpose +* ======= +* +* SQRT01P tests SGEQRFP, which computes the QR factorization of an m-by-n +* matrix A, and partially tests SORGQR which forms the m-by-m +* orthogonal matrix Q. +* +* SQRT01P compares R with Q'*A, and checks that Q is orthogonal. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0. +* +* A (input) REAL array, dimension (LDA,N) +* The m-by-n matrix A. +* +* AF (output) REAL array, dimension (LDA,N) +* Details of the QR factorization of A, as returned by SGEQRFP. +* See SGEQRFP for further details. +* +* Q (output) REAL array, dimension (LDA,M) +* The m-by-m orthogonal matrix Q. +* +* R (workspace) REAL array, dimension (LDA,max(M,N)) +* +* LDA (input) INTEGER +* The leading dimension of the arrays A, AF, Q and R. +* LDA >= max(M,N). +* +* TAU (output) REAL array, dimension (min(M,N)) +* The scalar factors of the elementary reflectors, as returned +* by SGEQRFP. +* +* WORK (workspace) REAL array, dimension (LWORK) +* +* LWORK (input) INTEGER +* The dimension of the array WORK. +* +* RWORK (workspace) REAL array, dimension (M) +* +* RESULT (output) REAL array, dimension (2) +* The test ratios: +* RESULT(1) = norm( R - Q'*A ) / ( M * norm(A) * EPS ) +* RESULT(2) = norm( I - Q'*Q ) / ( M * EPS ) +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) + REAL ROGUE + PARAMETER ( ROGUE = -1.0E+10 ) +* .. +* .. Local Scalars .. + INTEGER INFO, MINMN + REAL ANORM, EPS, RESID +* .. +* .. External Functions .. + REAL SLAMCH, SLANGE, SLANSY + EXTERNAL SLAMCH, SLANGE, SLANSY +* .. +* .. External Subroutines .. + EXTERNAL SGEMM, SGEQRFP, SLACPY, SLASET, SORGQR, SSYRK +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN, REAL +* .. +* .. Scalars in Common .. + CHARACTER*32 SRNAMT +* .. +* .. Common blocks .. + COMMON / SRNAMC / SRNAMT +* .. +* .. Executable Statements .. +* + MINMN = MIN( M, N ) + EPS = SLAMCH( 'Epsilon' ) +* +* Copy the matrix A to the array AF. +* + CALL SLACPY( 'Full', M, N, A, LDA, AF, LDA ) +* +* Factorize the matrix A in the array AF. +* + SRNAMT = 'SGEQRFP' + CALL SGEQRFP( M, N, AF, LDA, TAU, WORK, LWORK, INFO ) +* +* Copy details of Q +* + CALL SLASET( 'Full', M, M, ROGUE, ROGUE, Q, LDA ) + CALL SLACPY( 'Lower', M-1, N, AF( 2, 1 ), LDA, Q( 2, 1 ), LDA ) +* +* Generate the m-by-m matrix Q +* + SRNAMT = 'SORGQR' + CALL SORGQR( M, M, MINMN, Q, LDA, TAU, WORK, LWORK, INFO ) +* +* Copy R +* + CALL SLASET( 'Full', M, N, ZERO, ZERO, R, LDA ) + CALL SLACPY( 'Upper', M, N, AF, LDA, R, LDA ) +* +* Compute R - Q'*A +* + CALL SGEMM( 'Transpose', 'No transpose', M, N, M, -ONE, Q, LDA, A, + $ LDA, ONE, R, LDA ) +* +* Compute norm( R - Q'*A ) / ( M * norm(A) * EPS ) . +* + ANORM = SLANGE( '1', M, N, A, LDA, RWORK ) + RESID = SLANGE( '1', M, N, R, LDA, RWORK ) + IF( ANORM.GT.ZERO ) THEN + RESULT( 1 ) = ( ( RESID / REAL( MAX( 1, M ) ) ) / ANORM ) / EPS + ELSE + RESULT( 1 ) = ZERO + END IF +* +* Compute I - Q'*Q +* + CALL SLASET( 'Full', M, M, ZERO, ONE, R, LDA ) + CALL SSYRK( 'Upper', 'Transpose', M, M, -ONE, Q, LDA, ONE, R, + $ LDA ) +* +* Compute norm( I - Q'*Q ) / ( M * EPS ) . +* + RESID = SLANSY( '1', 'Upper', M, R, LDA, RWORK ) +* + RESULT( 2 ) = ( RESID / REAL( MAX( 1, M ) ) ) / EPS +* + RETURN +* +* End of SQRT01P +* + END diff --git a/TESTING/LIN/zchkqr.f b/TESTING/LIN/zchkqr.f index 126b6054..0800101e 100644 --- a/TESTING/LIN/zchkqr.f +++ b/TESTING/LIN/zchkqr.f @@ -103,7 +103,7 @@ * * .. Parameters .. INTEGER NTESTS - PARAMETER ( NTESTS = 7 ) + PARAMETER ( NTESTS = 9 ) INTEGER NTYPES PARAMETER ( NTYPES = 8 ) DOUBLE PRECISION ZERO @@ -121,10 +121,14 @@ INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) DOUBLE PRECISION RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL ZGENND + EXTERNAL ZGENND +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRQR, ZGEQRS, $ ZGET02, ZLACPY, ZLARHS, ZLATB4, ZLATMS, ZQRT01, - $ ZQRT02, ZQRT03 + $ ZQRT01P, ZQRT02, ZQRT03 * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -240,6 +244,15 @@ * CALL ZQRT01( M, N, A, AF, AQ, AR, LDA, TAU, $ WORK, LWORK, RWORK, RESULT( 1 ) ) +* +* Test ZGEQRFP +* + CALL ZQRT01P( M, N, A, AF, AQ, AR, LDA, TAU, + $ WORK, LWORK, RWORK, RESULT( 8 ) ) + + IF( .NOT. ZGENND( M, N, AF, LDA ) ) + $ RESULT( 9 ) = 2*THRESH + NT = NT + 1 ELSE IF( M.GE.N ) THEN * * Test ZUNGQR, using factorization diff --git a/TESTING/LIN/zerrqr.f b/TESTING/LIN/zerrqr.f index 1e938955..44cf2ae6 100644 --- a/TESTING/LIN/zerrqr.f +++ b/TESTING/LIN/zerrqr.f @@ -38,8 +38,8 @@ $ W( NMAX ), X( NMAX ) * .. * .. External Subroutines .. - EXTERNAL ALAESM, CHKXER, ZGEQR2, ZGEQRF, ZGEQRS, ZUNG2R, - $ ZUNGQR, ZUNM2R, ZUNMQR + EXTERNAL ALAESM, CHKXER, ZGEQR2, ZGEQR2P, ZGEQRF, ZGEQRFP, + $ ZGEQRS, ZUNG2R, UNGQR, ZUNM2R, ZUNMQR * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -91,6 +91,22 @@ CALL ZGEQRF( 1, 2, A, 1, B, W, 1, INFO ) CALL CHKXER( 'ZGEQRF', INFOT, NOUT, LERR, OK ) * +* ZGEQRFP +* + SRNAMT = 'ZGEQRFP' + INFOT = 1 + CALL ZGEQRFP( -1, 0, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'ZGEQRFP', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL ZGEQRFP( 0, -1, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'ZGEQRFP', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL ZGEQRFP( 2, 1, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'ZGEQRFP', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL ZGEQRFP( 1, 2, A, 1, B, W, 1, INFO ) + CALL CHKXER( 'ZGEQRFP', INFOT, NOUT, LERR, OK ) +* * ZGEQR2 * SRNAMT = 'ZGEQR2' @@ -104,6 +120,19 @@ CALL ZGEQR2( 2, 1, A, 1, B, W, INFO ) CALL CHKXER( 'ZGEQR2', INFOT, NOUT, LERR, OK ) * +* ZGEQR2P +* + SRNAMT = 'ZGEQR2P' + INFOT = 1 + CALL ZGEQR2P( -1, 0, A, 1, B, W, INFO ) + CALL CHKXER( 'ZGEQR2P', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL ZGEQR2P( 0, -1, A, 1, B, W, INFO ) + CALL CHKXER( 'ZGEQR2P', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL ZGEQR2P( 2, 1, A, 1, B, W, INFO ) + CALL CHKXER( 'ZGEQR2P', INFOT, NOUT, LERR, OK ) +* * ZGEQRS * SRNAMT = 'ZGEQRS' diff --git a/TESTING/LIN/zqrt01p.f b/TESTING/LIN/zqrt01p.f new file mode 100644 index 00000000..e27f6c7b --- /dev/null +++ b/TESTING/LIN/zqrt01p.f @@ -0,0 +1,157 @@ + SUBROUTINE ZQRT01P( M, N, A, AF, Q, R, LDA, TAU, WORK, LWORK, + $ RWORK, RESULT ) +* +* -- LAPACK test routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER LDA, LWORK, M, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION RESULT( * ), RWORK( * ) + COMPLEX*16 A( LDA, * ), AF( LDA, * ), Q( LDA, * ), + $ R( LDA, * ), TAU( * ), WORK( LWORK ) +* .. +* +* Purpose +* ======= +* +* ZQRT01P tests ZGEQRFP, which computes the QR factorization of an m-by-n +* matrix A, and partially tests ZUNGQR which forms the m-by-m +* orthogonal matrix Q. +* +* ZQRT01P compares R with Q'*A, and checks that Q is orthogonal. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0. +* +* A (input) COMPLEX*16 array, dimension (LDA,N) +* The m-by-n matrix A. +* +* AF (output) COMPLEX*16 array, dimension (LDA,N) +* Details of the QR factorization of A, as returned by ZGEQRFP. +* See ZGEQRFP for further details. +* +* Q (output) COMPLEX*16 array, dimension (LDA,M) +* The m-by-m orthogonal matrix Q. +* +* R (workspace) COMPLEX*16 array, dimension (LDA,max(M,N)) +* +* LDA (input) INTEGER +* The leading dimension of the arrays A, AF, Q and R. +* LDA >= max(M,N). +* +* TAU (output) COMPLEX*16 array, dimension (min(M,N)) +* The scalar factors of the elementary reflectors, as returned +* by ZGEQRFP. +* +* WORK (workspace) COMPLEX*16 array, dimension (LWORK) +* +* LWORK (input) INTEGER +* The dimension of the array WORK. +* +* RWORK (workspace) DOUBLE PRECISION array, dimension (M) +* +* RESULT (output) DOUBLE PRECISION array, dimension (2) +* The test ratios: +* RESULT(1) = norm( R - Q'*A ) / ( M * norm(A) * EPS ) +* RESULT(2) = norm( I - Q'*Q ) / ( M * EPS ) +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) + COMPLEX*16 ROGUE + PARAMETER ( ROGUE = ( -1.0D+10, -1.0D+10 ) ) +* .. +* .. Local Scalars .. + INTEGER INFO, MINMN + DOUBLE PRECISION ANORM, EPS, RESID +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSY + EXTERNAL DLAMCH, ZLANGE, ZLANSY +* .. +* .. External Subroutines .. + EXTERNAL ZGEMM, ZGEQRFP, ZHERK, ZLACPY, ZLASET, ZUNGQR +* .. +* .. Intrinsic Functions .. + INTRINSIC DBLE, DCMPLX, MAX, MIN +* .. +* .. Scalars in Common .. + CHARACTER*32 SRNAMT +* .. +* .. Common blocks .. + COMMON / SRNAMC / SRNAMT +* .. +* .. Executable Statements .. +* + MINMN = MIN( M, N ) + EPS = DLAMCH( 'Epsilon' ) +* +* Copy the matrix A to the array AF. +* + CALL ZLACPY( 'Full', M, N, A, LDA, AF, LDA ) +* +* Factorize the matrix A in the array AF. +* + SRNAMT = 'ZGEQRFP' + CALL ZGEQRFP( M, N, AF, LDA, TAU, WORK, LWORK, INFO ) +* +* Copy details of Q +* + CALL ZLASET( 'Full', M, M, ROGUE, ROGUE, Q, LDA ) + CALL ZLACPY( 'Lower', M-1, N, AF( 2, 1 ), LDA, Q( 2, 1 ), LDA ) +* +* Generate the m-by-m matrix Q +* + SRNAMT = 'ZUNGQR' + CALL ZUNGQR( M, M, MINMN, Q, LDA, TAU, WORK, LWORK, INFO ) +* +* Copy R +* + CALL ZLASET( 'Full', M, N, DCMPLX( ZERO ), DCMPLX( ZERO ), R, + $ LDA ) + CALL ZLACPY( 'Upper', M, N, AF, LDA, R, LDA ) +* +* Compute R - Q'*A +* + CALL ZGEMM( 'Conjugate transpose', 'No transpose', M, N, M, + $ DCMPLX( -ONE ), Q, LDA, A, LDA, DCMPLX( ONE ), R, + $ LDA ) +* +* Compute norm( R - Q'*A ) / ( M * norm(A) * EPS ) . +* + ANORM = ZLANGE( '1', M, N, A, LDA, RWORK ) + RESID = ZLANGE( '1', M, N, R, LDA, RWORK ) + IF( ANORM.GT.ZERO ) THEN + RESULT( 1 ) = ( ( RESID / DBLE( MAX( 1, M ) ) ) / ANORM ) / EPS + ELSE + RESULT( 1 ) = ZERO + END IF +* +* Compute I - Q'*Q +* + CALL ZLASET( 'Full', M, M, DCMPLX( ZERO ), DCMPLX( ONE ), R, LDA ) + CALL ZHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, Q, LDA, + $ ONE, R, LDA ) +* +* Compute norm( I - Q'*Q ) / ( M * EPS ) . +* + RESID = ZLANSY( '1', 'Upper', M, R, LDA, RWORK ) +* + RESULT( 2 ) = ( RESID / DBLE( MAX( 1, M ) ) ) / EPS +* + RETURN +* +* End of ZQRT01P +* + END -- 2.34.1