INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, & N4 ) ! ! -- LAPACK auxiliary routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! September 30, 1994 ! ! .. Scalar Arguments .. CHARACTER*( * ) NAME, OPTS INTEGER ISPEC, N1, N2, N3, N4 ! .. ! ! Purpose ! ======= ! ! ILAENV is called from the LAPACK routines to choose problem-dependent ! parameters for the local environment. See ISPEC for a description of ! the parameters. ! ! This version provides a set of parameters which should give good, ! but not optimal, performance on many of the currently available ! computers. Users are encouraged to modify this subroutine to set ! the tuning parameters for their particular machine using the option ! and problem size information in the arguments. ! ! This routine will not function correctly if it is converted to all ! lower case. Converting it to all upper case is allowed. ! ! Arguments ! ========= ! ! ISPEC (input) INTEGER ! Specifies the parameter to be returned as the value of ! ILAENV. ! = 1: the optimal blocksize; if this value is 1, an unblocked ! algorithm will give the best performance. ! = 2: the minimum block size for which the block routine ! should be used; if the usable block size is less than ! this value, an unblocked routine should be used. ! = 3: the crossover point (in a block routine, for N less ! than this value, an unblocked routine should be used) ! = 4: the number of shifts, used in the nonsymmetric ! eigenvalue routines ! = 5: the minimum column dimension for blocking to be used; ! rectangular blocks must have dimension at least k by m, ! where k is given by ILAENV(2,...) and m by ILAENV(5,...) ! = 6: the crossover point for the SVD (when reducing an m by n ! matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds ! this value, a QR factorization is used first to reduce ! the matrix to a triangular form.) ! = 7: the number of processors ! = 8: the crossover point for the multishift QR and QZ methods ! for nonsymmetric eigenvalue problems. ! ! NAME (input) CHARACTER*(*) ! The name of the calling subroutine, in either upper case or ! lower case. ! ! OPTS (input) CHARACTER*(*) ! The character options to the subroutine NAME, concatenated ! into a single character string. For example, UPLO = 'U', ! TRANS = 'T', and DIAG = 'N' for a triangular routine would ! be specified as OPTS = 'UTN'. ! ! N1 (input) INTEGER ! N2 (input) INTEGER ! N3 (input) INTEGER ! N4 (input) INTEGER ! Problem dimensions for the subroutine NAME; these may not all ! be required. ! ! (ILAENV) (output) INTEGER ! >= 0: the value of the parameter specified by ISPEC ! < 0: if ILAENV = -k, the k-th argument had an illegal value. ! ! Further Details ! =============== ! ! The following conventions have been used when calling ILAENV from the ! LAPACK routines: ! 1) OPTS is a concatenation of all of the character options to ! subroutine NAME, in the same order that they appear in the ! argument list for NAME, even if they are not used in determining ! the value of the parameter specified by ISPEC. ! 2) The problem dimensions N1, N2, N3, N4 are specified in the order ! that they appear in the argument list for NAME. N1 is used ! first, N2 second, and so on, and unused problem dimensions are ! passed a value of -1. ! 3) The parameter value returned by ILAENV is checked for validity in ! the calling subroutine. For example, ILAENV is used to retrieve ! the optimal blocksize for STRTRI as follows: ! ! NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) ! IF( NB<=1 ) NB = MAX( 1, N ) ! ! ===================================================================== ! ! .. Local Scalars .. LOGICAL CNAME, SNAME CHARACTER*1 C1 CHARACTER*2 C2, C4 CHARACTER*3 C3 CHARACTER*6 SUBNAM INTEGER I, IC, IZ, NB, NBMIN, NX ! .. ! .. Intrinsic Functions .. INTRINSIC CHAR, ICHAR, INT, MIN, REAL ! .. ! .. Executable Statements .. ! GO TO ( 100, 100, 100, 400, 500, 600, 700, 800 ) ISPEC ! ! Invalid value for ISPEC ! ILAENV = -1 RETURN ! 100 CONTINUE ! ! Convert NAME to upper case if the first character is lower case. ! ILAENV = 1 SUBNAM = NAME IC = ICHAR( SUBNAM( 1:1 ) ) IZ = ICHAR( 'Z' ) IF( IZ==90 .OR. IZ==122 ) THEN ! ! ASCII character set ! IF( IC>=97 .AND. IC<=122 ) THEN SUBNAM( 1:1 ) = CHAR( IC-32 ) DO 10 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( IC>=97 .AND. IC<=122 ) & SUBNAM( I:I ) = CHAR( IC-32 ) 10 CONTINUE END IF ! ELSE IF( IZ==233 .OR. IZ==169 ) THEN ! ! EBCDIC character set ! IF( ( IC>=129 .AND. IC<=137 ) .OR. & ( IC>=145 .AND. IC<=153 ) .OR. & ( IC>=162 .AND. IC<=169 ) ) THEN SUBNAM( 1:1 ) = CHAR( IC+64 ) DO 20 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( ( IC>=129 .AND. IC<=137 ) .OR. & ( IC>=145 .AND. IC<=153 ) .OR. & ( IC>=162 .AND. IC<=169 ) ) & SUBNAM( I:I ) = CHAR( IC+64 ) 20 CONTINUE END IF ! ELSE IF( IZ==218 .OR. IZ==250 ) THEN ! ! Prime machines: ASCII+128 ! IF( IC>=225 .AND. IC<=250 ) THEN SUBNAM( 1:1 ) = CHAR( IC-32 ) DO 30 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( IC>=225 .AND. IC<=250 ) & SUBNAM( I:I ) = CHAR( IC-32 ) 30 CONTINUE END IF END IF ! C1 = SUBNAM( 1:1 ) SNAME = C1=='S' .OR. C1=='D' CNAME = C1=='C' .OR. C1=='Z' IF( .NOT.( CNAME .OR. SNAME ) ) & RETURN C2 = SUBNAM( 2:3 ) C3 = SUBNAM( 4:6 ) C4 = C3( 2:3 ) ! GO TO ( 110, 200, 300 ) ISPEC ! 110 CONTINUE ! ! ISPEC = 1: block size ! ! In these examples, separate code is provided for setting NB for ! real and complex. We assume that NB will take the same value in ! single or double precision. ! NB = 1 ! IF( C2=='GE' ) THEN IF( C3=='TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF ELSE IF( C3=='QRF' .OR. C3=='RQF' .OR. C3=='LQF' .OR. & C3=='QLF' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3=='HRD' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3=='BRD' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3=='TRI' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2=='PO' ) THEN IF( C3=='TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2=='SY' ) THEN IF( C3=='TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF ELSE IF( SNAME .AND. C3=='TRD' ) THEN NB = 1 ELSE IF( SNAME .AND. C3=='GST' ) THEN NB = 64 END IF ELSE IF( CNAME .AND. C2=='HE' ) THEN IF( C3=='TRF' ) THEN NB = 64 ELSE IF( C3=='TRD' ) THEN NB = 1 ELSE IF( C3=='GST' ) THEN NB = 64 END IF ELSE IF( SNAME .AND. C2=='OR' ) THEN IF( C3( 1:1 )=='G' ) THEN IF( C4=='QR' .OR. C4=='RQ' .OR. C4=='LQ' .OR. & C4=='QL' .OR. C4=='HR' .OR. C4=='TR' .OR. & C4=='BR' ) THEN NB = 32 END IF ELSE IF( C3( 1:1 )=='M' ) THEN IF( C4=='QR' .OR. C4=='RQ' .OR. C4=='LQ' .OR. & C4=='QL' .OR. C4=='HR' .OR. C4=='TR' .OR. & C4=='BR' ) THEN NB = 32 END IF END IF ELSE IF( CNAME .AND. C2=='UN' ) THEN IF( C3( 1:1 )=='G' ) THEN IF( C4=='QR' .OR. C4=='RQ' .OR. C4=='LQ' .OR. & C4=='QL' .OR. C4=='HR' .OR. C4=='TR' .OR. & C4=='BR' ) THEN NB = 32 END IF ELSE IF( C3( 1:1 )=='M' ) THEN IF( C4=='QR' .OR. C4=='RQ' .OR. C4=='LQ' .OR. & C4=='QL' .OR. C4=='HR' .OR. C4=='TR' .OR. & C4=='BR' ) THEN NB = 32 END IF END IF ELSE IF( C2=='GB' ) THEN IF( C3=='TRF' ) THEN IF( SNAME ) THEN IF( N4<=64 ) THEN NB = 1 ELSE NB = 32 END IF ELSE IF( N4<=64 ) THEN NB = 1 ELSE NB = 32 END IF END IF END IF ELSE IF( C2=='PB' ) THEN IF( C3=='TRF' ) THEN IF( SNAME ) THEN IF( N2<=64 ) THEN NB = 1 ELSE NB = 32 END IF ELSE IF( N2<=64 ) THEN NB = 1 ELSE NB = 32 END IF END IF END IF ELSE IF( C2=='TR' ) THEN IF( C3=='TRI' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2=='LA' ) THEN IF( C3=='UUM' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( SNAME .AND. C2=='ST' ) THEN IF( C3=='EBZ' ) THEN NB = 1 END IF END IF ILAENV = NB RETURN ! 200 CONTINUE ! ! ISPEC = 2: minimum block size ! NBMIN = 2 IF( C2=='GE' ) THEN IF( C3=='QRF' .OR. C3=='RQF' .OR. C3=='LQF' .OR. & C3=='QLF' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3=='HRD' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3=='BRD' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3=='TRI' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF END IF ELSE IF( C2=='SY' ) THEN IF( C3=='TRF' ) THEN IF( SNAME ) THEN NBMIN = 8 ELSE NBMIN = 8 END IF ELSE IF( SNAME .AND. C3=='TRD' ) THEN NBMIN = 2 END IF ELSE IF( CNAME .AND. C2=='HE' ) THEN IF( C3=='TRD' ) THEN NBMIN = 2 END IF ELSE IF( SNAME .AND. C2=='OR' ) THEN IF( C3( 1:1 )=='G' ) THEN IF( C4=='QR' .OR. C4=='RQ' .OR. C4=='LQ' .OR. & C4=='QL' .OR. C4=='HR' .OR. C4=='TR' .OR. & C4=='BR' ) THEN NBMIN = 2 END IF ELSE IF( C3( 1:1 )=='M' ) THEN IF( C4=='QR' .OR. C4=='RQ' .OR. C4=='LQ' .OR. & C4=='QL' .OR. C4=='HR' .OR. C4=='TR' .OR. & C4=='BR' ) THEN NBMIN = 2 END IF END IF ELSE IF( CNAME .AND. C2=='UN' ) THEN IF( C3( 1:1 )=='G' ) THEN IF( C4=='QR' .OR. C4=='RQ' .OR. C4=='LQ' .OR. & C4=='QL' .OR. C4=='HR' .OR. C4=='TR' .OR. & C4=='BR' ) THEN NBMIN = 2 END IF ELSE IF( C3( 1:1 )=='M' ) THEN IF( C4=='QR' .OR. C4=='RQ' .OR. C4=='LQ' .OR. & C4=='QL' .OR. C4=='HR' .OR. C4=='TR' .OR. & C4=='BR' ) THEN NBMIN = 2 END IF END IF END IF ILAENV = NBMIN RETURN ! 300 CONTINUE ! ! ISPEC = 3: crossover point ! NX = 0 IF( C2=='GE' ) THEN IF( C3=='QRF' .OR. C3=='RQF' .OR. C3=='LQF' .OR. & C3=='QLF' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF ELSE IF( C3=='HRD' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF ELSE IF( C3=='BRD' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF END IF ELSE IF( C2=='SY' ) THEN IF( SNAME .AND. C3=='TRD' ) THEN NX = 1 END IF ELSE IF( CNAME .AND. C2=='HE' ) THEN IF( C3=='TRD' ) THEN NX = 1 END IF ELSE IF( SNAME .AND. C2=='OR' ) THEN IF( C3( 1:1 )=='G' ) THEN IF( C4=='QR' .OR. C4=='RQ' .OR. C4=='LQ' .OR. & C4=='QL' .OR. C4=='HR' .OR. C4=='TR' .OR. & C4=='BR' ) THEN NX = 128 END IF END IF ELSE IF( CNAME .AND. C2=='UN' ) THEN IF( C3( 1:1 )=='G' ) THEN IF( C4=='QR' .OR. C4=='RQ' .OR. C4=='LQ' .OR. & C4=='QL' .OR. C4=='HR' .OR. C4=='TR' .OR. & C4=='BR' ) THEN NX = 128 END IF END IF END IF ILAENV = NX RETURN ! 400 CONTINUE ! ! ISPEC = 4: number of shifts (used by xHSEQR) ! ILAENV = 6 RETURN ! 500 CONTINUE ! ! ISPEC = 5: minimum column dimension (not used) ! ILAENV = 2 RETURN ! 600 CONTINUE ! ! ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) ! ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 ) RETURN ! 700 CONTINUE ! ! ISPEC = 7: number of processors (not used) ! ILAENV = 1 RETURN ! 800 CONTINUE ! ! ISPEC = 8: crossover point for multishift (used by xHSEQR) ! ILAENV = 50 RETURN ! ! End of ILAENV ! END LOGICAL FUNCTION LSAME( CA, CB ) ! ! -- LAPACK auxiliary routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! September 30, 1994 ! ! .. Scalar Arguments .. CHARACTER CA, CB ! .. ! ! Purpose ! ======= ! ! LSAME returns .TRUE. if CA is the same letter as CB regardless of ! case. ! ! Arguments ! ========= ! ! CA (input) CHARACTER*1 ! CB (input) CHARACTER*1 ! CA and CB specify the single characters to be compared. ! ! ===================================================================== ! ! .. Intrinsic Functions .. INTRINSIC ICHAR ! .. ! .. Local Scalars .. INTEGER INTA, INTB, ZCODE ! .. ! .. Executable Statements .. ! ! Test if the characters are equal ! LSAME = CA==CB IF( LSAME ) & RETURN ! ! Now test for equivalence if both characters are alphabetic. ! ZCODE = ICHAR( 'Z' ) ! ! Use 'Z' rather than 'A' so that ASCII can be detected on Prime ! machines, on which ICHAR returns a value with bit 8 set. ! ICHAR('A') on Prime machines returns 193 which is the same as ! ICHAR('A') on an EBCDIC machine. ! INTA = ICHAR( CA ) INTB = ICHAR( CB ) ! IF( ZCODE==90 .OR. ZCODE==122 ) THEN ! ! ASCII is assumed - ZCODE is the ASCII code of either lower or ! upper case 'Z'. ! IF( INTA>=97 .AND. INTA<=122 ) INTA = INTA - 32 IF( INTB>=97 .AND. INTB<=122 ) INTB = INTB - 32 ! ELSE IF( ZCODE==233 .OR. ZCODE==169 ) THEN ! ! EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or ! upper case 'Z'. ! IF( INTA>=129 .AND. INTA<=137 .OR. & INTA>=145 .AND. INTA<=153 .OR. & INTA>=162 .AND. INTA<=169 ) INTA = INTA + 64 IF( INTB>=129 .AND. INTB<=137 .OR. & INTB>=145 .AND. INTB<=153 .OR. & INTB>=162 .AND. INTB<=169 ) INTB = INTB + 64 ! ELSE IF( ZCODE==218 .OR. ZCODE==250 ) THEN ! ! ASCII is assumed, on Prime machines - ZCODE is the ASCII code ! plus 128 of either lower or upper case 'Z'. ! IF( INTA>=225 .AND. INTA<=250 ) INTA = INTA - 32 IF( INTB>=225 .AND. INTB<=250 ) INTB = INTB - 32 END IF LSAME = INTA==INTB ! ! RETURN ! ! End of LSAME ! END LOGICAL FUNCTION LSAMEN( N, CA, CB ) ! ! -- LAPACK auxiliary routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! September 30, 1994 ! ! .. Scalar Arguments .. CHARACTER*( * ) CA, CB INTEGER N ! .. ! ! Purpose ! ======= ! ! LSAMEN tests if the first N letters of CA are the same as the ! first N letters of CB, regardless of case. ! LSAMEN returns .TRUE. if CA and CB are equivalent except for case ! and .FALSE. otherwise. LSAMEN also returns .FALSE. if LEN( CA ) ! or LEN( CB ) is less than N. ! ! Arguments ! ========= ! ! N (input) INTEGER ! The number of characters in CA and CB to be compared. ! ! CA (input) CHARACTER*(*) ! CB (input) CHARACTER*(*) ! CA and CB specify two character strings of length at least N. ! Only the first N characters of each string will be accessed. ! ! ===================================================================== ! ! .. Local Scalars .. INTEGER I ! .. ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. Intrinsic Functions .. INTRINSIC LEN ! .. ! .. Executable Statements .. ! LSAMEN = .FALSE. IF( LEN( CA ) 0. ! ! ===================================================================== ! ! .. Local Scalars .. INTEGER I, NINCX REAL STEMP ! .. ! .. Intrinsic Functions .. INTRINSIC ABS ! .. ! .. Executable Statements .. ! SCSUM1 = 0.0E0 STEMP = 0.0E0 IF( N<=0 ) & RETURN IF( INCX==1 ) & GO TO 20 ! ! CODE FOR INCREMENT NOT EQUAL TO 1 ! NINCX = N*INCX DO 10 I = 1, NINCX, INCX ! ! NEXT LINE MODIFIED. ! STEMP = STEMP + ABS( CX( I ) ) 10 CONTINUE SCSUM1 = STEMP RETURN ! ! CODE FOR INCREMENT EQUAL TO 1 ! 20 CONTINUE DO 30 I = 1, N ! ! NEXT LINE MODIFIED. ! STEMP = STEMP + ABS( CX( I ) ) 30 CONTINUE SCSUM1 = STEMP RETURN ! ! End of SCSUM1 ! END SUBROUTINE SGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) ! ! -- LAPACK routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! February 29, 1992 ! ! .. Scalar Arguments .. INTEGER INFO, KL, KU, LDAB, M, N ! .. ! .. Array Arguments .. INTEGER IPIV( * ) REAL AB( LDAB, * ) ! .. ! ! Purpose ! ======= ! ! SGBTF2 computes an LU factorization of a real m-by-n band matrix A ! using partial pivoting with row interchanges. ! ! This is the unblocked version of the algorithm, calling Level 2 BLAS. ! ! Arguments ! ========= ! ! M (input) INTEGER ! The number of rows of the matrix A. M >= 0. ! ! N (input) INTEGER ! The number of columns of the matrix A. N >= 0. ! ! KL (input) INTEGER ! The number of subdiagonals within the band of A. KL >= 0. ! ! KU (input) INTEGER ! The number of superdiagonals within the band of A. KU >= 0. ! ! AB (input/output) REAL array, dimension (LDAB,N) ! On entry, the matrix A in band storage, in rows KL+1 to ! 2*KL+KU+1; rows 1 to KL of the array need not be set. ! The j-th column of A is stored in the j-th column of the ! array AB as follows: ! AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl) ! ! On exit, details of the factorization: U is stored as an ! upper triangular band matrix with KL+KU superdiagonals in ! rows 1 to KL+KU+1, and the multipliers used during the ! factorization are stored in rows KL+KU+2 to 2*KL+KU+1. ! See below for further details. ! ! LDAB (input) INTEGER ! The leading dimension of the array AB. LDAB >= 2*KL+KU+1. ! ! IPIV (output) INTEGER array, dimension (min(M,N)) ! The pivot indices; for 1 <= i <= min(M,N), row i of the ! matrix was interchanged with row IPIV(i). ! ! INFO (output) INTEGER ! = 0: successful exit ! < 0: if INFO = -i, the i-th argument had an illegal value ! > 0: if INFO = +i, U(i,i) is exactly zero. The factorization ! has been completed, but the factor U is exactly ! singular, and division by zero will occur if it is used ! to solve a system of equations. ! ! Further Details ! =============== ! ! The band storage scheme is illustrated by the following example, when ! M = N = 6, KL = 2, KU = 1: ! ! On entry: On exit: ! ! * * * + + + * * * u14 u25 u36 ! * * + + + + * * u13 u24 u35 u46 ! * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 ! a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 ! a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 * ! a31 a42 a53 a64 * * m31 m42 m53 m64 * * ! ! Array elements marked * are not used by the routine; elements marked ! + need not be set on entry, but are required by the routine to store ! elements of U, because of fill-in resulting from the row ! interchanges. ! ! ===================================================================== ! ! .. Parameters .. REAL ONE, ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) ! .. ! .. Local Scalars .. INTEGER I, J, JP, JU, KM, KV ! .. ! .. External Functions .. INTEGER ISAMAX EXTERNAL ISAMAX ! .. ! .. External Subroutines .. EXTERNAL SGER, SSCAL, SSWAP, XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC MAX, MIN ! .. ! .. Executable Statements .. ! ! KV is the number of superdiagonals in the factor U, allowing for ! fill-in. ! KV = KU + KL ! ! Test the input parameters. ! INFO = 0 IF( M<0 ) THEN INFO = -1 ELSE IF( N<0 ) THEN INFO = -2 ELSE IF( KL<0 ) THEN INFO = -3 ELSE IF( KU<0 ) THEN INFO = -4 ELSE IF( LDAB0 ) THEN ! ! Compute multipliers. ! CALL SSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 ) ! ! Update trailing submatrix within the band. ! IF( JU>J ) & CALL SGER( KM, JU-J, -ONE, AB( KV+2, J ), 1, & AB( KV, J+1 ), LDAB-1, AB( KV+1, J+1 ), & LDAB-1 ) END IF ELSE ! ! If pivot is zero, set INFO to the index of the pivot ! unless a zero pivot has already been found. ! IF( INFO==0 ) & INFO = J END IF 40 CONTINUE RETURN ! ! End of SGBTF2 ! END SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) ! ! -- LAPACK routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! February 29, 1992 ! ! .. Scalar Arguments .. INTEGER INFO, KL, KU, LDAB, M, N ! .. ! .. Array Arguments .. INTEGER IPIV( * ) REAL AB( LDAB, * ) ! .. ! ! Purpose ! ======= ! ! SGBTRF computes an LU factorization of a real m-by-n band matrix A ! using partial pivoting with row interchanges. ! ! This is the blocked version of the algorithm, calling Level 3 BLAS. ! ! Arguments ! ========= ! ! M (input) INTEGER ! The number of rows of the matrix A. M >= 0. ! ! N (input) INTEGER ! The number of columns of the matrix A. N >= 0. ! ! KL (input) INTEGER ! The number of subdiagonals within the band of A. KL >= 0. ! ! KU (input) INTEGER ! The number of superdiagonals within the band of A. KU >= 0. ! ! AB (input/output) REAL array, dimension (LDAB,N) ! On entry, the matrix A in band storage, in rows KL+1 to ! 2*KL+KU+1; rows 1 to KL of the array need not be set. ! The j-th column of A is stored in the j-th column of the ! array AB as follows: ! AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl) ! ! On exit, details of the factorization: U is stored as an ! upper triangular band matrix with KL+KU superdiagonals in ! rows 1 to KL+KU+1, and the multipliers used during the ! factorization are stored in rows KL+KU+2 to 2*KL+KU+1. ! See below for further details. ! ! LDAB (input) INTEGER ! The leading dimension of the array AB. LDAB >= 2*KL+KU+1. ! ! IPIV (output) INTEGER array, dimension (min(M,N)) ! The pivot indices; for 1 <= i <= min(M,N), row i of the ! matrix was interchanged with row IPIV(i). ! ! INFO (output) INTEGER ! = 0: successful exit ! < 0: if INFO = -i, the i-th argument had an illegal value ! > 0: if INFO = +i, U(i,i) is exactly zero. The factorization ! has been completed, but the factor U is exactly ! singular, and division by zero will occur if it is used ! to solve a system of equations. ! ! Further Details ! =============== ! ! The band storage scheme is illustrated by the following example, when ! M = N = 6, KL = 2, KU = 1: ! ! On entry: On exit: ! ! * * * + + + * * * u14 u25 u36 ! * * + + + + * * u13 u24 u35 u46 ! * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 ! a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 ! a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 * ! a31 a42 a53 a64 * * m31 m42 m53 m64 * * ! ! Array elements marked * are not used by the routine; elements marked ! + need not be set on entry, but are required by the routine to store ! elements of U because of fill-in resulting from the row interchanges. ! ! ===================================================================== ! ! .. Parameters .. REAL ONE, ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) INTEGER NBMAX, LDWORK PARAMETER ( NBMAX = 64, LDWORK = NBMAX+1 ) ! .. ! .. Local Scalars .. INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP, & JU, K2, KM, KV, NB, NW REAL TEMP ! .. ! .. Local Arrays .. REAL WORK13( LDWORK, NBMAX ), & WORK31( LDWORK, NBMAX ) ! .. ! .. External Functions .. INTEGER ILAENV, ISAMAX EXTERNAL ILAENV, ISAMAX ! .. ! .. External Subroutines .. EXTERNAL SCOPY, SGBTF2, SGEMM, SGER, SLASWP, SSCAL, & SSWAP, STRSM, XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC MAX, MIN ! .. ! .. Executable Statements .. ! ! KV is the number of superdiagonals in the factor U, allowing for ! fill-in ! KV = KU + KL ! ! Test the input parameters. ! INFO = 0 IF( M<0 ) THEN INFO = -1 ELSE IF( N<0 ) THEN INFO = -2 ELSE IF( KL<0 ) THEN INFO = -3 ELSE IF( KU<0 ) THEN INFO = -4 ELSE IF( LDABKL ) THEN ! ! Use unblocked code ! CALL SGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) ELSE ! ! Use blocked code ! ! Zero the superdiagonal elements of the work array WORK13 ! DO 20 J = 1, NB DO 10 I = 1, J - 1 WORK13( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ! ! Zero the subdiagonal elements of the work array WORK31 ! DO 40 J = 1, NB DO 30 I = J + 1, NB WORK31( I, J ) = ZERO 30 CONTINUE 40 CONTINUE ! ! Gaussian elimination with partial pivoting ! ! Set fill-in elements in columns KU+2 to KV to zero ! DO 60 J = KU + 2, MIN( KV, N ) DO 50 I = KV - J + 2, KL AB( I, J ) = ZERO 50 CONTINUE 60 CONTINUE ! ! JU is the index of the last column affected by the current ! stage of the factorization ! JU = 1 ! DO 180 J = 1, MIN( M, N ), NB JB = MIN( NB, MIN( M, N )-J+1 ) ! ! The active part of the matrix is partitioned ! ! A11 A12 A13 ! A21 A22 A23 ! A31 A32 A33 ! ! Here A11, A21 and A31 denote the current block of JB columns ! which is about to be factorized. The number of rows in the ! partitioning are JB, I2, I3 respectively, and the numbers ! of columns are JB, J2, J3. The superdiagonal elements of A13 ! and the subdiagonal elements of A31 lie outside the band. ! I2 = MIN( KL-JB, M-J-JB+1 ) I3 = MIN( JB, M-J-KL+1 ) ! ! J2 and J3 are computed after JU has been updated. ! ! Factorize the current block of JB columns ! DO 80 JJ = J, J + JB - 1 ! ! Set fill-in elements in column JJ+KV to zero ! IF( JJ+KV<=N ) THEN DO 70 I = 1, KL AB( I, JJ+KV ) = ZERO 70 CONTINUE END IF ! ! Find pivot and test for singularity. KM is the number of ! subdiagonal elements in the current column. ! KM = MIN( KL, M-JJ ) JP = ISAMAX( KM+1, AB( KV+1, JJ ), 1 ) IPIV( JJ ) = JP + JJ - J IF( AB( KV+JP, JJ )/=ZERO ) THEN JU = MAX( JU, MIN( JJ+KU+JP-1, N ) ) IF( JP/=1 ) THEN ! ! Apply interchange to columns J to J+JB-1 ! IF( JP+JJ-1JJ ) & CALL SGER( KM, JM-JJ, -ONE, AB( KV+2, JJ ), 1, & AB( KV, JJ+1 ), LDAB-1, & AB( KV+1, JJ+1 ), LDAB-1 ) ELSE ! ! If pivot is zero, set INFO to the index of the pivot ! unless a zero pivot has already been found. ! IF( INFO==0 ) & INFO = JJ END IF ! ! Copy current column of A31 into the work array WORK31 ! NW = MIN( JJ-J+1, I3 ) IF( NW>0 ) & CALL SCOPY( NW, AB( KV+KL+1-JJ+J, JJ ), 1, & WORK31( 1, JJ-J+1 ), 1 ) 80 CONTINUE IF( J+JB<=N ) THEN ! ! Apply the row interchanges to the other blocks. ! J2 = MIN( JU-J+1, KV ) - JB J3 = MAX( 0, JU-J-KV+1 ) ! ! Use SLASWP to apply the row interchanges to A12, A22, and ! A32. ! CALL SLASWP( J2, AB( KV+1-JB, J+JB ), LDAB-1, 1, JB, & IPIV( J ), 1 ) ! ! Adjust the pivot indices. ! DO 90 I = J, J + JB - 1 IPIV( I ) = IPIV( I ) + J - 1 90 CONTINUE ! ! Apply the row interchanges to A13, A23, and A33 ! columnwise. ! K2 = J - 1 + JB + J2 DO 110 I = 1, J3 JJ = K2 + I DO 100 II = J + I - 1, J + JB - 1 IP = IPIV( II ) IF( IP/=II ) THEN TEMP = AB( KV+1+II-JJ, JJ ) AB( KV+1+II-JJ, JJ ) = AB( KV+1+IP-JJ, JJ ) AB( KV+1+IP-JJ, JJ ) = TEMP END IF 100 CONTINUE 110 CONTINUE ! ! Update the relevant part of the trailing submatrix ! IF( J2>0 ) THEN ! ! Update A12 ! CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', & JB, J2, ONE, AB( KV+1, J ), LDAB-1, & AB( KV+1-JB, J+JB ), LDAB-1 ) ! IF( I2>0 ) THEN ! ! Update A22 ! CALL SGEMM( 'No transpose', 'No transpose', I2, J2, & JB, -ONE, AB( KV+1+JB, J ), LDAB-1, & AB( KV+1-JB, J+JB ), LDAB-1, ONE, & AB( KV+1, J+JB ), LDAB-1 ) END IF ! IF( I3>0 ) THEN ! ! Update A32 ! CALL SGEMM( 'No transpose', 'No transpose', I3, J2, & JB, -ONE, WORK31, LDWORK, & AB( KV+1-JB, J+JB ), LDAB-1, ONE, & AB( KV+KL+1-JB, J+JB ), LDAB-1 ) END IF END IF ! IF( J3>0 ) THEN ! ! Copy the lower triangle of A13 into the work array ! WORK13 ! DO 130 JJ = 1, J3 DO 120 II = JJ, JB WORK13( II, JJ ) = AB( II-JJ+1, JJ+J+KV-1 ) 120 CONTINUE 130 CONTINUE ! ! Update A13 in the work array ! CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', & JB, J3, ONE, AB( KV+1, J ), LDAB-1, & WORK13, LDWORK ) ! IF( I2>0 ) THEN ! ! Update A23 ! CALL SGEMM( 'No transpose', 'No transpose', I2, J3, & JB, -ONE, AB( KV+1+JB, J ), LDAB-1, & WORK13, LDWORK, ONE, AB( 1+JB, J+KV ), & LDAB-1 ) END IF ! IF( I3>0 ) THEN ! ! Update A33 ! CALL SGEMM( 'No transpose', 'No transpose', I3, J3, & JB, -ONE, WORK31, LDWORK, WORK13, & LDWORK, ONE, AB( 1+KL, J+KV ), LDAB-1 ) END IF ! ! Copy the lower triangle of A13 back into place ! DO 150 JJ = 1, J3 DO 140 II = JJ, JB AB( II-JJ+1, JJ+J+KV-1 ) = WORK13( II, JJ ) 140 CONTINUE 150 CONTINUE END IF ELSE ! ! Adjust the pivot indices. ! DO 160 I = J, J + JB - 1 IPIV( I ) = IPIV( I ) + J - 1 160 CONTINUE END IF ! ! Partially undo the interchanges in the current block to ! restore the upper triangular form of A31 and copy the upper ! triangle of A31 back into place ! DO 170 JJ = J + JB - 1, J, -1 JP = IPIV( JJ ) - JJ + 1 IF( JP/=1 ) THEN ! ! Apply interchange to columns J to JJ-1 ! IF( JP+JJ-10 ) & CALL SCOPY( NW, WORK31( 1, JJ-J+1 ), 1, & AB( KV+KL+1-JJ+J, JJ ), 1 ) 170 CONTINUE 180 CONTINUE END IF ! RETURN ! ! End of SGBTRF ! END SUBROUTINE SGBTRS( TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, & INFO ) ! ! -- LAPACK routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! March 31, 1993 ! ! .. Scalar Arguments .. CHARACTER TRANS INTEGER INFO, KL, KU, LDAB, LDB, N, NRHS ! .. ! .. Array Arguments .. INTEGER IPIV( * ) REAL AB( LDAB, * ), B( LDB, * ) ! .. ! ! Purpose ! ======= ! ! SGBTRS solves a system of linear equations ! A * X = B or A' * X = B ! with a general band matrix A using the LU factorization computed ! by SGBTRF. ! ! Arguments ! ========= ! ! TRANS (input) CHARACTER*1 ! Specifies the form of the system of equations. ! = 'N': A * X = B (No transpose) ! = 'T': A'* X = B (Transpose) ! = 'C': A'* X = B (Conjugate transpose = Transpose) ! ! N (input) INTEGER ! The order of the matrix A. N >= 0. ! ! KL (input) INTEGER ! The number of subdiagonals within the band of A. KL >= 0. ! ! KU (input) INTEGER ! The number of superdiagonals within the band of A. KU >= 0. ! ! NRHS (input) INTEGER ! The number of right hand sides, i.e., the number of columns ! of the matrix B. NRHS >= 0. ! ! AB (input) REAL array, dimension (LDAB,N) ! Details of the LU factorization of the band matrix A, as ! computed by SGBTRF. U is stored as an upper triangular band ! matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and ! the multipliers used during the factorization are stored in ! rows KL+KU+2 to 2*KL+KU+1. ! ! LDAB (input) INTEGER ! The leading dimension of the array AB. LDAB >= 2*KL+KU+1. ! ! IPIV (input) INTEGER array, dimension (N) ! The pivot indices; for 1 <= i <= N, row i of the matrix was ! interchanged with row IPIV(i). ! ! B (input/output) REAL array, dimension (LDB,NRHS) ! On entry, the right hand side matrix B. ! On exit, the solution matrix X. ! ! LDB (input) INTEGER ! The leading dimension of the array B. LDB >= max(1,N). ! ! INFO (output) INTEGER ! = 0: successful exit ! < 0: if INFO = -i, the i-th argument had an illegal value ! ! ===================================================================== ! ! .. Parameters .. REAL ONE PARAMETER ( ONE = 1.0E+0 ) ! .. ! .. Local Scalars .. LOGICAL LNOTI, NOTRAN INTEGER I, J, KD, L, LM ! .. ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. External Subroutines .. EXTERNAL SGEMV, SGER, SSWAP, STBSV, XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC MAX, MIN ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 NOTRAN = LSAME( TRANS, 'N' ) IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. & LSAME( TRANS, 'C' ) ) THEN INFO = -1 ELSE IF( N<0 ) THEN INFO = -2 ELSE IF( KL<0 ) THEN INFO = -3 ELSE IF( KU<0 ) THEN INFO = -4 ELSE IF( NRHS<0 ) THEN INFO = -5 ELSE IF( LDAB<( 2*KL+KU+1 ) ) THEN INFO = -7 ELSE IF( LDB0 ! IF( NOTRAN ) THEN ! ! Solve A*X = B. ! ! Solve L*X = B, overwriting B with X. ! ! L is represented as a product of permutations and unit lower ! triangular matrices L = P(1) * L(1) * ... * P(n-1) * L(n-1), ! where each transformation L(i) is a rank-one modification of ! the identity matrix. ! IF( LNOTI ) THEN DO 10 J = 1, N - 1 LM = MIN( KL, N-J ) L = IPIV( J ) IF( L/=J ) & CALL SSWAP( NRHS, B( L, 1 ), LDB, B( J, 1 ), LDB ) CALL SGER( LM, NRHS, -ONE, AB( KD+1, J ), 1, B( J, 1 ), & LDB, B( J+1, 1 ), LDB ) 10 CONTINUE END IF ! DO 20 I = 1, NRHS ! ! Solve U*X = B, overwriting B with X. ! CALL STBSV( 'Upper', 'No transpose', 'Non-unit', N, KL+KU, & AB, LDAB, B( 1, I ), 1 ) 20 CONTINUE ! ELSE ! ! Solve A'*X = B. ! DO 30 I = 1, NRHS ! ! Solve U'*X = B, overwriting B with X. ! CALL STBSV( 'Upper', 'Transpose', 'Non-unit', N, KL+KU, AB, & LDAB, B( 1, I ), 1 ) 30 CONTINUE ! ! Solve L'*X = B, overwriting B with X. ! IF( LNOTI ) THEN DO 40 J = N - 1, 1, -1 LM = MIN( KL, N-J ) CALL SGEMV( 'Transpose', LM, NRHS, -ONE, B( J+1, 1 ), & LDB, AB( KD+1, J ), 1, ONE, B( J, 1 ), LDB ) L = IPIV( J ) IF( L/=J ) & CALL SSWAP( NRHS, B( L, 1 ), LDB, B( J, 1 ), LDB ) 40 CONTINUE END IF END IF RETURN ! ! End of SGBTRS ! END SUBROUTINE SGEQR2( M, N, A, LDA, TAU, WORK, INFO ) ! ! -- LAPACK routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! February 29, 1992 ! ! .. Scalar Arguments .. INTEGER INFO, LDA, M, N ! .. ! .. Array Arguments .. REAL A( LDA, * ), TAU( * ), WORK( * ) ! .. ! ! Purpose ! ======= ! ! SGEQR2 computes a QR factorization of a real m by n matrix A: ! A = Q * R. ! ! Arguments ! ========= ! ! M (input) INTEGER ! The number of rows of the matrix A. M >= 0. ! ! N (input) INTEGER ! The number of columns of the matrix A. N >= 0. ! ! A (input/output) REAL array, dimension (LDA,N) ! On entry, the m by n matrix A. ! On exit, the elements on and above the diagonal of the array ! contain the min(m,n) by n upper trapezoidal matrix R (R is ! upper triangular if m >= n); the elements below the diagonal, ! with the array TAU, represent the orthogonal matrix Q as a ! product of elementary reflectors (see Further Details). ! ! LDA (input) INTEGER ! The leading dimension of the array A. LDA >= max(1,M). ! ! TAU (output) REAL array, dimension (min(M,N)) ! The scalar factors of the elementary reflectors (see Further ! Details). ! ! WORK (workspace) REAL array, dimension (N) ! ! INFO (output) INTEGER ! = 0: successful exit ! < 0: if INFO = -i, the i-th argument had an illegal value ! ! Further Details ! =============== ! ! The matrix Q is represented as a product of elementary reflectors ! ! Q = H(1) H(2) . . . H(k), where k = min(m,n). ! ! Each H(i) has the form ! ! H(i) = I - tau * v * v' ! ! where tau is a real scalar, and v is a real vector with ! v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), ! and tau in TAU(i). ! ! ===================================================================== ! ! .. Parameters .. REAL ONE PARAMETER ( ONE = 1.0E+0 ) ! .. ! .. Local Scalars .. INTEGER I, K REAL AII ! .. ! .. External Subroutines .. EXTERNAL SLARF, SLARFG, XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC MAX, MIN ! .. ! .. Executable Statements .. ! ! Test the input arguments ! INFO = 0 IF( M<0 ) THEN INFO = -1 ELSE IF( N<0 ) THEN INFO = -2 ELSE IF( LDA= 0. ! ! DL (input/output) REAL array, dimension (N-1) ! On entry, DL must contain the (n-1) subdiagonal elements of ! A. ! On exit, DL is overwritten by the (n-1) multipliers that ! define the matrix L from the LU factorization of A. ! ! D (input/output) REAL array, dimension (N) ! On entry, D must contain the diagonal elements of A. ! On exit, D is overwritten by the n diagonal elements of the ! upper triangular matrix U from the LU factorization of A. ! ! DU (input/output) REAL array, dimension (N-1) ! On entry, DU must contain the (n-1) superdiagonal elements ! of A. ! On exit, DU is overwritten by the (n-1) elements of the first ! superdiagonal of U. ! ! DU2 (output) REAL array, dimension (N-2) ! On exit, DU2 is overwritten by the (n-2) elements of the ! second superdiagonal of U. ! ! IPIV (output) INTEGER array, dimension (N) ! The pivot indices; for 1 <= i <= n, row i of the matrix was ! interchanged with row IPIV(i). IPIV(i) will always be either ! i or i+1; IPIV(i) = i indicates a row interchange was not ! required. ! ! INFO (output) INTEGER ! = 0: successful exit ! < 0: if INFO = -i, the i-th argument had an illegal value ! > 0: if INFO = i, U(i,i) is exactly zero. The factorization ! has been completed, but the factor U is exactly ! singular, and division by zero will occur if it is used ! to solve a system of equations. ! ! ===================================================================== ! ! .. Local Scalars .. INTEGER I REAL FACT, TEMP ! .. ! .. Intrinsic Functions .. INTRINSIC ABS ! .. ! .. External Subroutines .. EXTERNAL XERBLA ! .. ! .. Parameters .. REAL ZERO PARAMETER ( ZERO = 0.0E+0 ) ! .. ! .. Executable Statements .. ! INFO = 0 IF( N<0 ) THEN INFO = -1 CALL XERBLA( 'SGTTRF', -INFO ) RETURN END IF ! ! Quick return if possible ! IF( N==0 ) & RETURN ! ! Initialize IPIV(i) = i ! DO 10 I = 1, N IPIV( I ) = I 10 CONTINUE ! DO 20 I = 1, N - 1 IF( DL( I )==ZERO ) THEN ! ! Subdiagonal is zero, no elimination is required. ! IF( D( I )==ZERO .AND. INFO==0 ) & INFO = I IF( I=ABS( DL( I ) ) ) THEN ! ! No row interchange required, eliminate DL(I) ! FACT = DL( I ) / D( I ) DL( I ) = FACT D( I+1 ) = D( I+1 ) - FACT*DU( I ) IF( I= 0. ! ! NRHS (input) INTEGER ! The number of right hand sides, i.e., the number of columns ! of the matrix B. NRHS >= 0. ! ! DL (input) REAL array, dimension (N-1) ! The (n-1) multipliers that define the matrix L from the ! LU factorization of A. ! ! D (input) REAL array, dimension (N) ! The n diagonal elements of the upper triangular matrix U from ! the LU factorization of A. ! ! DU (input) REAL array, dimension (N-1) ! The (n-1) elements of the first superdiagonal of U. ! ! DU2 (input) REAL array, dimension (N-2) ! The (n-2) elements of the second superdiagonal of U. ! ! IPIV (input) INTEGER array, dimension (N) ! The pivot indices; for 1 <= i <= n, row i of the matrix was ! interchanged with row IPIV(i). IPIV(i) will always be either ! i or i+1; IPIV(i) = i indicates a row interchange was not ! required. ! ! B (input/output) REAL array, dimension (LDB,NRHS) ! On entry, the right hand side matrix B. ! On exit, B is overwritten by the solution matrix X. ! ! LDB (input) INTEGER ! The leading dimension of the array B. LDB >= max(1,N). ! ! INFO (output) INTEGER ! = 0: successful exit ! < 0: if INFO = -i, the i-th argument had an illegal value ! ! ===================================================================== ! ! .. Local Scalars .. LOGICAL NOTRAN INTEGER I, J REAL TEMP ! .. ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. External Subroutines .. EXTERNAL XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Executable Statements .. ! INFO = 0 NOTRAN = LSAME( TRANS, 'N' ) IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. & LSAME( TRANS, 'C' ) ) THEN INFO = -1 ELSE IF( N<0 ) THEN INFO = -2 ELSE IF( NRHS<0 ) THEN INFO = -3 ELSE IF( LDB1 ) & B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / & D( N-1 ) DO 20 I = N - 2, 1, -1 B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )* & B( I+2, J ) ) / D( I ) 20 CONTINUE 30 CONTINUE ELSE ! ! Solve A' * X = B. ! DO 60 J = 1, NRHS ! ! Solve U'*x = b. ! B( 1, J ) = B( 1, J ) / D( 1 ) IF( N>1 ) & B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 ) DO 40 I = 3, N B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-DU2( I-2 )* & B( I-2, J ) ) / D( I ) 40 CONTINUE ! ! Solve L'*x = b. ! DO 50 I = N - 1, 1, -1 IF( IPIV( I )==I ) THEN B( I, J ) = B( I, J ) - DL( I )*B( I+1, J ) ELSE TEMP = B( I+1, J ) B( I+1, J ) = B( I, J ) - DL( I )*TEMP B( I, J ) = TEMP END IF 50 CONTINUE 60 CONTINUE END IF ! ! End of SGTTRS ! END SUBROUTINE SLABAD( SMALL, LARGE ) ! ! -- LAPACK auxiliary routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! October 31, 1992 ! ! .. Scalar Arguments .. REAL LARGE, SMALL ! .. ! ! Purpose ! ======= ! ! SLABAD takes as input the values computed by SLAMCH for underflow and ! overflow, and returns the square root of each of these values if the ! log of LARGE is sufficiently large. This subroutine is intended to ! identify machines with a large exponent range, such as the Crays, and ! redefine the underflow and overflow limits to be the square roots of ! the values computed by SLAMCH. This subroutine is needed because ! SLAMCH does not compensate for poor arithmetic in the upper half of ! the exponent range, as is found on a Cray. ! ! Arguments ! ========= ! ! SMALL (input/output) REAL ! On entry, the underflow threshold as computed by SLAMCH. ! On exit, if LOG10(LARGE) is sufficiently large, the square ! root of SMALL, otherwise unchanged. ! ! LARGE (input/output) REAL ! On entry, the overflow threshold as computed by SLAMCH. ! On exit, if LOG10(LARGE) is sufficiently large, the square ! root of LARGE, otherwise unchanged. ! ! ===================================================================== ! ! .. Intrinsic Functions .. INTRINSIC LOG10, SQRT ! .. ! .. Executable Statements .. ! ! If it looks like we're on a Cray, take the square root of ! SMALL and LARGE to avoid overflow and underflow problems. ! IF( LOG10( LARGE )>2000. ) THEN SMALL = SQRT( SMALL ) LARGE = SQRT( LARGE ) END IF ! RETURN ! ! End of SLABAD ! END SUBROUTINE SLACON( N, V, X, ISGN, EST, KASE ) ! ! -- LAPACK auxiliary routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! February 29, 1992 ! ! .. Scalar Arguments .. INTEGER KASE, N REAL EST ! .. ! .. Array Arguments .. INTEGER ISGN( * ) REAL V( * ), X( * ) ! .. ! ! Purpose ! ======= ! ! SLACON estimates the 1-norm of a square, real matrix A. ! Reverse communication is used for evaluating matrix-vector products. ! ! Arguments ! ========= ! ! N (input) INTEGER ! The order of the matrix. N >= 1. ! ! V (workspace) REAL array, dimension (N) ! On the final return, V = A*W, where EST = norm(V)/norm(W) ! (W is not returned). ! ! X (input/output) REAL array, dimension (N) ! On an intermediate return, X should be overwritten by ! A * X, if KASE=1, ! A' * X, if KASE=2, ! and SLACON must be re-called with all the other parameters ! unchanged. ! ! ISGN (workspace) INTEGER array, dimension (N) ! ! EST (output) REAL ! An estimate (a lower bound) for norm(A). ! ! KASE (input/output) INTEGER ! On the initial call to SLACON, KASE should be 0. ! On an intermediate return, KASE will be 1 or 2, indicating ! whether X should be overwritten by A * X or A' * X. ! On the final return from SLACON, KASE will again be 0. ! ! Further Details ! ======= ======= ! ! Contributed by Nick Higham, University of Manchester. ! Originally named SONEST, dated March 16, 1988. ! ! Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of ! a real or complex matrix, with applications to condition estimation", ! ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. ! ! ===================================================================== ! ! .. Parameters .. INTEGER ITMAX PARAMETER ( ITMAX = 5 ) REAL ZERO, ONE, TWO PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 ) ! .. ! .. Local Scalars .. INTEGER I, ITER, J, JLAST, JUMP REAL ALTSGN, ESTOLD, TEMP ! .. ! .. External Functions .. INTEGER ISAMAX REAL SASUM EXTERNAL ISAMAX, SASUM ! .. ! .. External Subroutines .. EXTERNAL SCOPY ! .. ! .. Intrinsic Functions .. INTRINSIC ABS, NINT, REAL, SIGN ! .. ! .. Save statement .. SAVE ! .. ! .. Executable Statements .. ! IF( KASE==0 ) THEN DO 10 I = 1, N X( I ) = ONE / REAL( N ) 10 CONTINUE KASE = 1 JUMP = 1 RETURN END IF ! GO TO ( 20, 40, 70, 110, 140 )JUMP ! ! ................ ENTRY (JUMP = 1) ! FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. ! 20 CONTINUE IF( N==1 ) THEN V( 1 ) = X( 1 ) EST = ABS( V( 1 ) ) ! ... QUIT GO TO 150 END IF EST = SASUM( N, X, 1 ) ! DO 30 I = 1, N X( I ) = SIGN( ONE, X( I ) ) ISGN( I ) = NINT( X( I ) ) 30 CONTINUE KASE = 2 JUMP = 2 RETURN ! ! ................ ENTRY (JUMP = 2) ! FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. ! 40 CONTINUE J = ISAMAX( N, X, 1 ) ITER = 2 ! ! MAIN LOOP - ITERATIONS 2,3,...,ITMAX. ! 50 CONTINUE DO 60 I = 1, N X( I ) = ZERO 60 CONTINUE X( J ) = ONE KASE = 1 JUMP = 3 RETURN ! ! ................ ENTRY (JUMP = 3) ! X HAS BEEN OVERWRITTEN BY A*X. ! 70 CONTINUE CALL SCOPY( N, X, 1, V, 1 ) ESTOLD = EST EST = SASUM( N, V, 1 ) DO 80 I = 1, N IF( NINT( SIGN( ONE, X( I ) ) )/=ISGN( I ) ) & GO TO 90 80 CONTINUE ! REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. GO TO 120 ! 90 CONTINUE ! TEST FOR CYCLING. IF( EST<=ESTOLD ) & GO TO 120 ! DO 100 I = 1, N X( I ) = SIGN( ONE, X( I ) ) ISGN( I ) = NINT( X( I ) ) 100 CONTINUE KASE = 2 JUMP = 4 RETURN ! ! ................ ENTRY (JUMP = 4) ! X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. ! 110 CONTINUE JLAST = J J = ISAMAX( N, X, 1 ) IF( ( X( JLAST )/=ABS( X( J ) ) ) .AND. ( ITEREST ) THEN CALL SCOPY( N, X, 1, V, 1 ) EST = TEMP END IF ! 150 CONTINUE KASE = 0 RETURN ! ! End of SLACON ! END SUBROUTINE SLACPY( UPLO, M, N, A, LDA, B, LDB ) ! ! -- LAPACK auxiliary routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! February 29, 1992 ! ! .. Scalar Arguments .. CHARACTER UPLO INTEGER LDA, LDB, M, N ! .. ! .. Array Arguments .. REAL A( LDA, * ), B( LDB, * ) ! .. ! ! Purpose ! ======= ! ! SLACPY copies all or part of a two-dimensional matrix A to another ! matrix B. ! ! Arguments ! ========= ! ! UPLO (input) CHARACTER*1 ! Specifies the part of the matrix A to be copied to B. ! = 'U': Upper triangular part ! = 'L': Lower triangular part ! Otherwise: All of the matrix A ! ! M (input) INTEGER ! The number of rows of the matrix A. M >= 0. ! ! N (input) INTEGER ! The number of columns of the matrix A. N >= 0. ! ! A (input) REAL array, dimension (LDA,N) ! The m by n matrix A. If UPLO = 'U', only the upper triangle ! or trapezoid is accessed; if UPLO = 'L', only the lower ! triangle or trapezoid is accessed. ! ! LDA (input) INTEGER ! The leading dimension of the array A. LDA >= max(1,M). ! ! B (output) REAL array, dimension (LDB,N) ! On exit, B = A in the locations specified by UPLO. ! ! LDB (input) INTEGER ! The leading dimension of the array B. LDB >= max(1,M). ! ! ===================================================================== ! ! .. Local Scalars .. INTEGER I, J ! .. ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. Intrinsic Functions .. INTRINSIC MIN ! .. ! .. Executable Statements .. ! IF( LSAME( UPLO, 'U' ) ) THEN DO 20 J = 1, N DO 10 I = 1, MIN( J, M ) B( I, J ) = A( I, J ) 10 CONTINUE 20 CONTINUE ELSE IF( LSAME( UPLO, 'L' ) ) THEN DO 40 J = 1, N DO 30 I = J, M B( I, J ) = A( I, J ) 30 CONTINUE 40 CONTINUE ELSE DO 60 J = 1, N DO 50 I = 1, M B( I, J ) = A( I, J ) 50 CONTINUE 60 CONTINUE END IF RETURN ! ! End of SLACPY ! END SUBROUTINE SLADIV( A, B, C, D, P, Q ) ! ! -- LAPACK auxiliary routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! October 31, 1992 ! ! .. Scalar Arguments .. REAL A, B, C, D, P, Q ! .. ! ! Purpose ! ======= ! ! SLADIV performs complex division in real arithmetic ! ! a + i*b ! p + i*q = --------- ! c + i*d ! ! The algorithm is due to Robert L. Smith and can be found ! in D. Knuth, The art of Computer Programming, Vol.2, p.195 ! ! Arguments ! ========= ! ! A (input) REAL ! B (input) REAL ! C (input) REAL ! D (input) REAL ! The scalars a, b, c, and d in the above expression. ! ! P (output) REAL ! Q (output) REAL ! The scalars p and q in the above expression. ! ! ===================================================================== ! ! .. Local Scalars .. REAL E, F ! .. ! .. Intrinsic Functions .. INTRINSIC ABS ! .. ! .. Executable Statements .. ! IF( ABS( D )ABS( C ) ) THEN ACMX = A ACMN = C ELSE ACMX = C ACMN = A END IF IF( ADF>AB ) THEN RT = ADF*SQRT( ONE+( AB / ADF )**2 ) ELSE IF( ADFZERO ) THEN RT1 = HALF*( SM+RT ) ! ! Order of execution important. ! To get fully accurate smaller eigenvalue, ! next line needs to be executed in higher precision. ! RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B ELSE ! ! Includes case RT1 = RT2 = 0 ! RT1 = HALF*RT RT2 = -HALF*RT END IF RETURN ! ! End of SLAE2 ! END SUBROUTINE SLAEV2( A, B, C, RT1, RT2, CS1, SN1 ) ! ! -- LAPACK auxiliary routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! October 31, 1992 ! ! .. Scalar Arguments .. REAL A, B, C, CS1, RT1, RT2, SN1 ! .. ! ! Purpose ! ======= ! ! SLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix ! [ A B ] ! [ B C ]. ! On return, RT1 is the eigenvalue of larger absolute value, RT2 is the ! eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right ! eigenvector for RT1, giving the decomposition ! ! [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ] ! [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ]. ! ! Arguments ! ========= ! ! A (input) REAL ! The (1,1) element of the 2-by-2 matrix. ! ! B (input) REAL ! The (1,2) element and the conjugate of the (2,1) element of ! the 2-by-2 matrix. ! ! C (input) REAL ! The (2,2) element of the 2-by-2 matrix. ! ! RT1 (output) REAL ! The eigenvalue of larger absolute value. ! ! RT2 (output) REAL ! The eigenvalue of smaller absolute value. ! ! CS1 (output) REAL ! SN1 (output) REAL ! The vector (CS1, SN1) is a unit right eigenvector for RT1. ! ! Further Details ! =============== ! ! RT1 is accurate to a few ulps barring over/underflow. ! ! RT2 may be inaccurate if there is massive cancellation in the ! determinant A*C-B*B; higher precision or correctly rounded or ! correctly truncated arithmetic would be needed to compute RT2 ! accurately in all cases. ! ! CS1 and SN1 are accurate to a few ulps barring over/underflow. ! ! Overflow is possible only if RT1 is within a factor of 5 of overflow. ! Underflow is harmless if the input data is 0 or exceeds ! underflow_threshold / macheps. ! ! ===================================================================== ! ! .. Parameters .. REAL ONE PARAMETER ( ONE = 1.0E0 ) REAL TWO PARAMETER ( TWO = 2.0E0 ) REAL ZERO PARAMETER ( ZERO = 0.0E0 ) REAL HALF PARAMETER ( HALF = 0.5E0 ) ! .. ! .. Local Scalars .. INTEGER SGN1, SGN2 REAL AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM, & TB, TN ! .. ! .. Intrinsic Functions .. INTRINSIC ABS, SQRT ! .. ! .. Executable Statements .. ! ! Compute the eigenvalues ! SM = A + C DF = A - C ADF = ABS( DF ) TB = B + B AB = ABS( TB ) IF( ABS( A )>ABS( C ) ) THEN ACMX = A ACMN = C ELSE ACMX = C ACMN = A END IF IF( ADF>AB ) THEN RT = ADF*SQRT( ONE+( AB / ADF )**2 ) ELSE IF( ADFZERO ) THEN RT1 = HALF*( SM+RT ) SGN1 = 1 ! ! Order of execution important. ! To get fully accurate smaller eigenvalue, ! next line needs to be executed in higher precision. ! RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B ELSE ! ! Includes case RT1 = RT2 = 0 ! RT1 = HALF*RT RT2 = -HALF*RT SGN1 = 1 END IF ! ! Compute the eigenvector ! IF( DF>=ZERO ) THEN CS = DF + RT SGN2 = 1 ELSE CS = DF - RT SGN2 = -1 END IF ACS = ABS( CS ) IF( ACS>AB ) THEN CT = -TB / CS SN1 = ONE / SQRT( ONE+CT*CT ) CS1 = CT*SN1 ELSE IF( AB==ZERO ) THEN CS1 = ONE SN1 = ZERO ELSE TN = -CS / TB CS1 = ONE / SQRT( ONE+TN*TN ) SN1 = TN*CS1 END IF END IF IF( SGN1==SGN2 ) THEN TN = CS1 CS1 = -SN1 SN1 = TN END IF RETURN ! ! End of SLAEV2 ! END SUBROUTINE SLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK, & INFO ) ! ! -- LAPACK auxiliary routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! February 29, 1992 ! ! .. Scalar Arguments .. LOGICAL WANTQ INTEGER INFO, J1, LDQ, LDT, N, N1, N2 ! .. ! .. Array Arguments .. REAL Q( LDQ, * ), T( LDT, * ), WORK( * ) ! .. ! ! Purpose ! ======= ! ! SLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in ! an upper quasi-triangular matrix T by an orthogonal similarity ! transformation. ! ! T must be in Schur canonical form, that is, block upper triangular ! with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block ! has its diagonal elemnts equal and its off-diagonal elements of ! opposite sign. ! ! Arguments ! ========= ! ! WANTQ (input) LOGICAL ! = .TRUE. : accumulate the transformation in the matrix Q; ! = .FALSE.: do not accumulate the transformation. ! ! N (input) INTEGER ! The order of the matrix T. N >= 0. ! ! T (input/output) REAL array, dimension (LDT,N) ! On entry, the upper quasi-triangular matrix T, in Schur ! canonical form. ! On exit, the updated matrix T, again in Schur canonical form. ! ! LDT (input) INTEGER ! The leading dimension of the array T. LDT >= max(1,N). ! ! Q (input/output) REAL array, dimension (LDQ,N) ! On entry, if WANTQ is .TRUE., the orthogonal matrix Q. ! On exit, if WANTQ is .TRUE., the updated matrix Q. ! If WANTQ is .FALSE., Q is not referenced. ! ! LDQ (input) INTEGER ! The leading dimension of the array Q. ! LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N. ! ! J1 (input) INTEGER ! The index of the first row of the first block T11. ! ! N1 (input) INTEGER ! The order of the first block T11. N1 = 0, 1 or 2. ! ! N2 (input) INTEGER ! The order of the second block T22. N2 = 0, 1 or 2. ! ! WORK (workspace) REAL array, dimension (N) ! ! INFO (output) INTEGER ! = 0: successful exit ! = 1: the transformed matrix T would be too far from Schur ! form; the blocks are not swapped and T and Q are ! unchanged. ! ! ===================================================================== ! ! .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) REAL TEN PARAMETER ( TEN = 1.0E+1 ) INTEGER LDD, LDX PARAMETER ( LDD = 4, LDX = 2 ) ! .. ! .. Local Scalars .. INTEGER IERR, J2, J3, J4, K, ND REAL CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22, & T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2, & WR1, WR2, XNORM ! .. ! .. Local Arrays .. REAL D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ), & X( LDX, 2 ) ! .. ! .. External Functions .. REAL SLAMCH, SLANGE EXTERNAL SLAMCH, SLANGE ! .. ! .. External Subroutines .. EXTERNAL SLACPY, SLANV2, SLARFG, SLARFX, SLARTG, SLASY2, & SROT ! .. ! .. Intrinsic Functions .. INTRINSIC ABS, MAX ! .. ! .. Executable Statements .. ! INFO = 0 ! ! Quick return if possible ! IF( N==0 .OR. N1==0 .OR. N2==0 ) & RETURN IF( J1+N1>N ) & RETURN ! J2 = J1 + 1 J3 = J1 + 2 J4 = J1 + 3 ! IF( N1==1 .AND. N2==1 ) THEN ! ! Swap two 1-by-1 blocks. ! T11 = T( J1, J1 ) T22 = T( J2, J2 ) ! ! Determine the transformation to perform the interchange. ! CALL SLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP ) ! ! Apply transformation to the matrix T. ! IF( J3<=N ) & CALL SROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS, & SN ) CALL SROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN ) ! T( J1, J1 ) = T22 T( J2, J2 ) = T11 ! IF( WANTQ ) THEN ! ! Accumulate transformation in the matrix Q. ! CALL SROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN ) END IF ! ELSE ! ! Swapping involves at least one 2-by-2 block. ! ! Copy the diagonal block of order N1+N2 to the local array D ! and compute its norm. ! ND = N1 + N2 CALL SLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD ) DNORM = SLANGE( 'Max', ND, ND, D, LDD, WORK ) ! ! Compute machine-dependent threshold for test for accepting ! swap. ! EPS = SLAMCH( 'P' ) SMLNUM = SLAMCH( 'S' ) / EPS THRESH = MAX( TEN*EPS*DNORM, SMLNUM ) ! ! Solve T11*X - X*T22 = scale*T12 for X. ! CALL SLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD, & D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X, & LDX, XNORM, IERR ) ! ! Swap the adjacent diagonal blocks. ! K = N1 + N1 + N2 - 3 GO TO ( 10, 20, 30 )K ! 10 CONTINUE ! ! N1 = 1, N2 = 2: generate elementary reflector H so that: ! ! ( scale, X11, X12 ) H = ( 0, 0, * ) ! U( 1 ) = SCALE U( 2 ) = X( 1, 1 ) U( 3 ) = X( 1, 2 ) CALL SLARFG( 3, U( 3 ), U, 1, TAU ) U( 3 ) = ONE T11 = T( J1, J1 ) ! ! Perform swap provisionally on diagonal block in D. ! CALL SLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK ) CALL SLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK ) ! ! Test whether to reject swap. ! IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3, & 3 )-T11 ) )>THRESH )GO TO 50 ! ! Accept swap: apply transformation to the entire matrix T. ! CALL SLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK ) CALL SLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK ) ! T( J3, J1 ) = ZERO T( J3, J2 ) = ZERO T( J3, J3 ) = T11 ! IF( WANTQ ) THEN ! ! Accumulate transformation in the matrix Q. ! CALL SLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK ) END IF GO TO 40 ! 20 CONTINUE ! ! N1 = 2, N2 = 1: generate elementary reflector H so that: ! ! H ( -X11 ) = ( * ) ! ( -X21 ) = ( 0 ) ! ( scale ) = ( 0 ) ! U( 1 ) = -X( 1, 1 ) U( 2 ) = -X( 2, 1 ) U( 3 ) = SCALE CALL SLARFG( 3, U( 1 ), U( 2 ), 1, TAU ) U( 1 ) = ONE T33 = T( J3, J3 ) ! ! Perform swap provisionally on diagonal block in D. ! CALL SLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK ) CALL SLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK ) ! ! Test whether to reject swap. ! IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1, & 1 )-T33 ) )>THRESH )GO TO 50 ! ! Accept swap: apply transformation to the entire matrix T. ! CALL SLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK ) CALL SLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK ) ! T( J1, J1 ) = T33 T( J2, J1 ) = ZERO T( J3, J1 ) = ZERO ! IF( WANTQ ) THEN ! ! Accumulate transformation in the matrix Q. ! CALL SLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK ) END IF GO TO 40 ! 30 CONTINUE ! ! N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so ! that: ! ! H(2) H(1) ( -X11 -X12 ) = ( * * ) ! ( -X21 -X22 ) ( 0 * ) ! ( scale 0 ) ( 0 0 ) ! ( 0 scale ) ( 0 0 ) ! U1( 1 ) = -X( 1, 1 ) U1( 2 ) = -X( 2, 1 ) U1( 3 ) = SCALE CALL SLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 ) U1( 1 ) = ONE ! TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) ) U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 ) U2( 2 ) = -TEMP*U1( 3 ) U2( 3 ) = SCALE CALL SLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 ) U2( 1 ) = ONE ! ! Perform swap provisionally on diagonal block in D. ! CALL SLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK ) CALL SLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK ) CALL SLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK ) CALL SLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK ) ! ! Test whether to reject swap. ! IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ), & ABS( D( 4, 2 ) ) )>THRESH )GO TO 50 ! ! Accept swap: apply transformation to the entire matrix T. ! CALL SLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK ) CALL SLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK ) CALL SLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK ) CALL SLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK ) ! T( J3, J1 ) = ZERO T( J3, J2 ) = ZERO T( J4, J1 ) = ZERO T( J4, J2 ) = ZERO ! IF( WANTQ ) THEN ! ! Accumulate transformation in the matrix Q. ! CALL SLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK ) CALL SLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK ) END IF ! 40 CONTINUE ! IF( N2==2 ) THEN ! ! Standardize new 2-by-2 block T11 ! CALL SLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ), & T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN ) CALL SROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT, & CS, SN ) CALL SROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN ) IF( WANTQ ) & CALL SROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN ) END IF ! IF( N1==2 ) THEN ! ! Standardize new 2-by-2 block T22 ! J3 = J1 + N2 J4 = J3 + 1 CALL SLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ), & T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN ) IF( J3+2<=N ) & CALL SROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ), & LDT, CS, SN ) CALL SROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN ) IF( WANTQ ) & CALL SROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN ) END IF ! END IF RETURN ! ! Exit with INFO = 1 if swap was rejected. ! 50 INFO = 1 RETURN ! ! End of SLAEXC ! END SUBROUTINE SLAGTM( TRANS, N, NRHS, ALPHA, DL, D, DU, X, LDX, BETA, & B, LDB ) ! ! -- LAPACK auxiliary routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! October 31, 1992 ! ! .. Scalar Arguments .. CHARACTER TRANS INTEGER LDB, LDX, N, NRHS REAL ALPHA, BETA ! .. ! .. Array Arguments .. REAL B( LDB, * ), D( * ), DL( * ), DU( * ), & X( LDX, * ) ! .. ! ! Purpose ! ======= ! ! SLAGTM performs a matrix-vector product of the form ! ! B := alpha * A * X + beta * B ! ! where A is a tridiagonal matrix of order N, B and X are N by NRHS ! matrices, and alpha and beta are real scalars, each of which may be ! 0., 1., or -1. ! ! Arguments ! ========= ! ! TRANS (input) CHARACTER ! Specifies the operation applied to A. ! = 'N': No transpose, B := alpha * A * X + beta * B ! = 'T': Transpose, B := alpha * A'* X + beta * B ! = 'C': Conjugate transpose = Transpose ! ! N (input) INTEGER ! The order of the matrix A. N >= 0. ! ! NRHS (input) INTEGER ! The number of right hand sides, i.e., the number of columns ! of the matrices X and B. ! ! ALPHA (input) REAL ! The scalar alpha. ALPHA must be 0., 1., or -1.; otherwise, ! it is assumed to be 0. ! ! DL (input) REAL array, dimension (N-1) ! The (n-1) sub-diagonal elements of T. ! ! D (input) REAL array, dimension (N) ! The diagonal elements of T. ! ! DU (input) REAL array, dimension (N-1) ! The (n-1) super-diagonal elements of T. ! ! X (input) REAL array, dimension (LDX,NRHS) ! The N by NRHS matrix X. ! LDX (input) INTEGER ! The leading dimension of the array X. LDX >= max(N,1). ! ! BETA (input) REAL ! The scalar beta. BETA must be 0., 1., or -1.; otherwise, ! it is assumed to be 1. ! ! B (input/output) REAL array, dimension (LDB,NRHS) ! On entry, the N by NRHS matrix B. ! On exit, B is overwritten by the matrix expression ! B := alpha * A * X + beta * B. ! ! LDB (input) INTEGER ! The leading dimension of the array B. LDB >= max(N,1). ! ! ===================================================================== ! ! .. Parameters .. REAL ONE, ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) ! .. ! .. Local Scalars .. INTEGER I, J ! .. ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. Executable Statements .. ! IF( N==0 ) & RETURN ! ! Multiply B by BETA if BETA/=1. ! IF( BETA==ZERO ) THEN DO 20 J = 1, NRHS DO 10 I = 1, N B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE IF( BETA==-ONE ) THEN DO 40 J = 1, NRHS DO 30 I = 1, N B( I, J ) = -B( I, J ) 30 CONTINUE 40 CONTINUE END IF ! IF( ALPHA==ONE ) THEN IF( LSAME( TRANS, 'N' ) ) THEN ! ! Compute B := B + A*X ! DO 60 J = 1, NRHS IF( N==1 ) THEN B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J ) ELSE B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J ) + & DU( 1 )*X( 2, J ) B( N, J ) = B( N, J ) + DL( N-1 )*X( N-1, J ) + & D( N )*X( N, J ) DO 50 I = 2, N - 1 B( I, J ) = B( I, J ) + DL( I-1 )*X( I-1, J ) + & D( I )*X( I, J ) + DU( I )*X( I+1, J ) 50 CONTINUE END IF 60 CONTINUE ELSE ! ! Compute B := B + A'*X ! DO 80 J = 1, NRHS IF( N==1 ) THEN B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J ) ELSE B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J ) + & DL( 1 )*X( 2, J ) B( N, J ) = B( N, J ) + DU( N-1 )*X( N-1, J ) + & D( N )*X( N, J ) DO 70 I = 2, N - 1 B( I, J ) = B( I, J ) + DU( I-1 )*X( I-1, J ) + & D( I )*X( I, J ) + DL( I )*X( I+1, J ) 70 CONTINUE END IF 80 CONTINUE END IF ELSE IF( ALPHA==-ONE ) THEN IF( LSAME( TRANS, 'N' ) ) THEN ! ! Compute B := B - A*X ! DO 100 J = 1, NRHS IF( N==1 ) THEN B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J ) ELSE B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J ) - & DU( 1 )*X( 2, J ) B( N, J ) = B( N, J ) - DL( N-1 )*X( N-1, J ) - & D( N )*X( N, J ) DO 90 I = 2, N - 1 B( I, J ) = B( I, J ) - DL( I-1 )*X( I-1, J ) - & D( I )*X( I, J ) - DU( I )*X( I+1, J ) 90 CONTINUE END IF 100 CONTINUE ELSE ! ! Compute B := B - A'*X ! DO 120 J = 1, NRHS IF( N==1 ) THEN B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J ) ELSE B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J ) - & DL( 1 )*X( 2, J ) B( N, J ) = B( N, J ) - DU( N-1 )*X( N-1, J ) - & D( N )*X( N, J ) DO 110 I = 2, N - 1 B( I, J ) = B( I, J ) - DU( I-1 )*X( I-1, J ) - & D( I )*X( I, J ) - DL( I )*X( I+1, J ) 110 CONTINUE END IF 120 CONTINUE END IF END IF RETURN ! ! End of SLAGTM ! END SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, & ILOZ, IHIZ, Z, LDZ, INFO ) ! ! -- LAPACK auxiliary routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! October 31, 1992 ! ! .. Scalar Arguments .. LOGICAL WANTT, WANTZ INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N ! .. ! .. Array Arguments .. REAL H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * ) ! .. ! ! Purpose ! ======= ! ! SLAHQR is an auxiliary routine called by SHSEQR to update the ! eigenvalues and Schur decomposition already computed by SHSEQR, by ! dealing with the Hessenberg submatrix in rows and columns ILO to IHI. ! ! Arguments ! ========= ! ! WANTT (input) LOGICAL ! = .TRUE. : the full Schur form T is required; ! = .FALSE.: only eigenvalues are required. ! ! WANTZ (input) LOGICAL ! = .TRUE. : the matrix of Schur vectors Z is required; ! = .FALSE.: Schur vectors are not required. ! ! N (input) INTEGER ! The order of the matrix H. N >= 0. ! ! ILO (input) INTEGER ! IHI (input) INTEGER ! It is assumed that H is already upper quasi-triangular in ! rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ! ILO = 1). SLAHQR works primarily with the Hessenberg ! submatrix in rows and columns ILO to IHI, but applies ! transformations to all of H if WANTT is .TRUE.. ! 1 <= ILO <= max(1,IHI); IHI <= N. ! ! H (input/output) REAL array, dimension (LDH,N) ! On entry, the upper Hessenberg matrix H. ! On exit, if WANTT is .TRUE., H is upper quasi-triangular in ! rows and columns ILO:IHI, with any 2-by-2 diagonal blocks in ! standard form. If WANTT is .FALSE., the contents of H are ! unspecified on exit. ! ! LDH (input) INTEGER ! The leading dimension of the array H. LDH >= max(1,N). ! ! WR (output) REAL array, dimension (N) ! WI (output) REAL array, dimension (N) ! The real and imaginary parts, respectively, of the computed ! eigenvalues ILO to IHI are stored in the corresponding ! elements of WR and WI. If two eigenvalues are computed as a ! complex conjugate pair, they are stored in consecutive ! elements of WR and WI, say the i-th and (i+1)th, with ! WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the ! eigenvalues are stored in the same order as on the 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 block, ! WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i). ! ! ILOZ (input) INTEGER ! IHIZ (input) INTEGER ! Specify the rows of Z to which transformations must be ! applied if WANTZ is .TRUE.. ! 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. ! ! Z (input/output) REAL array, dimension (LDZ,N) ! If WANTZ is .TRUE., on entry Z must contain the current ! matrix Z of transformations accumulated by SHSEQR, and on ! exit Z has been updated; transformations are applied only to ! the submatrix Z(ILOZ:IHIZ,ILO:IHI). ! If WANTZ is .FALSE., Z is not referenced. ! ! LDZ (input) INTEGER ! The leading dimension of the array Z. LDZ >= max(1,N). ! ! INFO (output) INTEGER ! = 0: successful exit ! > 0: SLAHQR failed to compute all the eigenvalues ILO to IHI ! in a total of 30*(IHI-ILO+1) iterations; if INFO = i, ! elements i+1:ihi of WR and WI contain those eigenvalues ! which have been successfully computed. ! ! ===================================================================== ! ! .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) REAL DAT1, DAT2 PARAMETER ( DAT1 = 0.75E+0, DAT2 = -0.4375E+0 ) ! .. ! .. Local Scalars .. INTEGER I, I1, I2, ITN, ITS, J, K, L, M, NH, NR, NZ REAL CS, H00, H10, H11, H12, H21, H22, H33, H33S, & H43H34, H44, H44S, OVFL, S, SMLNUM, SN, SUM, & T1, T2, T3, TST1, ULP, UNFL, V1, V2, V3 ! .. ! .. Local Arrays .. REAL V( 3 ), WORK( 1 ) ! .. ! .. External Functions .. REAL SLAMCH, SLANHS EXTERNAL SLAMCH, SLANHS ! .. ! .. External Subroutines .. EXTERNAL SCOPY, SLABAD, SLANV2, SLARFG, SROT ! .. ! .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN ! .. ! .. Executable Statements .. ! INFO = 0 ! ! Quick return if possible ! IF( N==0 ) & RETURN IF( ILO==IHI ) THEN WR( ILO ) = H( ILO, ILO ) WI( ILO ) = ZERO RETURN END IF ! NH = IHI - ILO + 1 NZ = IHIZ - ILOZ + 1 ! ! Set machine-dependent constants for the stopping criterion. ! If norm(H) <= sqrt(OVFL), overflow should not occur. ! UNFL = SLAMCH( 'Safe minimum' ) OVFL = ONE / UNFL CALL SLABAD( UNFL, OVFL ) ULP = SLAMCH( 'Precision' ) SMLNUM = UNFL*( NH / ULP ) ! ! I1 and I2 are the indices of the first row and last column of H ! to which transformations must be applied. If eigenvalues only are ! being computed, I1 and I2 are set inside the main loop. ! IF( WANTT ) THEN I1 = 1 I2 = N END IF ! ! ITN is the total number of QR iterations allowed. ! ITN = 30*NH ! ! The main loop begins here. I is the loop index and decreases from ! IHI to ILO in steps of 1 or 2. Each iteration of the loop works ! with the active submatrix in rows and columns L to I. ! Eigenvalues I+1 to IHI have already converged. Either L = ILO or ! H(L,L-1) is negligible so that the matrix splits. ! I = IHI 10 CONTINUE L = ILO IF( IILO ) THEN ! ! H(L,L-1) is negligible ! H( L, L-1 ) = ZERO END IF ! ! Exit from loop if a submatrix of order 1 or 2 has split off. ! IF( L>=I-1 ) & GO TO 140 ! ! Now the active submatrix is in rows and columns L to I. If ! eigenvalues only are being computed, only the active submatrix ! need be transformed. ! IF( .NOT.WANTT ) THEN I1 = L I2 = I END IF ! IF( ITS==10 .OR. ITS==20 ) THEN ! ! Exceptional shift. ! S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) ) H44 = DAT1*S H33 = H44 H43H34 = DAT2*S*S ELSE ! ! Prepare to use Wilkinson's double shift ! H44 = H( I, I ) H33 = H( I-1, I-1 ) H43H34 = H( I, I-1 )*H( I-1, I ) END IF ! ! Look for two consecutive small subdiagonal elements. ! DO 40 M = I - 2, L, -1 ! ! Determine the effect of starting the double-shift QR ! iteration at row M, and see if this would make H(M,M-1) ! negligible. ! H11 = H( M, M ) H22 = H( M+1, M+1 ) H21 = H( M+1, M ) H12 = H( M, M+1 ) H44S = H44 - H11 H33S = H33 - H11 V1 = ( H33S*H44S-H43H34 ) / H21 + H12 V2 = H22 - H11 - H33S - H44S V3 = H( M+2, M+1 ) S = ABS( V1 ) + ABS( V2 ) + ABS( V3 ) V1 = V1 / S V2 = V2 / S V3 = V3 / S V( 1 ) = V1 V( 2 ) = V2 V( 3 ) = V3 IF( M==L ) & GO TO 50 H00 = H( M-1, M-1 ) H10 = H( M, M-1 ) TST1 = ABS( V1 )*( ABS( H00 )+ABS( H11 )+ABS( H22 ) ) IF( ABS( H10 )*( ABS( V2 )+ABS( V3 ) )<=ULP*TST1 ) & GO TO 50 40 CONTINUE 50 CONTINUE ! ! Double-shift QR step ! DO 120 K = M, I - 1 ! ! The first iteration of this loop determines a reflection G ! from the vector V and applies it from left and right to H, ! thus creating a nonzero bulge below the subdiagonal. ! ! Each subsequent iteration determines a reflection G to ! restore the Hessenberg form in the (K-1)th column, and thus ! chases the bulge one step toward the bottom of the active ! submatrix. NR is the order of G. ! NR = MIN( 3, I-K+1 ) IF( K>M ) & CALL SCOPY( NR, H( K, K-1 ), 1, V, 1 ) CALL SLARFG( NR, V( 1 ), V( 2 ), 1, T1 ) IF( K>M ) THEN H( K, K-1 ) = V( 1 ) H( K+1, K-1 ) = ZERO IF( KL ) THEN H( K, K-1 ) = -H( K, K-1 ) END IF V2 = V( 2 ) T2 = T1*V2 IF( NR==3 ) THEN V3 = V( 3 ) T3 = T1*V3 ! ! Apply G from the left to transform the rows of the matrix ! in columns K to I2. ! DO 60 J = K, I2 SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J ) H( K, J ) = H( K, J ) - SUM*T1 H( K+1, J ) = H( K+1, J ) - SUM*T2 H( K+2, J ) = H( K+2, J ) - SUM*T3 60 CONTINUE ! ! Apply G from the right to transform the columns of the ! matrix in rows I1 to min(K+3,I). ! DO 70 J = I1, MIN( K+3, I ) SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 ) H( J, K ) = H( J, K ) - SUM*T1 H( J, K+1 ) = H( J, K+1 ) - SUM*T2 H( J, K+2 ) = H( J, K+2 ) - SUM*T3 70 CONTINUE ! IF( WANTZ ) THEN ! ! Accumulate transformations in the matrix Z ! DO 80 J = ILOZ, IHIZ SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 ) Z( J, K ) = Z( J, K ) - SUM*T1 Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2 Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3 80 CONTINUE END IF ELSE IF( NR==2 ) THEN ! ! Apply G from the left to transform the rows of the matrix ! in columns K to I2. ! DO 90 J = K, I2 SUM = H( K, J ) + V2*H( K+1, J ) H( K, J ) = H( K, J ) - SUM*T1 H( K+1, J ) = H( K+1, J ) - SUM*T2 90 CONTINUE ! ! Apply G from the right to transform the columns of the ! matrix in rows I1 to min(K+3,I). ! DO 100 J = I1, I SUM = H( J, K ) + V2*H( J, K+1 ) H( J, K ) = H( J, K ) - SUM*T1 H( J, K+1 ) = H( J, K+1 ) - SUM*T2 100 CONTINUE ! IF( WANTZ ) THEN ! ! Accumulate transformations in the matrix Z ! DO 110 J = ILOZ, IHIZ SUM = Z( J, K ) + V2*Z( J, K+1 ) Z( J, K ) = Z( J, K ) - SUM*T1 Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2 110 CONTINUE END IF END IF 120 CONTINUE ! 130 CONTINUE ! ! Failure to converge in remaining number of iterations ! INFO = I RETURN ! 140 CONTINUE ! IF( L==I ) THEN ! ! H(I,I-1) is negligible: one eigenvalue has converged. ! WR( I ) = H( I, I ) WI( I ) = ZERO ELSE IF( L==I-1 ) THEN ! ! H(I-1,I-2) is negligible: a pair of eigenvalues have converged. ! ! Transform the 2-by-2 submatrix to standard Schur form, ! and compute and store the eigenvalues. ! CALL SLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ), & H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ), & CS, SN ) ! IF( WANTT ) THEN ! ! Apply the transformation to the rest of H. ! IF( I2>I ) & CALL SROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH, & CS, SN ) CALL SROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN ) END IF IF( WANTZ ) THEN ! ! Apply the transformation to Z. ! CALL SROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN ) END IF END IF ! ! Decrement number of remaining iterations, and return to start of ! the main loop with new value of I. ! ITN = ITN - ITS I = L - 1 GO TO 10 ! 150 CONTINUE RETURN ! ! End of SLAHQR ! END SUBROUTINE SLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, & LDB, WR, WI, X, LDX, SCALE, XNORM, INFO ) ! ! -- LAPACK auxiliary routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! October 31, 1992 ! ! .. Scalar Arguments .. LOGICAL LTRANS INTEGER INFO, LDA, LDB, LDX, NA, NW REAL CA, D1, D2, SCALE, SMIN, WI, WR, XNORM ! .. ! .. Array Arguments .. REAL A( LDA, * ), B( LDB, * ), X( LDX, * ) ! .. ! ! Purpose ! ======= ! ! SLALN2 solves a system of the form (ca A - w D ) X = s B ! or (ca A' - w D) X = s B with possible scaling ("s") and ! perturbation of A. (A' means A-transpose.) ! ! A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA ! real diagonal matrix, w is a real or complex value, and X and B are ! NA x 1 matrices -- real if w is real, complex if w is complex. NA ! may be 1 or 2. ! ! If w is complex, X and B are represented as NA x 2 matrices, ! the first column of each being the real part and the second ! being the imaginary part. ! ! "s" is a scaling factor (<= 1), computed by SLALN2, which is ! so chosen that X can be computed without overflow. X is further ! scaled if necessary to assure that norm(ca A - w D)*norm(X) is less ! than overflow. ! ! If both singular values of (ca A - w D) are less than SMIN, ! SMIN*identity will be used instead of (ca A - w D). If only one ! singular value is less than SMIN, one element of (ca A - w D) will be ! perturbed enough to make the smallest singular value roughly SMIN. ! If both singular values are at least SMIN, (ca A - w D) will not be ! perturbed. In any case, the perturbation will be at most some small ! multiple of max( SMIN, ulp*norm(ca A - w D) ). The singular values ! are computed by infinity-norm approximations, and thus will only be ! correct to a factor of 2 or so. ! ! Note: all input quantities are assumed to be smaller than overflow ! by a reasonable factor. (See BIGNUM.) ! ! Arguments ! ========== ! ! LTRANS (input) LOGICAL ! =.TRUE.: A-transpose will be used. ! =.FALSE.: A will be used (not transposed.) ! ! NA (input) INTEGER ! The size of the matrix A. It may (only) be 1 or 2. ! ! NW (input) INTEGER ! 1 if "w" is real, 2 if "w" is complex. It may only be 1 ! or 2. ! ! SMIN (input) REAL ! The desired lower bound on the singular values of A. This ! should be a safe distance away from underflow or overflow, ! say, between (underflow/machine precision) and (machine ! precision * overflow ). (See BIGNUM and ULP.) ! ! CA (input) REAL ! The coefficient c, which A is multiplied by. ! ! A (input) REAL array, dimension (LDA,NA) ! The NA x NA matrix A. ! ! LDA (input) INTEGER ! The leading dimension of A. It must be at least NA. ! ! D1 (input) REAL ! The 1,1 element in the diagonal matrix D. ! ! D2 (input) REAL ! The 2,2 element in the diagonal matrix D. Not used if NW=1. ! ! B (input) REAL array, dimension (LDB,NW) ! The NA x NW matrix B (right-hand side). If NW=2 ("w" is ! complex), column 1 contains the real part of B and column 2 ! contains the imaginary part. ! ! LDB (input) INTEGER ! The leading dimension of B. It must be at least NA. ! ! WR (input) REAL ! The real part of the scalar "w". ! ! WI (input) REAL ! The imaginary part of the scalar "w". Not used if NW=1. ! ! X (output) REAL array, dimension (LDX,NW) ! The NA x NW matrix X (unknowns), as computed by SLALN2. ! If NW=2 ("w" is complex), on exit, column 1 will contain ! the real part of X and column 2 will contain the imaginary ! part. ! ! LDX (input) INTEGER ! The leading dimension of X. It must be at least NA. ! ! SCALE (output) REAL ! The scale factor that B must be multiplied by to insure ! that overflow does not occur when computing X. Thus, ! (ca A - w D) X will be SCALE*B, not B (ignoring ! perturbations of A.) It will be at most 1. ! ! XNORM (output) REAL ! The infinity-norm of X, when X is regarded as an NA x NW ! real matrix. ! ! INFO (output) INTEGER ! An error flag. It will be set to zero if no error occurs, ! a negative number if an argument is in error, or a positive ! number if ca A - w D had to be perturbed. ! The possible values are: ! = 0: No error occurred, and (ca A - w D) did not have to be ! perturbed. ! = 1: (ca A - w D) had to be perturbed to make its smallest ! (or only) singular value greater than SMIN. ! NOTE: In the interests of speed, this routine does not ! check the inputs for errors. ! ! ===================================================================== ! ! .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) REAL TWO PARAMETER ( TWO = 2.0E0 ) ! .. ! .. Local Scalars .. INTEGER ICMAX, J REAL BBND, BI1, BI2, BIGNUM, BNORM, BR1, BR2, CI21, & CI22, CMAX, CNORM, CR21, CR22, CSI, CSR, LI21, & LR21, SMINI, SMLNUM, TEMP, U22ABS, UI11, UI11R, & UI12, UI12S, UI22, UR11, UR11R, UR12, UR12S, & UR22, XI1, XI2, XR1, XR2 ! .. ! .. Local Arrays .. LOGICAL CSWAP( 4 ), RSWAP( 4 ) INTEGER IPIVOT( 4, 4 ) REAL CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 ) ! .. ! .. External Functions .. REAL SLAMCH EXTERNAL SLAMCH ! .. ! .. External Subroutines .. EXTERNAL SLADIV ! .. ! .. Intrinsic Functions .. INTRINSIC ABS, MAX ! .. ! .. Equivalences .. EQUIVALENCE ( CI( 1, 1 ), CIV( 1 ) ), & ( CR( 1, 1 ), CRV( 1 ) ) ! .. ! .. Data statements .. DATA CSWAP / .FALSE., .FALSE., .TRUE., .TRUE. / DATA RSWAP / .FALSE., .TRUE., .FALSE., .TRUE. / DATA IPIVOT / 1, 2, 3, 4, 2, 1, 4, 3, 3, 4, 1, 2, 4, & 3, 2, 1 / ! .. ! .. Executable Statements .. ! ! Compute BIGNUM ! SMLNUM = TWO*SLAMCH( 'Safe minimum' ) BIGNUM = ONE / SMLNUM SMINI = MAX( SMIN, SMLNUM ) ! ! Don't check for input errors ! INFO = 0 ! ! Standard Initializations ! SCALE = ONE ! IF( NA==1 ) THEN ! ! 1 x 1 (i.e., scalar) system C X = B ! IF( NW==1 ) THEN ! ! Real 1x1 system. ! ! C = ca A - w D ! CSR = CA*A( 1, 1 ) - WR*D1 CNORM = ABS( CSR ) ! ! If | C | < SMINI, use C = SMINI ! IF( CNORMONE ) THEN IF( BNORM>BIGNUM*CNORM ) & SCALE = ONE / BNORM END IF ! ! Compute X ! X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / CSR XNORM = ABS( X( 1, 1 ) ) ELSE ! ! Complex 1x1 system (w is complex) ! ! C = ca A - w D ! CSR = CA*A( 1, 1 ) - WR*D1 CSI = -WI*D1 CNORM = ABS( CSR ) + ABS( CSI ) ! ! If | C | < SMINI, use C = SMINI ! IF( CNORMONE ) THEN IF( BNORM>BIGNUM*CNORM ) & SCALE = ONE / BNORM END IF ! ! Compute X ! CALL SLADIV( SCALE*B( 1, 1 ), SCALE*B( 1, 2 ), CSR, CSI, & X( 1, 1 ), X( 1, 2 ) ) XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) ) END IF ! ELSE ! ! 2x2 System ! ! Compute the real part of C = ca A - w D (or ca A' - w D ) ! CR( 1, 1 ) = CA*A( 1, 1 ) - WR*D1 CR( 2, 2 ) = CA*A( 2, 2 ) - WR*D2 IF( LTRANS ) THEN CR( 1, 2 ) = CA*A( 2, 1 ) CR( 2, 1 ) = CA*A( 1, 2 ) ELSE CR( 2, 1 ) = CA*A( 2, 1 ) CR( 1, 2 ) = CA*A( 1, 2 ) END IF ! IF( NW==1 ) THEN ! ! Real 2x2 system (w is real) ! ! Find the largest element in C ! CMAX = ZERO ICMAX = 0 ! DO 10 J = 1, 4 IF( ABS( CRV( J ) )>CMAX ) THEN CMAX = ABS( CRV( J ) ) ICMAX = J END IF 10 CONTINUE ! ! If norm(C) < SMINI, use SMINI*identity. ! IF( CMAXONE ) THEN IF( BNORM>BIGNUM*SMINI ) & SCALE = ONE / BNORM END IF TEMP = SCALE / SMINI X( 1, 1 ) = TEMP*B( 1, 1 ) X( 2, 1 ) = TEMP*B( 2, 1 ) XNORM = TEMP*BNORM INFO = 1 RETURN END IF ! ! Gaussian elimination with complete pivoting. ! UR11 = CRV( ICMAX ) CR21 = CRV( IPIVOT( 2, ICMAX ) ) UR12 = CRV( IPIVOT( 3, ICMAX ) ) CR22 = CRV( IPIVOT( 4, ICMAX ) ) UR11R = ONE / UR11 LR21 = UR11R*CR21 UR22 = CR22 - UR12*LR21 ! ! If smaller pivot < SMINI, use SMINI ! IF( ABS( UR22 )ONE .AND. ABS( UR22 )=BIGNUM*ABS( UR22 ) ) & SCALE = ONE / BBND END IF ! XR2 = ( BR2*SCALE ) / UR22 XR1 = ( SCALE*BR1 )*UR11R - XR2*( UR11R*UR12 ) IF( CSWAP( ICMAX ) ) THEN X( 1, 1 ) = XR2 X( 2, 1 ) = XR1 ELSE X( 1, 1 ) = XR1 X( 2, 1 ) = XR2 END IF XNORM = MAX( ABS( XR1 ), ABS( XR2 ) ) ! ! Further scaling if norm(A) norm(X) > overflow ! IF( XNORM>ONE .AND. CMAX>ONE ) THEN IF( XNORM>BIGNUM / CMAX ) THEN TEMP = CMAX / BIGNUM X( 1, 1 ) = TEMP*X( 1, 1 ) X( 2, 1 ) = TEMP*X( 2, 1 ) XNORM = TEMP*XNORM SCALE = TEMP*SCALE END IF END IF ELSE ! ! Complex 2x2 system (w is complex) ! ! Find the largest element in C ! CI( 1, 1 ) = -WI*D1 CI( 2, 1 ) = ZERO CI( 1, 2 ) = ZERO CI( 2, 2 ) = -WI*D2 CMAX = ZERO ICMAX = 0 ! DO 20 J = 1, 4 IF( ABS( CRV( J ) )+ABS( CIV( J ) )>CMAX ) THEN CMAX = ABS( CRV( J ) ) + ABS( CIV( J ) ) ICMAX = J END IF 20 CONTINUE ! ! If norm(C) < SMINI, use SMINI*identity. ! IF( CMAXONE ) THEN IF( BNORM>BIGNUM*SMINI ) & SCALE = ONE / BNORM END IF TEMP = SCALE / SMINI X( 1, 1 ) = TEMP*B( 1, 1 ) X( 2, 1 ) = TEMP*B( 2, 1 ) X( 1, 2 ) = TEMP*B( 1, 2 ) X( 2, 2 ) = TEMP*B( 2, 2 ) XNORM = TEMP*BNORM INFO = 1 RETURN END IF ! ! Gaussian elimination with complete pivoting. ! UR11 = CRV( ICMAX ) UI11 = CIV( ICMAX ) CR21 = CRV( IPIVOT( 2, ICMAX ) ) CI21 = CIV( IPIVOT( 2, ICMAX ) ) UR12 = CRV( IPIVOT( 3, ICMAX ) ) UI12 = CIV( IPIVOT( 3, ICMAX ) ) CR22 = CRV( IPIVOT( 4, ICMAX ) ) CI22 = CIV( IPIVOT( 4, ICMAX ) ) IF( ICMAX==1 .OR. ICMAX==4 ) THEN ! ! Code when off-diagonals of pivoted C are real ! IF( ABS( UR11 )>ABS( UI11 ) ) THEN TEMP = UI11 / UR11 UR11R = ONE / ( UR11*( ONE+TEMP**2 ) ) UI11R = -TEMP*UR11R ELSE TEMP = UR11 / UI11 UI11R = -ONE / ( UI11*( ONE+TEMP**2 ) ) UR11R = -TEMP*UI11R END IF LR21 = CR21*UR11R LI21 = CR21*UI11R UR12S = UR12*UR11R UI12S = UR12*UI11R UR22 = CR22 - UR12*LR21 UI22 = CI22 - UR12*LI21 ELSE ! ! Code when diagonals of pivoted C are real ! UR11R = ONE / UR11 UI11R = ZERO LR21 = CR21*UR11R LI21 = CI21*UR11R UR12S = UR12*UR11R UI12S = UI12*UR11R UR22 = CR22 - UR12*LR21 + UI12*LI21 UI22 = -UR12*LI21 - UI12*LR21 END IF U22ABS = ABS( UR22 ) + ABS( UI22 ) ! ! If smaller pivot < SMINI, use SMINI ! IF( U22ABSONE .AND. U22ABS=BIGNUM*U22ABS ) THEN SCALE = ONE / BBND BR1 = SCALE*BR1 BI1 = SCALE*BI1 BR2 = SCALE*BR2 BI2 = SCALE*BI2 END IF END IF ! CALL SLADIV( BR2, BI2, UR22, UI22, XR2, XI2 ) XR1 = UR11R*BR1 - UI11R*BI1 - UR12S*XR2 + UI12S*XI2 XI1 = UI11R*BR1 + UR11R*BI1 - UI12S*XR2 - UR12S*XI2 IF( CSWAP( ICMAX ) ) THEN X( 1, 1 ) = XR2 X( 2, 1 ) = XR1 X( 1, 2 ) = XI2 X( 2, 2 ) = XI1 ELSE X( 1, 1 ) = XR1 X( 2, 1 ) = XR2 X( 1, 2 ) = XI1 X( 2, 2 ) = XI2 END IF XNORM = MAX( ABS( XR1 )+ABS( XI1 ), ABS( XR2 )+ABS( XI2 ) ) ! ! Further scaling if norm(A) norm(X) > overflow ! IF( XNORM>ONE .AND. CMAX>ONE ) THEN IF( XNORM>BIGNUM / CMAX ) THEN TEMP = CMAX / BIGNUM X( 1, 1 ) = TEMP*X( 1, 1 ) X( 2, 1 ) = TEMP*X( 2, 1 ) X( 1, 2 ) = TEMP*X( 1, 2 ) X( 2, 2 ) = TEMP*X( 2, 2 ) XNORM = TEMP*XNORM SCALE = TEMP*SCALE END IF END IF END IF END IF ! RETURN ! ! End of SLALN2 ! END REAL FUNCTION SLAMCH( CMACH ) ! ! -- LAPACK auxiliary routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! October 31, 1992 ! ! .. Scalar Arguments .. CHARACTER CMACH ! .. ! ! Purpose ! ======= ! ! SLAMCH determines single precision machine parameters. ! ! Arguments ! ========= ! ! CMACH (input) CHARACTER*1 ! Specifies the value to be returned by SLAMCH: ! = 'E' or 'e', SLAMCH := eps ! = 'S' or 's , SLAMCH := sfmin ! = 'B' or 'b', SLAMCH := base ! = 'P' or 'p', SLAMCH := eps*base ! = 'N' or 'n', SLAMCH := t ! = 'R' or 'r', SLAMCH := rnd ! = 'M' or 'm', SLAMCH := emin ! = 'U' or 'u', SLAMCH := rmin ! = 'L' or 'l', SLAMCH := emax ! = 'O' or 'o', SLAMCH := rmax ! ! where ! ! eps = relative machine precision ! sfmin = safe minimum, such that 1/sfmin does not overflow ! base = base of the machine ! prec = eps*base ! t = number of (base) digits in the mantissa ! rnd = 1.0 when rounding occurs in addition, 0.0 otherwise ! emin = minimum exponent before (gradual) underflow ! rmin = underflow threshold - base**(emin-1) ! emax = largest exponent before overflow ! rmax = overflow threshold - (base**emax)*(1-eps) ! ! ===================================================================== ! ! .. Parameters .. REAL ONE, ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) ! .. ! .. Local Scalars .. LOGICAL FIRST, LRND INTEGER BETA, IMAX, IMIN, IT REAL BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN, & RND, SFMIN, SMALL, T ! .. ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. External Subroutines .. EXTERNAL SLAMC2 ! .. ! .. Save statement .. SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN, & EMAX, RMAX, PREC ! .. ! .. Data statements .. DATA FIRST / .TRUE. / ! .. ! .. Executable Statements .. ! IF( FIRST ) THEN FIRST = .FALSE. CALL SLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX ) BASE = BETA T = IT IF( LRND ) THEN RND = ONE EPS = ( BASE**( 1-IT ) ) / 2 ELSE RND = ZERO EPS = BASE**( 1-IT ) END IF PREC = EPS*BASE EMIN = IMIN EMAX = IMAX SFMIN = RMIN SMALL = ONE / RMAX IF( SMALL>=SFMIN ) THEN ! ! Use SMALL plus a bit, to avoid the possibility of rounding ! causing overflow when computing 1/sfmin. ! SFMIN = SMALL*( ONE+EPS ) END IF END IF ! IF( LSAME( CMACH, 'E' ) ) THEN RMACH = EPS ELSE IF( LSAME( CMACH, 'S' ) ) THEN RMACH = SFMIN ELSE IF( LSAME( CMACH, 'B' ) ) THEN RMACH = BASE ELSE IF( LSAME( CMACH, 'P' ) ) THEN RMACH = PREC ELSE IF( LSAME( CMACH, 'N' ) ) THEN RMACH = T ELSE IF( LSAME( CMACH, 'R' ) ) THEN RMACH = RND ELSE IF( LSAME( CMACH, 'M' ) ) THEN RMACH = EMIN ELSE IF( LSAME( CMACH, 'U' ) ) THEN RMACH = RMIN ELSE IF( LSAME( CMACH, 'L' ) ) THEN RMACH = EMAX ELSE IF( LSAME( CMACH, 'O' ) ) THEN RMACH = RMAX END IF ! SLAMCH = RMACH RETURN END SUBROUTINE SLAMC1( BETA, T, RND, IEEE1 ) ! ! -- LAPACK auxiliary routine (version 2.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! October 31, 1992 ! ! .. Scalar Arguments .. LOGICAL IEEE1, RND INTEGER BETA, T ! .. ! ! Purpose ! ======= ! ! SLAMC1 determines the machine parameters given by BETA, T, RND, and ! IEEE1. ! ! Arguments ! ========= ! ! BETA (output) INTEGER ! The base of the machine. ! ! T (output) INTEGER ! The number of ( BETA ) digits in the mantissa. ! ! RND (output) LOGICAL ! Specifies whether proper rounding ( RND = .TRUE. ) or ! chopping ( RND = .FALSE. ) occurs in addition. This may not ! be a reliable guide to the way in which the machine performs ! its arithmetic. ! ! IEEE1 (output) LOGICAL ! Specifies whether rounding appears to be done in the IEEE ! 'round to nearest' style. ! ! Further Details ! =============== ! ! The routine is based on the routine ENVRON by Malcolm and ! incorporates suggestions by Gentleman and Marovich. See ! ! Malcolm M. A. (1972) Algorithms to reveal properties of ! floating-point ar