phynexis-deps/lapack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fix ?GETSLS

Browse Source
This commit is contained in:
eugene.chereshnev
2016-12-14 11:30:51 -08:00
committed by eugene.chereshnev
parent 3136b938aa
commit ca017fe2b6
4 changed files with 285 additions and 253 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SRC/cgetsls.f
+72 -65
View File
@@ -81,7 +81,7 @@
*> On entry, the M-by-N matrix A.
*> On exit,
*> A is overwritten by details of its QR or LQ
*> factorization as returned by CGEQR or CGEQL.
*> factorization as returned by CGEQR or CGELQ.
*> \endverbatim
*>
*> \param[in] LDA
@@ -152,18 +152,18 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*> \date November 2016
*
*> \ingroup complexGEsolve
*
* =====================================================================
SUBROUTINE CGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
$ WORK, LWORK, INFO )
$ WORK, LWORK, INFO )
*
* -- LAPACK driver routine (version 3.4.0) --
* -- LAPACK driver routine (version 3.7.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2011
* November 2016
*
* .. Scalar Arguments ..
CHARACTER TRANS
@@ -177,18 +177,18 @@
* =====================================================================
*
* .. Parameters ..
REAL ZERO, ONE
REAL ZERO, ONE
PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
COMPLEX CZERO
COMPLEX CZERO
PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) )
* ..
* .. Local Scalars ..
LOGICAL LQUERY, TRAN
INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
$ SCLLEN, MNK, TSZ, TSZM, LW, LWM, LW1, LW2,
$ INFO2
$ SCLLEN, MNK, TSZO, TSZM, LWO, LWM, LW1, LW2,
$ WSIZEO, WSIZEM, INFO2
REAL ANRM, BIGNUM, BNRM, SMLNUM
COMPLEX WRK1(5), WRK2
COMPLEX TQ( 5 ), WORKQ
* ..
* .. External Functions ..
LOGICAL LSAME
@@ -207,10 +207,10 @@
*
* Test the input arguments.
*
INFO=0
INFO = 0
MINMN = MIN( M, N )
MAXMN = MAX( M, N )
MNK = MAX(MINMN,NRHS)
MNK = MAX( MINMN, NRHS )
TRAN = LSAME( TRANS, 'C' )
*
LQUERY = ( LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
@@ -229,65 +229,71 @@
INFO = -8
END IF
*
IF( INFO.EQ.0 ) THEN
IF( INFO.EQ.0 ) THEN
*
* Determine the block size and minimum LWORK
*
IF ( M.GE.N ) THEN
CALL CGEQR( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
TSZ = INT( WRK1(1) )
LW = INT( WRK2 )
CALL CGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
$ TSZ, B, LDB, WRK2, -1 , INFO2 )
LW = MAX( LW, INT(WRK2) )
CALL CGEQR( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
TSZM = INT( WRK1(1) )
LWM = INT( WRK2 )
CALL CGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
$ TSZM, B, LDB, WRK2, -2, INFO2 )
LWM = MAX( LWM, INT(WRK2) )
IF( M.GE.N ) THEN
CALL CGEQR( M, N, A, LDA, TQ, -1, WORKQ, -1, INFO2 )
TSZO = INT( TQ( 1 ) )
LWO = INT( WORKQ )
CALL CGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, TQ,
$ TSZO, B, LDB, WORKQ, -1, INFO2 )
LWO = MAX( LWO, INT( WORKQ ) )
CALL CGEQR( M, N, A, LDA, TQ, -2, WORKQ, -2, INFO2 )
TSZM = INT( TQ( 1 ) )
LWM = INT( WORKQ )
CALL CGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, TQ,
$ TSZM, B, LDB, WORKQ, -1, INFO2 )
LWM = MAX( LWM, INT( WORKQ ) )
WSIZEO = TSZO + LWO
WSIZEM = TSZM + LWM
ELSE
CALL CGELQ( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
TSZ = INT( WRK1(1) )
LW = INT( WRK2 )
CALL CGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
$ TSZ, B, LDB, WRK2, -1, INFO2 )
LW = MAX( LW, INT(WRK2) )
CALL CGELQ( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
TSZM = INT( WRK1(1) )
LWM = INT( WRK2 )
CALL CGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
$ TSZ, B, LDB, WRK2, -2, INFO2 )
LWM = MAX( LWM, INT(WRK2) )
CALL CGELQ( M, N, A, LDA, TQ, -1, WORKQ, -1, INFO2 )
TSZO = INT( TQ( 1 ) )
LWO = INT( WORKQ )
CALL CGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, TQ,
$ TSZO, B, LDB, WORKQ, -1, INFO2 )
LWO = MAX( LWO, INT( WORKQ ) )
CALL CGELQ( M, N, A, LDA, TQ, -2, WORKQ, -2, INFO2 )
TSZM = INT( TQ( 1 ) )
LWM = INT( WORKQ )
CALL CGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, TQ,
$ TSZO, B, LDB, WORKQ, -1, INFO2 )
LWM = MAX( LWM, INT( WORKQ ) )
WSIZEO = TSZO + LWO
WSIZEM = TSZM + LWM
END IF
*
IF( ( LWORK.LT.( TSZM + LWM ) ) .AND. ( .NOT.LQUERY ) ) THEN
INFO=-10
IF( ( LWORK.LT.WSIZEM ).AND.( .NOT.LQUERY ) ) THEN
INFO = -10
END IF
*
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'CGETSLS', -INFO )
WORK( 1 ) = REAL( TSZ + LW )
WORK( 1 ) = REAL( WSIZEO )
RETURN
END IF
IF ( LQUERY ) THEN
WORK( 1 ) = REAL( TSZ + LW )
IF( LQUERY ) THEN
IF( LWORK.EQ.-1 ) WORK( 1 ) = REAL( WSIZEO )
IF( LWORK.EQ.-2 ) WORK( 1 ) = REAL( WSIZEM )
RETURN
END IF
IF( LWORK.LT.( TSZ + LW ) ) THEN
IF( LWORK.LT.WSIZEO ) THEN
LW1 = TSZM
LW2 = LWM
ELSE
LW1 = TSZ
LW2 = LW
LW1 = TSZO
LW2 = LWO
END IF
*
* Quick return if possible
*
IF( MIN( M, N, NRHS ).EQ.0 ) THEN
CALL CLASET( 'FULL', MAX( M, N ), NRHS, CZERO, CZERO,
$ B, LDB )
$ B, LDB )
RETURN
END IF
*
@@ -343,12 +349,12 @@
IBSCL = 2
END IF
*
IF ( M.GE.N) THEN
IF ( M.GE.N ) THEN
*
* compute QR factorization of A
*
CALL CGEQR( M, N, A, LDA, WORK(LW2+1), LW1,
$ WORK(1), LW2, INFO )
CALL CGEQR( M, N, A, LDA, WORK( LW2+1 ), LW1,
$ WORK( 1 ), LW2, INFO )
IF ( .NOT.TRAN ) THEN
*
* Least-Squares Problem min || A * X - B ||
@@ -356,13 +362,14 @@
* B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
*
CALL CGEMQR( 'L' , 'C', M, NRHS, N, A, LDA,
$ WORK(LW2+1), LW1, B, LDB, WORK(1), LW2, INFO )
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
* B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
*
CALL CTRTRS( 'U', 'N', 'N', N, NRHS,
$ A, LDA, B, LDB, INFO )
IF(INFO.GT.0) THEN
$ A, LDA, B, LDB, INFO )
IF( INFO.GT.0 ) THEN
RETURN
END IF
SCLLEN = N
@@ -390,7 +397,7 @@
* B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
*
CALL CGEMQR( 'L', 'N', M, NRHS, N, A, LDA,
$ WORK( LW2+1), LW1, B, LDB, WORK( 1 ), LW2,
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
SCLLEN = M
@@ -401,8 +408,8 @@
*
* Compute LQ factorization of A
*
CALL CGELQ( M, N, A, LDA, WORK(LW2+1), LW1,
$ WORK(1), LW2, INFO )
CALL CGELQ( M, N, A, LDA, WORK( LW2+1 ), LW1,
$ WORK( 1 ), LW2, INFO )
*
* workspace at least M, optimally M*NB.
*
@@ -430,7 +437,7 @@
* B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
*
CALL CGEMLQ( 'L', 'C', N, NRHS, M, A, LDA,
$ WORK( LW2+1), LW1, B, LDB, WORK( 1 ), LW2,
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
* workspace at least NRHS, optimally NRHS*NB
@@ -444,7 +451,7 @@
* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
*
CALL CGEMLQ( 'L', 'N', N, NRHS, M, A, LDA,
$ WORK( LW2+1), LW1, B, LDB, WORK( 1 ), LW2,
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
* workspace at least NRHS, optimally NRHS*NB
@@ -468,21 +475,21 @@
*
IF( IASCL.EQ.1 ) THEN
CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
$ INFO )
$ INFO )
ELSE IF( IASCL.EQ.2 ) THEN
CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
$ INFO )
$ INFO )
END IF
IF( IBSCL.EQ.1 ) THEN
CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
ELSE IF( IBSCL.EQ.2 ) THEN
CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
END IF
*
50 CONTINUE
WORK( 1 ) = REAL( TSZ + LW )
WORK( 1 ) = REAL( TSZO + LWO )
RETURN
*
* End of ZGETSLS
SRC/dgetsls.f
+66 -59
View File
@@ -81,7 +81,7 @@
*> On entry, the M-by-N matrix A.
*> On exit,
*> A is overwritten by details of its QR or LQ
*> factorization as returned by DGEQR or DGEQL.
*> factorization as returned by DGEQR or DGELQ.
*> \endverbatim
*>
*> \param[in] LDA
@@ -152,18 +152,18 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*> \date November 2016
*
*> \ingroup doubleGEsolve
*
* =====================================================================
SUBROUTINE DGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
$ WORK, LWORK, INFO )
$ WORK, LWORK, INFO )
*
* -- LAPACK driver routine (version 3.4.0) --
* -- LAPACK driver routine (version 3.7.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2011
* November 2016
*
* .. Scalar Arguments ..
CHARACTER TRANS
@@ -183,9 +183,9 @@
* .. Local Scalars ..
LOGICAL LQUERY, TRAN
INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
$ SCLLEN, MNK, TSZ, TSZM, LW, LWM, LW1, LW2
$ INFO2
DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM, WRK1(5), WRK2
$ SCLLEN, MNK, TSZO, TSZM, LWO, LWM, LW1, LW2,
$ WSIZEO, WSIZEM, INFO2
DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM, TQ( 5 ), WORKQ
* ..
* .. External Functions ..
LOGICAL LSAME
@@ -204,7 +204,7 @@
*
* Test the input arguments.
*
INFO=0
INFO = 0
MINMN = MIN( M, N )
MAXMN = MAX( M, N )
MNK = MAX( MINMN, NRHS )
@@ -226,65 +226,71 @@
INFO = -8
END IF
*
IF( INFO.EQ.0 ) THEN
IF( INFO.EQ.0 ) THEN
*
* Determine the block size and minimum LWORK
*
IF ( M.GE.N ) THEN
CALL DGEQR( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
TSZ = INT( WRK1(1) )
LW = INT( WRK2 )
CALL DGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
$ TSZ, B, LDB, WRK2, -1 , INFO2 )
LW = MAX( LW, INT(WRK2) )
CALL DGEQR( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
TSZM = INT( WRK1(1) )
LWM = INT( WRK2 )
CALL DGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
$ TSZM, B, LDB, WRK2, -2, INFO2 )
LWM = MAX( LWM, INT(WRK2) )
IF( M.GE.N ) THEN
CALL DGEQR( M, N, A, LDA, TQ, -1, WORKQ, -1, INFO2 )
TSZO = INT( TQ( 1 ) )
LWO = INT( WORKQ )
CALL DGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, TQ,
$ TSZO, B, LDB, WORKQ, -1, INFO2 )
LWO = MAX( LWO, INT( WORKQ ) )
CALL DGEQR( M, N, A, LDA, TQ, -2, WORKQ, -2, INFO2 )
TSZM = INT( TQ( 1 ) )
LWM = INT( WORKQ )
CALL DGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, TQ,
$ TSZM, B, LDB, WORKQ, -1, INFO2 )
LWM = MAX( LWM, INT( WORKQ ) )
WSIZEO = TSZO + LWO
WSIZEM = TSZM + LWM
ELSE
CALL DGELQ( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
TSZ = INT( WRK1(1) )
LW = INT( WRK2 )
CALL DGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
$ TSZ, B, LDB, WRK2, -1, INFO2 )
LW = MAX( LW, INT(WRK2) )
CALL DGELQ( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
TSZM = INT( WRK1(1) )
LWM = INT( WRK2 )
CALL DGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
$ TSZ, B, LDB, WRK2, -2, INFO2 )
LWM = MAX( LWM, INT(WRK2) )
CALL DGELQ( M, N, A, LDA, TQ, -1, WORKQ, -1, INFO2 )
TSZO = INT( TQ( 1 ) )
LWO = INT( WORKQ )
CALL DGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, TQ,
$ TSZO, B, LDB, WORKQ, -1, INFO2 )
LWO = MAX( LWO, INT( WORKQ ) )
CALL DGELQ( M, N, A, LDA, TQ, -2, WORKQ, -2, INFO2 )
TSZM = INT( TQ( 1 ) )
LWM = INT( WORKQ )
CALL DGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, TQ,
$ TSZO, B, LDB, WORKQ, -1, INFO2 )
LWM = MAX( LWM, INT( WORKQ ) )
WSIZEO = TSZO + LWO
WSIZEM = TSZM + LWM
END IF
*
IF( ( LWORK.LT.( TSZM + LWM ) ) .AND. (.NOT.LQUERY)) THEN
INFO=-10
IF( ( LWORK.LT.WSIZEM ).AND.( .NOT.LQUERY ) ) THEN
INFO = -10
END IF
*
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGETSLS', -INFO )
WORK( 1 ) = DBLE( TSZ + LW )
WORK( 1 ) = DBLE( WSIZEO )
RETURN
END IF
IF ( LQUERY ) THEN
WORK( 1 ) = DBLE( TSZ + LW )
IF( LQUERY ) THEN
IF( LWORK.EQ.-1 ) WORK( 1 ) = REAL( WSIZEO )
IF( LWORK.EQ.-2 ) WORK( 1 ) = REAL( WSIZEM )
RETURN
END IF
IF( LWORK.LT.( TSZ + LW ) ) THEN
IF( LWORK.LT.WSIZEO ) THEN
LW1 = TSZM
LW2 = LWM
ELSE
LW1 = TSZ
LW2 = LW
LW1 = TSZO
LW2 = LWO
END IF
*
* Quick return if possible
*
IF( MIN( M, N, NRHS ).EQ.0 ) THEN
CALL DLASET( 'FULL', MAX( M, N ), NRHS, ZERO, ZERO,
$ B, LDB )
$ B, LDB )
RETURN
END IF
*
@@ -344,8 +350,8 @@
*
* compute QR factorization of A
*
CALL DGEQR( M, N, A, LDA, WORK(LW2+1), LW1,
$ WORK(1), LW2, INFO )
CALL DGEQR( M, N, A, LDA, WORK( LW2+1 ), LW1,
$ WORK( 1 ), LW2, INFO )
IF ( .NOT.TRAN ) THEN
*
* Least-Squares Problem min || A * X - B ||
@@ -353,13 +359,14 @@
* B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
*
CALL DGEMQR( 'L' , 'T', M, NRHS, N, A, LDA,
$ WORK(LW2+1), LW1, B, LDB, WORK(1), LW2, INFO )
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
* B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
*
CALL DTRTRS( 'U', 'N', 'N', N, NRHS,
$ A, LDA, B, LDB, INFO )
IF(INFO.GT.0) THEN
$ A, LDA, B, LDB, INFO )
IF( INFO.GT.0 ) THEN
RETURN
END IF
SCLLEN = N
@@ -387,7 +394,7 @@
* B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
*
CALL DGEMQR( 'L', 'N', M, NRHS, N, A, LDA,
$ WORK( LW2+1), LW1, B, LDB, WORK( 1 ), LW2,
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
SCLLEN = M
@@ -398,8 +405,8 @@
*
* Compute LQ factorization of A
*
CALL DGELQ( M, N, A, LDA, WORK(LW2+1), LW1,
$ WORK(1), LW2, INFO )
CALL DGELQ( M, N, A, LDA, WORK( LW2+1 ), LW1,
$ WORK( 1 ), LW2, INFO )
*
* workspace at least M, optimally M*NB.
*
@@ -427,7 +434,7 @@
* B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
*
CALL DGEMLQ( 'L', 'T', N, NRHS, M, A, LDA,
$ WORK( LW2+1), LW1, B, LDB, WORK( 1 ), LW2,
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
* workspace at least NRHS, optimally NRHS*NB
@@ -441,7 +448,7 @@
* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
*
CALL DGEMLQ( 'L', 'N', N, NRHS, M, A, LDA,
$ WORK( LW2+1), LW1, B, LDB, WORK( 1 ), LW2,
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
* workspace at least NRHS, optimally NRHS*NB
@@ -465,21 +472,21 @@
*
IF( IASCL.EQ.1 ) THEN
CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
$ INFO )
$ INFO )
ELSE IF( IASCL.EQ.2 ) THEN
CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
$ INFO )
$ INFO )
END IF
IF( IBSCL.EQ.1 ) THEN
CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
$ INFO )
ELSE IF( IBSCL.EQ.2 ) THEN
CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
$ INFO )
END IF
*
50 CONTINUE
WORK( 1 ) = DBLE( TSZ + LW )
WORK( 1 ) = DBLE( TSZO + LWO )
RETURN
*
* End of DGETSLS
SRC/sgetsls.f
+77 -66
View File
@@ -29,16 +29,21 @@
*> 1. If TRANS = 'N' and m >= n: find the least squares solution of
*> an overdetermined system, i.e., solve the least squares problem
*> minimize || B - A*X ||.
*>
*> 2. If TRANS = 'N' and m < n: find the minimum norm solution of
*> an underdetermined system A * X = B.
*>
*> 3. If TRANS = 'T' and m >= n: find the minimum norm solution of
*> an undetermined system A**T * X = B.
*>
*> 4. If TRANS = 'T' and m < n: find the least squares solution of
*> an overdetermined system, i.e., solve the least squares problem
*> minimize || B - A**T * X ||.
*>
*> Several right hand side vectors b and solution vectors x can be
*> handled in a single call; they are stored as the columns of the
*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
*> matrix X.
*> \endverbatim
*
* Arguments:
@@ -76,7 +81,7 @@
*> On entry, the M-by-N matrix A.
*> On exit,
*> A is overwritten by details of its QR or LQ
*> factorization as returned by SGEQR or SGEQL.
*> factorization as returned by SGEQR or SGELQ.
*> \endverbatim
*>
*> \param[in] LDA
@@ -147,18 +152,18 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*> \date November 2016
*
*> \ingroup doubleGEsolve
*
* =====================================================================
SUBROUTINE SGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
$ WORK, LWORK, INFO )
$ WORK, LWORK, INFO )
*
* -- LAPACK driver routine (version 3.4.0) --
* -- LAPACK driver routine (version 3.7.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2011
* November 2016
*
* .. Scalar Arguments ..
CHARACTER TRANS
@@ -178,10 +183,9 @@
* .. Local Scalars ..
LOGICAL LQUERY, TRAN
INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
$ SCLLEN, MNK, TSZ, TSZM, LW, LWM, LW1, LW2
$ INFO2
REAL ANRM, BIGNUM, BNRM, SMLNUM
COMPLEX WRK1(5), WRK2
$ SCLLEN, MNK, TSZO, TSZM, LWO, LWM, LW1, LW2,
$ WSIZEO, WSIZEM, INFO2
REAL ANRM, BIGNUM, BNRM, SMLNUM, TQ( 5 ), WORKQ
* ..
* .. External Functions ..
LOGICAL LSAME
@@ -200,7 +204,7 @@
*
* Test the input arguments.
*
INFO=0
INFO = 0
MINMN = MIN( M, N )
MAXMN = MAX( M, N )
MNK = MAX( MINMN, NRHS )
@@ -222,65 +226,71 @@
INFO = -8
END IF
*
IF( INFO.EQ.0 ) THEN
IF( INFO.EQ.0 ) THEN
*
* Determine the block size and minimum LWORK
*
IF ( M.GE.N ) THEN
CALL SGEQR( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
TSZ = INT( WRK1(1) )
LW = INT( WRK2 )
CALL SGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
$ TSZ, B, LDB, WRK2, -1 , INFO2 )
LW = MAX( LW, INT(WRK2) )
CALL SGEQR( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
TSZM = INT( WRK1(1) )
LWM = INT( WRK2 )
CALL SGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
$ TSZM, B, LDB, WRK2, -2, INFO2 )
LWM = MAX( LWM, INT(WRK2) )
IF( M.GE.N ) THEN
CALL SGEQR( M, N, A, LDA, TQ, -1, WORKQ, -1, INFO2 )
TSZO = INT( TQ( 1 ) )
LWO = INT( WORKQ )
CALL SGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, TQ,
$ TSZO, B, LDB, WORKQ, -1, INFO2 )
LWO = MAX( LWO, INT( WORKQ ) )
CALL SGEQR( M, N, A, LDA, TQ, -2, WORKQ, -2, INFO2 )
TSZM = INT( TQ( 1 ) )
LWM = INT( WORKQ )
CALL SGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, TQ,
$ TSZM, B, LDB, WORKQ, -1, INFO2 )
LWM = MAX( LWM, INT( WORKQ ) )
WSIZEO = TSZO + LWO
WSIZEM = TSZM + LWM
ELSE
CALL SGELQ( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
TSZ = INT( WRK1(1) )
LW = INT( WRK2 )
CALL SGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
$ TSZ, B, LDB, WRK2, -1, INFO2 )
LW = MAX( LW, INT(WRK2) )
CALL SGELQ( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
TSZM = INT( WRK1(1) )
LWM = INT( WRK2 )
CALL SGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
$ TSZ, B, LDB, WRK2, -2, INFO2 )
LWM = MAX( LWM, INT(WRK2) )
CALL SGELQ( M, N, A, LDA, TQ, -1, WORKQ, -1, INFO2 )
TSZO = INT( TQ( 1 ) )
LWO = INT( WORKQ )
CALL SGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, TQ,
$ TSZO, B, LDB, WORKQ, -1, INFO2 )
LWO = MAX( LWO, INT( WORKQ ) )
CALL SGELQ( M, N, A, LDA, TQ, -2, WORKQ, -2, INFO2 )
TSZM = INT( TQ( 1 ) )
LWM = INT( WORKQ )
CALL SGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, TQ,
$ TSZO, B, LDB, WORKQ, -1, INFO2 )
LWM = MAX( LWM, INT( WORKQ ) )
WSIZEO = TSZO + LWO
WSIZEM = TSZM + LWM
END IF
*
IF( ( LWORK.LT.( TSZM + LWM ) ) .AND. (.NOT.LQUERY)) THEN
INFO=-10
IF( ( LWORK.LT.WSIZEM ).AND.( .NOT.LQUERY ) ) THEN
INFO = -10
END IF
*
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'SGETSLS', -INFO )
WORK( 1 ) = REAL( TSZ + LW )
WORK( 1 ) = REAL( WSIZEO )
RETURN
END IF
IF ( LQUERY ) THEN
WORK( 1 ) = REAL( TSZ + LW )
IF( LQUERY ) THEN
IF( LWORK.EQ.-1 ) WORK( 1 ) = REAL( WSIZEO )
IF( LWORK.EQ.-2 ) WORK( 1 ) = REAL( WSIZEM )
RETURN
END IF
IF( LWORK.LT.( TSZ + LW ) ) THEN
IF( LWORK.LT.WSIZEO ) THEN
LW1 = TSZM
LW2 = LWM
ELSE
LW1 = TSZ
LW2 = LW
LW1 = TSZO
LW2 = LWO
END IF
*
* Quick return if possible
*
IF( MIN( M, N, NRHS ).EQ.0 ) THEN
CALL SLASET( 'FULL', MAX( M, N ), NRHS, ZERO, ZERO,
$ B, LDB )
$ B, LDB )
RETURN
END IF
*
@@ -336,12 +346,12 @@
IBSCL = 2
END IF
*
IF ( M.GE.N) THEN
IF ( M.GE.N ) THEN
*
* compute QR factorization of A
*
CALL SGEQR( M, N, A, LDA, WORK(LW2+1), LW1,
$ WORK(1), LW2, INFO )
CALL SGEQR( M, N, A, LDA, WORK( LW2+1 ), LW1,
$ WORK( 1 ), LW2, INFO )
IF ( .NOT.TRAN ) THEN
*
* Least-Squares Problem min || A * X - B ||
@@ -349,13 +359,14 @@
* B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
*
CALL SGEMQR( 'L' , 'T', M, NRHS, N, A, LDA,
$ WORK(LW2+1), LW1, B, LDB, WORK(1), LW2, INFO )
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
* B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
*
CALL STRTRS( 'U', 'N', 'N', N, NRHS,
$ A, LDA, B, LDB, INFO )
IF(INFO.GT.0) THEN
$ A, LDA, B, LDB, INFO )
IF( INFO.GT.0 ) THEN
RETURN
END IF
SCLLEN = N
@@ -383,7 +394,7 @@
* B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
*
CALL SGEMQR( 'L', 'N', M, NRHS, N, A, LDA,
$ WORK( LW2+1), LW1, B, LDB, WORK( 1 ), LW2,
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
SCLLEN = M
@@ -394,8 +405,8 @@
*
* Compute LQ factorization of A
*
CALL SGELQ( M, N, A, LDA, WORK(LW2+1), LW1,
$ WORK(1), LW2, INFO )
CALL SGELQ( M, N, A, LDA, WORK( LW2+1 ), LW1,
$ WORK( 1 ), LW2, INFO )
*
* workspace at least M, optimally M*NB.
*
@@ -423,7 +434,7 @@
* B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
*
CALL SGEMLQ( 'L', 'T', N, NRHS, M, A, LDA,
$ WORK( LW2+1), LW1, B, LDB, WORK( 1 ), LW2,
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
* workspace at least NRHS, optimally NRHS*NB
@@ -437,7 +448,7 @@
* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
*
CALL SGEMLQ( 'L', 'N', N, NRHS, M, A, LDA,
$ WORK( LW2+1), LW1, B, LDB, WORK( 1 ), LW2,
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
* workspace at least NRHS, optimally NRHS*NB
@@ -461,21 +472,21 @@
*
IF( IASCL.EQ.1 ) THEN
CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
$ INFO )
$ INFO )
ELSE IF( IASCL.EQ.2 ) THEN
CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
$ INFO )
$ INFO )
END IF
IF( IBSCL.EQ.1 ) THEN
CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
ELSE IF( IBSCL.EQ.2 ) THEN
CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
END IF
*
50 CONTINUE
WORK( 1 ) = REAL( TSZ + LW )
WORK( 1 ) = REAL( TSZO + LWO )
RETURN
*
* End of SGETSLS
SRC/zgetsls.f
+70 -63
View File
@@ -81,7 +81,7 @@
*> On entry, the M-by-N matrix A.
*> On exit,
*> A is overwritten by details of its QR or LQ
*> factorization as returned by ZGEQR or ZGEQL.
*> factorization as returned by ZGEQR or ZGELQ.
*> \endverbatim
*>
*> \param[in] LDA
@@ -152,18 +152,18 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*> \date November 2016
*
*> \ingroup complex16GEsolve
*
* =====================================================================
SUBROUTINE ZGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
$ WORK, LWORK, INFO )
$ WORK, LWORK, INFO )
*
* -- LAPACK driver routine (version 3.4.0) --
* -- LAPACK driver routine (version 3.7.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2011
* November 2016
*
* .. Scalar Arguments ..
CHARACTER TRANS
@@ -185,10 +185,10 @@
* .. Local Scalars ..
LOGICAL LQUERY, TRAN
INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
$ SCLLEN, MNK, TSZ, TSZM, LW, LWM, LW1, LW2
$ INFO2
$ SCLLEN, MNK, TSZO, TSZM, LWO, LWM, LW1, LW2,
$ WSIZEO, WSIZEM, INFO2
DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
COMPLEX*16 WRK1(5), WRK2
COMPLEX*16 TQ( 5 ), WORKQ
* ..
* .. External Functions ..
LOGICAL LSAME
@@ -207,10 +207,10 @@
*
* Test the input arguments.
*
INFO=0
INFO = 0
MINMN = MIN( M, N )
MAXMN = MAX( M, N )
MNK = MAX(MINMN,NRHS)
MNK = MAX( MINMN, NRHS )
TRAN = LSAME( TRANS, 'C' )
*
LQUERY = ( LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
@@ -229,65 +229,71 @@
INFO = -8
END IF
*
IF( INFO.EQ.0 ) THEN
IF( INFO.EQ.0 ) THEN
*
* Determine the block size and minimum LWORK
*
IF ( M.GE.N ) THEN
CALL ZGEQR( M, N, A, LDA, WRK1, -1, WRK2, -1, INFO2 )
TSZ = INT( WRK1(1) )
LW = INT( WRK2 )
CALL ZGEQR( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
TSZM = INT( WRK1(1) )
LWM = INT( WRK2 )
CALL ZGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
$ TSZ, B, LDB, WRK2, -1 , INFO2 )
LW = MAX( LW, INT(WRK2) )
CALL ZGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
$ TSZM, B, LDB, WRK2, -2, INFO2 )
LWM = MAX( LWM, INT(WRK2) )
IF( M.GE.N ) THEN
CALL ZGEQR( M, N, A, LDA, TQ, -1, WORKQ, -1, INFO2 )
TSZO = INT( TQ( 1 ) )
LWO = INT( WORKQ )
CALL ZGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, TQ,
$ TSZO, B, LDB, WORKQ, -1, INFO2 )
LWO = MAX( LWO, INT( WORKQ ) )
CALL ZGEQR( M, N, A, LDA, TQ, -2, WORKQ, -2, INFO2 )
TSZM = INT( TQ( 1 ) )
LWM = INT( WORKQ )
CALL ZGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, TQ,
$ TSZM, B, LDB, WORKQ, -1, INFO2 )
LWM = MAX( LWM, INT( WORKQ ) )
WSIZEO = TSZO + LWO
WSIZEM = TSZM + LWM
ELSE
CALL ZGELQ( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
TSZ = INT( WRK1(1) )
LW = INT( WRK2 )
CALL ZGELQ( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
TSZM = INT( WRK1(1) )
LWM = INT( WRK2 )
CALL ZGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
$ TSZ, B, LDB, WRK2, -1, INFO2 )
LW = MAX( LW, INT(WRK2) )
CALL ZGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
$ TSZ, B, LDB, WRK2, -2, INFO2 )
LWM = MAX( LWM, INT(WRK2) )
CALL ZGELQ( M, N, A, LDA, TQ, -1, WORKQ, -1, INFO2 )
TSZO = INT( TQ( 1 ) )
LWO = INT( WORKQ )
CALL ZGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, TQ,
$ TSZO, B, LDB, WORKQ, -1, INFO2 )
LWO = MAX( LWO, INT( WORKQ ) )
CALL ZGELQ( M, N, A, LDA, TQ, -2, WORKQ, -2, INFO2 )
TSZM = INT( TQ( 1 ) )
LWM = INT( WORKQ )
CALL ZGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, TQ,
$ TSZO, B, LDB, WORKQ, -1, INFO2 )
LWM = MAX( LWM, INT( WORKQ ) )
WSIZEO = TSZO + LWO
WSIZEM = TSZM + LWM
END IF
*
IF( ( LWORK.LT.( TSZM + LWM ) ) .AND. (.NOT.LQUERY)) THEN
INFO=-10
IF( ( LWORK.LT.WSIZEM ).AND.( .NOT.LQUERY ) ) THEN
INFO = -10
END IF
*
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'ZGETSLS', -INFO )
WORK( 1 ) = DBLE( TSZ + LW )
WORK( 1 ) = DBLE( WSIZEO )
RETURN
END IF
IF ( LQUERY ) THEN
WORK( 1 ) = DBLE( TSZ + LW )
IF( LQUERY ) THEN
IF( LWORK.EQ.-1 ) WORK( 1 ) = REAL( WSIZEO )
IF( LWORK.EQ.-2 ) WORK( 1 ) = REAL( WSIZEM )
RETURN
END IF
IF( LWORK.LT.( TSZ + LW ) ) THEN
IF( LWORK.LT.WSIZEO ) THEN
LW1 = TSZM
LW2 = LWM
ELSE
LW1 = TSZ
LW2 = LW
LW1 = TSZO
LW2 = LWO
END IF
*
* Quick return if possible
*
IF( MIN( M, N, NRHS ).EQ.0 ) THEN
CALL ZLASET( 'FULL', MAX( M, N ), NRHS, CZERO, CZERO,
$ B, LDB )
$ B, LDB )
RETURN
END IF
*
@@ -343,12 +349,12 @@
IBSCL = 2
END IF
*
IF ( M.GE.N) THEN
IF ( M.GE.N ) THEN
*
* compute QR factorization of A
*
CALL ZGEQR( M, N, A, LDA, WORK(LW2+1), LW1,
$ WORK(1), LW2, INFO )
CALL ZGEQR( M, N, A, LDA, WORK( LW2+1 ), LW1,
$ WORK( 1 ), LW2, INFO )
IF ( .NOT.TRAN ) THEN
*
* Least-Squares Problem min || A * X - B ||
@@ -356,13 +362,14 @@
* B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
*
CALL ZGEMQR( 'L' , 'C', M, NRHS, N, A, LDA,
$ WORK(LW2+1), LW1, B, LDB, WORK(1), LW2, INFO )
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
* B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
*
CALL ZTRTRS( 'U', 'N', 'N', N, NRHS,
$ A, LDA, B, LDB, INFO )
IF(INFO.GT.0) THEN
$ A, LDA, B, LDB, INFO )
IF( INFO.GT.0 ) THEN
RETURN
END IF
SCLLEN = N
@@ -390,7 +397,7 @@
* B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
*
CALL ZGEMQR( 'L', 'N', M, NRHS, N, A, LDA,
$ WORK( LW2+1), LW1, B, LDB, WORK( 1 ), LW2,
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
SCLLEN = M
@@ -401,8 +408,8 @@
*
* Compute LQ factorization of A
*
CALL ZGELQ( M, N, A, LDA, WORK(LW2+1), LW1,
$ WORK(1), LW2, INFO )
CALL ZGELQ( M, N, A, LDA, WORK( LW2+1 ), LW1,
$ WORK( 1 ), LW2, INFO )
*
* workspace at least M, optimally M*NB.
*
@@ -430,7 +437,7 @@
* B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
*
CALL ZGEMLQ( 'L', 'C', N, NRHS, M, A, LDA,
$ WORK( LW2+1), LW1, B, LDB, WORK( 1 ), LW2,
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
* workspace at least NRHS, optimally NRHS*NB
@@ -444,7 +451,7 @@
* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
*
CALL ZGEMLQ( 'L', 'N', N, NRHS, M, A, LDA,
$ WORK( LW2+1), LW1, B, LDB, WORK( 1 ), LW2,
$ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
$ INFO )
*
* workspace at least NRHS, optimally NRHS*NB
@@ -468,21 +475,21 @@
*
IF( IASCL.EQ.1 ) THEN
CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
$ INFO )
$ INFO )
ELSE IF( IASCL.EQ.2 ) THEN
CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
$ INFO )
$ NFO )
END IF
IF( IBSCL.EQ.1 ) THEN
CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
ELSE IF( IBSCL.EQ.2 ) THEN
CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
END IF
*
50 CONTINUE
WORK( 1 ) = DBLE( TSZ + LW )
WORK( 1 ) = DBLE( TSZO + LWO )
RETURN
*
* End of ZGETSLS