Update LAPACK to 3.9.0

This commit is contained in:
Martin Kroeker 2019-12-29 23:43:51 +01:00 committed by GitHub
parent ef93c62eb7
commit ff45990663
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
51 changed files with 1115 additions and 306 deletions

View File

@ -127,7 +127,7 @@
*> \param[in,out] AUXV *> \param[in,out] AUXV
*> \verbatim *> \verbatim
*> AUXV is REAL array, dimension (NB) *> AUXV is REAL array, dimension (NB)
*> Auxiliar vector. *> Auxiliary vector.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] F *> \param[in,out] F

View File

@ -67,7 +67,7 @@
*> \param[in] N *> \param[in] N
*> \verbatim *> \verbatim
*> N is INTEGER *> N is INTEGER
*> The order of the matrix H. N .GE. 0. *> The order of the matrix H. N >= 0.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] ILO *> \param[in] ILO
@ -79,12 +79,12 @@
*> \verbatim *> \verbatim
*> IHI is INTEGER *> IHI is INTEGER
*> It is assumed that H is already upper triangular in rows *> It is assumed that H is already upper triangular in rows
*> and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1, *> and columns 1:ILO-1 and IHI+1:N and, if ILO > 1,
*> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a *> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
*> previous call to SGEBAL, and then passed to SGEHRD when the *> previous call to SGEBAL, and then passed to SGEHRD when the
*> matrix output by SGEBAL is reduced to Hessenberg form. *> matrix output by SGEBAL is reduced to Hessenberg form.
*> Otherwise, ILO and IHI should be set to 1 and N, *> Otherwise, ILO and IHI should be set to 1 and N,
*> respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N. *> respectively. If N > 0, then 1 <= ILO <= IHI <= N.
*> If N = 0, then ILO = 1 and IHI = 0. *> If N = 0, then ILO = 1 and IHI = 0.
*> \endverbatim *> \endverbatim
*> *>
@ -97,19 +97,19 @@
*> decomposition (the Schur form); 2-by-2 diagonal blocks *> decomposition (the Schur form); 2-by-2 diagonal blocks
*> (corresponding to complex conjugate pairs of eigenvalues) *> (corresponding to complex conjugate pairs of eigenvalues)
*> are returned in standard form, with H(i,i) = H(i+1,i+1) *> are returned in standard form, with H(i,i) = H(i+1,i+1)
*> and H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and WANTT is *> and H(i+1,i)*H(i,i+1) < 0. If INFO = 0 and WANTT is
*> .FALSE., then the contents of H are unspecified on exit. *> .FALSE., then the contents of H are unspecified on exit.
*> (The output value of H when INFO.GT.0 is given under the *> (The output value of H when INFO > 0 is given under the
*> description of INFO below.) *> description of INFO below.)
*> *>
*> This subroutine may explicitly set H(i,j) = 0 for i.GT.j and *> This subroutine may explicitly set H(i,j) = 0 for i > j and
*> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. *> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] LDH *> \param[in] LDH
*> \verbatim *> \verbatim
*> LDH is INTEGER *> LDH is INTEGER
*> The leading dimension of the array H. LDH .GE. max(1,N). *> The leading dimension of the array H. LDH >= max(1,N).
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WR *> \param[out] WR
@ -125,7 +125,7 @@
*> and WI(ILO:IHI). If two eigenvalues are computed as a *> and WI(ILO:IHI). If two eigenvalues are computed as a
*> complex conjugate pair, they are stored in consecutive *> complex conjugate pair, they are stored in consecutive
*> elements of WR and WI, say the i-th and (i+1)th, with *> elements of WR and WI, say the i-th and (i+1)th, with
*> WI(i) .GT. 0 and WI(i+1) .LT. 0. If WANTT is .TRUE., then *> WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., then
*> the eigenvalues are stored in the same order as on the *> the eigenvalues are stored in the same order as on the
*> diagonal of the Schur form returned in H, with *> diagonal of the Schur form returned in H, with
*> WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal *> WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal
@ -143,7 +143,7 @@
*> IHIZ is INTEGER *> IHIZ is INTEGER
*> Specify the rows of Z to which transformations must be *> Specify the rows of Z to which transformations must be
*> applied if WANTZ is .TRUE.. *> applied if WANTZ is .TRUE..
*> 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N. *> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] Z *> \param[in,out] Z
@ -153,7 +153,7 @@
*> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is *> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
*> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the *> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
*> orthogonal Schur factor of H(ILO:IHI,ILO:IHI). *> orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
*> (The output value of Z when INFO.GT.0 is given under *> (The output value of Z when INFO > 0 is given under
*> the description of INFO below.) *> the description of INFO below.)
*> \endverbatim *> \endverbatim
*> *>
@ -161,7 +161,7 @@
*> \verbatim *> \verbatim
*> LDZ is INTEGER *> LDZ is INTEGER
*> The leading dimension of the array Z. if WANTZ is .TRUE. *> The leading dimension of the array Z. if WANTZ is .TRUE.
*> then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1. *> then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WORK *> \param[out] WORK
@ -174,7 +174,7 @@
*> \param[in] LWORK *> \param[in] LWORK
*> \verbatim *> \verbatim
*> LWORK is INTEGER *> LWORK is INTEGER
*> The dimension of the array WORK. LWORK .GE. max(1,N) *> The dimension of the array WORK. LWORK >= max(1,N)
*> is sufficient, but LWORK typically as large as 6*N may *> is sufficient, but LWORK typically as large as 6*N may
*> be required for optimal performance. A workspace query *> be required for optimal performance. A workspace query
*> to determine the optimal workspace size is recommended. *> to determine the optimal workspace size is recommended.
@ -191,18 +191,18 @@
*> \verbatim *> \verbatim
*> INFO is INTEGER *> INFO is INTEGER
*> = 0: successful exit *> = 0: successful exit
*> .GT. 0: if INFO = i, SLAQR0 failed to compute all of *> > 0: if INFO = i, SLAQR0 failed to compute all of
*> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR *> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
*> and WI contain those eigenvalues which have been *> and WI contain those eigenvalues which have been
*> successfully computed. (Failures are rare.) *> successfully computed. (Failures are rare.)
*> *>
*> If INFO .GT. 0 and WANT is .FALSE., then on exit, *> If INFO > 0 and WANT is .FALSE., then on exit,
*> the remaining unconverged eigenvalues are the eigen- *> the remaining unconverged eigenvalues are the eigen-
*> values of the upper Hessenberg matrix rows and *> values of the upper Hessenberg matrix rows and
*> columns ILO through INFO of the final, output *> columns ILO through INFO of the final, output
*> value of H. *> value of H.
*> *>
*> If INFO .GT. 0 and WANTT is .TRUE., then on exit *> If INFO > 0 and WANTT is .TRUE., then on exit
*> *>
*> (*) (initial value of H)*U = U*(final value of H) *> (*) (initial value of H)*U = U*(final value of H)
*> *>
@ -210,7 +210,7 @@
*> value of H is upper Hessenberg and quasi-triangular *> value of H is upper Hessenberg and quasi-triangular
*> in rows and columns INFO+1 through IHI. *> in rows and columns INFO+1 through IHI.
*> *>
*> If INFO .GT. 0 and WANTZ is .TRUE., then on exit *> If INFO > 0 and WANTZ is .TRUE., then on exit
*> *>
*> (final value of Z(ILO:IHI,ILOZ:IHIZ) *> (final value of Z(ILO:IHI,ILOZ:IHIZ)
*> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U *> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
@ -218,7 +218,7 @@
*> where U is the orthogonal matrix in (*) (regard- *> where U is the orthogonal matrix in (*) (regard-
*> less of the value of WANTT.) *> less of the value of WANTT.)
*> *>
*> If INFO .GT. 0 and WANTZ is .FALSE., then Z is not *> If INFO > 0 and WANTZ is .FALSE., then Z is not
*> accessed. *> accessed.
*> \endverbatim *> \endverbatim
* *
@ -677,7 +677,7 @@
END IF END IF
END IF END IF
* *
* ==== Use up to NS of the the smallest magnatiude * ==== Use up to NS of the the smallest magnitude
* . shifts. If there aren't NS shifts available, * . shifts. If there aren't NS shifts available,
* . then use them all, possibly dropping one to * . then use them all, possibly dropping one to
* . make the number of shifts even. ==== * . make the number of shifts even. ====

View File

@ -69,7 +69,7 @@
*> \verbatim *> \verbatim
*> LDH is INTEGER *> LDH is INTEGER
*> The leading dimension of H as declared in *> The leading dimension of H as declared in
*> the calling procedure. LDH.GE.N *> the calling procedure. LDH >= N
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] SR1 *> \param[in] SR1

View File

@ -103,7 +103,7 @@
*> \param[in] NW *> \param[in] NW
*> \verbatim *> \verbatim
*> NW is INTEGER *> NW is INTEGER
*> Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1). *> Deflation window size. 1 <= NW <= (KBOT-KTOP+1).
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] H *> \param[in,out] H
@ -121,7 +121,7 @@
*> \verbatim *> \verbatim
*> LDH is INTEGER *> LDH is INTEGER
*> Leading dimension of H just as declared in the calling *> Leading dimension of H just as declared in the calling
*> subroutine. N .LE. LDH *> subroutine. N <= LDH
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] ILOZ *> \param[in] ILOZ
@ -133,7 +133,7 @@
*> \verbatim *> \verbatim
*> IHIZ is INTEGER *> IHIZ is INTEGER
*> Specify the rows of Z to which transformations must be *> Specify the rows of Z to which transformations must be
*> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N. *> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] Z *> \param[in,out] Z
@ -149,7 +149,7 @@
*> \verbatim *> \verbatim
*> LDZ is INTEGER *> LDZ is INTEGER
*> The leading dimension of Z just as declared in the *> The leading dimension of Z just as declared in the
*> calling subroutine. 1 .LE. LDZ. *> calling subroutine. 1 <= LDZ.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] NS *> \param[out] NS
@ -194,13 +194,13 @@
*> \verbatim *> \verbatim
*> LDV is INTEGER *> LDV is INTEGER
*> The leading dimension of V just as declared in the *> The leading dimension of V just as declared in the
*> calling subroutine. NW .LE. LDV *> calling subroutine. NW <= LDV
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] NH *> \param[in] NH
*> \verbatim *> \verbatim
*> NH is INTEGER *> NH is INTEGER
*> The number of columns of T. NH.GE.NW. *> The number of columns of T. NH >= NW.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] T *> \param[out] T
@ -212,14 +212,14 @@
*> \verbatim *> \verbatim
*> LDT is INTEGER *> LDT is INTEGER
*> The leading dimension of T just as declared in the *> The leading dimension of T just as declared in the
*> calling subroutine. NW .LE. LDT *> calling subroutine. NW <= LDT
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] NV *> \param[in] NV
*> \verbatim *> \verbatim
*> NV is INTEGER *> NV is INTEGER
*> The number of rows of work array WV available for *> The number of rows of work array WV available for
*> workspace. NV.GE.NW. *> workspace. NV >= NW.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WV *> \param[out] WV
@ -231,7 +231,7 @@
*> \verbatim *> \verbatim
*> LDWV is INTEGER *> LDWV is INTEGER
*> The leading dimension of W just as declared in the *> The leading dimension of W just as declared in the
*> calling subroutine. NW .LE. LDV *> calling subroutine. NW <= LDV
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WORK *> \param[out] WORK

View File

@ -100,7 +100,7 @@
*> \param[in] NW *> \param[in] NW
*> \verbatim *> \verbatim
*> NW is INTEGER *> NW is INTEGER
*> Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1). *> Deflation window size. 1 <= NW <= (KBOT-KTOP+1).
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] H *> \param[in,out] H
@ -118,7 +118,7 @@
*> \verbatim *> \verbatim
*> LDH is INTEGER *> LDH is INTEGER
*> Leading dimension of H just as declared in the calling *> Leading dimension of H just as declared in the calling
*> subroutine. N .LE. LDH *> subroutine. N <= LDH
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] ILOZ *> \param[in] ILOZ
@ -130,7 +130,7 @@
*> \verbatim *> \verbatim
*> IHIZ is INTEGER *> IHIZ is INTEGER
*> Specify the rows of Z to which transformations must be *> Specify the rows of Z to which transformations must be
*> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N. *> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] Z *> \param[in,out] Z
@ -146,7 +146,7 @@
*> \verbatim *> \verbatim
*> LDZ is INTEGER *> LDZ is INTEGER
*> The leading dimension of Z just as declared in the *> The leading dimension of Z just as declared in the
*> calling subroutine. 1 .LE. LDZ. *> calling subroutine. 1 <= LDZ.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] NS *> \param[out] NS
@ -191,13 +191,13 @@
*> \verbatim *> \verbatim
*> LDV is INTEGER *> LDV is INTEGER
*> The leading dimension of V just as declared in the *> The leading dimension of V just as declared in the
*> calling subroutine. NW .LE. LDV *> calling subroutine. NW <= LDV
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] NH *> \param[in] NH
*> \verbatim *> \verbatim
*> NH is INTEGER *> NH is INTEGER
*> The number of columns of T. NH.GE.NW. *> The number of columns of T. NH >= NW.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] T *> \param[out] T
@ -209,14 +209,14 @@
*> \verbatim *> \verbatim
*> LDT is INTEGER *> LDT is INTEGER
*> The leading dimension of T just as declared in the *> The leading dimension of T just as declared in the
*> calling subroutine. NW .LE. LDT *> calling subroutine. NW <= LDT
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] NV *> \param[in] NV
*> \verbatim *> \verbatim
*> NV is INTEGER *> NV is INTEGER
*> The number of rows of work array WV available for *> The number of rows of work array WV available for
*> workspace. NV.GE.NW. *> workspace. NV >= NW.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WV *> \param[out] WV
@ -228,7 +228,7 @@
*> \verbatim *> \verbatim
*> LDWV is INTEGER *> LDWV is INTEGER
*> The leading dimension of W just as declared in the *> The leading dimension of W just as declared in the
*> calling subroutine. NW .LE. LDV *> calling subroutine. NW <= LDV
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WORK *> \param[out] WORK

View File

@ -74,7 +74,7 @@
*> \param[in] N *> \param[in] N
*> \verbatim *> \verbatim
*> N is INTEGER *> N is INTEGER
*> The order of the matrix H. N .GE. 0. *> The order of the matrix H. N >= 0.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] ILO *> \param[in] ILO
@ -86,12 +86,12 @@
*> \verbatim *> \verbatim
*> IHI is INTEGER *> IHI is INTEGER
*> It is assumed that H is already upper triangular in rows *> It is assumed that H is already upper triangular in rows
*> and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1, *> and columns 1:ILO-1 and IHI+1:N and, if ILO > 1,
*> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a *> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
*> previous call to SGEBAL, and then passed to SGEHRD when the *> previous call to SGEBAL, and then passed to SGEHRD when the
*> matrix output by SGEBAL is reduced to Hessenberg form. *> matrix output by SGEBAL is reduced to Hessenberg form.
*> Otherwise, ILO and IHI should be set to 1 and N, *> Otherwise, ILO and IHI should be set to 1 and N,
*> respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N. *> respectively. If N > 0, then 1 <= ILO <= IHI <= N.
*> If N = 0, then ILO = 1 and IHI = 0. *> If N = 0, then ILO = 1 and IHI = 0.
*> \endverbatim *> \endverbatim
*> *>
@ -104,19 +104,19 @@
*> decomposition (the Schur form); 2-by-2 diagonal blocks *> decomposition (the Schur form); 2-by-2 diagonal blocks
*> (corresponding to complex conjugate pairs of eigenvalues) *> (corresponding to complex conjugate pairs of eigenvalues)
*> are returned in standard form, with H(i,i) = H(i+1,i+1) *> are returned in standard form, with H(i,i) = H(i+1,i+1)
*> and H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and WANTT is *> and H(i+1,i)*H(i,i+1) < 0. If INFO = 0 and WANTT is
*> .FALSE., then the contents of H are unspecified on exit. *> .FALSE., then the contents of H are unspecified on exit.
*> (The output value of H when INFO.GT.0 is given under the *> (The output value of H when INFO > 0 is given under the
*> description of INFO below.) *> description of INFO below.)
*> *>
*> This subroutine may explicitly set H(i,j) = 0 for i.GT.j and *> This subroutine may explicitly set H(i,j) = 0 for i > j and
*> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. *> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] LDH *> \param[in] LDH
*> \verbatim *> \verbatim
*> LDH is INTEGER *> LDH is INTEGER
*> The leading dimension of the array H. LDH .GE. max(1,N). *> The leading dimension of the array H. LDH >= max(1,N).
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WR *> \param[out] WR
@ -132,7 +132,7 @@
*> and WI(ILO:IHI). If two eigenvalues are computed as a *> and WI(ILO:IHI). If two eigenvalues are computed as a
*> complex conjugate pair, they are stored in consecutive *> complex conjugate pair, they are stored in consecutive
*> elements of WR and WI, say the i-th and (i+1)th, with *> elements of WR and WI, say the i-th and (i+1)th, with
*> WI(i) .GT. 0 and WI(i+1) .LT. 0. If WANTT is .TRUE., then *> WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., then
*> the eigenvalues are stored in the same order as on the *> the eigenvalues are stored in the same order as on the
*> diagonal of the Schur form returned in H, with *> diagonal of the Schur form returned in H, with
*> WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal *> WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal
@ -150,7 +150,7 @@
*> IHIZ is INTEGER *> IHIZ is INTEGER
*> Specify the rows of Z to which transformations must be *> Specify the rows of Z to which transformations must be
*> applied if WANTZ is .TRUE.. *> applied if WANTZ is .TRUE..
*> 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N. *> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] Z *> \param[in,out] Z
@ -160,7 +160,7 @@
*> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is *> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
*> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the *> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
*> orthogonal Schur factor of H(ILO:IHI,ILO:IHI). *> orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
*> (The output value of Z when INFO.GT.0 is given under *> (The output value of Z when INFO > 0 is given under
*> the description of INFO below.) *> the description of INFO below.)
*> \endverbatim *> \endverbatim
*> *>
@ -168,7 +168,7 @@
*> \verbatim *> \verbatim
*> LDZ is INTEGER *> LDZ is INTEGER
*> The leading dimension of the array Z. if WANTZ is .TRUE. *> The leading dimension of the array Z. if WANTZ is .TRUE.
*> then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1. *> then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WORK *> \param[out] WORK
@ -181,7 +181,7 @@
*> \param[in] LWORK *> \param[in] LWORK
*> \verbatim *> \verbatim
*> LWORK is INTEGER *> LWORK is INTEGER
*> The dimension of the array WORK. LWORK .GE. max(1,N) *> The dimension of the array WORK. LWORK >= max(1,N)
*> is sufficient, but LWORK typically as large as 6*N may *> is sufficient, but LWORK typically as large as 6*N may
*> be required for optimal performance. A workspace query *> be required for optimal performance. A workspace query
*> to determine the optimal workspace size is recommended. *> to determine the optimal workspace size is recommended.
@ -200,18 +200,18 @@
*> \verbatim *> \verbatim
*> INFO is INTEGER *> INFO is INTEGER
*> = 0: successful exit *> = 0: successful exit
*> .GT. 0: if INFO = i, SLAQR4 failed to compute all of *> > 0: if INFO = i, SLAQR4 failed to compute all of
*> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR *> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
*> and WI contain those eigenvalues which have been *> and WI contain those eigenvalues which have been
*> successfully computed. (Failures are rare.) *> successfully computed. (Failures are rare.)
*> *>
*> If INFO .GT. 0 and WANT is .FALSE., then on exit, *> If INFO > 0 and WANT is .FALSE., then on exit,
*> the remaining unconverged eigenvalues are the eigen- *> the remaining unconverged eigenvalues are the eigen-
*> values of the upper Hessenberg matrix rows and *> values of the upper Hessenberg matrix rows and
*> columns ILO through INFO of the final, output *> columns ILO through INFO of the final, output
*> value of H. *> value of H.
*> *>
*> If INFO .GT. 0 and WANTT is .TRUE., then on exit *> If INFO > 0 and WANTT is .TRUE., then on exit
*> *>
*> (*) (initial value of H)*U = U*(final value of H) *> (*) (initial value of H)*U = U*(final value of H)
*> *>
@ -219,7 +219,7 @@
*> value of H is upper Hessenberg and triangular in *> value of H is upper Hessenberg and triangular in
*> rows and columns INFO+1 through IHI. *> rows and columns INFO+1 through IHI.
*> *>
*> If INFO .GT. 0 and WANTZ is .TRUE., then on exit *> If INFO > 0 and WANTZ is .TRUE., then on exit
*> *>
*> (final value of Z(ILO:IHI,ILOZ:IHIZ) *> (final value of Z(ILO:IHI,ILOZ:IHIZ)
*> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U *> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
@ -227,7 +227,7 @@
*> where U is the orthogonal matrix in (*) (regard- *> where U is the orthogonal matrix in (*) (regard-
*> less of the value of WANTT.) *> less of the value of WANTT.)
*> *>
*> If INFO .GT. 0 and WANTZ is .FALSE., then Z is not *> If INFO > 0 and WANTZ is .FALSE., then Z is not
*> accessed. *> accessed.
*> \endverbatim *> \endverbatim
* *
@ -680,7 +680,7 @@
END IF END IF
END IF END IF
* *
* ==== Use up to NS of the the smallest magnatiude * ==== Use up to NS of the the smallest magnitude
* . shifts. If there aren't NS shifts available, * . shifts. If there aren't NS shifts available,
* . then use them all, possibly dropping one to * . then use them all, possibly dropping one to
* . make the number of shifts even. ==== * . make the number of shifts even. ====

View File

@ -133,7 +133,7 @@
*> \verbatim *> \verbatim
*> LDH is INTEGER *> LDH is INTEGER
*> LDH is the leading dimension of H just as declared in the *> LDH is the leading dimension of H just as declared in the
*> calling procedure. LDH.GE.MAX(1,N). *> calling procedure. LDH >= MAX(1,N).
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] ILOZ *> \param[in] ILOZ
@ -145,7 +145,7 @@
*> \verbatim *> \verbatim
*> IHIZ is INTEGER *> IHIZ is INTEGER
*> Specify the rows of Z to which transformations must be *> Specify the rows of Z to which transformations must be
*> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N *> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] Z *> \param[in,out] Z
@ -161,7 +161,7 @@
*> \verbatim *> \verbatim
*> LDZ is INTEGER *> LDZ is INTEGER
*> LDA is the leading dimension of Z just as declared in *> LDA is the leading dimension of Z just as declared in
*> the calling procedure. LDZ.GE.N. *> the calling procedure. LDZ >= N.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] V *> \param[out] V
@ -173,7 +173,7 @@
*> \verbatim *> \verbatim
*> LDV is INTEGER *> LDV is INTEGER
*> LDV is the leading dimension of V as declared in the *> LDV is the leading dimension of V as declared in the
*> calling procedure. LDV.GE.3. *> calling procedure. LDV >= 3.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] U *> \param[out] U
@ -185,33 +185,14 @@
*> \verbatim *> \verbatim
*> LDU is INTEGER *> LDU is INTEGER
*> LDU is the leading dimension of U just as declared in the *> LDU is the leading dimension of U just as declared in the
*> in the calling subroutine. LDU.GE.3*NSHFTS-3. *> in the calling subroutine. LDU >= 3*NSHFTS-3.
*> \endverbatim
*>
*> \param[in] NH
*> \verbatim
*> NH is INTEGER
*> NH is the number of columns in array WH available for
*> workspace. NH.GE.1.
*> \endverbatim
*>
*> \param[out] WH
*> \verbatim
*> WH is REAL array, dimension (LDWH,NH)
*> \endverbatim
*>
*> \param[in] LDWH
*> \verbatim
*> LDWH is INTEGER
*> Leading dimension of WH just as declared in the
*> calling procedure. LDWH.GE.3*NSHFTS-3.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] NV *> \param[in] NV
*> \verbatim *> \verbatim
*> NV is INTEGER *> NV is INTEGER
*> NV is the number of rows in WV agailable for workspace. *> NV is the number of rows in WV agailable for workspace.
*> NV.GE.1. *> NV >= 1.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WV *> \param[out] WV
@ -223,9 +204,28 @@
*> \verbatim *> \verbatim
*> LDWV is INTEGER *> LDWV is INTEGER
*> LDWV is the leading dimension of WV as declared in the *> LDWV is the leading dimension of WV as declared in the
*> in the calling subroutine. LDWV.GE.NV. *> in the calling subroutine. LDWV >= NV.
*> \endverbatim *> \endverbatim
* *
*> \param[in] NH
*> \verbatim
*> NH is INTEGER
*> NH is the number of columns in array WH available for
*> workspace. NH >= 1.
*> \endverbatim
*>
*> \param[out] WH
*> \verbatim
*> WH is REAL array, dimension (LDWH,NH)
*> \endverbatim
*>
*> \param[in] LDWH
*> \verbatim
*> LDWH is INTEGER
*> Leading dimension of WH just as declared in the
*> calling procedure. LDWH >= 3*NSHFTS-3.
*> \endverbatim
*>
* Authors: * Authors:
* ======== * ========
* *

View File

@ -92,6 +92,8 @@
*> K is INTEGER *> K is INTEGER
*> The order of the matrix T (= the number of elementary *> The order of the matrix T (= the number of elementary
*> reflectors whose product defines the block reflector). *> reflectors whose product defines the block reflector).
*> If SIDE = 'L', M >= K >= 0;
*> if SIDE = 'R', N >= K >= 0.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] V *> \param[in] V

View File

@ -94,7 +94,7 @@
*> \param[in] LDC *> \param[in] LDC
*> \verbatim *> \verbatim
*> LDC is INTEGER *> LDC is INTEGER
*> The leading dimension of the array C. LDA >= (1,M). *> The leading dimension of the array C. LDC >= (1,M).
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WORK *> \param[out] WORK

View File

@ -103,7 +103,7 @@
* *
*> \date December 2016 *> \date December 2016
* *
*> \ingroup single_eig *> \ingroup realOTHERauxiliary
* *
* ===================================================================== * =====================================================================
SUBROUTINE SLARFY( UPLO, N, V, INCV, TAU, C, LDC, WORK ) SUBROUTINE SLARFY( UPLO, N, V, INCV, TAU, C, LDC, WORK )

View File

@ -91,7 +91,7 @@
*> RTOL2 is REAL *> RTOL2 is REAL
*> Tolerance for the convergence of the bisection intervals. *> Tolerance for the convergence of the bisection intervals.
*> An interval [LEFT,RIGHT] has converged if *> An interval [LEFT,RIGHT] has converged if
*> RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) *> RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
*> where GAP is the (estimated) distance to the nearest *> where GAP is the (estimated) distance to the nearest
*> eigenvalue. *> eigenvalue.
*> \endverbatim *> \endverbatim
@ -117,7 +117,7 @@
*> WGAP is REAL array, dimension (N-1) *> WGAP is REAL array, dimension (N-1)
*> On input, the (estimated) gaps between consecutive *> On input, the (estimated) gaps between consecutive
*> eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between *> eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between
*> eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST *> eigenvalues I and I+1. Note that if IFIRST = ILAST
*> then WGAP(IFIRST-OFFSET) must be set to ZERO. *> then WGAP(IFIRST-OFFSET) must be set to ZERO.
*> On output, these gaps are refined. *> On output, these gaps are refined.
*> \endverbatim *> \endverbatim

View File

@ -150,7 +150,7 @@
*> RTOL2 is REAL *> RTOL2 is REAL
*> Parameters for bisection. *> Parameters for bisection.
*> An interval [LEFT,RIGHT] has converged if *> An interval [LEFT,RIGHT] has converged if
*> RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) *> RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] SPLTOL *> \param[in] SPLTOL

View File

@ -85,7 +85,7 @@
*> RTOL is REAL *> RTOL is REAL
*> Tolerance for the convergence of the bisection intervals. *> Tolerance for the convergence of the bisection intervals.
*> An interval [LEFT,RIGHT] has converged if *> An interval [LEFT,RIGHT] has converged if
*> RIGHT-LEFT.LT.RTOL*MAX(|LEFT|,|RIGHT|). *> RIGHT-LEFT < RTOL*MAX(|LEFT|,|RIGHT|).
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] OFFSET *> \param[in] OFFSET

View File

@ -149,7 +149,7 @@
*> RTOL2 is REAL *> RTOL2 is REAL
*> Parameters for bisection. *> Parameters for bisection.
*> An interval [LEFT,RIGHT] has converged if *> An interval [LEFT,RIGHT] has converged if
*> RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) *> RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] W *> \param[in,out] W

View File

@ -400,7 +400,7 @@
VL( I ) = VLW( IDXI ) VL( I ) = VLW( IDXI )
50 CONTINUE 50 CONTINUE
* *
* Calculate the allowable deflation tolerence * Calculate the allowable deflation tolerance
* *
EPS = SLAMCH( 'Epsilon' ) EPS = SLAMCH( 'Epsilon' )
TOL = MAX( ABS( ALPHA ), ABS( BETA ) ) TOL = MAX( ABS( ALPHA ), ABS( BETA ) )

View File

@ -60,7 +60,7 @@
*> *>
*> \param[in] X *> \param[in] X
*> \verbatim *> \verbatim
*> X is REAL array, dimension (N) *> X is REAL array, dimension (1+(N-1)*INCX)
*> The vector for which a scaled sum of squares is computed. *> The vector for which a scaled sum of squares is computed.
*> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. *> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n.
*> \endverbatim *> \endverbatim

View File

@ -1,3 +1,4 @@
*> \brief \b SLASWLQ
* *
* Definition: * Definition:
* =========== * ===========
@ -18,9 +19,20 @@
*> *>
*> \verbatim *> \verbatim
*> *>
*> SLASWLQ computes a blocked Short-Wide LQ factorization of a *> SLASWLQ computes a blocked Tall-Skinny LQ factorization of
*> M-by-N matrix A, where N >= M: *> a real M-by-N matrix A for M <= N:
*> A = L * Q *>
*> A = ( L 0 ) * Q,
*>
*> where:
*>
*> Q is a n-by-N orthogonal matrix, stored on exit in an implicit
*> form in the elements above the digonal of the array A and in
*> the elemenst of the array T;
*> L is an lower-triangular M-by-M matrix stored on exit in
*> the elements on and below the diagonal of the array A.
*> 0 is a M-by-(N-M) zero matrix, if M < N, and is not stored.
*>
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -150,10 +162,10 @@
SUBROUTINE SLASWLQ( M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK, SUBROUTINE SLASWLQ( M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK,
$ INFO) $ INFO)
* *
* -- LAPACK computational routine (version 3.8.0) -- * -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
* November 2017 * November 2019
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
INTEGER INFO, LDA, M, N, MB, NB, LWORK, LDT INTEGER INFO, LDA, M, N, MB, NB, LWORK, LDT

View File

@ -284,7 +284,8 @@
* *
* Swap A(I1, I2+1:M) with A(I2, I2+1:M) * Swap A(I1, I2+1:M) with A(I2, I2+1:M)
* *
CALL SSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA, IF( I2.LT.M )
$ CALL SSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
$ A( J1+I2-1, I2+1 ), LDA ) $ A( J1+I2-1, I2+1 ), LDA )
* *
* Swap A(I1, I1) with A(I2,I2) * Swap A(I1, I1) with A(I2,I2)
@ -325,6 +326,7 @@
* Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1), * Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
* where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1) * where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
* *
IF( J.LT.(M-1) ) THEN
IF( A( K, J+1 ).NE.ZERO ) THEN IF( A( K, J+1 ).NE.ZERO ) THEN
ALPHA = ONE / A( K, J+1 ) ALPHA = ONE / A( K, J+1 )
CALL SCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA ) CALL SCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA )
@ -334,6 +336,7 @@
$ A( K, J+2 ), LDA) $ A( K, J+2 ), LDA)
END IF END IF
END IF END IF
END IF
J = J + 1 J = J + 1
GO TO 10 GO TO 10
20 CONTINUE 20 CONTINUE
@ -432,7 +435,8 @@
* *
* Swap A(I2+1:M, I1) with A(I2+1:M, I2) * Swap A(I2+1:M, I1) with A(I2+1:M, I2)
* *
CALL SSWAP( M-I2, A( I2+1, J1+I1-1 ), 1, IF( I2.LT.M )
$ CALL SSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
$ A( I2+1, J1+I2-1 ), 1 ) $ A( I2+1, J1+I2-1 ), 1 )
* *
* Swap A(I1, I1) with A(I2, I2) * Swap A(I1, I1) with A(I2, I2)
@ -473,6 +477,7 @@
* Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1), * Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
* where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1) * where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
* *
IF( J.LT.(M-1) ) THEN
IF( A( J+1, K ).NE.ZERO ) THEN IF( A( J+1, K ).NE.ZERO ) THEN
ALPHA = ONE / A( J+1, K ) ALPHA = ONE / A( J+1, K )
CALL SCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 ) CALL SCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 )
@ -482,6 +487,7 @@
$ A( J+2, K ), LDA ) $ A( J+2, K ), LDA )
END IF END IF
END IF END IF
END IF
J = J + 1 J = J + 1
GO TO 30 GO TO 30
40 CONTINUE 40 CONTINUE

View File

@ -321,7 +321,7 @@
* of A and working backwards, and compute the matrix W = U12*D * of A and working backwards, and compute the matrix W = U12*D
* for use in updating A11 * for use in updating A11
* *
* Initilize the first entry of array E, where superdiagonal * Initialize the first entry of array E, where superdiagonal
* elements of D are stored * elements of D are stored
* *
E( 1 ) = ZERO E( 1 ) = ZERO
@ -649,7 +649,7 @@
* of A and working forwards, and compute the matrix W = L21*D * of A and working forwards, and compute the matrix W = L21*D
* for use in updating A22 * for use in updating A22
* *
* Initilize the unused last entry of the subdiagonal array E. * Initialize the unused last entry of the subdiagonal array E.
* *
E( N ) = ZERO E( N ) = ZERO
* *

View File

@ -85,7 +85,7 @@
*> RHS is REAL array, dimension N. *> RHS is REAL array, dimension N.
*> On entry, RHS contains contributions from other subsystems. *> On entry, RHS contains contributions from other subsystems.
*> On exit, RHS contains the solution of the subsystem with *> On exit, RHS contains the solution of the subsystem with
*> entries acoording to the value of IJOB (see above). *> entries according to the value of IJOB (see above).
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] RDSUM *> \param[in,out] RDSUM
@ -260,7 +260,7 @@
* *
* Solve for U-part, look-ahead for RHS(N) = +-1. This is not done * Solve for U-part, look-ahead for RHS(N) = +-1. This is not done
* in BSOLVE and will hopefully give us a better estimate because * in BSOLVE and will hopefully give us a better estimate because
* any ill-conditioning of the original matrix is transfered to U * any ill-conditioning of the original matrix is transferred to U
* and not to L. U(N, N) is an approximation to sigma_min(LU). * and not to L. U(N, N) is an approximation to sigma_min(LU).
* *
CALL SCOPY( N-1, RHS, 1, XP, 1 ) CALL SCOPY( N-1, RHS, 1, XP, 1 )

View File

@ -1,3 +1,4 @@
*> \brief \b SLATSQR
* *
* Definition: * Definition:
* =========== * ===========
@ -19,8 +20,22 @@
*> \verbatim *> \verbatim
*> *>
*> SLATSQR computes a blocked Tall-Skinny QR factorization of *> SLATSQR computes a blocked Tall-Skinny QR factorization of
*> an M-by-N matrix A, where M >= N: *> a real M-by-N matrix A for M >= N:
*> A = Q * R . *>
*> A = Q * ( R ),
*> ( 0 )
*>
*> where:
*>
*> Q is a M-by-M orthogonal matrix, stored on exit in an implicit
*> form in the elements below the digonal of the array A and in
*> the elemenst of the array T;
*>
*> R is an upper-triangular N-by-N matrix, stored on exit in
*> the elements on and above the diagonal of the array A.
*>
*> 0 is a (M-N)-by-N zero matrix, and is not stored.
*>
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -149,10 +164,10 @@
SUBROUTINE SLATSQR( M, N, MB, NB, A, LDA, T, LDT, WORK, SUBROUTINE SLATSQR( M, N, MB, NB, A, LDA, T, LDT, WORK,
$ LWORK, INFO) $ LWORK, INFO)
* *
* -- LAPACK computational routine (version 3.7.0) -- * -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
* December 2016 * November 2019
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
INTEGER INFO, LDA, M, N, MB, NB, LDT, LWORK INTEGER INFO, LDA, M, N, MB, NB, LDT, LWORK

View File

@ -0,0 +1,306 @@
*> \brief \b SORGTSQR
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download SORGTSQR + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorgtsqr.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorgtsqr.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorgtsqr.f">
*> [TXT]</a>
*>
* Definition:
* ===========
*
* SUBROUTINE SORGTSQR( M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK,
* $ INFO )
*
* .. Scalar Arguments ..
* INTEGER INFO, LDA, LDT, LWORK, M, N, MB, NB
* ..
* .. Array Arguments ..
* REAL A( LDA, * ), T( LDT, * ), WORK( * )
* ..
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> SORGTSQR generates an M-by-N real matrix Q_out with orthonormal columns,
*> which are the first N columns of a product of real orthogonal
*> matrices of order M which are returned by SLATSQR
*>
*> Q_out = first_N_columns_of( Q(1)_in * Q(2)_in * ... * Q(k)_in ).
*>
*> See the documentation for SLATSQR.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> The number of rows of the matrix A. M >= 0.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The number of columns of the matrix A. M >= N >= 0.
*> \endverbatim
*>
*> \param[in] MB
*> \verbatim
*> MB is INTEGER
*> The row block size used by SLATSQR to return
*> arrays A and T. MB > N.
*> (Note that if MB > M, then M is used instead of MB
*> as the row block size).
*> \endverbatim
*>
*> \param[in] NB
*> \verbatim
*> NB is INTEGER
*> The column block size used by SLATSQR to return
*> arrays A and T. NB >= 1.
*> (Note that if NB > N, then N is used instead of NB
*> as the column block size).
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*> A is REAL array, dimension (LDA,N)
*>
*> On entry:
*>
*> The elements on and above the diagonal are not accessed.
*> The elements below the diagonal represent the unit
*> lower-trapezoidal blocked matrix V computed by SLATSQR
*> that defines the input matrices Q_in(k) (ones on the
*> diagonal are not stored) (same format as the output A
*> below the diagonal in SLATSQR).
*>
*> On exit:
*>
*> The array A contains an M-by-N orthonormal matrix Q_out,
*> i.e the columns of A are orthogonal unit vectors.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> The leading dimension of the array A. LDA >= max(1,M).
*> \endverbatim
*>
*> \param[in] T
*> \verbatim
*> T is REAL array,
*> dimension (LDT, N * NIRB)
*> where NIRB = Number_of_input_row_blocks
*> = MAX( 1, CEIL((M-N)/(MB-N)) )
*> Let NICB = Number_of_input_col_blocks
*> = CEIL(N/NB)
*>
*> The upper-triangular block reflectors used to define the
*> input matrices Q_in(k), k=(1:NIRB*NICB). The block
*> reflectors are stored in compact form in NIRB block
*> reflector sequences. Each of NIRB block reflector sequences
*> is stored in a larger NB-by-N column block of T and consists
*> of NICB smaller NB-by-NB upper-triangular column blocks.
*> (same format as the output T in SLATSQR).
*> \endverbatim
*>
*> \param[in] LDT
*> \verbatim
*> LDT is INTEGER
*> The leading dimension of the array T.
*> LDT >= max(1,min(NB1,N)).
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*> (workspace) REAL array, dimension (MAX(2,LWORK))
*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*> The dimension of the array WORK. LWORK >= (M+NB)*N.
*> 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.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
*> = 0: successful exit
*> < 0: if INFO = -i, the i-th argument had an illegal value
*> \endverbatim
*>
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2019
*
*> \ingroup singleOTHERcomputational
*
*> \par Contributors:
* ==================
*>
*> \verbatim
*>
*> November 2019, Igor Kozachenko,
*> Computer Science Division,
*> University of California, Berkeley
*>
*> \endverbatim
*
* =====================================================================
SUBROUTINE SORGTSQR( M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK,
$ INFO )
IMPLICIT NONE
*
* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LDT, LWORK, M, N, MB, NB
* ..
* .. Array Arguments ..
REAL A( LDA, * ), T( LDT, * ), WORK( * )
* ..
*
* =====================================================================
*
* .. Parameters ..
REAL ONE, ZERO
PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
* ..
* .. Local Scalars ..
LOGICAL LQUERY
INTEGER IINFO, LDC, LWORKOPT, LC, LW, NBLOCAL, J
* ..
* .. External Subroutines ..
EXTERNAL SCOPY, SLAMTSQR, SLASET, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC REAL, MAX, MIN
* ..
* .. Executable Statements ..
*
* Test the input parameters
*
LQUERY = LWORK.EQ.-1
INFO = 0
IF( M.LT.0 ) THEN
INFO = -1
ELSE IF( N.LT.0 .OR. M.LT.N ) THEN
INFO = -2
ELSE IF( MB.LE.N ) THEN
INFO = -3
ELSE IF( NB.LT.1 ) THEN
INFO = -4
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -6
ELSE IF( LDT.LT.MAX( 1, MIN( NB, N ) ) ) THEN
INFO = -8
ELSE
*
* Test the input LWORK for the dimension of the array WORK.
* This workspace is used to store array C(LDC, N) and WORK(LWORK)
* in the call to DLAMTSQR. See the documentation for DLAMTSQR.
*
IF( LWORK.LT.2 .AND. (.NOT.LQUERY) ) THEN
INFO = -10
ELSE
*
* Set block size for column blocks
*
NBLOCAL = MIN( NB, N )
*
* LWORK = -1, then set the size for the array C(LDC,N)
* in DLAMTSQR call and set the optimal size of the work array
* WORK(LWORK) in DLAMTSQR call.
*
LDC = M
LC = LDC*N
LW = N * NBLOCAL
*
LWORKOPT = LC+LW
*
IF( ( LWORK.LT.MAX( 1, LWORKOPT ) ).AND.(.NOT.LQUERY) ) THEN
INFO = -10
END IF
END IF
*
END IF
*
* Handle error in the input parameters and return workspace query.
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'SORGTSQR', -INFO )
RETURN
ELSE IF ( LQUERY ) THEN
WORK( 1 ) = REAL( LWORKOPT )
RETURN
END IF
*
* Quick return if possible
*
IF( MIN( M, N ).EQ.0 ) THEN
WORK( 1 ) = REAL( LWORKOPT )
RETURN
END IF
*
* (1) Form explicitly the tall-skinny M-by-N left submatrix Q1_in
* of M-by-M orthogonal matrix Q_in, which is implicitly stored in
* the subdiagonal part of input array A and in the input array T.
* Perform by the following operation using the routine DLAMTSQR.
*
* Q1_in = Q_in * ( I ), where I is a N-by-N identity matrix,
* ( 0 ) 0 is a (M-N)-by-N zero matrix.
*
* (1a) Form M-by-N matrix in the array WORK(1:LDC*N) with ones
* on the diagonal and zeros elsewhere.
*
CALL SLASET( 'F', M, N, ZERO, ONE, WORK, LDC )
*
* (1b) On input, WORK(1:LDC*N) stores ( I );
* ( 0 )
*
* On output, WORK(1:LDC*N) stores Q1_in.
*
CALL SLAMTSQR( 'L', 'N', M, N, N, MB, NBLOCAL, A, LDA, T, LDT,
$ WORK, LDC, WORK( LC+1 ), LW, IINFO )
*
* (2) Copy the result from the part of the work array (1:M,1:N)
* with the leading dimension LDC that starts at WORK(1) into
* the output array A(1:M,1:N) column-by-column.
*
DO J = 1, N
CALL SCOPY( M, WORK( (J-1)*LDC + 1 ), 1, A( 1, J ), 1 )
END DO
*
WORK( 1 ) = REAL( LWORKOPT )
RETURN
*
* End of SORGTSQR
*
END

View File

@ -0,0 +1,439 @@
*> \brief \b SORHR_COL
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download SORHR_COL + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorhr_col.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorhr_col.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorhr_col.f">
*> [TXT]</a>
*>
* Definition:
* ===========
*
* SUBROUTINE SORHR_COL( M, N, NB, A, LDA, T, LDT, D, INFO )
*
* .. Scalar Arguments ..
* INTEGER INFO, LDA, LDT, M, N, NB
* ..
* .. Array Arguments ..
* REAL A( LDA, * ), D( * ), T( LDT, * )
* ..
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> SORHR_COL takes an M-by-N real matrix Q_in with orthonormal columns
*> as input, stored in A, and performs Householder Reconstruction (HR),
*> i.e. reconstructs Householder vectors V(i) implicitly representing
*> another M-by-N matrix Q_out, with the property that Q_in = Q_out*S,
*> where S is an N-by-N diagonal matrix with diagonal entries
*> equal to +1 or -1. The Householder vectors (columns V(i) of V) are
*> stored in A on output, and the diagonal entries of S are stored in D.
*> Block reflectors are also returned in T
*> (same output format as SGEQRT).
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> The number of rows of the matrix A. M >= 0.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The number of columns of the matrix A. M >= N >= 0.
*> \endverbatim
*>
*> \param[in] NB
*> \verbatim
*> NB is INTEGER
*> The column block size to be used in the reconstruction
*> of Householder column vector blocks in the array A and
*> corresponding block reflectors in the array T. NB >= 1.
*> (Note that if NB > N, then N is used instead of NB
*> as the column block size.)
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*> A is REAL array, dimension (LDA,N)
*>
*> On entry:
*>
*> The array A contains an M-by-N orthonormal matrix Q_in,
*> i.e the columns of A are orthogonal unit vectors.
*>
*> On exit:
*>
*> The elements below the diagonal of A represent the unit
*> lower-trapezoidal matrix V of Householder column vectors
*> V(i). The unit diagonal entries of V are not stored
*> (same format as the output below the diagonal in A from
*> SGEQRT). The matrix T and the matrix V stored on output
*> in A implicitly define Q_out.
*>
*> The elements above the diagonal contain the factor U
*> of the "modified" LU-decomposition:
*> Q_in - ( S ) = V * U
*> ( 0 )
*> where 0 is a (M-N)-by-(M-N) zero matrix.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> The leading dimension of the array A. LDA >= max(1,M).
*> \endverbatim
*>
*> \param[out] T
*> \verbatim
*> T is REAL array,
*> dimension (LDT, N)
*>
*> Let NOCB = Number_of_output_col_blocks
*> = CEIL(N/NB)
*>
*> On exit, T(1:NB, 1:N) contains NOCB upper-triangular
*> block reflectors used to define Q_out stored in compact
*> form as a sequence of upper-triangular NB-by-NB column
*> blocks (same format as the output T in SGEQRT).
*> The matrix T and the matrix V stored on output in A
*> implicitly define Q_out. NOTE: The lower triangles
*> below the upper-triangular blcoks will be filled with
*> zeros. See Further Details.
*> \endverbatim
*>
*> \param[in] LDT
*> \verbatim
*> LDT is INTEGER
*> The leading dimension of the array T.
*> LDT >= max(1,min(NB,N)).
*> \endverbatim
*>
*> \param[out] D
*> \verbatim
*> D is REAL array, dimension min(M,N).
*> The elements can be only plus or minus one.
*>
*> D(i) is constructed as D(i) = -SIGN(Q_in_i(i,i)), where
*> 1 <= i <= min(M,N), and Q_in_i is Q_in after performing
*> i-1 steps of “modified” Gaussian elimination.
*> See Further Details.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
*> = 0: successful exit
*> < 0: if INFO = -i, the i-th argument had an illegal value
*> \endverbatim
*>
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> The computed M-by-M orthogonal factor Q_out is defined implicitly as
*> a product of orthogonal matrices Q_out(i). Each Q_out(i) is stored in
*> the compact WY-representation format in the corresponding blocks of
*> matrices V (stored in A) and T.
*>
*> The M-by-N unit lower-trapezoidal matrix V stored in the M-by-N
*> matrix A contains the column vectors V(i) in NB-size column
*> blocks VB(j). For example, VB(1) contains the columns
*> V(1), V(2), ... V(NB). NOTE: The unit entries on
*> the diagonal of Y are not stored in A.
*>
*> The number of column blocks is
*>
*> NOCB = Number_of_output_col_blocks = CEIL(N/NB)
*>
*> where each block is of order NB except for the last block, which
*> is of order LAST_NB = N - (NOCB-1)*NB.
*>
*> For example, if M=6, N=5 and NB=2, the matrix V is
*>
*>
*> V = ( VB(1), VB(2), VB(3) ) =
*>
*> = ( 1 )
*> ( v21 1 )
*> ( v31 v32 1 )
*> ( v41 v42 v43 1 )
*> ( v51 v52 v53 v54 1 )
*> ( v61 v62 v63 v54 v65 )
*>
*>
*> For each of the column blocks VB(i), an upper-triangular block
*> reflector TB(i) is computed. These blocks are stored as
*> a sequence of upper-triangular column blocks in the NB-by-N
*> matrix T. The size of each TB(i) block is NB-by-NB, except
*> for the last block, whose size is LAST_NB-by-LAST_NB.
*>
*> For example, if M=6, N=5 and NB=2, the matrix T is
*>
*> T = ( TB(1), TB(2), TB(3) ) =
*>
*> = ( t11 t12 t13 t14 t15 )
*> ( t22 t24 )
*>
*>
*> The M-by-M factor Q_out is given as a product of NOCB
*> orthogonal M-by-M matrices Q_out(i).
*>
*> Q_out = Q_out(1) * Q_out(2) * ... * Q_out(NOCB),
*>
*> where each matrix Q_out(i) is given by the WY-representation
*> using corresponding blocks from the matrices V and T:
*>
*> Q_out(i) = I - VB(i) * TB(i) * (VB(i))**T,
*>
*> where I is the identity matrix. Here is the formula with matrix
*> dimensions:
*>
*> Q(i){M-by-M} = I{M-by-M} -
*> VB(i){M-by-INB} * TB(i){INB-by-INB} * (VB(i))**T {INB-by-M},
*>
*> where INB = NB, except for the last block NOCB
*> for which INB=LAST_NB.
*>
*> =====
*> NOTE:
*> =====
*>
*> If Q_in is the result of doing a QR factorization
*> B = Q_in * R_in, then:
*>
*> B = (Q_out*S) * R_in = Q_out * (S * R_in) = O_out * R_out.
*>
*> So if one wants to interpret Q_out as the result
*> of the QR factorization of B, then corresponding R_out
*> should be obtained by R_out = S * R_in, i.e. some rows of R_in
*> should be multiplied by -1.
*>
*> For the details of the algorithm, see [1].
*>
*> [1] "Reconstructing Householder vectors from tall-skinny QR",
*> G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
*> E. Solomonik, J. Parallel Distrib. Comput.,
*> vol. 85, pp. 3-31, 2015.
*> \endverbatim
*>
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2019
*
*> \ingroup singleOTHERcomputational
*
*> \par Contributors:
* ==================
*>
*> \verbatim
*>
*> November 2019, Igor Kozachenko,
*> Computer Science Division,
*> University of California, Berkeley
*>
*> \endverbatim
*
* =====================================================================
SUBROUTINE SORHR_COL( M, N, NB, A, LDA, T, LDT, D, INFO )
IMPLICIT NONE
*
* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LDT, M, N, NB
* ..
* .. Array Arguments ..
REAL A( LDA, * ), D( * ), T( LDT, * )
* ..
*
* =====================================================================
*
* .. Parameters ..
REAL ONE, ZERO
PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
* ..
* .. Local Scalars ..
INTEGER I, IINFO, J, JB, JBTEMP1, JBTEMP2, JNB,
$ NPLUSONE
* ..
* .. External Subroutines ..
EXTERNAL SCOPY, SLAORHR_COL_GETRFNP, SSCAL, STRSM, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
* ..
* .. Executable Statements ..
*
* Test the input parameters
*
INFO = 0
IF( M.LT.0 ) THEN
INFO = -1
ELSE IF( N.LT.0 .OR. N.GT.M ) THEN
INFO = -2
ELSE IF( NB.LT.1 ) THEN
INFO = -3
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -5
ELSE IF( LDT.LT.MAX( 1, MIN( NB, N ) ) ) THEN
INFO = -7
END IF
*
* Handle error in the input parameters.
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'SORHR_COL', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( MIN( M, N ).EQ.0 ) THEN
RETURN
END IF
*
* On input, the M-by-N matrix A contains the orthogonal
* M-by-N matrix Q_in.
*
* (1) Compute the unit lower-trapezoidal V (ones on the diagonal
* are not stored) by performing the "modified" LU-decomposition.
*
* Q_in - ( S ) = V * U = ( V1 ) * U,
* ( 0 ) ( V2 )
*
* where 0 is an (M-N)-by-N zero matrix.
*
* (1-1) Factor V1 and U.
CALL SLAORHR_COL_GETRFNP( N, N, A, LDA, D, IINFO )
*
* (1-2) Solve for V2.
*
IF( M.GT.N ) THEN
CALL STRSM( 'R', 'U', 'N', 'N', M-N, N, ONE, A, LDA,
$ A( N+1, 1 ), LDA )
END IF
*
* (2) Reconstruct the block reflector T stored in T(1:NB, 1:N)
* as a sequence of upper-triangular blocks with NB-size column
* blocking.
*
* Loop over the column blocks of size NB of the array A(1:M,1:N)
* and the array T(1:NB,1:N), JB is the column index of a column
* block, JNB is the column block size at each step JB.
*
NPLUSONE = N + 1
DO JB = 1, N, NB
*
* (2-0) Determine the column block size JNB.
*
JNB = MIN( NPLUSONE-JB, NB )
*
* (2-1) Copy the upper-triangular part of the current JNB-by-JNB
* diagonal block U(JB) (of the N-by-N matrix U) stored
* in A(JB:JB+JNB-1,JB:JB+JNB-1) into the upper-triangular part
* of the current JNB-by-JNB block T(1:JNB,JB:JB+JNB-1)
* column-by-column, total JNB*(JNB+1)/2 elements.
*
JBTEMP1 = JB - 1
DO J = JB, JB+JNB-1
CALL SCOPY( J-JBTEMP1, A( JB, J ), 1, T( 1, J ), 1 )
END DO
*
* (2-2) Perform on the upper-triangular part of the current
* JNB-by-JNB diagonal block U(JB) (of the N-by-N matrix U) stored
* in T(1:JNB,JB:JB+JNB-1) the following operation in place:
* (-1)*U(JB)*S(JB), i.e the result will be stored in the upper-
* triangular part of T(1:JNB,JB:JB+JNB-1). This multiplication
* of the JNB-by-JNB diagonal block U(JB) by the JNB-by-JNB
* diagonal block S(JB) of the N-by-N sign matrix S from the
* right means changing the sign of each J-th column of the block
* U(JB) according to the sign of the diagonal element of the block
* S(JB), i.e. S(J,J) that is stored in the array element D(J).
*
DO J = JB, JB+JNB-1
IF( D( J ).EQ.ONE ) THEN
CALL SSCAL( J-JBTEMP1, -ONE, T( 1, J ), 1 )
END IF
END DO
*
* (2-3) Perform the triangular solve for the current block
* matrix X(JB):
*
* X(JB) * (A(JB)**T) = B(JB), where:
*
* A(JB)**T is a JNB-by-JNB unit upper-triangular
* coefficient block, and A(JB)=V1(JB), which
* is a JNB-by-JNB unit lower-triangular block
* stored in A(JB:JB+JNB-1,JB:JB+JNB-1).
* The N-by-N matrix V1 is the upper part
* of the M-by-N lower-trapezoidal matrix V
* stored in A(1:M,1:N);
*
* B(JB) is a JNB-by-JNB upper-triangular right-hand
* side block, B(JB) = (-1)*U(JB)*S(JB), and
* B(JB) is stored in T(1:JNB,JB:JB+JNB-1);
*
* X(JB) is a JNB-by-JNB upper-triangular solution
* block, X(JB) is the upper-triangular block
* reflector T(JB), and X(JB) is stored
* in T(1:JNB,JB:JB+JNB-1).
*
* In other words, we perform the triangular solve for the
* upper-triangular block T(JB):
*
* T(JB) * (V1(JB)**T) = (-1)*U(JB)*S(JB).
*
* Even though the blocks X(JB) and B(JB) are upper-
* triangular, the routine STRSM will access all JNB**2
* elements of the square T(1:JNB,JB:JB+JNB-1). Therefore,
* we need to set to zero the elements of the block
* T(1:JNB,JB:JB+JNB-1) below the diagonal before the call
* to STRSM.
*
* (2-3a) Set the elements to zero.
*
JBTEMP2 = JB - 2
DO J = JB, JB+JNB-2
DO I = J-JBTEMP2, NB
T( I, J ) = ZERO
END DO
END DO
*
* (2-3b) Perform the triangular solve.
*
CALL STRSM( 'R', 'L', 'T', 'U', JNB, JNB, ONE,
$ A( JB, JB ), LDA, T( 1, JB ), LDT )
*
END DO
*
RETURN
*
* End of SORHR_COL
*
END

View File

@ -135,7 +135,7 @@
*> \param[in,out] S *> \param[in,out] S
*> \verbatim *> \verbatim
*> S is REAL array, dimension (N) *> S is REAL array, dimension (N)
*> The row scale factors for A. If EQUED = 'Y', A is multiplied on *> The scale factors for A. If EQUED = 'Y', A is multiplied on
*> the left and right by diag(S). S is an input argument if FACT = *> the left and right by diag(S). S is an input argument if FACT =
*> 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED *> 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED
*> = 'Y', each element of S must be positive. If S is output, each *> = 'Y', each element of S must be positive. If S is output, each
@ -263,7 +263,7 @@
*> information as described below. There currently are up to three *> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If *> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned. *> the first (:,N_ERR_BNDS) entries are returned.
*> *>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
@ -299,14 +299,14 @@
*> \param[in] NPARAMS *> \param[in] NPARAMS
*> \verbatim *> \verbatim
*> NPARAMS is INTEGER *> NPARAMS is INTEGER
*> Specifies the number of parameters set in PARAMS. If .LE. 0, the *> Specifies the number of parameters set in PARAMS. If <= 0, the
*> PARAMS array is never referenced and default values are used. *> PARAMS array is never referenced and default values are used.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] PARAMS *> \param[in,out] PARAMS
*> \verbatim *> \verbatim
*> PARAMS is REAL array, dimension NPARAMS *> PARAMS is REAL array, dimension NPARAMS
*> Specifies algorithm parameters. If an entry is .LT. 0.0, then *> Specifies algorithm parameters. If an entry is < 0.0, then
*> that entry will be filled with default value used for that *> that entry will be filled with default value used for that
*> parameter. Only positions up to NPARAMS are accessed; defaults *> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters. *> are used for higher-numbered parameters.

View File

@ -366,7 +366,7 @@
*> information as described below. There currently are up to three *> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If *> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned. *> the first (:,N_ERR_BNDS) entries are returned.
*> *>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
@ -402,14 +402,14 @@
*> \param[in] NPARAMS *> \param[in] NPARAMS
*> \verbatim *> \verbatim
*> NPARAMS is INTEGER *> NPARAMS is INTEGER
*> Specifies the number of parameters set in PARAMS. If .LE. 0, the *> Specifies the number of parameters set in PARAMS. If <= 0, the
*> PARAMS array is never referenced and default values are used. *> PARAMS array is never referenced and default values are used.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] PARAMS *> \param[in,out] PARAMS
*> \verbatim *> \verbatim
*> PARAMS is REAL array, dimension NPARAMS *> PARAMS is REAL array, dimension NPARAMS
*> Specifies algorithm parameters. If an entry is .LT. 0.0, then *> Specifies algorithm parameters. If an entry is < 0.0, then
*> that entry will be filled with default value used for that *> that entry will be filled with default value used for that
*> parameter. Only positions up to NPARAMS are accessed; defaults *> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters. *> are used for higher-numbered parameters.

View File

@ -124,7 +124,7 @@
*> LDVT is INTEGER. *> LDVT is INTEGER.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is REAL array. Workspace of size nb. *> WORK is REAL array. Workspace of size nb.
*> \endverbatim *> \endverbatim

View File

@ -233,13 +233,13 @@
*> \param[in,out] TRYRAC *> \param[in,out] TRYRAC
*> \verbatim *> \verbatim
*> TRYRAC is LOGICAL *> TRYRAC is LOGICAL
*> If TRYRAC.EQ..TRUE., indicates that the code should check whether *> If TRYRAC = .TRUE., indicates that the code should check whether
*> the tridiagonal matrix defines its eigenvalues to high relative *> the tridiagonal matrix defines its eigenvalues to high relative
*> accuracy. If so, the code uses relative-accuracy preserving *> accuracy. If so, the code uses relative-accuracy preserving
*> algorithms that might be (a bit) slower depending on the matrix. *> algorithms that might be (a bit) slower depending on the matrix.
*> If the matrix does not define its eigenvalues to high relative *> If the matrix does not define its eigenvalues to high relative
*> accuracy, the code can uses possibly faster algorithms. *> accuracy, the code can uses possibly faster algorithms.
*> If TRYRAC.EQ..FALSE., the code is not required to guarantee *> If TRYRAC = .FALSE., the code is not required to guarantee
*> relatively accurate eigenvalues and can use the fastest possible *> relatively accurate eigenvalues and can use the fastest possible
*> techniques. *> techniques.
*> On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix *> On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix

View File

@ -291,7 +291,7 @@
* *
* Convert PERMUTATIONS and IPIV * Convert PERMUTATIONS and IPIV
* *
* Apply permutaions to submatrices of upper part of A * Apply permutations to submatrices of upper part of A
* in factorization order where i decreases from N to 1 * in factorization order where i decreases from N to 1
* *
I = N I = N
@ -344,7 +344,7 @@
* *
* Revert PERMUTATIONS and IPIV * Revert PERMUTATIONS and IPIV
* *
* Apply permutaions to submatrices of upper part of A * Apply permutations to submatrices of upper part of A
* in reverse factorization order where i increases from 1 to N * in reverse factorization order where i increases from 1 to N
* *
I = 1 I = 1
@ -435,7 +435,7 @@
* *
* Convert PERMUTATIONS and IPIV * Convert PERMUTATIONS and IPIV
* *
* Apply permutaions to submatrices of lower part of A * Apply permutations to submatrices of lower part of A
* in factorization order where k increases from 1 to N * in factorization order where k increases from 1 to N
* *
I = 1 I = 1
@ -488,7 +488,7 @@
* *
* Revert PERMUTATIONS and IPIV * Revert PERMUTATIONS and IPIV
* *
* Apply permutaions to submatrices of lower part of A * Apply permutations to submatrices of lower part of A
* in reverse factorization order where i decreases from N to 1 * in reverse factorization order where i decreases from N to 1
* *
I = N I = N

View File

@ -282,7 +282,7 @@
* *
* Convert PERMUTATIONS * Convert PERMUTATIONS
* *
* Apply permutaions to submatrices of upper part of A * Apply permutations to submatrices of upper part of A
* in factorization order where i decreases from N to 1 * in factorization order where i decreases from N to 1
* *
I = N I = N
@ -333,7 +333,7 @@
* *
* Revert PERMUTATIONS * Revert PERMUTATIONS
* *
* Apply permutaions to submatrices of upper part of A * Apply permutations to submatrices of upper part of A
* in reverse factorization order where i increases from 1 to N * in reverse factorization order where i increases from 1 to N
* *
I = 1 I = 1
@ -423,7 +423,7 @@
* *
* Convert PERMUTATIONS * Convert PERMUTATIONS
* *
* Apply permutaions to submatrices of lower part of A * Apply permutations to submatrices of lower part of A
* in factorization order where i increases from 1 to N * in factorization order where i increases from 1 to N
* *
I = 1 I = 1
@ -474,7 +474,7 @@
* *
* Revert PERMUTATIONS * Revert PERMUTATIONS
* *
* Apply permutaions to submatrices of lower part of A * Apply permutations to submatrices of lower part of A
* in reverse factorization order where i decreases from N to 1 * in reverse factorization order where i decreases from N to 1
* *
I = N I = N

View File

@ -317,7 +317,7 @@
IF( .NOT.WANTZ ) THEN IF( .NOT.WANTZ ) THEN
CALL SSTERF( N, W, WORK( INDE ), INFO ) CALL SSTERF( N, W, WORK( INDE ), INFO )
ELSE ELSE
* Not available in this release, and agrument checking should not * Not available in this release, and argument checking should not
* let it getting here * let it getting here
RETURN RETURN
CALL SORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ), CALL SORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),

View File

@ -385,7 +385,7 @@
IF( .NOT.WANTZ ) THEN IF( .NOT.WANTZ ) THEN
CALL SSTERF( N, W, WORK( INDE ), INFO ) CALL SSTERF( N, W, WORK( INDE ), INFO )
ELSE ELSE
* Not available in this release, and agrument checking should not * Not available in this release, and argument checking should not
* let it getting here * let it getting here
RETURN RETURN
CALL SSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N, CALL SSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,

View File

@ -271,7 +271,7 @@
*> information as described below. There currently are up to three *> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If *> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned. *> the first (:,N_ERR_BNDS) entries are returned.
*> *>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
@ -307,14 +307,14 @@
*> \param[in] NPARAMS *> \param[in] NPARAMS
*> \verbatim *> \verbatim
*> NPARAMS is INTEGER *> NPARAMS is INTEGER
*> Specifies the number of parameters set in PARAMS. If .LE. 0, the *> Specifies the number of parameters set in PARAMS. If <= 0, the
*> PARAMS array is never referenced and default values are used. *> PARAMS array is never referenced and default values are used.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] PARAMS *> \param[in,out] PARAMS
*> \verbatim *> \verbatim
*> PARAMS is REAL array, dimension NPARAMS *> PARAMS is REAL array, dimension NPARAMS
*> Specifies algorithm parameters. If an entry is .LT. 0.0, then *> Specifies algorithm parameters. If an entry is < 0.0, then
*> that entry will be filled with default value used for that *> that entry will be filled with default value used for that
*> parameter. Only positions up to NPARAMS are accessed; defaults *> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters. *> are used for higher-numbered parameters.

View File

@ -42,7 +42,7 @@
*> matrices. *> matrices.
*> *>
*> Aasen's algorithm is used to factor A as *> Aasen's algorithm is used to factor A as
*> A = U * T * U**T, if UPLO = 'U', or *> A = U**T * T * U, if UPLO = 'U', or
*> A = L * T * L**T, if UPLO = 'L', *> A = L * T * L**T, if UPLO = 'L',
*> where U (or L) is a product of permutation and unit upper (lower) *> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is symmetric tridiagonal. The factored *> triangular matrices, and T is symmetric tridiagonal. The factored
@ -86,7 +86,7 @@
*> *>
*> On exit, if INFO = 0, the tridiagonal matrix T and the *> On exit, if INFO = 0, the tridiagonal matrix T and the
*> multipliers used to obtain the factor U or L from the *> multipliers used to obtain the factor U or L from the
*> factorization A = U*T*U**T or A = L*T*L**T as computed by *> factorization A = U**T*T*U or A = L*T*L**T as computed by
*> SSYTRF. *> SSYTRF.
*> \endverbatim *> \endverbatim
*> *>
@ -229,7 +229,7 @@
RETURN RETURN
END IF END IF
* *
* Compute the factorization A = U*T*U**T or A = L*T*L**T. * Compute the factorization A = U**T*T*U or A = L*T*L**T.
* *
CALL SSYTRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) CALL SSYTRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
IF( INFO.EQ.0 ) THEN IF( INFO.EQ.0 ) THEN

View File

@ -44,7 +44,7 @@
*> matrices. *> matrices.
*> *>
*> Aasen's 2-stage algorithm is used to factor A as *> Aasen's 2-stage algorithm is used to factor A as
*> A = U * T * U**T, if UPLO = 'U', or *> A = U**T * T * U, if UPLO = 'U', or
*> A = L * T * L**T, if UPLO = 'L', *> A = L * T * L**T, if UPLO = 'L',
*> where U (or L) is a product of permutation and unit upper (lower) *> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is symmetric and band. The matrix T is *> triangular matrices, and T is symmetric and band. The matrix T is
@ -258,7 +258,7 @@
END IF END IF
* *
* *
* Compute the factorization A = U*T*U**T or A = L*T*L**T. * Compute the factorization A = U**T*T*U or A = L*T*L**T.
* *
CALL SSYTRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2, CALL SSYTRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2,
$ WORK, LWORK, INFO ) $ WORK, LWORK, INFO )

View File

@ -377,7 +377,7 @@
*> information as described below. There currently are up to three *> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If *> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned. *> the first (:,N_ERR_BNDS) entries are returned.
*> *>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
@ -413,14 +413,14 @@
*> \param[in] NPARAMS *> \param[in] NPARAMS
*> \verbatim *> \verbatim
*> NPARAMS is INTEGER *> NPARAMS is INTEGER
*> Specifies the number of parameters set in PARAMS. If .LE. 0, the *> Specifies the number of parameters set in PARAMS. If <= 0, the
*> PARAMS array is never referenced and default values are used. *> PARAMS array is never referenced and default values are used.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] PARAMS *> \param[in,out] PARAMS
*> \verbatim *> \verbatim
*> PARAMS is REAL array, dimension NPARAMS *> PARAMS is REAL array, dimension NPARAMS
*> Specifies algorithm parameters. If an entry is .LT. 0.0, then *> Specifies algorithm parameters. If an entry is < 0.0, then
*> that entry will be filled with default value used for that *> that entry will be filled with default value used for that
*> parameter. Only positions up to NPARAMS are accessed; defaults *> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters. *> are used for higher-numbered parameters.

View File

@ -312,7 +312,7 @@
* *
* Factorize A as U*D*U**T using the upper triangle of A * Factorize A as U*D*U**T using the upper triangle of A
* *
* Initilize the first entry of array E, where superdiagonal * Initialize the first entry of array E, where superdiagonal
* elements of D are stored * elements of D are stored
* *
E( 1 ) = ZERO E( 1 ) = ZERO
@ -623,7 +623,7 @@
* *
* Factorize A as L*D*L**T using the lower triangle of A * Factorize A as L*D*L**T using the lower triangle of A
* *
* Initilize the unused last entry of the subdiagonal array E. * Initialize the unused last entry of the subdiagonal array E.
* *
E( N ) = ZERO E( N ) = ZERO
* *

View File

@ -123,23 +123,22 @@
*> *>
*> \param[out] HOUS2 *> \param[out] HOUS2
*> \verbatim *> \verbatim
*> HOUS2 is REAL array, dimension LHOUS2, that *> HOUS2 is REAL array, dimension (LHOUS2)
*> store the Householder representation of the stage2 *> Stores the Householder representation of the stage2
*> band to tridiagonal. *> band to tridiagonal.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] LHOUS2 *> \param[in] LHOUS2
*> \verbatim *> \verbatim
*> LHOUS2 is INTEGER *> LHOUS2 is INTEGER
*> The dimension of the array HOUS2. LHOUS2 = MAX(1, dimension) *> The dimension of the array HOUS2.
*> If LWORK = -1, or LHOUS2 = -1, *> If LWORK = -1, or LHOUS2 = -1,
*> then a query is assumed; the routine *> then a query is assumed; the routine
*> only calculates the optimal size of the HOUS2 array, returns *> only calculates the optimal size of the HOUS2 array, returns
*> this value as the first entry of the HOUS2 array, and no error *> this value as the first entry of the HOUS2 array, and no error
*> message related to LHOUS2 is issued by XERBLA. *> message related to LHOUS2 is issued by XERBLA.
*> LHOUS2 = MAX(1, dimension) where *> If VECT='N', LHOUS2 = max(1, 4*n);
*> dimension = 4*N if VECT='N' *> if VECT='V', option not yet available.
*> not available now if VECT='H'
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WORK *> \param[out] WORK

View File

@ -50,9 +50,9 @@
* Arguments: * Arguments:
* ========== * ==========
* *
*> \param[in] STAGE *> \param[in] STAGE1
*> \verbatim *> \verbatim
*> STAGE is CHARACTER*1 *> STAGE1 is CHARACTER*1
*> = 'N': "No": to mention that the stage 1 of the reduction *> = 'N': "No": to mention that the stage 1 of the reduction
*> from dense to band using the ssytrd_sy2sb routine *> from dense to band using the ssytrd_sy2sb routine
*> was not called before this routine to reproduce AB. *> was not called before this routine to reproduce AB.
@ -481,7 +481,7 @@
* *
* Call the kernel * Call the kernel
* *
#if defined(_OPENMP) && _OPENMP >= 201307 #if defined(_OPENMP)
IF( TTYPE.NE.1 ) THEN IF( TTYPE.NE.1 ) THEN
!$OMP TASK DEPEND(in:WORK(MYID+SHIFT-1)) !$OMP TASK DEPEND(in:WORK(MYID+SHIFT-1))
!$OMP$ DEPEND(in:WORK(MYID-1)) !$OMP$ DEPEND(in:WORK(MYID-1))

View File

@ -39,7 +39,7 @@
*> the Bunch-Kaufman diagonal pivoting method. The form of the *> the Bunch-Kaufman diagonal pivoting method. The form of the
*> factorization is *> factorization is
*> *>
*> A = U*D*U**T or A = L*D*L**T *> A = U**T*D*U or A = L*D*L**T
*> *>
*> where U (or L) is a product of permutation and unit upper (lower) *> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and D is symmetric and block diagonal with *> triangular matrices, and D is symmetric and block diagonal with
@ -144,7 +144,7 @@
*> *>
*> \verbatim *> \verbatim
*> *>
*> If UPLO = 'U', then A = U*D*U**T, where *> If UPLO = 'U', then A = U**T*D*U, where
*> U = P(n)*U(n)* ... *P(k)U(k)* ..., *> U = P(n)*U(n)* ... *P(k)U(k)* ...,
*> i.e., U is a product of terms P(k)*U(k), where k decreases from n to *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to
*> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
@ -262,7 +262,7 @@
* *
IF( UPPER ) THEN IF( UPPER ) THEN
* *
* Factorize A as U*D*U**T using the upper triangle of A * Factorize A as U**T*D*U using the upper triangle of A
* *
* K is the main loop index, decreasing from N to 1 in steps of * K is the main loop index, decreasing from N to 1 in steps of
* KB, where KB is the number of columns factorized by SLASYF; * KB, where KB is the number of columns factorized by SLASYF;

View File

@ -37,7 +37,7 @@
*> SSYTRF_AA computes the factorization of a real symmetric matrix A *> SSYTRF_AA computes the factorization of a real symmetric matrix A
*> using the Aasen's algorithm. The form of the factorization is *> using the Aasen's algorithm. The form of the factorization is
*> *>
*> A = U*T*U**T or A = L*T*L**T *> A = U**T*T*U or A = L*T*L**T
*> *>
*> where U (or L) is a product of permutation and unit upper (lower) *> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is a symmetric tridiagonal matrix. *> triangular matrices, and T is a symmetric tridiagonal matrix.
@ -223,7 +223,7 @@
IF( UPPER ) THEN IF( UPPER ) THEN
* *
* ..................................................... * .....................................................
* Factorize A as L*D*L**T using the upper triangle of A * Factorize A as U**T*D*U using the upper triangle of A
* ..................................................... * .....................................................
* *
* Copy first row A(1, 1:N) into H(1:n) (stored in WORK(1:N)) * Copy first row A(1, 1:N) into H(1:n) (stored in WORK(1:N))
@ -256,7 +256,7 @@
$ A( MAX(1, J), J+1 ), LDA, $ A( MAX(1, J), J+1 ), LDA,
$ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ) ) $ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ) )
* *
* Ajust IPIV and apply it back (J-th step picks (J+1)-th pivot) * Adjust IPIV and apply it back (J-th step picks (J+1)-th pivot)
* *
DO J2 = J+2, MIN(N, J+JB+1) DO J2 = J+2, MIN(N, J+JB+1)
IPIV( J2 ) = IPIV( J2 ) + J IPIV( J2 ) = IPIV( J2 ) + J
@ -375,7 +375,7 @@
$ A( J+1, MAX(1, J) ), LDA, $ A( J+1, MAX(1, J) ), LDA,
$ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ) ) $ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ) )
* *
* Ajust IPIV and apply it back (J-th step picks (J+1)-th pivot) * Adjust IPIV and apply it back (J-th step picks (J+1)-th pivot)
* *
DO J2 = J+2, MIN(N, J+JB+1) DO J2 = J+2, MIN(N, J+JB+1)
IPIV( J2 ) = IPIV( J2 ) + J IPIV( J2 ) = IPIV( J2 ) + J

View File

@ -38,7 +38,7 @@
*> SSYTRF_AA_2STAGE computes the factorization of a real symmetric matrix A *> SSYTRF_AA_2STAGE computes the factorization of a real symmetric matrix A
*> using the Aasen's algorithm. The form of the factorization is *> using the Aasen's algorithm. The form of the factorization is
*> *>
*> A = U*T*U**T or A = L*T*L**T *> A = U**T*T*U or A = L*T*L**T
*> *>
*> where U (or L) is a product of permutation and unit upper (lower) *> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is a symmetric band matrix with the *> triangular matrices, and T is a symmetric band matrix with the
@ -275,7 +275,7 @@
IF( UPPER ) THEN IF( UPPER ) THEN
* *
* ..................................................... * .....................................................
* Factorize A as L*D*L**T using the upper triangle of A * Factorize A as U**T*D*U using the upper triangle of A
* ..................................................... * .....................................................
* *
DO J = 0, NT-1 DO J = 0, NT-1
@ -443,10 +443,12 @@ c END IF
CALL SSWAP( K-1, A( (J+1)*NB+1, I1 ), 1, CALL SSWAP( K-1, A( (J+1)*NB+1, I1 ), 1,
$ A( (J+1)*NB+1, I2 ), 1 ) $ A( (J+1)*NB+1, I2 ), 1 )
* > Swap A(I1+1:M, I1) with A(I2, I1+1:M) * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
CALL SSWAP( I2-I1-1, A( I1, I1+1 ), LDA, IF( I2.GT.(I1+1) )
$ CALL SSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
$ A( I1+1, I2 ), 1 ) $ A( I1+1, I2 ), 1 )
* > Swap A(I2+1:M, I1) with A(I2+1:M, I2) * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
CALL SSWAP( N-I2, A( I1, I2+1 ), LDA, IF( I2.LT.N )
$ CALL SSWAP( N-I2, A( I1, I2+1 ), LDA,
$ A( I2, I2+1 ), LDA ) $ A( I2, I2+1 ), LDA )
* > Swap A(I1, I1) with A(I2, I2) * > Swap A(I1, I1) with A(I2, I2)
PIV = A( I1, I1 ) PIV = A( I1, I1 )
@ -616,10 +618,12 @@ c END IF
CALL SSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA, CALL SSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA,
$ A( I2, (J+1)*NB+1 ), LDA ) $ A( I2, (J+1)*NB+1 ), LDA )
* > Swap A(I1+1:M, I1) with A(I2, I1+1:M) * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
CALL SSWAP( I2-I1-1, A( I1+1, I1 ), 1, IF( I2.GT.(I1+1) )
$ CALL SSWAP( I2-I1-1, A( I1+1, I1 ), 1,
$ A( I2, I1+1 ), LDA ) $ A( I2, I1+1 ), LDA )
* > Swap A(I2+1:M, I1) with A(I2+1:M, I2) * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
CALL SSWAP( N-I2, A( I2+1, I1 ), 1, IF( I2.LT.N )
$ CALL SSWAP( N-I2, A( I2+1, I1 ), 1,
$ A( I2+1, I2 ), 1 ) $ A( I2+1, I2 ), 1 )
* > Swap A(I1, I1) with A(I2, I2) * > Swap A(I1, I1) with A(I2, I2)
PIV = A( I1, I1 ) PIV = A( I1, I1 )

View File

@ -62,7 +62,7 @@
*> \param[in,out] A *> \param[in,out] A
*> \verbatim *> \verbatim
*> A is REAL array, dimension (LDA,N) *> A is REAL array, dimension (LDA,N)
*> On entry, the NB diagonal matrix D and the multipliers *> On entry, the block diagonal matrix D and the multipliers
*> used to obtain the factor U or L as computed by SSYTRF. *> used to obtain the factor U or L as computed by SSYTRF.
*> *>
*> On exit, if INFO = 0, the (symmetric) inverse of the original *> On exit, if INFO = 0, the (symmetric) inverse of the original
@ -82,7 +82,7 @@
*> \param[in] IPIV *> \param[in] IPIV
*> \verbatim *> \verbatim
*> IPIV is INTEGER array, dimension (N) *> IPIV is INTEGER array, dimension (N)
*> Details of the interchanges and the NB structure of D *> Details of the interchanges and the block structure of D
*> as determined by SSYTRF. *> as determined by SSYTRF.
*> \endverbatim *> \endverbatim
*> *>

View File

@ -37,7 +37,7 @@
*> \verbatim *> \verbatim
*> *>
*> SSYTRS_AA solves a system of linear equations A*X = B with a real *> SSYTRS_AA solves a system of linear equations A*X = B with a real
*> symmetric matrix A using the factorization A = U*T*U**T or *> symmetric matrix A using the factorization A = U**T*T*U or
*> A = L*T*L**T computed by SSYTRF_AA. *> A = L*T*L**T computed by SSYTRF_AA.
*> \endverbatim *> \endverbatim
* *
@ -49,7 +49,7 @@
*> UPLO is CHARACTER*1 *> UPLO is CHARACTER*1
*> Specifies whether the details of the factorization are stored *> Specifies whether the details of the factorization are stored
*> as an upper or lower triangular matrix. *> as an upper or lower triangular matrix.
*> = 'U': Upper triangular, form is A = U*T*U**T; *> = 'U': Upper triangular, form is A = U**T*T*U;
*> = 'L': Lower triangular, form is A = L*T*L**T. *> = 'L': Lower triangular, form is A = L*T*L**T.
*> \endverbatim *> \endverbatim
*> *>
@ -97,14 +97,16 @@
*> The leading dimension of the array B. LDB >= max(1,N). *> The leading dimension of the array B. LDB >= max(1,N).
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is DOUBLE array, dimension (MAX(1,LWORK)) *> WORK is REAL array, dimension (MAX(1,LWORK))
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] LWORK *> \param[in] LWORK
*> \verbatim *> \verbatim
*> LWORK is INTEGER, LWORK >= MAX(1,3*N-2). *> LWORK is INTEGER
*> The dimension of the array WORK. LWORK >= max(1,3*N-2).
*> \endverbatim
*> *>
*> \param[out] INFO *> \param[out] INFO
*> \verbatim *> \verbatim
@ -198,9 +200,13 @@
* *
IF( UPPER ) THEN IF( UPPER ) THEN
* *
* Solve A*X = B, where A = U*T*U**T. * Solve A*X = B, where A = U**T*T*U.
* *
* Pivot, P**T * B * 1) Forward substitution with U**T
*
IF( N.GT.1 ) THEN
*
* Pivot, P**T * B -> B
* *
K = 1 K = 1
DO WHILE ( K.LE.N ) DO WHILE ( K.LE.N )
@ -210,12 +216,15 @@
K = K + 1 K = K + 1
END DO END DO
* *
* Compute (U \P**T * B) -> B [ (U \P**T * B) ] * Compute U**T \ B -> B [ (U**T \P**T * B) ]
* *
CALL STRSM('L', 'U', 'T', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA, CALL STRSM( 'L', 'U', 'T', 'U', N-1, NRHS, ONE, A( 1, 2 ),
$ B( 2, 1 ), LDB) $ LDA, B( 2, 1 ), LDB)
END IF
* *
* Compute T \ B -> B [ T \ (U \P**T * B) ] * 2) Solve with triangular matrix T
*
* Compute T \ B -> B [ T \ (U**T \P**T * B) ]
* *
CALL SLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1) CALL SLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1)
IF( N.GT.1 ) THEN IF( N.GT.1 ) THEN
@ -225,13 +234,17 @@
CALL SGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB, CALL SGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB,
$ INFO) $ INFO)
* *
* 3) Backward substitution with U
* *
* Compute (U**T \ B) -> B [ U**T \ (T \ (U \P**T * B) ) ] IF( N.GT.1 ) THEN
* *
CALL STRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA,
$ B(2, 1), LDB)
* *
* Pivot, P * B [ P * (U**T \ (T \ (U \P**T * B) )) ] * Compute U \ B -> B [ U \ (T \ (U**T \P**T * B) ) ]
*
CALL STRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ),
$ LDA, B(2, 1), LDB)
*
* Pivot, P * B -> B [ P * (U \ (T \ (U**T \P**T * B) )) ]
* *
K = N K = N
DO WHILE ( K.GE.1 ) DO WHILE ( K.GE.1 )
@ -240,12 +253,17 @@
$ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
K = K - 1 K = K - 1
END DO END DO
END IF
* *
ELSE ELSE
* *
* Solve A*X = B, where A = L*T*L**T. * Solve A*X = B, where A = L*T*L**T.
* *
* Pivot, P**T * B * 1) Forward substitution with L
*
IF( N.GT.1 ) THEN
*
* Pivot, P**T * B -> B
* *
K = 1 K = 1
DO WHILE ( K.LE.N ) DO WHILE ( K.LE.N )
@ -255,10 +273,13 @@
K = K + 1 K = K + 1
END DO END DO
* *
* Compute (L \P**T * B) -> B [ (L \P**T * B) ] * Compute L \ B -> B [ (L \P**T * B) ]
* *
CALL STRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1), LDA, CALL STRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1),
$ B(2, 1), LDB) $ LDA, B(2, 1), LDB)
END IF
*
* 2) Solve with triangular matrix T
* *
* Compute T \ B -> B [ T \ (L \P**T * B) ] * Compute T \ B -> B [ T \ (L \P**T * B) ]
* *
@ -270,12 +291,16 @@
CALL SGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB, CALL SGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB,
$ INFO) $ INFO)
* *
* Compute (L**T \ B) -> B [ L**T \ (T \ (L \P**T * B) ) ] * 3) Backward substitution with L**T
* *
CALL STRSM( 'L', 'L', 'T', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA, IF( N.GT.1 ) THEN
$ B( 2, 1 ), LDB)
* *
* Pivot, P * B [ P * (L**T \ (T \ (L \P**T * B) )) ] * Compute L**T \ B -> B [ L**T \ (T \ (L \P**T * B) ) ]
*
CALL STRSM( 'L', 'L', 'T', 'U', N-1, NRHS, ONE, A( 2, 1 ),
$ LDA, B( 2, 1 ), LDB)
*
* Pivot, P * B -> B [ P * (L**T \ (T \ (L \P**T * B) )) ]
* *
K = N K = N
DO WHILE ( K.GE.1 ) DO WHILE ( K.GE.1 )
@ -284,6 +309,7 @@
$ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
K = K - 1 K = K - 1
END DO END DO
END IF
* *
END IF END IF
* *

View File

@ -36,7 +36,7 @@
*> \verbatim *> \verbatim
*> *>
*> SSYTRS_AA_2STAGE solves a system of linear equations A*X = B with a real *> SSYTRS_AA_2STAGE solves a system of linear equations A*X = B with a real
*> symmetric matrix A using the factorization A = U*T*U**T or *> symmetric matrix A using the factorization A = U**T*T*U or
*> A = L*T*L**T computed by SSYTRF_AA_2STAGE. *> A = L*T*L**T computed by SSYTRF_AA_2STAGE.
*> \endverbatim *> \endverbatim
* *
@ -48,7 +48,7 @@
*> UPLO is CHARACTER*1 *> UPLO is CHARACTER*1
*> Specifies whether the details of the factorization are stored *> Specifies whether the details of the factorization are stored
*> as an upper or lower triangular matrix. *> as an upper or lower triangular matrix.
*> = 'U': Upper triangular, form is A = U*T*U**T; *> = 'U': Upper triangular, form is A = U**T*T*U;
*> = 'L': Lower triangular, form is A = L*T*L**T. *> = 'L': Lower triangular, form is A = L*T*L**T.
*> \endverbatim *> \endverbatim
*> *>
@ -208,15 +208,15 @@
* *
IF( UPPER ) THEN IF( UPPER ) THEN
* *
* Solve A*X = B, where A = U*T*U**T. * Solve A*X = B, where A = U**T*T*U.
* *
IF( N.GT.NB ) THEN IF( N.GT.NB ) THEN
* *
* Pivot, P**T * B * Pivot, P**T * B -> B
* *
CALL SLASWP( NRHS, B, LDB, NB+1, N, IPIV, 1 ) CALL SLASWP( NRHS, B, LDB, NB+1, N, IPIV, 1 )
* *
* Compute (U**T \P**T * B) -> B [ (U**T \P**T * B) ] * Compute (U**T \ B) -> B [ (U**T \P**T * B) ]
* *
CALL STRSM( 'L', 'U', 'T', 'U', N-NB, NRHS, ONE, A(1, NB+1), CALL STRSM( 'L', 'U', 'T', 'U', N-NB, NRHS, ONE, A(1, NB+1),
$ LDA, B(NB+1, 1), LDB) $ LDA, B(NB+1, 1), LDB)
@ -234,7 +234,7 @@
CALL STRSM( 'L', 'U', 'N', 'U', N-NB, NRHS, ONE, A(1, NB+1), CALL STRSM( 'L', 'U', 'N', 'U', N-NB, NRHS, ONE, A(1, NB+1),
$ LDA, B(NB+1, 1), LDB) $ LDA, B(NB+1, 1), LDB)
* *
* Pivot, P * B [ P * (U \ (T \ (U**T \P**T * B) )) ] * Pivot, P * B -> B [ P * (U \ (T \ (U**T \P**T * B) )) ]
* *
CALL SLASWP( NRHS, B, LDB, NB+1, N, IPIV, -1 ) CALL SLASWP( NRHS, B, LDB, NB+1, N, IPIV, -1 )
* *
@ -246,11 +246,11 @@
* *
IF( N.GT.NB ) THEN IF( N.GT.NB ) THEN
* *
* Pivot, P**T * B * Pivot, P**T * B -> B
* *
CALL SLASWP( NRHS, B, LDB, NB+1, N, IPIV, 1 ) CALL SLASWP( NRHS, B, LDB, NB+1, N, IPIV, 1 )
* *
* Compute (L \P**T * B) -> B [ (L \P**T * B) ] * Compute (L \ B) -> B [ (L \P**T * B) ]
* *
CALL STRSM( 'L', 'L', 'N', 'U', N-NB, NRHS, ONE, A(NB+1, 1), CALL STRSM( 'L', 'L', 'N', 'U', N-NB, NRHS, ONE, A(NB+1, 1),
$ LDA, B(NB+1, 1), LDB) $ LDA, B(NB+1, 1), LDB)
@ -268,7 +268,7 @@
CALL STRSM( 'L', 'L', 'T', 'U', N-NB, NRHS, ONE, A(NB+1, 1), CALL STRSM( 'L', 'L', 'T', 'U', N-NB, NRHS, ONE, A(NB+1, 1),
$ LDA, B(NB+1, 1), LDB) $ LDA, B(NB+1, 1), LDB)
* *
* Pivot, P * B [ P * (L**T \ (T \ (L \P**T * B) )) ] * Pivot, P * B -> B [ P * (L**T \ (T \ (L \P**T * B) )) ]
* *
CALL SLASWP( NRHS, B, LDB, NB+1, N, IPIV, -1 ) CALL SLASWP( NRHS, B, LDB, NB+1, N, IPIV, -1 )
* *

View File

@ -71,7 +71,7 @@
*> R * B**T + L * E**T = scale * -F *> R * B**T + L * E**T = scale * -F
*> *>
*> This case is used to compute an estimate of Dif[(A, D), (B, E)] = *> This case is used to compute an estimate of Dif[(A, D), (B, E)] =
*> sigma_min(Z) using reverse communicaton with SLACON. *> sigma_min(Z) using reverse communication with SLACON.
*> *>
*> STGSY2 also (IJOB >= 1) contributes to the computation in STGSYL *> STGSY2 also (IJOB >= 1) contributes to the computation in STGSYL
*> of an upper bound on the separation between to matrix pairs. Then *> of an upper bound on the separation between to matrix pairs. Then
@ -85,7 +85,7 @@
*> \param[in] TRANS *> \param[in] TRANS
*> \verbatim *> \verbatim
*> TRANS is CHARACTER*1 *> TRANS is CHARACTER*1
*> = 'N', solve the generalized Sylvester equation (1). *> = 'N': solve the generalized Sylvester equation (1).
*> = 'T': solve the 'transposed' system (3). *> = 'T': solve the 'transposed' system (3).
*> \endverbatim *> \endverbatim
*> *>

View File

@ -88,8 +88,8 @@
*> \param[in] TRANS *> \param[in] TRANS
*> \verbatim *> \verbatim
*> TRANS is CHARACTER*1 *> TRANS is CHARACTER*1
*> = 'N', solve the generalized Sylvester equation (1). *> = 'N': solve the generalized Sylvester equation (1).
*> = 'T', solve the 'transposed' system (3). *> = 'T': solve the 'transposed' system (3).
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] IJOB *> \param[in] IJOB

View File

@ -94,7 +94,7 @@
*> *>
*> \param[in] V *> \param[in] V
*> \verbatim *> \verbatim
*> V is REAL array, dimension (LDA,K) *> V is REAL array, dimension (LDV,K)
*> The i-th row must contain the vector which defines the *> The i-th row must contain the vector which defines the
*> elementary reflector H(i), for i = 1,2,...,k, as returned by *> elementary reflector H(i), for i = 1,2,...,k, as returned by
*> DTPLQT in B. See Further Details. *> DTPLQT in B. See Further Details.

View File

@ -94,7 +94,7 @@
*> *>
*> \param[in] V *> \param[in] V
*> \verbatim *> \verbatim
*> V is REAL array, dimension (LDA,K) *> V is REAL array, dimension (LDV,K)
*> The i-th column must contain the vector which defines the *> The i-th column must contain the vector which defines the
*> elementary reflector H(i), for i = 1,2,...,k, as returned by *> elementary reflector H(i), for i = 1,2,...,k, as returned by
*> CTPQRT in B. See Further Details. *> CTPQRT in B. See Further Details.

View File

@ -152,8 +152,8 @@
*> \verbatim *> \verbatim
*> LDA is INTEGER *> LDA is INTEGER
*> The leading dimension of the array A. *> The leading dimension of the array A.
*> If SIDE = 'L', LDC >= max(1,K); *> If SIDE = 'L', LDA >= max(1,K);
*> If SIDE = 'R', LDC >= max(1,M). *> If SIDE = 'R', LDA >= max(1,M).
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] B *> \param[in,out] B